Is it possible to set an upper limit on the step size used by fsolve?

I'm trying to use fsolve to solve a small nonlinear system (3 equations) for which I have constraints on the unknowns. I know I can't explicitely enforce constraints using fsolve, but I always have access to a good initial guess, which generally allows the solver to find the solution I need without looking outside the constraints boundaries.

However, changing the equation parameters (which I need to be able to do) sometimes results in fsolve trying to find a solution outside the boundaries, which is not acceptable.

I've tried using fmincon with the same initial guesses as fsolve without success: without applying constraints, it tends to search outside the constraints boundaries and make the solution crash faster than fsolve. With the constraints, it's unable to find the solution because it gets stuck on a boundary.

I'm thinking that if I could add some kind of relaxation in the fsolve algorithm, it could help me stay inside my constraints (without guarantees, of course). In other words, I'm trying to find if there's a way to limit the step size taken by fsolve from one iteration to another. I've found the following options in the documentation, but they do not do what I need:

- FiniteDifferenceStepSize: this allows to set the step size used to estimate the gradient inside a given fsolve iteration, but offers no control on the step taken on the solution vector from one iteration to another

- StepTolerance: this is a lower bound on the step size from one iteration to another (what I'm looking is an upper bound on that value)

Basically, I'd like fsolve to estimate the gradient, but let me decide how far along that gradient I want to go. Is there any other option, or even other solver I didn't think of that could do something like that?

Thank you

Answer 1

[Matt J]

I've tried using fmincon with the same initial guesses as fsolve without success:

A simpler option than fmincon would be to use lsqnonlin, which is quite similar to fsolve except that it allows you to stipulate simple upper and lower bounds.

With the constraints, it's unable to find the solution because it gets stuck on a boundary.

Based on your description, I'm skeptical that the solver is "getting stuck". It sounds to me like the unconstrained solution really is outside the boundaries, and the constrained solution really is at the boundary. If I'm correct, tweaking the search mechanism will not help. Your equations themselves are the problem. There is some error or physical information lacking in the equations that would lead to a solution in the expected region.

One way to verify this would be to sample your equation function on a grid covering the bounded region. This is something that should be computationally tractable since you have only 3 unknowns. This will let you inspect directly where the solution lies (to within the resolution of the sampling grid).

[Jean-Philippe Gauthier]

Thank you for your input @Matt J. Here is some additional info regarding my problem.

A simpler option than fmincon would be to use lsqnonlin, which is quite similar to fsolve except that it allows you to stipulate simple upper and lower bounds.

I've actually also looked into lsqnonlin, but my constraints are more complex than lower and upper bounds.

Based on your description, I'm skeptical that the solver is "getting stuck". It sounds to me like the unconstrained solution really is outside the boundaries, and the constrained solution really is at the boundary. If I'm correct, tweaking the search mechanism will not help. Your equations themselves are the problem. There is some error or physical information lacking in the equations that would lead to a solution in the expected region.

One way to verify this would be to sample your equation function on a grid covering the bounded region. This is something that should be computationally tractable since you have only 3 unknowns. This will let you inspect directly where the solution lies (to within the resolution of the sampling grid).

I never used fmincon before, so the first thing I did is test it on the following toy problem:

Here the code I wrote for that (inspired from https://www.mathworks.com/help/optim/ug/nonlinear-systems-with-constraints.html#NonlinearSystemsWithConstraintsExample-5)

% Initial guess and bounds

x0 = [0.6,0.2];

%x0 = [0.1,0.1];

% fsolve solution

opts = optimoptions(@fsolve,'Display','off','FunctionTolerance',1e-12);

x_fsolve = fsolve(@fbnd,x0,opts);

% fmincon solution

A = [1,1];

b = 1;

lb = [0,0];

opts = optimoptions(@fmincon,'Algorithm','interior-point',...

'Display','off',...

'OptimalityTolerance',1e-12);

x_fmincon = fmincon(@(x)0,x0,A,b,[],[],lb,[],@fminconstr,opts);

% Plot

x = linspace(-1,2,1001);

y1 = 0.1./(x.^2);

y2 = sin(1.5*pi*x);

figure(1)

clf(1)

hold on

plot(x,y1)

plot(x,y2)

plot([0,1,0,0],[0,0,1,0],'k')

plot(x0(1),x0(2),'ok','MarkerFaceColor','k')

plot(x_fsolve(1),x_fsolve(2),'xk','MarkerSize',10)

plot(x_fmincon(1),x_fmincon(2),'+k','MarkerSize',10)

axis([-0.1,1.3,-0.1,1.1])

xlabel('x_1')

ylabel('x_2')

legend('x^2_1 x_2 - 0.1 = 0', 'sin(1.5\pix_1) - x_2 = 0', ...

'Constraint', 'Initial guess', 'fsolve', 'fmincon')

% Objective functions

function F = fbnd(x)

F(1) = x(1)^2*x(2) - 0.1;

F(2) = sin(1.5*pi*x(1)) - x(2);

end

% Objectives functions as nonlinear equality constraints for fmincon

function [c,ceq] = fminconstr(x)

c = [];%constraints(x); % constraint is fmincon nonlinear inequality constraints

ceq = fbnd(x); % fsolve objective is fmincon nonlinear equality constraints

end

There is exactly one solution inside the constraints boundaries, which both fsolve and fmincon are able find if I start close enough:

With another initial guess, I get different results, which is not surpring. There actually is a solution inside the boundaries, but it's still trying to find the one which is outside by travelling along the solution of the 2nd equation (red curve) up the constraint. This is what I meant by fmincon getting stuck at the boundary: it converges to an infeasible point. I still tried fmincon on my actual problem to see what would happen but got the same result.

I suppose this is a limitation of a local optimization algorith such as fmincon? Are there any global constrained algorithms available?

Back to my actual problem, I can elaborate a little: it is part a of numerical simulation code to model tilting-pad thrust bearing. Without going into details, I iteratively resolve PDEs for the pressure and temperature distributions inside a thin oil film between two rigid surfaces, one of which (the pad) can freely tilt on a point pivot, until global convergence is achieved for all equations.

The film geometry, on which the pression distribution is highly dependent, is described by 3 parameters which are the thickness at the pivot location and the two tilt angles of the pad. During each global iteration of the solution, I need to optimize the film geometry and find the equilibrium point where the pressure distribution matches the applied axial load/force (model input) and generate zero moments on the pad around both tilting axes. This is my nonlinear system of 3 equations. The constraint is that the local film thickness must be positive everywhere (otherwise the pad would penetrate the opposite surface, which doesn't make physical sense).

After solving this system, the temperature is updated for the new film geometry, which will in turn affect the pressure field and so a new global iteration is started. To aid with convergence, the 3 geometric parameters from the previous global iteration are used as an initial guess for solving the nonlinear system at the current global iteration.

As I've mentionned in my initial post, this works fine with fsolve in most cases (different bearing designs).

One curious observation:

For a given case that is problematic, fsolve works well at the beginning of the simulation: for the first ~20 global iterations, it always find the correct solution without looking outside the constraints (which are not enforced, obviously). But at some point, when everything is converging smoothly, it starts to take very large steps and searches for a solution where the film thickness becomes locally negatives. This is why I was looking for a way to limit the step size.

What is curious is that for this same given case, if I always reset the initial guess for fsolve at the values used on iteration 0, fsolve never tries to take large steps and always remain inside the constraints, and the global simulation converges. Of course this requires more solutions of the pressure PDE as fsolve always start further from the desired solution.

I'm not sure why this would happen, but is it possible that starting fsolve with an initial guess that is very close to the final solution would lead to problems? This is what seems to be happening when I update the initial guess at each global iteration: as the global solution converges, changes are becoming very small from one global iteration to the next. It's really the only difference between a diverging and converged solution (for this particular case at least).

In any case, as you suggested, I will try to investigate more thoroughly my constrained parameter space to see if I can find something strange, but I'm not sure I will. The physics of the problem should allow exactly one equilibrium state to exist, or zero if you have a bad bearing design. But then again, I might be surprised.

Thanks again for your help, and I'll gladly look into any other suggestion you might have.

[Matt J]

But at some point, when everything is converging smoothly, it starts to take very large steps and searches for a solution where the film thickness becomes locally negatives.

Sounds like some discontinuous change occurred in the capture regions of your optimization problem. In other words, you are hoping that the solution to the previous optimization will give you an initial point lying in the capture region of the global solution to the new problem. This may be true much of the time, but there really is no gaurantee.

One solution may be to come up with a more independent way to generate an initial guess. For your toy problem, it is very simple to do so, as below. I don't know if something analogous is possible in your true problem.

% Initial guess and bounds

x1=linspace(0,1,1000);
x2=sin(1.5*pi*x1);

loc=abs(x1.^2.*x2 - 0.1)<=1e-4 & x1>=0 & x2>=0 & x1+x2<=1;

x0=[x1(loc), x2(loc)]

x0 = 1×2

0.6086 0.2702

% fsolve solution
opts = optimoptions(@fsolve,'FunctionTolerance',1e-12,'Display','off');
x_fsolve = fsolve(@fbnd,x0,opts)

x_fsolve = 1×2

0.6087 0.2699

% fmincon solution
A = [1,1];
b = 1;
lb = [0,0];
opts = optimoptions(@fmincon,'Algorithm','interior-point','Display','off',...
                             'OptimalityTolerance',1e-12);
x_fmincon = fmincon(@(x)norm(fbnd(x)).^2,x0,A,b,[],[],lb,[],[],opts)

x_fmincon = 1×2

0.6087 0.2699

% Objective functions
function F = fbnd(x)
F(1) = x(1)^2*x(2) - 0.1;
F(2) = sin(1.5*pi*x(1)) - x(2);
end
% Objectives functions as nonlinear equality constraints for fmincon

if I always reset the initial guess for fsolve at the values used on iteration 0, fsolve never tries to take large steps and always remain inside the constraints, and the global simulation converges. Of course this requires more solutions of the pressure PDE as fsolve always start further from the desired solution.

Maybe you can use this as a fallback initialization strategy only in the case where the optimization fails. Maybe it will be necessary only in rare instances, when discontinuous changes in the capture regions occur.

[Jean-Philippe Gauthier]

One solution may be to come up with a more independent way to generate an initial guess. For your toy problem, it is very simple to do so, as below. I don't know if something analogous is possible in your true problem.

I've tried doing something similar by evaluating a few points inside the constrained parameter space and keeping the best point as the initial guess for fsolve, but it's not very robust. I don't have explicit expressions for my equations, which also change between global iteration as other parts of the solution evolve. The best option I've found so far really is to use the same initial guess every time.

Maybe you can use this as a fallback initialization strategy only in the case where the optimization fails. Maybe it will be necessary only in rare instances, when discontinuous changes in the capture regions occur.

I think this is the best approach right now. The code performs very well for most cases as it is, so I'll just include special treatment when this problem, which should indeed be rare, arises.

Answer 2

[Matt J]

You could try playing with the InitDampling option parameter used by the Levenberg-Marquardt method.