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.



How to fit the parameters of differential equations with known data?
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.](http://www.sciencedirect.com/science/article/pii/S0005109810004279)" 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](http://sigpromu.org/idtoolbox/) and [this example](http://www.mathworks.com/help/ident/examples/represent-nonlinear-dynamics-using-matlab-file-for-grey-box-estimation.html) from the system identification toolbox of MATLAB.


**Update (08/23/2014):**


Following this [example](http://www.mathworks.com/help/ident/examples/represent-nonlinear-dynamics-using-matlab-file-for-grey-box-estimation.html) 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.


