MATLAB: Does fsolve depend on the stability of the system

fsolveMATLABOptimization Toolboxstabilitysystem of nonlinear equationstrust-region algorithms

Hello all,

I am trying to solve a system of non linear equations using fsolve; lets say

F(x;lambda) = 0, where lambda is a vector of parameters, and x the vector I want to solve for.

I have 2 values of the parameter lambda, that I want to solve the system for. For the one value of lambda I get a solution, which seems alright.

For the other value of lambda I get a solution again (matlab exits with a flag of 1. However I know this is not an actual solution For example I know that some of the dimensions of x have to be equal to each other, and this is not the case in the solution I get from fsolve.

I have tried both trust-region and the levenberg-marquardt algorithm, and I am not getting any better results. (explicitly enforcing those x's to be the same, still seems to give solutions that are not consistent with what I would be expecting from the properties of the system)

My question is: do the algorithms used by fsolve depend on any kind of stability of the system? Could it be that changing the parameter lambda in the second case I mention above, I make the system unstable, and could that make fsolve unable to solve it correctly?

Thank you, George

Best Answer

Well, the first formulation in your Comment to Alan definitely looks like it could be unstable if the g_i(x)--->0 as x becomes large. That would also mean that f_i(x)--->0 with large x and thus would be insensitive to changes in x in that region.

You could try running fsolve with your second formulation, but that formulation is non-differentiable and undefined for g_i(x)<0, so you may need to move to a constrained solver.

Related Solutions

MATLAB: How to evaluate a unknown variable

Using fsolve:

alpha=0.9;
g=[3, 2, 1];
N=10;
P=5;
f = @(lambda) sum(( lambda*(2-alpha)-sqrt(lambda*alpha*(lambda*alpha+g.*4*(1-alpha))) ) ...
    ./ ( g.*2*lambda*(alpha-1) )) - N*P;
lambda0 = 1;
options = optimoptions('fsolve','Display','iter','Algorithm','levenberg-marquardt');
lambda = fsolve(f, lambda0, options);

Learn more at https://www.mathworks.com/help/optim/ug/fsolve.html

MATLAB: Does the objective function of the optimization problem “drift away” from an optimal solution in the Optimization Toolbox

There are several reasons why your objective function might suddenly increase.

What may be happening is that FMINCON is getting near the solution and is trying to model the function as a quadratic (i.e. the SQP sub-problem) but cannot do so successfully. This means that the quadratic model doesn't fit the objective function well near the solution: perhaps it is not convex, or there is a loss of continuity, or it is trying to meet the constraints but is unable to do so.

There may also be spikes in the computed gradient which make the magnitude of the directional derivative and the 1st order optimality incorrect. That would fool FMINCON into thinking that it wasn't at a stationary (flat) region and that it should continue.

This is also something that can commonly happen if you reduce the tolerances to unnecessarily small values. Basically, at such low tolerances, there is some numerical "fuzziness" near the optimal solution. FMINCON is near the optimum, but because of the very low exit tolerances it is unable to converge; at the next iteration, FMINCON is getting near the solution and is trying to model the function as a quadratic: this is the SQP iteration, (i.e. the cost function needs to be modelled as a quadratic in the neighborhood of the current iteration); it is trying to approximate the Hessian (i.e. the matrix of second derivatives) as a positive definite matrix but cannot do so successfully. This indicates that the quadratic model doesn't fit the your function well near the solution (i.e. the quasi-Newton algorithm fails). In certain cases, this could lead to the SQP iteration maxing out and the Hessian being modified twice. The message "Hessian modified twice; MaxSQPIter" then appears, indicating the SQP failure.

From this iteration onward the 'search directions' are not be optimal, as the quasi-Newton method failed. This would explain the divergence from the optimal solution.

If the algorithm drifts away from the optimal solution at very tight tolerances, there are a few options. The most obvious workaround is to loosen the tolerances. If you wish to maintain low values for 'TolFun' and 'TolX', however, there may be other things you can do, depending on your specific case.

The round-off errors in the function evaluations when requiring such strict tolerances are tantamount to introducing discontinuities in the objective function causing the gradient to vary greatly: this explains the large value of the directional derivative at the "bad" iteration, causing the next step to be far away from the optimum point. One possibility is therefore to adjust the finite-difefrence perturbation level 'DiffMinChange' in the gradient evaluation.

The last option, if it is applicable to your situation, is to work around the problems with finite-differences altogether and specify the gradient and Hessian analytically.

Related Question