Using Multistart in an unconstrained optimization

Hello, suprisingly working with this model I have imposed nonlinear and upper/lower bounds in a separate script and have come up with a better fit. I am surprised to find that the unconstrained case, the code I have here, is not giving me correct solutions. I tried messing with ms.XTolerance and the functionTolerance but nothing helped (since I have only seen an exit flag of 2). It does not seem like multistart is minimizing my objective function, as the error for a solution is a massive number.

Is there anything you all see wrong, or an explanation as to why this is happening? (I would preferably use lsqnonlin, not lsqcurvefit just to stay consistent)

 D = readmatrix('Treloar_Data.xlsx');
    x = D(:,1);
    y = D(:,2);

  
% Define Ogden Function
    ogden_funct1 = @(c) (c(1).*(x.^(c(4)-1)-x.^(-1/2*c(4)-1)) + ...
    c(2).*(x.^(c(5)-1)-x.^(-1/2*c(5)-1)) + ...
    c(3).*(x.^(c(6)-1)-x.^(-1/2*c(6)-1))) - y;


% Initial_Guess
    Initial_Guess1 = [1 1 1 1 1 1];
    lb =[];
    ub=[];


%Run multistart 

    problem1 = createOptimProblem('lsqnonlin', 'objective', ogden_funct1, 'x0', Initial_Guess1, 'lb', lb ,'ub', ub);
    ms = MultiStart; ms.XTolerance =1*10.^(-13); ms.FunctionTolerance =1*10.^(-13);
    [xmultinonlin,errormultinonlin] =run(ms,problem1, 1000)   

    
% Plot Our Data    
  ogden_plot_func =@(c) c(1).*(x.^(c(4)-1)-x.^((-1)/2*c(4)-1)) + ...
    c(2).*(x.^(c(5)-1)-x.^(-1/2*c(5)-1)) + ...
    c(3).*(x.^(c(6)-1)-x.^((-1)/2*c(6)-1)); 
plot(x, y, 'ko', x, ...
        ogden_plot_func(xmultinonlin), '-b')
    legend('Data', 'lsqnonlin-Unconstrained')
    title('Ogden Model Unconstrained')
    xlabel('Stretch')
    ylabel('Stress')

Answer 1

[Alan Weiss]

This type of problem is well-known to be numerically touchy, and to have multiple local minima. You really should give bounds, the tightest ones you can. Otherwise, MultiStart has way too large a space to search. And let MultiStart run for more iterations than I show (I gave 500 as a limit because of the limitations of online publishing, 50 second time limit).

x = [1.0000
    1.0100
    1.1200
    1.2400
    1.3900
    1.6100
    1.8900
    2.1700
    2.4200
    3.0100
    3.5800
    4.0300
    4.7600
    5.3600
    5.7600
    6.1600
    6.4000
    6.6200
    6.8700
    7.0500
    7.1600
    7.2700
    7.4300
    7.5000
    7.6100];

y = [         0
    0.0300
    0.1400
    0.2300
    0.3200
    0.4100
    0.5000
    0.5800
    0.6700
    0.8500
    1.0400
    1.2100
    1.5800
    1.9400
    2.2900
    2.6700
    3.0200
    3.3900
    3.7500
    4.1200
    4.4700
    4.8500
    5.2100
    5.5700
    6.3000];
  
% Define Ogden Function
    ogden_funct1 = @(c) (c(1).*(x.^(c(4)-1)-x.^(-1/2*c(4)-1)) + ...
    c(2).*(x.^(c(5)-1)-x.^(-1/2*c(5)-1)) + ...
    c(3).*(x.^(c(6)-1)-x.^(-1/2*c(6)-1))) - y;


% Initial_Guess
    Initial_Guess1 = [1 1 1 1 1 1];
    lb =-5*Initial_Guess1;
    ub=5*Initial_Guess1;

rng(1)
%Run multistart 

    problem1 = createOptimProblem('lsqnonlin', 'objective', ogden_funct1, 'x0', Initial_Guess1, 'lb', lb ,'ub', ub);
    ms = MultiStart; 
    % ms.XTolerance =1*10.^(-13); ms.FunctionTolerance =1*10.^(-13);
    [xmultinonlin,errormultinonlin] =run(ms,problem1, 500)   

MultiStart completed some of the runs from the start points.

135 out of 500 local solver runs converged with a positive local solver exit flag.

xmultinonlin = 1×6

0.9277 -0.2801 -0.9546 4.5049 -5.0000 4.4908

errormultinonlin = 0.3153

    
% Plot Our Data    
  ogden_plot_func =@(c) c(1).*(x.^(c(4)-1)-x.^((-1)/2*c(4)-1)) + ...
    c(2).*(x.^(c(5)-1)-x.^(-1/2*c(5)-1)) + ...
    c(3).*(x.^(c(6)-1)-x.^((-1)/2*c(6)-1)); 
plot(x, y, 'ko', x, ...
        ogden_plot_func(xmultinonlin), '-b')
    legend('Data', 'lsqnonlin+Bounds')
    title('Ogden Model With Bounds')
    xlabel('Stretch')
    ylabel('Stress')

Alan Weiss

