5
$\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.

$\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$ – Shamisen Aug 9 '14 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 Aug 10 '14 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 Aug 10 '14 at 0:03
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.

$\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 Aug 10 '14 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$ – Shamisen Aug 23 '14 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 Aug 24 '14 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 Jul 14 '20 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$ – Shamisen Jul 14 '20 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.

$\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 Aug 15 '14 at 23:10
  • $\begingroup$ Are the differential equation correct? Can you give me the actual stoichometric reactionsequations? $\endgroup$ – Samir Aug 19 '14 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 Aug 20 '14 at 13:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.