fsolve with bound constraints - transformation method

I would like to solve a system of nonlinear equations with bound constraints, i.e. lb<=x<=ub. Matlab suggests to convert the problem to a minimization and use lsqnonlin or fmincon. I want instead to use a transformation with fsolve. As an example, I picked this one:

https://www.mathworks.com/help/optim/ug/nonlinear-systems-with-constraints.html

This is a system of 2 nonlinear equations in two unknowns, with 4 different solutions. I impose the bounds x>=0 and only one solution satisfies the bounds, namely x*=(10,20). To solve the system

F(x)=0 s.t. x>=0, where F maps R^2 to R^2

I call fsolve with F(z.^2) and I find the solution z* and transform it back to x=z.^2. The square transformation makes sure that the variable x never becomes negative. As a starting point I choose a feasible value, avoiding x0=(0,0). However, the transformation method does not work. I also tried lsqnonlin with bounds and it does not work as well. It works only if I choose an initial condition very close to the true solution (10,20), but suppose in a more complicated problem I don't know what the solution is.

Q: Should I use other transformations, such as x=exp(z)? Does the transformation have to be a one-to-one function (the square is obviously not)?

Can someone help me on this? Thanks a lot!

%% fsolve with bound constraints
% Reference: https://www.mathworks.com/help/optim/ug/nonlinear-systems-with-constraints.html
% The system of equations have 4 different solutions but only one obeys the
% nonnegativity constraint x(1)>=0, x(2)>=0
% (-1,-2) NOT FEASIBLE
% (10,-2) NOT FEASIBLE
% (-1,20) NOT FEASIBLE
% (10,20) NOT FEASIBLE
% Transformation: instead of solving F(x)=0 where x is unbounded, solve for
% F(z.^2)=0 where z is unbounded but of course z.^2>=0. 
clear
clc
close 
x0 = [5,10];
%opts = optimoptions(@fsolve,'Display','off');
%% Transformation method
x0 = sqrt(x0);
[x_star,~,flag] = fsolve(@fbnd_trans,x0);

No solution found.

fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the value of the function tolerance. 

x_star = x_star.^2;
if flag<=0
    warning('Equations not solved!')
end

Warning: Equations not solved!

disp('Solution x* = ')

Solution x* = 

disp(x_star)

    0.0000   19.9996

disp('Residuals = ')

Residuals = 

disp(fbnd(x_star))

    9.5249    0.0097

%% lsqnonlin
lb = [0,0]';
ub = [inf,inf]';
[x_star,~,~,flag] = lsqnonlin(@fbnd,x0,lb,ub);

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

if flag<=0
    warning('Equations not solved!')
end
disp('Solution x* = ')

Solution x* = 

disp(x_star)

    0.6933    0.0000

disp('Residuals = ')

Residuals = 

disp(fbnd(x_star))

   15.7586   27.2438

%------------------------- SUBFUNCTIONS ----------------------------------% 
function F = fbnd(x)
F(1) = (x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
F(2) = (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));
end
function F = fbnd_trans(x)
% Impose x>=0 bound
x = x.^2;
F(1) = (x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
F(2) = (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));
end

Answer 1

[John D'Errico]

The transformation need not be one to one. For example, a common and reasonably good way to perform an optimization subject to bound constraints is a trig function map. Yes, it does it result in multiple solutions. SO WHAT? All of those multiple solutions are equivalent, so infinitely many solutions that are all equivalent is irrelevant.

In fact, this is what my fminsearchbnd function (on the File Exchange) does. It transforms the doubly bounded problem L <= x <= U into an unbounded one, using a sine transformation. Thus if theta is unbounded, then the transformation

x = (sin(theta)+1)*(U-L)/2 + L

will result in a variable x that is bounded on the desired finite interval. Of course, you can also go the other way. Given a starting value x0, then you can get the corresponding starting value in terms of theta.

theta0 = asin(2*(x0 - L)/(U - L) - 1)

Does this absolutely insure you will find a solution to the problem? Of course not. It merely insures that any solutions stay bounded inside your chosen domain.

Answer 2

[Torsten]

Transformation and initial values are two independent problems that both have to be solved. You cannot compensate "bad" initial guesses by a "good" transformation.

[Alessandro Maria Marco]

Thanks for your answer. Regarding the issue of transformation, does the transformation have to be a one-to-one function? Whenver I wish to impose a non-negativity constraint like here, is it better to use x=z^2 or x=exp(z)?

[Torsten]

The main point is that the map is differentiable and that the image under the map is your feasible region. Thus if you can exclude that x = 0 is a solution, both transformations are ok.