How to implement MultiStart to find a good fit of two curves at the same time with real data such as our curves are governed by an ODE system?

% Set up the problem for MultiStart
problem = createOptimProblem('lsqcurvefit','x0',k0,'objective',@simulatedhs,...
    'lb',zeros(size(k0)),'xdata',tforward ,'ydata',[Hdata,HSdata] );
ms = MultiStart;
[k,fval,Exitflag,Output,Solutions] = run(ms,problem,50);
simulated_data = simulatedhs(k,tforward);
X=simulated_data(:,1);
Z=simulated_data(:,2);
%plot the result
figure(5)
plot(tforward,Hdata ,tforward,X)
plot(tforward,Hdata,tforward,Y)
legend('Data','Fitted result')

[Torsten]

Every step is explained in detail here:

https://uk.mathworks.com/help/gads/multistart-using-lsqcurvefit-or-lsqnonlin.html

So what problem do you encounter ?

[Khadija]

I followed the instructions cited in the site that you recommended to me, the problem is that for 'ydata' there is a vector, that's why; the execution is not done for my problem and the MultiStart function no longer works!!

[Torsten]

There are problems with ode15s, but MultiStart seems to work in principle. Do you use an older MATLAB version ?

% 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:1:2019; 
%tmeasure = [ 1:100:1001]'; 
% initial values
gamma = 1.5;
phi_S = 0.0006;    % transmission prob
phi_H = 0.000051;    % trans proba
c1=3;
c2=1.5;
theta1 = 100;      % djustment parameters for syph
theta2 = 4;      %djustment parameters for hi
alpha = 0.6;      % progression rate
beta = 0.2;       %Complications rate
rho1 = 0.4;       % adjustment parameters
rho2 = 1.5;       %  adjustment parameters
rho3 = 1.5;       % adjustment parameters
k0 = [phi_S phi_H c1 c2 theta1 theta2 alpha beta rho1 rho2 rho3 ];
% solve equ with initial value of parameters
[t, Y] = ode15s(@(t, y)modelhs(t, y, k0), tforward, [ 5000.0 20.0 2.0  0.0  6.0 2.0  0.0 ],odeset('RelTol',1e-10,'AbsTol',1e-10));
yintM = Y(:,1);
yintS = Y(:,2);
yintH1 = Y(:,3);
yintH2 = Y(:,4);
yintH1S=Y(:,5);
yintH2S=Y(:,6);
Hh=phi_H * ( yintH1 + c1 *  yintH1S ).*yintM + rho1 * gamma * yintH1S + alpha* yintH2; % new cases from monoH
HShs=theta1*phi_S * ( yintS + c2 *  yintH1S ).*yintH1 + theta2*phi_H * ( yintH1 + c1 *  yintH1S ).*yintS+ rho2*alpha*yintH1S;   %new case of co h-s
%H2q = Y(:,4);% assignts the y-coordinates of .
% Plotting specific data and solutions
% Display the results
figure(1)
%subplot(1,2,1);
plot(tdata, Hdata, 'r*');
hold on
plot(tdata, Hh, 'b-');
xlabel('time in days');
ylabel('Number of monhiv cases');

%axis([2009 2019 0 500]);
figure(2)
%subplot(1,2,1);
plot(tdata, HSdata, 'r*');
hold on
plot(tdata, HShs, 'b-');
xlabel('time in days');
ylabel('Number of Coinfection cases');

%axis([2009 2019 0 500]);
% Set up the problem for MultiStart
problem = createOptimProblem('lsqcurvefit','x0',k0,'objective',@simulatedhs,...
                              'lb',zeros(size(k0)),'xdata',tforward ,'ydata',[Hdata,HSdata] );
ms = MultiStart;
[k,fval,Exitflag,Output,Solutions] = run(ms,problem,5);

Warning: Failure at t=2.009000e+03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.637979e-12) at time t.

Warning: Failure at t=2.009000e+03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.637979e-12) at time t.

Warning: Failure at t=2.009000e+03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.637979e-12) at time t.

Warning: Failure at t=2.009000e+03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.637979e-12) at time t.

MultiStart completed some of the runs from the start points. 

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

simulated_data = simulatedhs(k,tforward);
X=simulated_data(:,1);
Z=simulated_data(:,2);
%plot the result
figure(5)
plot(tforward,Hdata ,tforward,X)

plot(tforward,HSdata,tforward,Z)
legend('Data','Fitted result')

