how to get global optimum with Multistart?

I used lsqcurvfit to give an estimate to the parameters, but certainly the initial condition does not give a good fit, I implemented Multistart, and here is the error that appears:

Error using lsqcurvefit (line 289)

Function value and YDATA sizes are not equal.

Error in globaloptim.multistart.fmultistart

Error in MultiStart/run (line 257)

globaloptim.multistart.fmultistart(problem,startPointSets,msoptions);

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Error in rectifModel (line 51)

[k,fval,Exitflag,Output,Solutions] = run(ms,problem,5);

^^^^^^^^^^^^^^^^^

Caused by:

Failure in call to the local solver with supplied problem structure.

>>

[Star Strider]

I do not know what you are doing, so I cannot comment with respect to your code.

This can occur when a differential equation solver encounters a singularity and stops the integration before finishing it, so the data and model vectors do not have the same number of elemenets, or if the model results are the transpose of the data to be fitted, even though the number of elements in the two vectors (or matrices) are the same.

[Khadija]

Hi Mr star,

Here is the code I am trying to develop to have a good fit, the implementation of MultiStart will allow me to have a good result for the most part, even if I change the data or the initial condition of my system, please take a look!!

% Données spécifiques
specific_data = [
    2009 2 8;
    2010   10 22;
    2011   30 45;
    2012    111 75;
    2013 125 96;
    2014   255 192;
    2015  379 227;
    2016  384 238
    2017   360 279;
    2018   399 229;
    2019   235 128
    ];
% Utilisez les données spécifiques
tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
HSdata = specific_data(:, 3);
tforward = 2009:0.1:2019; 
%tmeasure = [ 1:100:1001]'; 
% initial values
gamma = 5.642;
alpha = 0.08;% progression rate
%fitted parameters
phi_S =0.0000034;    
phi_H = 0.000000028;    
c1=3;
c2=1.5;
theta1 =19;      
theta2 =10;      
beta = 0.6;      
rho1 =0.4;       
rho2 =1.2;      
rho3 = 2;       
rho4 =0.3;
rho5 =1.2;
k0 = [phi_S phi_H c1 c2 theta1 theta2  beta rho1 rho2 rho3 rho4 rho5 ];

%kstart = k0.*(1+0.01*rand(size(k0)));   % Change k by a little amount
% Set up the problem for MultiStart
problem = createOptimProblem('lsqcurvefit','x0',k0,'objective',@simulatedhs,...
    'lb',[0 0 1 1 1 1 0 0 1 1 1 1],'ub',[1 1 inf inf inf inf 1 1 inf inf inf inf],'xdata',tforward ,'ydata',[Hdata,HSdata] );
ms = MultiStart;
[k,fval,Exitflag,Output,Solutions] = run(ms,problem,5);
display(k);
simulated_data = simulatedhs(k,tforward);
H = simulated_data(:,1);
HS = simulated_data(:,2);
figure(1)
plot (tdata.', HSdata,'r*')
hold on
plot(tforward, HS, 'k--','linewidth',2)
legend('H-S coinfection Data', 'Model with fit');
axis([2009 2019 0 1000]);
figure(2)
plot(tforward, H, 'b--','linewidth',2)
hold on 
plot (tdata.', Hdata,'ro')
legend('Mono-Hiv Data', 'Model with fit');
axis([2009 2019 0 1000]);
function dy=modelhs(~,y,k)
gamma = 5.642;
alpha = 0.08; 
delta = 0.0056773;  
delta_S = 4.7889*10^(-8);     
delta_H = 9.538*10^(-6);   
Lambda =111130.48; 
phi_S =k(1);
phi_H =k(2);
c1 = k(3);
c2=k(4) ;
theta1 =k(5) ;
theta2 = k(6);
beta = k(7);
rho1 = k(8);
rho2= k(9);
rho3=k(10);
rho4=k(11);
rho5=k(12);
dy = zeros(7,1);
dy(1) = Lambda + gamma * y(2) - (phi_S * ( y(2) + c2 * y(5)) + phi_H * (y(3) + c1 * y(5)) + delta ) * y(1) ;%M
dy(2)= phi_S * ( y(2) + c2  *  y(5) ) * y(1) - ( gamma + theta2 * phi_H * ( y(3) + c1 * y(5) ) + delta_S ) * y(2) ;%S
dy(3) = phi_H * ( y(3) + c1 * y(5) ) * y(1) + rho1 * gamma * y(5) - (theta1 * phi_S * ( y(2) + c2 * y(5)) + delta + delta_H + alpha+rho4*beta) * y(3) ;%H1
dy(4) = alpha * y(3) - (beta  + delta + delta_H) * y(4) ;%H2
dy(5) = theta2 * phi_H * ( y(3) + c1 * y(5) ) * y(2) + ( theta1 * phi_S * ( y(2) + c2 * y(5) ) ) * y(3) - ( rho1 * gamma + rho2 * alpha + delta_H + delta_S + delta + rho5*beta) * y(5) ;%H1S
dy(6)= rho2 * alpha * y(5)  - ( rho3 * beta + rho1 * gamma + delta_S + delta_H + delta ) * y(6) ;%H2S
dy(7)= rho4 * beta * y(3)+rho5 * beta * y(5)+beta * y(4) + rho3 * beta * y(6) - ( delta + delta_H ) * y(7) ;%C
end
function simulated_data = simulatedhs(k,tdata)
gamma = 5.642;
alpha = 0.08; 
delta = 0.0056773; 
delta_S = 4.7889*10^(-8);     
delta_H = 9.538*10^(-6);  
Lambda =111130.48;
phi_S =k(1);
phi_H =k(2);
c1 = k(3);
c2=k(4) ;
theta1 =k(5) ;
theta2 = k(6);
beta = k(7);
rho1 = k(8);
rho2= k(9);
rho3=k(10);
rho4=k(11);
rho5=k(12);
[T, Y] = ode15s(@(t,y)modelhs(t,y,k),tdata,[ 12280000.0 100.0 2.0  0.0  6.0 2.0  0.0 ],odeset('RelTol',1e-10,'AbsTol',1e-10));
%[T, Y] = ode15s(@(t,y)modelhs(t,y,k),tdata,[ 12280000.0 100.0 2.0  0.0  6.0 2.0  0.0 ]);
M = Y(:,1);
S = Y(:,2);
H1 = Y(:,3);
H2 = Y(:,4);
H1S=Y(:,5);
H2S=Y(:,6);
H=phi_H * ( H1 + c1 *  H1S ).*M + rho1 * gamma * H1S + alpha.* H2+rho4*beta.* H1;
HS=theta1*phi_S .* ( S + c2 .*  H1S ).*H1 + theta2*phi_H .* ( H1 + c1 .*  H1S ).*S+ rho2*alpha.*H1S+rho5*beta.*H1S;   
simulated_data = [H,HS];
end

Answer 1

[Walter Roberson]

tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
HSdata = specific_data(:, 3);

tdata is an 11 x 1 column vector formed from the first column of the input. Similarly for HData and HSdata

tforward = 2009:0.1:2019; 

That is a 1 x 101 row vector

problem = createOptimProblem('lsqcurvefit','x0',k0,'objective',@simulatedhs,...
    'lb',[0 0 1 1 1 1 0 0 1 1 1 1],'ub',[1 1 inf inf inf inf 1 1 inf inf inf inf],'xdata',tforward ,'ydata',[Hdata,HSdata] );

So the lsqcurvefit is to pass in a 1 x 101 vector (tforward) as the xdata. And according to the problem creation, it expects back ydata which is [11 x 1, 11x1] for an overall solution of 11 x 2.

function simulated_data = simulatedhs(k,tdata)

The tdata received here will be the same as xdata -- a 1 x 101 vector.

[T, Y] = ode15s(@(t,y)modelhs(t,y,k),tdata,[ 12280000.0 100.0 2.0  0.0  6.0 2.0  0.0 ],odeset('RelTol',1e-10,'AbsTol',1e-10));

The tdata gets passed to ode15s as time data. It is a vector with length greater than 2, so the output of ode15s will normally have length(tdata) time entries (unless ode15s aborts early.)

M = Y(:,1);
S = Y(:,2);
H1 = Y(:,3);
H2 = Y(:,4);
H1S=Y(:,5);
H2S=Y(:,6);

Those will each be length(tdata) x 1.

H=phi_H * ( H1 + c1 *  H1S ).*M + rho1 * gamma * H1S + alpha.* H2+rho4*beta.* H1;
HS=theta1*phi_S .* ( S + c2 .*  H1S ).*H1 + theta2*phi_H .* ( H1 + c1 .*  H1S ).*S+ rho2*alpha.*H1S+rho5*beta.*H1S;   

H and HS are calculated in terms of length(tdata) x 1 vectors, so they will each be length(tdata) x 1 results.

simulated_data = [H,HS];

so simulated_data will be length(tdata) x 2.

But this does not match ydata, which is 11 x 2.

The fix is to use 'xdata', tdata and get rid of tforward

[Khadija]

Hi,

I changed forward to tdata, the problem is solved. but the lsqcurvefit function with Multistart does not give a good fit!!

[Khadija]

I think that MultiStart doesn't work in my program, is it because its implementation is not correct?

[Torsten]

If for some parameter vector during the MultiStart procedure, the ODE integrator gives up before reaching "tend", the size of "simulated_data" is too short because outputs are missing compared to the size of "[Hdata,HSdata]".

As a consequence, "lsqcurvefit" will complain that "Function value and YDATA sizes are not equal.".

Maybe it is this problem that you encounter - we don't know what you else mean by "I think that MultiStart doesn't work in my program".

[Star Strider]

@Torsten — Exactly what I mentioned in my earlier Comment.

[Torsten]

@Star Strider

Yes, I read your comment, but since @Khadija didn't go into it, I thought it might help to repeat it here in the "active" part of the discussion.

[Khadija]

Hi,

If I understand correctly, the problem is in my system?

[Torsten]

Are we still talking about the error message that you posted in your first question

Error using lsqcurvefit (line 289)

Function value and YDATA sizes are not equal.

Error in globaloptim.multistart.fmultistart

Error in MultiStart/run (line 257)

globaloptim.multistart.fmultistart(problem,startPointSets,msoptions);

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Error in rectifModel (line 51)

[k,fval,Exitflag,Output,Solutions] = run(ms,problem,5);

^^^^^^^^^^^^^^^^^

Caused by:

Failure in call to the local solver with supplied problem structure.

when you now write

"I think that MultiStart doesn't work in my program"

?

Or do you encounter a different problem now ?

[Star Strider]

@Torsten — I tried to help @Khadija with this problem in early July 2024 and ran iinto the same problem. (I attempted to use ga with it, because I thought ga might be able to discover the best parameters that would also keep the system from crashing, since ga is forgiving and resiliant in my experience, and has successfully done that in other situations.)

Obviously, the problem is iin the system of differential equations. However since we do not have the system being modeled, we cannot independently derive the equations or check the existing version for coding errors.

.

[Khadija]

Hi,

The system is analytically verified before being implemented. Could this behavior be due to the choice of the initial condition?

[Torsten]

The system is analytically verified before being implemented.

What does that mean ?

Look: you have 11 data points and you want to estimate 12 model parameters. This is simply not reasonable.

[Walter Roberson]

After the line

[T, Y] = ode15s(@(t,y)modelhs(t,y,k),tdata,[ 12280000.0 100.0 2.0  0.0  6.0 2.0  0.0 ],odeset('RelTol',1e-10,'AbsTol',1e-10));

add in a line

assert(numel(T) == numel(tdata), sprintf('ode15s returned %d values instead of %d', numel(T), numel(tdata)));

and run the code again.

[Khadija]

Hi,

I added the command you mentioned and nothing is displayed, I hope it's a good sign.

The existence of the solutions and their continuities is verified theoretically; the resolution process by discretizing, using ode15s solvers or other, escapes me. such problems I still get stuck!!

i modified the initial condition, and here is the message that is displayed:

" MultiStart completed some of the runs from the start points.

1 out of 5 local solver runs converged with a positive local solver exitflag."

for your remark Mr Torsten, of parameter reduction, if I do it, will the problem be solved?

[Walter Roberson]

I would suggest,

opts = optimoptions("lsqcurvefit",MaxFunctionEvaluations=1e5,MaxIterations=1e5);

problem = createOptimProblem('lsqcurvefit','x0',k0,'objective',@simulatedhs,...
    'lb',[0 0 1 1 1 1 0 0 1 1 1 1],'ub',[1 1 inf inf inf inf 1 1 inf inf inf inf],'xdata',tdata ,'ydata',[Hdata,HSdata], 'options', opts );

[Khadija]

Hi,

I added the "opts" option to the program, but no improvement in the result!!

[Khadija]

Hi,

i hope you are well, i am still struggling!!

i am still struggling to find a good fit.

i followed your advice Mr. Torsten, i reduced the number of parameters. As a result, my curve is close to the experimental data than before. I am looking for a good improvement of the fit. for a good visualization of the fit.

To minimize "fval", can genetic algorithm be a good approach to solve my problem?

[Walter Roberson]

You can always try genetic algorithm.

My personal experience is that ga() does not work especially well on most curve fitting tasks