You do not give the initial value used for ? in your first NDSolve. I assume there was some trial value and I picked zero.

You seem to be using { and } in an unconventional way to try to limit the scope of variable substitution. { and } are almost exclusively used for list construction. You might try using ( and ) instead.

For what seems to be the essential bit of your question, you might look at using

FindMinimum[d[t] - (?1[t]*(y1[t] /. s) + ?2[t]*(y2[t] /. s)),{t, 50, 0, 100}]

which should give you a non-negative value if your constraint is satisfied and a negative value otherwise. Look at the information hidden behind Details and Options on the documentation page for FindMinimum to see how this should restrict the search for the minimum to lie between 0 and 100. That is only guaranteed to find local minima. There are a few other functions in Mathmatica, but I don't think they may guarantee a global minima with several interpolating functions. A function that has an option of using random search might stumble onto a minima less than zero when others might be trapped in local minima greater than zero.

I was able to plot the function you desired by using

Plot[First[d[t] - (\[Rho]1[t]*(y1[t] /. s) + \[Rho]2[t]*(y2[t] /. s))], {t, 0, 100}]

but I can't be certain that the use of First has not missed essential behavior perhaps hiding behind other solutions.

Hopefully you will see how you might adapt this to serve your purposes. All this should work with version 9.0.