%kstart = k0.*(1+0.01*rand(size(k0)))   % Change k by a little amount
%[k,fval] = lsqcurvefit(@simulatedhs,kstart,tforward,[Hh,HShs],zeros(numel(k0),1),[],optimset('MaxFunEvals',10000,'MaxIter',10000));
%k(:)-k0(:)
%simulated_data = simulatedhs(k,tforward);
%H = simulated_data(:,1);
%HS = simulated_data(:,2);
%max(H-Hh)
%max(HS-HShs)
%figure(3)
%plot(tforward.',[H,Hh])
%figure(4)
%plot(tforward.',[HS,HShs])
function dy=modelhs(~,y,k)
delta = 0.01;  % Taux de mortalité
delta_S = 0.05;     % Taux de mort de Syphilis.
delta_H = 0.4;
Lambda =4.04 *100; 
gamma=1.5;
phi_S =k(1);
phi_H =k(2);
c1 = k(3);
c2=k(4) ;
theta1 =k(5) ;
theta2 = k(6);
alpha = k(7);
beta = k(8);
rho1 = k(9);
rho2= k(10);
rho3=k(11);
dy = zeros(7,1);
%lambda_s=phi_S * ( y(2) + c2 * y(5))
%lambda_H= phi_H * ( y(3) + c1 * y(5) )
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) * 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 ) * y(5) ;%H1S
dy(6)= rho2 * alpha * y(5)  - ( rho3 * beta + delta_S + delta_H + delta ) * y(6) ;%H2S
dy(7)= beta * y(4) + rho3 * beta * y(6) - ( delta + delta_H ) * y(7) ;%C
end

function simulated_data = simulatedhs(k,tdata)
% 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
    ];
gamma = 1.5;
phi_S =k(1);
phi_H =k(2);
c1 = k(3);
c2=k(4) ;
theta1 =k(5) ;
theta2 = k(6);
alpha = k(7);
beta = k(8);
rho1 = k(9);
rho2= k(10);
rho3=k(11);
[T, Y] = ode15s(@(t,y)modelhs(t,y,k),tdata,[ 5000.0 20.0 2.0  0.0  6.0 2.0  0.0 ],odeset('RelTol',1e-10,'AbsTol',1e-10));
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; % new cases from mono-HIV
HS=theta1*phi_S * ( S + c2 *  H1S ).*H1 + theta2*phi_H * ( H1 + c1 *  H1S ).*S+ rho2*alpha*H1S;   %new case of coinfection hiv+syphilis
simulated_data = [H,HS];
end

[Khadija]

when i'am running the code, this message appears:

MultiStart stopped without completing the runs from all start points.

1 out of 5 local solver runs exceeded the

iteration limit (problem.options.MaxIterations) or

the function evaluation limit (problem.options.MaxFunctionEvaluations).

None of the 5 local solver runs converged with a positive local solver exitflag.

>>

[Torsten]

This can happen. Either "ode15s" or "lsqcurvefit" didn't converge for the initial value vectors k chosen by "MultiStart".

You should control the parameter vectors being chosen by "MultiStart" so that they don't become nonsense:

https://uk.mathworks.com/help/gads/specify-multistart-start-points-problem.html

But as said: An ODE model with 11 free parameters and the only restriction that they are all positive is a challenge. And keep in mind that you only have 11 (ok: 22) data to fit 11 parameters - these are far too few.

[Khadija]

At the moment, I feel really lost!!

Is the problem due to the two data vectors (11*1,11*1)?

[Walter Roberson]

You should increase the iteration limit

opts = optimoptions('lsqcurvefit', 'MaxIterations', 1e5, 'MaxFunctionEvaluations', 1e5);
problem = createOptimProblem('lsqcurvefit','x0',k0,'objective',@simulatedhs,...
                              'lb',zeros(size(k0)),'xdata',tforward ,'ydata',[Hdata,HSdata], ...
                              'options', opts);

[Torsten]

Is the problem due to the two data vectors (11*1,11*1)?

No. The problem is due to 11 free parameters that can attain any positive value according to your settings, the fact that the simulated data come from a system of complicated ordinary differential equaltions and the fact that you supply far too few data for the fitting process. Isn't it possible to reduce the number of fitting parameters and to restrict them to certain lower and upper bounds via the arrays lb and ub ?

[Khadija]

I have eliminated three parameters that I will fix later. I have added the lower bound; but for the upper bound I cannot put it only to the three of these parameters, do I have to decompose my parameter vector into two vectors?

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

[Torsten]

do I have to decompose my parameter vector into two vectors?

No, you can use "Inf" in the list of upper bounds for the parameters where no upper bound is known.

[Khadija]

Is there any problem of this ligne?

[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@simulatedhs,k0,tforward,[Hdata,HSdata],...

'lb',[0 0 1 1 1 1 0 0 0 1 1],'ub',[1 1 inf inf inf inf 1 1 1 inf inf],...

optimset('MaxFunEvals',10000,'MaxIter',10000));

[Torsten]

lb = [0 0 1 1 1 1 0 0 0 1 1];
ub = [1 1 inf inf inf inf 1 1 1 inf inf];
options = optimset('MaxFunEvals',10000,'MaxIter',10000);
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@simulatedhs,k0,tforward,[Hdata,HSdata],lb,ub,options);

[Khadija]

for now, MultiStart does not give a better result, I have not been able to solve the problem. I keep the estimate without multistart which gives a result close to the suitable one. Thank you again for your patience!!