6
$\begingroup$

I have the following data from chemical kinetics research to fit the parameters of ordinary differential equations:

$$ \left[ \begin{array}{ccccccc} \text{No.}& t & y_1(t)&y_2(t) & y_3(t) & y_4(t) & y_5(t)\\ 1&30.0000 & 9.1300 & 0.0931 & 0.0899 & 0.1000 & 0.0000 \\ 2&60.0000 & 8.9300 & 0.1270 & 0.1230 & 0.2270 & 0.0049 \\ 3&90.0000 & 8.6000 & 0.1510 & 0.1390 & 0.4920 & 0.0153 \\ 4&120.0000 & 8.2800 & 0.1540 & 0.1490 & 0.7780 & 0.0249 \\ 5&150.0000 & 7.9700 & 0.1540 & 0.1570 & 1.0700 & 0.0329 \\ 6&180.0000 & 7.8600 & 0.1540 & 0.1600 & 1.1700 & 0.0348 \\ 7&210.0000 & 7.8100 & 0.1530 & 0.1530 & 1.2100 & 0.0404 \\ 8&240.0000 & 7.7700 & 0.1400 & 0.1420 & 1.2800 & 0.0432 \\ \end{array} \right] $$

The ordinary differential equations to fit have $k_1,k_2,k_3,k_4,k_5,k_6$ to be determined. $$ \left\{ \begin{array}{l} {y_1}'(t)=-{k_1} {y_1}(t)-{k_2} {y_1}(t),\\ {y_2}'(t)={k_2} {y_1}(t)-{k_3} {y_2}(t),\\ {y_3}'(t)={k_1} {y_1}(t)+{k_3} {y_2}(t)-{k_4} {y_3}(t),\\ {y_4}'(t)={k_4} {y_3}(t)-{k_5} {y_2}(t) {y_4}(t)+{k_6} {y_5}(t),\\ {y_5}'(t)={k_5} {y_2}(t) {y_4}(t)-{k_6} {y_5}(t)\\ \end{array} \right.$$

In order to solve it from conventional numerical optimization methods, my original thoughts are: first convert it into least square problems, then apply numerical optimization to it, but this requires symbolically solve a nonlinear system of ordinary differential equations into explicit solutions first, which seems difficult.

My questions are:

(1)Is it possible to determine the global (least square or similarly converted) solution to all the parameters $k_i(i=1,\cdots, 6)$ ?

(2)Is the solution unique?

(3) are there general approaches to solve such problems (globally if possible)?

Update: Further data from replicates:

$$ \begin{array}{ccccccc} \text{No}&t& y_1(t) &y_2(t) &y_3(t) & y_4(t) &y_5(t)\\ 9&30.0000 & 9.0400 & 0.1190 & 0.1040 & 0.1390 & 0.0044 \\ 10& 60.0000 & 8.8000 & 0.1640 & 0.1120 & 0.3210 & 0.0097 \\ 11&90.0000 & 8.5300 & 0.1640 & 0.1140 & 0.5630 & 0.0219 \\ 12&120.0000 & 8.1800 & 0.1600 & 0.1250 & 0.8730 & 0.0369 \\ 13&150.0000 & 7.9700 & 0.1550 & 0.1380 & 1.0600 & 0.0459 \\ 14&180.0000 & 7.7900 & 0.1580 & 0.1510 & 1.2000 & 0.0545 \\ 15&210.0000 & 7.4900 & 0.1480 & 0.1430 & 1.5000 & 0.0636 \\ 16&240.0000 & 7.0800 & 0.1390 & 0.1380 & 1.9100 & 0.0756 \\ \end{array} $$

Update 2:

I am still seeking for a working approach to obtaining solutions to the problem.

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Isn't there a typo in the equations for ${y_4}'$ and ${y_5}'$? As it is stated, they aren't linear. $\endgroup$
    – Tadashi
    Commented Aug 9, 2014 at 23:36
  • 1
    $\begingroup$ Is there a typo in your equation for $y_1'$? $y_1$ appears twice on the right-hand side. $\endgroup$
    – Paul Tupper
    Commented Aug 10, 2014 at 0:01
  • $\begingroup$ This time it is not a typo, but because $k_1, k_2$ are different dynamic parameters. thank you $\endgroup$
    – LCFactorization
    Commented Aug 10, 2014 at 0:03

2 Answers 2

Reset to default
4
$\begingroup$

This kind of problem is known in the literature as a nonlinear state-space system identification. Several algorithms have been proposed in the literature to solve these problems. I think a good starting point would be (1) and the references therein, in particular the works of L. Ljung.

As far I know, in general if you don't have a good initial estimate of the parameters $k_i$, there is no guarantee that the algorithm will converge to a global solution.

(1) Schön, Thomas B., Adrian Wills, and Brett Ninness. "System identification of nonlinear state-space models." Automatica 47.1 (2011): 39-49.

P.S. Many of these identification algorithms have been implemented in MATLAB (maybe there is an equivalent implementation in Octave). See for instance this freely available third party toolbox and this example from the system identification toolbox of MATLAB.

Update (08/23/2014):

Following this example and using the first data set I found the following solution: $$\begin{array}{cccc} k_1 &= &1.100350 &\pm& 8.3610 \\ k_2 &= &2.588940 &\pm& 9.6769 \\ k_3 &= &0.400037 &\pm& 14.8607 \\ k_4 &= &0.338143 &\pm& 7.8111 \\ k_5 &= &3.757100 &\pm& 168.1370 \\ k_6 &= &0.915842 &\pm& 11.7712 \end{array}$$

The algorithm used was trust-region reflective Newton method of nonlinear least-squares, where the cost is the sum of squares of errors between the measured and simulated outputs.

Notice that the high standard deviation must be due to the small amount of samples.

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ It seems I will need more time to learn it. Can you please help obtain one solution to the above problem and indicate the method used? thank you $\endgroup$
    – LCFactorization
    Commented Aug 10, 2014 at 0:59
  • 1
    $\begingroup$ @LCFactorization, sorry for the delay in answering. I updated the answer with one solution (which isn't so good, but I think it'd need a higher sample rate to get a better model). $\endgroup$
    – Tadashi
    Commented Aug 23, 2014 at 16:44
  • 1
    $\begingroup$ thank you very much! there are more data here: ilovematlab.cn/thread-151440-1-1.html (simplified Chinese version, original question raised by a Chinese student), together there was also a solution by 1stOpt. I am only curious whether there is already systematic approach to solve such problems. thank you $\endgroup$
    – LCFactorization
    Commented Aug 24, 2014 at 1:51
  • 1
    $\begingroup$ @Shamisen When you mention the references 30 and 31 in your answer, which are these references? in the reference number [1] the bibliography is not enumerated. Please advse $\endgroup$
    – JuanMuñoz
    Commented Jul 14, 2020 at 19:15
  • $\begingroup$ Hi @JuanMuñoz, thanks for noticing it. When I put these numbers I was reading the pre-print version of the paper: user.it.uu.se/~thosc112/pubpdf/schonwn2011-2.pdf $\endgroup$
    – Tadashi
    Commented Jul 14, 2020 at 21:37
2
$\begingroup$

If you have access to Maple and the Global Optimization Toolbox, you could try something like this.

data:=[[0.0000,9.1300,0.0931,0.0899,0.1000,0.0000],
[30.0000,8.9300,0.1270,0.1230,0.2270,0.0049],
[60.0000,8.6000,0.1510,0.1390,0.4920,0.0153],
[90.0000,8.2800,0.1540,0.1490,0.7780,0.0249],
[120.0000,7.9700,0.1540,0.1570,1.0700,0.0329],
[150.0000,7.8600,0.1540,0.1600,1.1700,0.0348],
[180.0000,7.8100,0.1530,0.1530,1.2100,0.0404],
[210.0000,7.7700,0.1400,0.1420,1.2800,0.0432]]:
des:=[diff(y1(t),t) = -k1*y1(t) - k2*y1(t),
diff(y2(t),t) = k2*y1(t) - k3*y2(t),
diff(y3(t),t) = k1*y1(t) +k3*y2(t) - k4*y3(t),
diff(y4(t),t) = k4*y3(t) - k5*y2(t)*y4(t) + k6*y5(t),
diff(y5(t),t) = k5*y2(t)*y4(t) - k6*y5(t)]:
ics := seq((y || i)(0) = data[1, i+1], i = 1 .. 5):
res := dsolve({ics, des[]}, numeric, parameters = [k1, k2, k3, k4, k5, k6]):
timeList := [0, 30, 60, 90, 120, 150, 180, 210];
sse := proc (k1, k2, k3, k4, k5, k6) 

   res(parameters = [k1, k2, k3, k4, k5, k6]); 

   add(
   (rhs(select(has, res(timeList[i]), y1)[])-data[i, 2])^2+
   (rhs(select(has, res(timeList[i]), y2)[])-data[i, 3])^2+
   (rhs(select(has, res(timeList[i]), y3)[])-data[i, 4])^2+
   (rhs(select(has, res(timeList[i]), y4)[])-data[i, 5])^2+
   (rhs(select(has, res(timeList[i]), y5)[])-data[i, 6])^2, 
   i = 2 .. 8);

end proc;
c := GlobalOptimization:-GlobalSolve('sse'(k1, k2, k3, k4, k5, k6), k1 = 0 .. 1, k2 = 0 .. 1, k3 = 0 .. 1, k4 = 0 .. 1, k5 = 0 .. 1, k6 = 0 .. 1, timelimit = 10);
   [0.219132447080011505, [k1 = 8.52482740113834e-4, 

     k2 = 5.2683998680924474e-5, k3 = 0.0, 

     k4 = 0.05113239298267808, 

     k5 = 0.004363021255887466, k6 = 0.0]]
res(parameters = c[2]):
p:=Array(1..5):
for n from 1 to 5 do
   p[n]:=plots:-display(plots:-odeplot(res, [t, (y || n)(t)], t = 0 .. 210),plots:-pointplot([seq([data[i, 1], data[i, n+1]], i = 1 .. 8)]));
end do;
plots:-display(p)

The code is a parameterized numeric solution of the differential equations (I shifted all the data down by 30s, so I could get the ICs at t=0) followed by a global optimization of the least-squares problem.

I haven't tested the code for correctness, and if you let the optimizer run for longer or change the optimization bounds, you may get better results.

Improve this answer
$\endgroup$
3
  • $\begingroup$ thank you! I have no access to the global solver you mentioned. -- It was said another optimization software 1stOpt 7d-soft.com supporting such fitting from its version 4.0 (latest 6.0). But I am interested in the mathematical principles not just the answer. $\endgroup$
    – LCFactorization
    Commented Aug 15, 2014 at 23:10
  • $\begingroup$ Are the differential equation correct? Can you give me the actual stoichometric reactionsequations? $\endgroup$
    – Samir
    Commented Aug 19, 2014 at 18:30
  • $\begingroup$ I saw the question with 1stOpt answers from this link: ilovematlab.cn/thread-151440-1-1.html, it seems like chemical reactions and side reactions being combined together. $\endgroup$
    – LCFactorization
    Commented Aug 20, 2014 at 13:59

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .