**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Parameter Estimation for a System of Differential Equations

365 views (last 30 days)

Show older comments

Hello. Ok, so I'm new to matlab and I've got a question regarding parameter estimation for a kinetic model. I have 4 different reactants and their concentrations are c1, c2, c3 and c4. I also have 4 differential equations, each one related to a concentration (c1, c2, c3 and c4, respectively -see below-) and experimental data for all these concentrations on 12 different times plus the initial condition. The k's are the rate coefficients. I want to solve this system of ODE's using ode45 and then use the output to compute the experimental data minus the observed data and use these results to estimate the values of k's using lsqnonlin, but apparently I can't solve these ODE's without numerical values for k -which is what I want to know-. Any help on how to set up the command to solve this?

function dcdt=batch(t,c,k)

dcdt=zeros(4,1);

dcdt(1)=-k(1)*c(1)-k(2)*c(1);

dcdt(2)= k(1)*c(1)+k(4)*c(3)-k(3)*c(2)-k(5)*c(2);

dcdt(3)= k(2)*c(1)+k(3)*c(2)-k(4)*c(3)+k(6)*c(4);

dcdt(4)= k(5)*c(2)-k(6)*c(4);

end

Data:

t c1 c2 c3 c4

0 1 0 0 0

0.1 0.902 0.06997 0.02463 0.00218

0.2 0.8072 0.1353 0.0482 0.008192

0.4 0.6757 0.2123 0.0864 0.0289

0.6 0.5569 0.2789 0.1063 0.06233

0.8 0.4297 0.3292 0.1476 0.09756

1 0.3774 0.3457 0.1485 0.1255

1.5 0.2149 0.3486 0.1821 0.2526

2 0.141 0.3254 0.194 0.3401

3 0.04921 0.2445 0.1742 0.5277

4 0.0178 0.1728 0.1732 0.6323

5 0.006431 0.1091 0.1137 0.7702

6 0.002595 0.08301 0.08224 0.835

Thanks in advance!

##### 7 Comments

Alex Sha
on 5 Sep 2019

the results below are much better:

k1 0.756713424357579

k2 0.246983152463589

k3 0.100637714178661

k4 0.24638788199029

k5 0.582535403727263

k6 -0.0269993611093679

Star Strider
on 5 Sep 2019

Torsten
on 5 Sep 2019

A negative reaction rate constant (k6) does not make sense.

You will have to constrain the ki to be non-negative.

Alex Sha
on 6 Sep 2019

if want to all parameters are positive， the results will be:

k1 0.757805051886583

k2 0.248942776988651

k3 0.182367452149268

k4 0.462511556999762

k5 0.620286121513934

k6 0

Star Strider has done this kind of job perfectly, I just used another package to obtain the result above, This result can be used to compare with Matlab results。

Star Strider
on 8 May 2020

Please do not post new problems to this thread unless they relate directly to it. Post them instead as new Questions.

Votes are always appreciated!

.

Budda
on 2 Jan 2021

Hi Star Srider, How would I change the code if only say the data for first variable ( equation) available in the Igor_Moura.m?

Thank you in advance

Star Strider
on 2 Jan 2021

Budda —

Make appropriate changes to the ‘kinetics’ function (specifically to the ‘C’ variable) to output only the variable you need. In that code, there are 4 columns, so address only the column you want returned in ‘C’.

### Accepted Answer

Star Strider
on 1 Dec 2016

##### 36 Comments

Igor Moura
on 1 Dec 2016

Hi, Star Rider, thanks for your help! Both problems are very alike, those techniques might be what I'm looking for, but I just can't get the code to work :( I created a kinetics.m file with the following code:

function C=kinetics(t,theta)

c0=[1;0;0;0];

[T,Cv]=ode45(@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(:,1);

end

and then I called it using:

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);

but I keep getting the following errors:

Error using odearguments (line 83)

The last entry in tspan must be different from the first entry.

%

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

%

Error in kinetics (line 3)

[T,Cv]=ode45(@DifEq,t,c0);

%

Error in lsqcurvefit (line 202)

initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

%

Error in cu (line 26)

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

%

Caused by:

Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

Star Strider
on 4 Dec 2016

I apologise for the delay. It’s been four years since I dealt with the Monod kinetics problem and I had to go back to remember what I’d done.

There were two problems:

- Your ‘kinetics’ function has to have the parameter variable (here ‘theta’) as the first argument, and ‘t’ as the second;
- You need to return all the values of your ‘Cv’ vector as ‘C’ since you are fitting all of them.

I made those changes and put all the corrected code, data, and plot calls in the attached ‘Igor_Moura.m’ file, so you can run it by typing Igor_Moura in your Command Window. I won’t post it here. Open it in your Editor to see the changes and my added code.

If you want to return ‘theta’, the concentrations, and other variables to your workspace, the easiest way is to add them as outputs to the ‘Igor_Moura’ function.

The Estimated Rate Constants:

Rate Constants:

Theta(1) = 0.89197

Theta(2) = 0.26069

Theta(3) = 0.24870

Theta(4) = 0.48183

Theta(5) = 0.60977

Theta(6) = -0.01301

The Plot:

Make appropriate changes in the legend code to describe the concentrations.

Khan Jan
on 3 Apr 2017

Star Strider
on 3 Apr 2017

‘How to plot them’

See the code in the ‘Igor_Moura.m’ file I attached to my previous Comment.

‘Is it possible to constrain the rate constant value between 0 and 0.80 using LSQCURVEFIT function, if yes how to do it?’

Yes. See the documentation for lsqcurvefit.

Taimur Akhtar
on 4 Apr 2017

Edited: Taimur Akhtar
on 4 Apr 2017

Khan Jan
on 4 Apr 2017

Taimur Akhtar
on 6 May 2017

I got plots after solving ODE'S. Now, I need to fit a power law line in the plot. Any idea.

Wouter Van Hecke
on 25 Aug 2017

Hi all,

Great post.

Do you have an idea on how to deal with incomplete data sets? Assume matrix c in the "Igor_Moura.m" file is incomplete and contains some unknown "NaN" values. In that case, lsqcurvefit gives the following error:

"Objective function is returning undefined values at initial point. lsqcurvefit cannot continue. " One option is to remove the entire row(s) with the unknown values, but aren't there more elegant approaches for this problem?

Many thanks for your suggestions!

Best regards, Wouter

Star Strider
on 25 Aug 2017

‘One option is to remove the entire row(s) with the unknown values, but aren't there more elegant approaches for this problem?’

Not that I’m aware of. The only other approach would be to use interpolation, and that has the extreme hazard of creating data that don’t actually exist, especially if the data are noisy or changing rapidly, making the parameter estimates essentially worthless. The Statistics and Machine Learning Toolbox functions ignore the rows with NaN values. I would follow that convention.

SHUBHAM PATEL
on 1 Oct 2019

Hello everyone, I had a very similar problem and used the above code by Star Strider for the fit. One of the problem i am facing is that, i am getting negative values for the parameters which for my specfic problem is not physical.

How do i constrain my parameters to give only postive values.

Thank you.

Star Strider
on 1 Oct 2019

Set the ‘lb’ (lower bound) argument in lsqcurvefit (or other optimization routines that allow bounded parameters) to a zeros vector the size of your parameter vector.

In the context of this code, that would be:

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));

↑

Gustavo Lunardon
on 24 Mar 2020

Edited: Gustavo Lunardon
on 24 Mar 2020

Could someone show me how to do the same minimization, (same data) but with fminsearch? I was trying to adapt the example for a while and tried to pass the variables and parameters in different ways but I keep getting :

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

My attempt is in attachment.

Star Strider
on 24 Mar 2020

Use this:

[theta,fval,exitflag] = fminsearch(@(theta)norm(c - kinetics(theta,t,c0)),theta0);

and:

Cfit = kinetics(theta, tv, c0);

I suggest optimising the initial conidtions as well, rather than using fixed initial conditions.

Gustavo Lunardon
on 24 Mar 2020

Thank you very much Star Strider!

I accompany many of your answers and they are always helpful, but answering this question in so little time exceeded my expectations!

Star Strider
on 24 Mar 2020

My pleasure!

Thank you!

I cannot always respond quickly, however I do my best!

Daphne Deidre Wong
on 27 Aug 2020

Hi Star Strider.

My the set of equations that i need to solve are very much similar to the file attached above 'Igor_Moura.m'. However I was wondering how should i implement the code if i want the reaction to occur over a range of temperature?

or would i need to put in the experimental results data at the certain temperature to get my theta value @ that temperature?

Kind regards, Daphne

Star Strider
on 27 Aug 2020

It depends.

If they are functions of temperature only and not of time, then use temperature as the independent variable and adjust the code accordingly.

If they are functions of time and temperature simultaneously, and both vary, that requires a completely different approach since the system is one of partial differential equations, not ordinary differential equations. The parameters would then not be scalars, but functions of temperature. That would be very difficult. One option would be to solve them with respect to time at each temperature (in a loop) as ordinary differential equations, and then later use appropriate parameter estimation techniques to characterise each parameter as a function of temperature. (That would likely be the approach I would use, at least at the outset, to determine if it was appropriate.)

Daniel Zhou
on 19 Nov 2020

Hi, Star.

I'm using the same code to fit data now. But I want to print this Rsdnrm value. How can I print that value?

Star Strider
on 19 Nov 2020

Star Strider
on 7 Dec 2020

ana ramirez de las heras’s Answer moved here —

Hi Star Strider,

Do I have to know c() values to use this method? I have a similar problem and I have data from the exercise statement but without knowing the values of theta() I don't know how to obtain c() values...

Thank you in advance!

Star Strider
on 7 Dec 2020

ana ramirez de las hera —

In order to estimate the parameters, the ‘c’ data (including the time vector associated with it) have to exist in your workspace. The code in my Answer estimates the kinetic parameters ‘theta’ vector from the data, so knowing them in advance is not necessary, however approximate estimates of their amplitudes are necessary, since the parameter estimation routines require an initial estimate of the parameters.

If you already have the parameters and are simply simulating the system, the ‘c’ values are not necessary.

If this Comment does not address your problems, please post a new Question with more detail as to what you are starting with and the result you want.

ana ramirez de las heras
on 7 Dec 2020

Yasin Islam
on 18 Jan 2021

Hi Star Strider,

thank you very much for the usefull solution to the problem! I adapted the code to fit for my 6 diff. eq. and 12 k-values. it works.

unfortunately i not only have unknown k-values but also unknown concentrations.

The problem in my case extends up to 18 concentrations with 18 differential equations (i do know time dependend concentration values for 6 of them)

How can I update the code to not only look for my k-values but also for the unknown concentrations in order to represent the experimental concentrations best? Although alot of the unknown concentrations will be equal to 0 and/or neglectable.

Star Strider
on 18 Jan 2021

Yasin Islam — You can obviously fit only the concentrations you have data for. If the other 12 concentrations (actually, differential equations describing them) can be expressed in terms of the parameters (‘k values’) you can estimate, simply simulate them outside the ‘kinetics’ function using the parameters and plot all of them.

If the other 12 differential equations involve parameters that you cannot estimate (or cannot derive from the parameters you can estimate), then you cannot simulate them. There is no way around that unfortunate situation.

si heon Seong
on 2 Feb 2021

Hi Star Strider,

first, thank you for introducing this page! It was great help to me

One thing I have a additional question about this example

What should I do to give each parameter the constraint I want? For instance, thetas must be bigger than 0 etc..

Thank you

Star Strider
on 2 Feb 2021

si heon Seong — The lsqcurvefit and many other optimisation functions allow parameter constraints, at least setting the upper and lower parameter bounds. Use those arguments.

Here, that would be:

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));

with the zeros call setting the lower limit. (The upper limit is not set here, since it is not necessary.)

si heon Seong
on 4 Feb 2021

Dear @Star Strider

I really appreciate your help! Can I ask the last one?

I have a model of 6 coupled differential equations for c1, c2, c3, c4, c5 and c6

but there are only two data set for c1 and c2

Can I still do the fitting job?

Always thank you

Star Strider
on 4 Feb 2021

si heon Seong — You can only fit ‘c1’ and ‘c2’ to your data, so in ‘kinetics’:

C = Cv(:,1:2);

(assuming that ‘c1’ and ‘c2’ are the first 2 differential equations in your system).

There have to be more data (more rows in the data matrix and time vector) than there are parameters to estimate.

My pleasure!

Khoi Ly
on 21 Feb 2021

Edited: Khoi Ly
on 21 Feb 2021

I was directed by your answer in another post regarding a similar problem.

In that post, I was trying to estimate the parameters of the differential equations using more than one time series data sets.

In the example above, given.

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];

we can estimate the parameters theta. Since the example above uses only one data set, what should I do if I have multiple data sets? Can I just concatenate the time-series data sets together and run the lsqcurvefit as usual?

Also, in my problem, each different time-series data set contains a different initial value. How do I modify the example to allow for different initial values depending on the chosen time-series data set?

Thank you

sumit kumar
on 28 Jun 2022

Edited: sumit kumar
on 28 Jun 2022

Star Strider
on 22 Sep 2022

UPDATE

Slightly revised code with some aethetic tweaks so that the lines and markers colours now match —

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 = rand(10,1);

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));

Local minimum found.
Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

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.76483
Theta( 2) = 0.23492
Theta( 3) = 0.20878
Theta( 4) = 0.49179
Theta( 5) = 0.62211
Theta( 6) = 0.00000
Theta( 7) = 0.90288
Theta( 8) = 0.07146
Theta( 9) = 0.02840
Theta(10) = 0.00000

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','best')

function C=kinetics(theta,t)

c0=theta(end-3:end);

[T,Cv]=ode45(@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

.

uki71319
on 24 Oct 2022

Edited: uki71319
on 24 Oct 2022

Dear @Star Strider, thanks for your detailed explanations. I've followed your codes, but i still got some problmes. Thus, i posted a new question under this link and hope to find out the solution for my lsqcurvefit problem. Would you mind taking a look and see where the issues are?

Any advice and help would be greatly appreciated. Thanks in advance!

Star Strider
on 24 Oct 2022

I looked at it just now, however it has already been answered.

Note that initial parameter estimates of zeros is never a good approach. Choose something nonzero and preferably not the same values for both parameters.

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)