Help needed on ODE solving coupled with fmincon

Hi, I need help to solve the problem with the code for my ODE solving. As I have many parameters and I have used lsqcurvefit to solve my problem, but there was an error message of:

"Solver stopped prematurely.

lsqcurvefit stopped because it exceeded the function evaluation limit,

options.MaxFunctionEvaluations = 800 (the default value)."

function ODE
% 2019 12 4

    function X=kinetics(k,t)
x0=[7.525;13.9;3;0;0;0];
[T,Xv]=ode45(@DifEq,t,x0);
%
    function dX=DifEq(t,x)
    dxdt=zeros(6,1);
    dxdt(1)= -(k(1)+k(6)).*x(1);
    dxdt(2)= -(k(2)+k(7)).*x(2);
    dxdt(3)= -k(8).*x(3);
    dxdt(4)= k(1).*x(1)-k(3).*x(4)-k(4).*x(4);
    dxdt(5)= k(2).*x(2)+k(3).*x(4)-k(5).*x(5);
    dxdt(6)= k(6).*x(1)+k(7).*x(2)+k(8).*x(3)+k(4).*x(4)+k(5).*x(5);
    dX=dxdt;
    end
X=Xv;
end

t=[0
5
10
15
20];

x=[7.525 13.9 3 0 0 0
5.585 9.51 8 1.44 2.89 2.0
5.275 8.27 7 1.65 3.88 2.35
4.925 7.93 6 1.90 3.97 2.7
3.885 6.00 5.5 2.89 5.70 2.95]; 
k0=[0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5];

[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,t,x);

fprintf(1,'\tRate Constants:\n')
for n1 = 1:length(k)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', n1, k(n1))
end

tv = linspace(min(t), max(t));
Xfit = kinetics(k, tv);

figure(1)
plot(t, x, 'p')
hold on
hlp = plot(tv, Xfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'X_1(t)', 'X_2(t)', 'X_3(t)', 'X_4(t)', 'X_5(t)', 'X_6(t)', 'Location','N')

end

Instead, I followed advice of another expert to use fmincon to substitute lsqcurvefit to evaluate the function to find k

k0=[0.25;0.25;0.25;0.25;0.25;0.25;0.25;0.25];
lb=[0;0;0;0;0;0;0;0];
ub=[0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5];
A=[];
b=[];
Aeq=[];
beq=[];

k = fmincon(@kinetics,k0,A,b,Aeq,beq,lb,ub);

But there was also an error message of:

"Not enough input arguments.

Error in kinetics

[T,Xv]=ode45(@DifEq,t,x0);

Error in fmincon

initVals.f = feval(funfcn{3},X,varargin{:});

Error in call_fmincon

k = fmincon(@kinetics,k0,A,b,Aeq,beq,lb,ub);

Caused by:

Failure in initial objective function evaluation. FMINCON cannot continue."

Answer 1

[Star Strider]

It is best to fit the initial conditions as well as the kinetics parameters. I would not put an upper bound on the kinetics parameters. A lower bound to be certain they are all is all that is necessary.

This is something of an improvement in my original code (I upgraded it later, including a change so the fitted lines and markers are the same colours):

    function X=kinetics(k,t)
% x0=[7.525;13.9;3;0;0;0];
x0 = k(9:14);                                                           % Fit The Initial Conditions As Well
[T,Xv]=ode45(@DifEq,t,x0);
%
    function dX=DifEq(t,x)
    dxdt=zeros(6,1);
    dxdt(1)= -(k(1)+k(6)).*x(1);
    dxdt(2)= -(k(2)+k(7)).*x(2);
    dxdt(3)= -k(8).*x(3);
    dxdt(4)= k(1).*x(1)-k(3).*x(4)-k(4).*x(4);
    dxdt(5)= k(2).*x(2)+k(3).*x(4)-k(5).*x(5);
    dxdt(6)= k(6).*x(1)+k(7).*x(2)+k(8).*x(3)+k(4).*x(4)+k(5).*x(5);
    dX=dxdt;
    end
