As I indicated last time you wrote, the best way to unravel these situations is to evaluate the objective function on a grid to make sure it is doing what you expect. You should also ensure that the various functions correctly handle out-of-domain errors (such as log(x) when x<=0). If there are domain issues, you need to add boundary constraints for your objective function.
The following statements plot the objective function. The objective function appears to maximize as x -> 0 from the right, but it is undefined for x <= 0, so that need to be added to the optimization.
proc iml;
start f1(x) global(a,b,y);
return(a*y-b##x);
finish;
start f2(x);
bounds = {0,1000};
return( froot( "f1", bounds) );
finish;
start f3(y) global(a,b);
z = f2(y);
return(z*a-y##2);
finish;
a=2; b=0.5;
y = do(0.05, 2, 0.05);
f = j(1, ncol(y), .);
do i = 1 to ncol(y);
f[i] = f3(y[i]);
end;
print (f[1:5]);
call series(y, f);
y0 = 1;
opt = {1,2};
blc = {1e-6, /* 1e-6 <= y < infinity */
.};
call nlpnra(rc, result, "f3", y0, opt) blc=blc;
Using your notation, I suggest the following:
1. You've set f1 to be a function of x, with a, b, y as parameters. That means that f2 is also a function of x, with a,b, and y as parameters.
2. It sounds like you want to optimize over the parameter y where a and b are fixed.
3. As I mentioned last time, I've written about similar examples. Since you didn't like my previous suggestions, try this one: "Optimizing a function that evaluates an integral"