(OPTIMIZATION) Initial point is a local minimum - PROBLEM!!

Greetings everybody,

My code is about fitting two curves (one from calculated values and the other from laboratory data) by optimizing one parameter. The thing is that the result is always the same:

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point

is less than the value of the optimality tolerance.

<stopping criteria details>

Optimization completed: The final point is the initial point.

The first-order optimality measure, 0.000000e+00, is less than

options.OptimalityTolerance = 1.000000e-05.

After that result, the parameter never changes and of course the optimization process never gets complete and the curves never fit, and I don't really know what to do. The summary of the code is on the next lines (I'm trying with Trust-Region-Reflective and levenberg-Marquardt, but the results are the same):

%% Optimization command

x0=0.20; %% Initial value of the parameters [u, X, br]

% Vectors of lower and upper bounds on the design variables in x. Then, lb ? x ? ub

% For Levenberg-Marquardt algorithm use lb=[-inf] and ub=[inf]

lb=-inf;

ub=inf;

% optimization options: The first one is levenberg-Marquardt, and the

% second, by default is Trust-Region-Reflective

% options=optimset('Algorithm','levenberg-Marquardt', 'LargeScale','off', 'DiffMaxChange',0.01, 'DiffMinChange',0.0001, 'TolFun',1e-5, 'TolX',0.001);

options=optimset('LargeScale','off', 'DiffMaxChange',0.01, 'DiffMinChange',0.0001, 'TolFun',1e-5, 'TolX',0.001);

[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@HSfun,x0,lb,ub,options);

%% Show Results and Optimization Details

dim=1; %make sure to use same values a fun file

X=x.*dim;

resnorm;

exitflag;

output;

HSfun basically takes the vectors (the mentioned curves) and calculates the difference in order to reduce that gradient:

%% Read Calculated raw data from CalculatedData.xls file

calc_sh_strain=xlsread('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\CalculatedData.xls','C46:C235')*100; % strain in percentage

calc_q=xlsread('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\CalculatedData.xls','D46:D235'); % calc q

%% Interpolate Calculated data to get values at same strain as observed data

sh_strain_transp=transpose(lab_sh_strain); % sames strain values as observed data

q_calculated=interp1(calc_sh_strain,calc_q,sh_strain_transp,'spline');

%% calculate the difference between observed and calculated response

F=lab_q-q_calculated';

I would really appreciate if somebody can help me, maybe I'm not using the best optimization tool, or stopping criteria are not correct.

Thanks in advance.

Answer 1

[Alan Weiss]

I suggest that you plot the square of the objective function over a reasonable range of the parameter and over a small range such as from 1 to 1.5. It is possible that the objective function is either flat or stair-stepped, meaning locally flat. In either case you will not get a good solution using lsqnonlin. If the objective function is stair-stepped, use fminbnd on the square of the function as the optimizer, recognizing that there are many local minima. If the objective function is flat, well, you'll have to debug it.

Alan Weiss

MATLAB mathematical toolbox documentation

[Daniel Valencia]

Greetings Mr Weiss,

Thanks a Lot for your help. About your suggestions I also have some questions (I'm a noob on these topics): 1) When You mention the reasonable range, are You talking about the lower and Upper bound limits? 2) what do you mean by "debug it"?

Again, thanks a Lot.

[Alan Weiss]

I mean your optimization variable, the parameter you are optimizing. Just plot the objective function over a range that you thtink is reasonable for your problem.

What I mean if you see that the objective function is a constant is you'll have to figure out what you did wrong. See Debugging and Analysis.

Alan Weiss

MATLAB mathematical toolbox documentation

[Daniel Valencia]

I understand Mr Weiss. One more thing... Is there any way to know if the stopping criteria such as Tol, diffmaxchange, and others, are suitable for the specific análisis i'm trying to run? That stopping criteria I posted up in the initial question was taken from someone else code (with a simmilar approach and engineering topic), and I'm using it in my code.

[Alan Weiss]

I see no real problems with the options, outdated as they are (the recommended option-setting function is now optimoptions, not optimset, and several of your options are outdated such as LargeScale and DiffMinChange and DiffMaxChange), but there should be no issue with using the ones you stated. The real problem lies elsewhere, as I tried to explain.

Alan Weiss

MATLAB mathematical toolbox documentation

[Daniel Valencia]

Thanks again Mr Weiss, you really helped me.

[Daniel Valencia]

Hello Again mr Weiss,

I just tried fminbnd as you suggested but the thing is that my function returns a vector not a scalar value, and that vector is the difference between the calculated and objective curves. I must say that the function is not an equation, it just brings data from excel sheets and one software (which makes the computing work). Which solver could be suitable for my problem, having that into account?

On the other hand, I updated the code with the last optimoptions instead of optimset, but result is still the same.

[Alan Weiss]

I'm sorry that you wasted your time changing to optimoptions. As I told you, that has nothing to do with your problem. And if you are trying to use fminbnd instead of lsqnonlin, you need to change your objective function to an explicit sum of squares. For an example, see Nonlinear Data-Fitting.

What was the result of plotting your objective function? Is it flat or stair-stepped as I suspect?

Alan Weiss

MATLAB mathematical toolbox documentation

[Daniel Valencia]

No no no, it was not a waste of time, it was something new I learned and now my script is updated, so thanks again.

I would not know how to convert a whole curve (a pair of vectors) to a single value and still get the optimization done. What my function actually does is to calculate F=(LaboratoryVector - CalculatedVector), and try to get to zero or at least reduce as much as possible that difference between values. I do not even know if that is a correct way to use the function. All examples I see in Matlab come from the programmed equations inside the script in order to be the function, which makes totally sense, but in my case programming what the software does is almost impossible and totally out of the scope of my research.

Answering your question, I'm trying to optimize not just one but three parameters, so it might be a little bit complicated, still, I'm showing an example of the curve I'm trying to fit. As I told you, that plot comes from the software results I have linked to the matlab script.

That is the general shape of the curve, and for the single value I'm refereing to, the variation varies from 0.05 to 0.35 in it's real physical range.

[Daniel Valencia]

Plase let me know if it would be useful to post the whole script or part of it, and thank in advance!

[Alan Weiss]

I have misunderstood several things about your problem. You now say that you have three parameters. When I first looked at your code I saw the lines

x0 = 0.20;
lb = -inf;
ub = inf;

All these indicate to me a single parameter to be optimized. If you really have three parameters then fminbnd is totally inappropriate, and I would not have suggested it.

It now sounds to me as if you have code that takes your parameter vector x and calculates F(x), which is a vector or matrix of response values. You then calculate F(x) - data, where data is measurements that you have. You are trying to change x to get the smallest deviation, where deviation is measured as the sum of squares of F(x) - data. I hope that I have a correct understanding now.

What I have been suggesting all along and I still have no idea whether you have tried this is for you to plot F(x) as x changes in some way. I don't care if x is a 3-D vector or not, you can change it in just one component or along a line. The question is whether F(x) changes at all (is it flat?) or whether it changes in a stair-stepping kind of way.

• If your response F(x) is flat, then you need to debug your code.
• If your response F(x) is stair-stepped, then you need to use a solver that is appropriate. I would usually suggest patternsearch to minimize the sum of squares of deviations for a nonsmooth multivariable function. But maybe you don't have Global Optimization Toolbox. In that case you might have to simply calculate F at a variety of points and look for a minimum manually. Or maybe a scaled version of fminsearch could work. We can discuss details later if necessary.
• If the response is smoothly differentiable then I am surprised, because your initial problem is that the first-order optimality measure (essentially, the gradient) is zero, meaning there is no change in the response for small changes in x.

Alan Weiss

MATLAB mathematical toolbox documentation

[Daniel Valencia]

Greetings mr Weiss.

I first would like to apologize about the description I first gave you. I first talked about only one parameter in order to make the explanation a little easier, but now I see it was actually an important fact. On the other hand, until today I had no access to the software again, so I could not get the curve you were asking for, but the following image summarizes the function's variation with respect to that parameter. Physically, the parameter cannot go over 0,35, so that is the plot thereshold.

"It now sounds to me as if you have code that takes your parameter vector x and calculates F(x), which is a vector or matrix of response values. You then calculate F(x) - data, where data is measurements that you have. You are trying to change x to get the smallest deviation, where deviation is measured as the sum of squares of F(x) - data. I hope that I have a correct understanding now" R//: Yes, essentialy what the script does is exactly that, but the F(x) is not considered as the calculated data but the subtraction result as you can see here:

Main script (example for only one parameter):

clc

clear

format short g

%% Variables to record the parameters during the iterations

global n count iter OptParameters

n=1; %% number of parameters to be optimized

count=0;

iter=0;

OptParameters=zeros(1,n); %% Matrix to store parameters

%% Optimization command

x0=0.20; %% Initial value of the parameters [u, X, br]

% Vectors of lower and upper bounds on the design variables in x. Then, lb ? x ? ub

% For Levenberg-Marquardt algorithm use lb=[-inf] and ub=[inf]

lb=-inf;

ub=inf;

% optimization options:

options=optimset('TypicalX', 0.1, 'ScaleProblem', 'Jacobian','FinDiffRelStep', 0.01, 'Display', 'iter-detailed', 'Algorithm','levenberg-marquardt', 'MaxIter',10, 'FinDiffType','forward', 'TolFun', 1e-3,'TolX' ,1e-3);

[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@HSfun,x0,lb,ub,options);

%% Show Results and Optimization Details

dim=1; %make sure to use same values a fun file

X=x.*dim;

resnorm;

exitflag;

output;

%% Read and plot observed and calculated data (q vs axial strain)

lab_sh_strain=xlsread('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\LabData.xlsx','A1:A2992');

lab_q=xlsread('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\LabData.xlsx','B1:B2992');

calc_sh_strain=xlsread('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\CalculatedData.xls','C46:C235')*100; % strain in percentage

calc_q=xlsread('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\CalculatedData.xls','D46:D235');

figure(1)

plot(lab_sh_strain,lab_q,'rx',calc_sh_strain,calc_q,'ks')

xlabel("Shear Strain [%]")

ylabel("Deviatoric Stress q [kPa]")

legend("Lab data","Calculated data")

Function script @HSfun:

function F=HSfun(x)

%% Record the parameter values during the optimization

global n count iter OptParameters

dim=1;

count=count+1;

if count==iter*n+iter+1

iter=iter+1;

OptParameters(iter,:)=x.*dim;

end

dlmwrite('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\OptParameters.txt',OptParameters,'delimiter','\t','newline','pc','precision','%15.10f');

%% Read observed data

lab_sh_strain=xlsread('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\LabData.xlsx','A1:A2992');

lab_q=xlsread('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\LabData.xlsx','B1:B2992');

%% Calculated values

X=x.*dim;

u=X;

dlmwrite('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\u.txt',u,'newline','pc','precision','%15.10f');

%%% Macro to write parameters into PLAXIS project

!C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\WriteParameters.exe

%%% Macro to run PLAXIS and extract calculated data into CalculatedData.xls file

!C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\RunPlaxis.exe

%% Read Calculated raw data from CalculatedData.xls file

calc_sh_strain=xlsread('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\CalculatedData.xls','C46:C235')*100; % strain in percentage

calc_q=xlsread('C:\Users\Afuribeh\Desktop\UNAL-MED-DFV\U TXC TEST - Only u\CalculatedData.xls','D46:D235'); % calc q

%% Interpolate Calculated data to get values at same strain as observed data

sh_strain_transp=transpose(lab_sh_strain); % sames strain values as observed data

q_calculated=interp1(calc_sh_strain,calc_q,sh_strain_transp,'spline');

%% calculate the difference between observed and calculated response

F=lab_q-q_calculated';

The macros you see called in the function script run the software and take the output data to Excel sheets in order to be read by MatLab, and repeat it again for each iteration.

Thanks again for all your support!

[Alan Weiss]

We are almost there. What I have been asking you to do is calculate

fun = @(x)sum(HSfun(x).^2,'all');

and then plot fun for a small range of x values. But we don't really need a plot to see if my suspicions are correct. Set x0 as your initial point, a 3-D vector. Then try the following.

delta = 1e-8*ones(size(x0)); % a small vector
D0 = norm(fun(x0+delta) - fun(x0))
delta2 = 1e7*delta; % a larger vector
D1 = norm(fun(x0+delta2) - fun(x0))

I am expecting that D0 = 0 and that D1 might equal zero or might not.

• If D0 = 0 then you are at least in the stair-stepping case I have been describing, where small changes in the parameter vector x lead to zero changes in the response.
• If D1 = 0 then you are in the case where the response is independent of the input variable x and you need to debug your program.

Alan Weiss

MATLAB mathematical toolbox documentation

[Daniel Valencia]

Hello again Alan,

Here is the curve:

For

delta = 1e-8*ones(size(x0)); % a small vector

D0 = norm(fun(x0+delta) - fun(x0))

delta2 = 1e7*delta; % a larger vector

D1 = norm(fun(x0+delta2) - fun(x0))

Since x0 is only one parameter, size of delta and delta2 will always be 1 and it will be a scalar, so I think the norm of the norm will be exactly equal to the value fun(x0+delta) - fun(x0) or fun(x0+delta2) - fun(x0). Please tell me if I'm missing something.

As you can see in the graph, sum of squares are large values because the function is not even close to fit. Please recall that evaluated function was F=calculateddata-measureddata.

Thanks in advance.

[Alan Weiss]

I'm sorry, but I cannot help you if you do not provide the information I request. Please run the code I asked you to run with your real x0 and real objective function.

Alan Weiss

MATLAB mathematical toolbox documentation

[Daniel Valencia]

Greetings mr Weiss,

The result of the evaluation you were asking for is D1=0 and D2 different from zero, so our function is stair-stepped.

How could this be interpreted?

[Alan Weiss]

Finally, we understand what is happening. Your objective function does not change at all when your optimization parameters change by a tiny amount. That is what D1 = 0 means.

The way that lsqnonlin works (in one dimension, higher dimensions are exactly similar) is that it starts at x0 and evaluates F(x0). It then adds a tiny amount delta and evaluates F(x0+delta) and computes (F(x0+delta) - F(x0))/delta, which is an estimate of the gradient of the function. If this estimate is exactly zero, which you have now established, then the solver cannot move. It knows that the function is flat near x0, meaning if the solver samples the function values near x0 they are all the same.

How can this be? It is a typical result of a simulation; many simulations are not sensitive to small changes in the evaluation point. It also arises in many other circumstances. I don't know why it is happening with your objective function, but it is what I suspected from the start because it happens in many circumstances.

The point is, what can you do about it? You can use a different solver such as patternsearch from Global Optimization Toolbox. In that case you will have to change your objective function to an explicit sum of squares:

fun = @(x)sum(HSfun(x).^2,'all');

Another thing you can try is to set the finite difference step size to a larger-than-default value, such as

options = optimoptions('lsqnonlin','FiniteDifferenceStepSize',1e-3);

For a discussion of all this, see Optimizing a Simulation or Ordinary Differential Equation, the first sections about problems in finite differences.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation