Optimisation of parameters in coupled ODE

Hiya, am massively in over my head with some code I have been trying to create.

Currently, I have solved 9 coupled ODE using 4th-order Runge Kutta method however I am trying to fit these solutions to experimental data by changing the kinetic parameters. I was hoping to be able to use the GA method for optimisation due to the large number of variables making curvefit not very effective.

The form of the differential currently is dn/dt = K*n (where K is the matrix of kinetic parameters and is 9 x 9 and n is the number of obects and is 9 x 1).

My experimental data is not quite compiled yet but will be a set of 9 x 1 vectors (measured values for n) for varying times.

I hope this isn't too vague and thanks in advance (have attached files in case my explanation is poor)

[William Rose]

@James,

Do I understand correctly that you have 81 independent parameters which you want to estimate? You said "the large number of variables making curvefit not very effective". Do you mean you are using the Curve Fitter app? I have not used that app, so I cannot comment on it. I have used the fmincon() function and I have generally been very happy with it. For example, I have used fmincon() to adjust 9 parameters of a system of 15 ODEs, in order to fit experimental data.

Here is a intersting paper describing the use of a micro-genetic algorithm (described in the paper) to optimize a 99-parameter model, which is comparable to your 81 parameters - if that is what you have.

Good luck!

Answer 1

[Star Strider]

I’m not certain what you’re doing (and 81 parameters seems excessive), however this is the approach I use with ga to solve kinetics problems —

t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];

c=[0.902	0.06997	0.02463	0.00218
0.8072	0.1353	0.0482	0.008192
0.6757	0.2123	0.0864	0.0289
0.5569	0.2789	0.1063	0.06233
0.4297	0.3292	0.1476	0.09756
0.3774	0.3457	0.1485	0.1255
0.2149	0.3486	0.1821	0.2526
0.141	0.3254	0.194	0.3401
0.04921	0.2445	0.1742	0.5277
0.0178	0.1728	0.1732	0.6323
0.006431	0.1091	0.1137	0.7702
0.002595	0.08301	0.08224	0.835]; 
% theta0=[1;1;1;1;1;1];

% % [theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 50;
Parms = 10;

optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10);                                                                  % Options Structure For 'Answers' Problems
% opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
% optshf = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',5E3, 'FunctionTolerance',1E-10, 'HybridFcn',@fmincon, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);   % With 'HybridFcn'

t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)

Start Time: 2022-12-17 04:45:30.6173

[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);

Optimization terminated: average change in the fitness value less than options.FunctionTolerance.

t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)

Stop Time: 2022-12-17 04:46:10.3813

GA_Time = etime(t1,t0)

GA_Time = 39.7640

DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)

Elapsed Time:   3.976399800000000E+01		00:00:39.763

fprintf('Fitness value at convergence =\t%.4f\nGenerations =\t\t\t%d\n\n',fval,output.generations)

Fitness value at convergence =	0.2431
Generations =			926

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end

		Theta( 1) =  0.04910
		Theta( 2) =  1.03806
		Theta( 3) =  6.22608
		Theta( 4) = 12.82898
		Theta( 5) =  1.05286
		Theta( 6) =  0.32294
		Theta( 7) =  1.01858
		Theta( 8) =  0.08014
		Theta( 9) =  0.03686
		Theta(10) =  0.00009

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','N')

function C=kinetics(theta,t)
c0 = theta(7:10);
% c0=[1;0;0;0];
[T,Cv]=ode15s(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(4,1);
    dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
    dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
    dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
    dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
    dC=dcdt;
    end
C=Cv;
end

I posted this code previously for similar problems. The first six parameters are the kietic constants to be extimated, and the last four parameters to be estimated are the initial conditions for the differential equations. (I decreased the population size so it would run in the 55 seconds allotted here, so these results are not the best that ga can provide. Increase it — and if necessary the number of generations — to allow ga to converge on a satisfactory solution.)

Also, if the parameters have widely-ranging amplitudes (more than three orders-of-magnitude), use ode15s or one of the other stiff solvers for best results, otherwise it could take days for ga to converge.

It should be fairly straightforward to adapt it to your problem.

.