fsolve with bound constraints - transformation method
12 views (last 30 days)
Show older comments
Alessandro Maria Marco
on 24 Feb 2024
Edited: John D'Errico
on 24 Feb 2024
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:
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);
x_star = x_star.^2;
if flag<=0
warning('Equations not solved!')
end
disp('Solution x* = ')
disp(x_star)
disp('Residuals = ')
disp(fbnd(x_star))
%% lsqnonlin
lb = [0,0]';
ub = [inf,inf]';
[x_star,~,~,flag] = lsqnonlin(@fbnd,x0,lb,ub);
if flag<=0
warning('Equations not solved!')
end
disp('Solution x* = ')
disp(x_star)
disp('Residuals = ')
disp(fbnd(x_star))
%------------------------- 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
0 Comments
Accepted Answer
John D'Errico
on 24 Feb 2024
Edited: John D'Errico
on 24 Feb 2024
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.
0 Comments
More Answers (1)
Torsten
on 24 Feb 2024
Moved: Torsten
on 24 Feb 2024
Transformation and initial values are two independent problems that both have to be solved. You cannot compensate "bad" initial guesses by a "good" transformation.
2 Comments
Torsten
on 24 Feb 2024
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.
See Also
Categories
Find more on Mathematics and Optimization in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!