fsolve with bound constraints - transformation method
10 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Alessandro Maria Marco
el 24 de Feb. de 2024
Editada: John D'Errico
el 24 de Feb. de 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 comentarios
Respuesta aceptada
John D'Errico
el 24 de Feb. de 2024
Editada: John D'Errico
el 24 de Feb. de 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 comentarios
Más respuestas (1)
Torsten
el 24 de Feb. de 2024
Movida: Torsten
el 24 de Feb. de 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 comentarios
Torsten
el 24 de Feb. de 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.
Ver también
Categorías
Más información sobre Solver Outputs and Iterative Display en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!