X=Xv;
end
t=[0
5
10
15
20];
x=[7.525 13.9 3 0 0 0
5.585 9.51 8 1.44 2.89 2.0
5.275 8.27 7 1.65 3.88 2.35
4.925 7.93 6 1.90 3.97 2.7
3.885 6.00 5.5 2.89 5.70 2.95]; 
k0=[0.5;0.5;0.5;0.5;0.5;0.5;0.5;0.5;rand(6,1)];
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,t,x,zeros(size(k0)));
fprintf(1,'\tRate Constants:\n')
for n1 = 1:length(k)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', n1, k(n1))
end
tv = linspace(min(t), max(t));
Xfit = kinetics(k, tv);
figure(1)

hd = plot(t, x, 'p');
for k1 = 1:size(x,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Xfit);
for k1 = 1:size(x,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'X_1(t)', 'X_2(t)', 'X_3(t)', 'X_4(t)', 'X_5(t)', 'X_6(t)', 'Location','N')

This takes a while to run, however it produces a better result.

If you have the Global Optimization Toolbox, I also wrote a version of this code using the ga (genetic algorithm) function that is generally much more robust at finding the optimal parameter set, although it can take a long time to run (as ga optiomisations usually do). I will post it as an edit to my Answer here if you request it.

[Bernard Chon Han Ho]

I do have the Global Optimization Toolbox. Could you show me the coding as well as I am interested to find a more optimized parameter set? Thank you for your help

[Star Strider]

Sure!

As always, my pleasure!

Here’s the complete revised version of that code. It also includes an 'OuptutFcn' argument that saves the best individuals from each generation in the ‘SaveBest.mat’ file, so if anything crashes or you want to interrupt the optimisation to go do something else for a while, you can use the ‘SaveBest.mat’ file contents as part of the 'InitialPopulationMatrix' to reconstitute the best previous individuals as part of the new optimisation. (I run all of this inside a function, so put everything here inside your ‘ODE’ function wrapper and it should run without problems.)

The Code —

function ODE

    function X=kinetics(k,t)
% x0=[7.525;13.9;3;0;0;0];
x0 = k(9:14);                                                           % Fit The Initial Conditions As Well
[T,Xv]=ode45(@DifEq,t,x0);
%
    function dX=DifEq(t,x)
    dxdt=zeros(6,1);
    dxdt(1)= -(k(1)+k(6)).*x(1);
    dxdt(2)= -(k(2)+k(7)).*x(2);
    dxdt(3)= -k(8).*x(3);
    dxdt(4)= k(1).*x(1)-k(3).*x(4)-k(4).*x(4);
    dxdt(5)= k(2).*x(2)+k(3).*x(4)-k(5).*x(5);
    dxdt(6)= k(6).*x(1)+k(7).*x(2)+k(8).*x(3)+k(4).*x(4)+k(5).*x(5);
    dX=dxdt;
    end
X=Xv;
end
t=[0
5
10
15
20];
x=[7.525 13.9 3 0 0 0
5.585 9.51 8 1.44 2.89 2.0
5.275 8.27 7 1.65 3.88 2.35
4.925 7.93 6 1.90 3.97 2.7
3.885 6.00 5.5 2.89 5.70 2.95]; 


ftns = @(k) norm(x-kinetics(k,t));
PopSz = 500;
Parms = 14;

opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1, 'OutputFcn',@SaveOut);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
% [B,fval,exitflag,output,population,scores] = ga(ftns, 53, [],[],[],[],zeros(53,1),Inf(53,1),[],[],opts)
[theta,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
fprintf('\nElapsed Time: %23.15E\t\t%02d:%02d:%02d.%03d\n', GA_Time, hmsdv(GA_Time))

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

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

figure
hd = plot(t, x, 'p');
for k1 = 1:size(x,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Xfit);
for k1 = 1:size(x,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'X_1(t)', 'X_2(t)', 'X_3(t)', 'X_4(t)', 'X_5(t)', 'X_6(t)', 'Location','N')

    function [state,options,changed,str] = SaveOut(options,state,flag)
        file_name = 'SaveBest.mat';                                     % Name File
        if strcmp(flag,'init')
            bestind = state.Population;
            save(file_name, 'bestind')                                  % Write ‘Best Individual’ To File
        elseif strcmp(flag,'iter')
            ibest = state.Best(end);
            ibest = find(state.Score == ibest,1,'last');
            bestx = state.Population(ibest,:);
            previous = load('SaveBest.mat');
            bestind = [previous.bestind; bestx];                        % Read Previous Results, Append New Value
            save(file_name, 'bestind')                                  % Write ‘Best Individual’ To File
        end
        changed = true;                                                 % Necessary For Code, Use Appropriate Value
    end


format long
BestInGen = load('SaveBest.mat')
LastFive = BestInGen.bestind(end-4:end,:)

end

That ran without error when I tested it just now.

The best fitness value in the three runs I made to test this (9.975) produced these parameters:

	Rate Constants:
		Theta(1) = 20.11896
		Theta(2) =  0.00673
		Theta(3) =  0.10235
		Theta(4) =  0.01484
		Theta(5) =  0.01843
		Theta(6) =  0.34963
		Theta(7) =  0.01023
		Theta(8) =  0.00242
		Theta(9) =  6.87681
		Theta(10) =  9.96659
		Theta(11) =  6.79948
		Theta(12) =  0.00019
		Theta(13) =  0.00213
		Theta(14) =  0.00519

I do not know how well they fit because I was still adapting my code to your function (to be certain the code runs) and the plot failed. Everything now works correctly.

[Bernard Chon Han Ho]

May I know what is this error in your code for general algorithm?

"Undefined function or variable 'hmsdv'.

Error in generic_algorithm_nanno (line 44)

fprintf('\nElapsed Time: %23.15E\t\t%02d:%02d:%02d.%03d\n', GA_Time, hmsdv(GA_Time))"

And how do I check whether the theta produced fits my initial condition since you also fit my initial condition as well?

"% x0=[7.525;13.9;3;0;0;0];

x0 = k(9:14); % Fit The Initial Conditions As Well"

And also, there is no plot of x against t although the coding had ask for a plot of x against t.

[Star Strider]

The ‘hmsdv’ call is to a utility function I wrote more than 20 years ago to convert the output of etime into hours, minutes, and seconds. I forgot to remove it from the code I posted.

Change that line to:

fprintf('\nElapsed Time: %23.15E seconds\n', GA_Time)

and the problem wil disappear.

[Bernard Chon Han Ho]

Hi, I am so sorry but the coding seems to always change its answer. Is it that the data I have provided had some problem or the code always runs differently at every run?

[Star Strider]

It could very easily return different results in every run. In most instances, the results should not change much between runs that converge on an acceptable fit. The ga function is more likely to converge than lsqcurvefit for this sort of problem.

However you are estimating 8 parameters with 5 data points. (There should be at least 2 more data than parameters.) You also need to be certain that the differential equations are an appropriate model of the system you are investigating. (I do not know the details, so I cannot review the model and its coding to be certain of that.)

So to fit the data more accurately, it would be necessary either to get more data, or use an appropriate model with fewer parameters. (I experimented with interpolating (interp1) to create more data using several interpolation methods. However the fit did not improve, leading me to believe that the model as coded may not be appropriate.)

I will help as much as I can. However this is the best I can do with the information I have.

[Bernard Chon Han Ho]

Thank you so much for your help. The Genetic Algorithm works. But I have a question on what is the meaning of this part of your code.

ftns = @(k) norm(x-kinetics(k,t));

Are you minimizing the objective function using norm? Why is it x - kinetic?

[Star Strider]

As always, my pleasure!

This function:

ftns = @(k) norm(x-kinetics(k,t));

is the fitness function that the ga function uses to assess the parameters (‘k’) for each individual. It calculates the norm of the difference between the data (‘x’) and the evaluated objective function (‘kinetics’) with that particular parameter vector.

[Bernard Chon Han Ho]

HI, can I ask any particular reason for the population size PopSz = 500 and the 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3 ?

[Star Strider]

Put simply, when I run kinetics optimisations, that sort of initial population matrix works for me better than the default option. Choose the options that best fit your requriements.