MATLAB mathematical toolbox documentation

[Reed]

Hey thanks @Alan Weiss... that is what I was starting to think, just wanted some reassurance. I really would love to see this unconstrained but it seems like there are just too many points that steer lsqnonlin in the 'wrong direction'. Appreciate it.

Answer 2

[Torsten]

D = [1.00	0.00
    1.01	0.03
    1.12	0.14
    1.24	0.23
    1.39	0.32
    1.61	0.41
    1.89	0.50
    2.17	0.58
    2.42	0.67
    3.01	0.85
    3.58	1.04
    4.03	1.21
    4.76	1.58
    5.36	1.94
    5.76	2.29
    6.16	2.67
    6.40	3.02
    6.62	3.39
    6.87	3.75
    7.05	4.12
    7.16	4.47
    7.27	4.85
    7.43	5.21
    7.50	5.57
    7.61	6.30];
x = D(:,1);
y = D(:,2);


% Define Ogden Function
ogden_funct1 = @(c) (c(1).*(x.^(c(4)-1)-x.^(-1/2*c(4)-1)) + ...
    c(2).*(x.^(c(5)-1)-x.^(-1/2*c(5)-1)) + ...
    c(3).*(x.^(c(6)-1)-x.^(-1/2*c(6)-1))) - y;
% Initial_Guess
Initial_Guess1 = [-1.7720e+01 4.7842e-02 6.7178e-01 -1.7202e-01 3.5402e+00 -3.7330e+00];
lb =[];
ub=[];
options = optimset('MaxFunEvals',10000);
c = lsqnonlin(ogden_funct1,Initial_Guess1,lb,ub,options)

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

c = 1×6

-20.8780 0.0216 16.5974 -1.3422 3.8170 -1.6466

plot(x, y, 'ko', x, ...
    ogden_funct1(c)+y, '-b')

[Reed]

I appreciate the comment, but I'm more concerned on why my code is not working. You are using a totally different X0 and not using multistart. The intersting point of my code is that as soon as I impose a lb of say -5 for each parameter, the answer returned has an fval of .08. The issue is quite puzzling.

[Torsten]

The crucial point is to give good initial values for the parameters. Even MultiStart cannot ensure that some of the chosen parameter vectors yield convergence.

To choose c(4) - c(6) as 1, e.g., is a bad idea since important parts of your equation cancel out right at the beginning.

Further I would guess that your problem is overfitted with 6 parameters.

Restricting parameters to reasonable ranges often eases convergence. That's how I would interpret that setting a lower bound yields better results.

[Alex Sha]

The results given above by Torsten is actually a local solution, not the optimal one; The global solution is hard to be obtained, but should be like below:

Sum Squared Error (SSE): 0.0263802341297132
Root of Mean Square Error (RMSE): 0.0324839862884549
Correlation Coef. (R): 0.999860924905598
R-Square: 0.999721869153078

Parameter	Best Estimate       
---------	-------------       
c1       	-1.91181986182648E-6
c2       	1.55279012410112E-82
c3       	-0.240651398862925  
c4       	-16.1524712665223   
c5       	93.4941437734185    
c6       	-4.3004835501308 

or

Sum Squared Error (SSE): 0.0263802481671276
Root of Mean Square Error (RMSE): 0.0324839949311212
Correlation Coef. (R): 0.999860924854578
R-Square: 0.999721869051053

Parameter	Best Estimate       
---------	-------------       
c1       	-0.240651587389618  
c2       	1.55282658451632E-82
c3       	1.91182680758376E-6 
c4       	-4.30048241101075   
c5       	93.4941322424665    
c6       	8.07623397445497

[Alex Sha]

One more solution (different order of parameters):

Sum Squared Error (SSE): 0.0263802481671275
Root of Mean Square Error (RMSE): 0.0324839949311211
Correlation Coef. (R): 0.999860924854462
R-Square: 0.99972186905082

Parameter	Best Estimate       
---------	-------------       
c1       	-0.240651586280672  
c2       	1.91182678462031E-6 
c3       	1.55282423193542E-82
c4       	-4.30048242034866   
c5       	8.07623397801575    
c6       	93.4941329971174

[Reed]

Hello @Alex Sha, I always see you pose those 'global solutions' and never know where you are getting them from. I am assuming this is not in matlab? Where could I get access to that tool?

[Alex Sha]

Hi, Reed, it is ture that the results provided above are not from Matlab.

Matlab is a great and popular math software product, including me, widely used around the world. However, in some areas, for example, the global optimization algorithms or the tool of global optimization, the improvements are relatively minor and no longer have an advantage over some other products, this makes some of Matlab function or toolboxes (lsqcurvefit, gatool, fsolve, vpasolve, cftool...) to be difficult to get optimal results for some problems, for example, the curve fitting problem in this topic, even the result given by Alan Weiss is not correct.

Because of forum rules, it is not appropriate to discuss software other than matlab here, so only giving the right answer, and hope Mathworks will keep to improve Matlab more powerful.

[Reed]

@Alex Sha hey alex... I know we are not allowed to discuss other software products on here but I was wondering if you could email me wink wink. I have a question about your answer [e-mail]