State and parameter estimation using unconstrained optimization

We present an efficient method for estimating variables and parameters of a given system of ordinary differential equations by adapting the model output to an observed time series from the (physical) process described by the model. The proposed method is based on (unconstrained) nonlinear optimization exploiting the particular structure of the relevant cost function. To illustrate the features and performance of the method, simulations are presented using chaotic time series generated by the Colpitts oscillator, the three-dimensional Hindmarsh-Rose neuron model, and a nine-dimensional extended Rössler system.

Figure 1

Example of a nonsmooth solution $y\left(t_m\right)$ (solid curve) oscillating close to the true solution $x\left(t_m\right)$ (dashed curve) due to an instability that leads to an alternating sequence of integration values.

Figure 2

Adaption of the Colpitts oscillator model (17) to the time series $\eta \left(t_n\right)$ [Eq. (19)]. (a) The model output given by the measurement function $h$ [Eq. (20); dark blue dots] was adapted to the time series $\eta \left(t_n\right)$ (bright green circles). (b), (d), and (f) The estimated variables $y_1$, $y_2$, and $y_3$ (blue dots) and original (“true”) variables $x_1$, $x_2$, and $x_3$ (orange dashed line; not directly observed) used to generate $\eta$. (c), (e), and (g) Residues [see Eq. (5)] $u_i$, $i=1,2,3$ representing deviations from the model ODE's. In addition to the model variables, the model parameters are estimated at $\mathbf{p}=(5.23,0.0777,0.677,6.28)$ whereas the parameters used to integrate (18) and to generate the time series $\eta \left(t_n\right)$ are $(5,0.08,0.7,6.3)$. The parameter of the measurement function is estimated at $q=1.29$ [while $1.3$ was used in Eq. (19) to generate $\eta \left(t_n\right)$]. For better visibility, the results are only shown from $t=0$ to 30 while the model time series (20) was adapted to the given data (19) from $t=0$ to 150.

Figure 3

Adaption of the Hindmarsh Rose model (21) to the time series $\eta \left(t_n\right)$ [Eq. (23)]. (a) The model output given by the measurement function $h$ [Eq. (24); dark blue dots] was adapted to the time series $\eta \left(t_n\right)$ (bright green circles). (b), (d), and (f) The estimated variables $y_1$, $y_2$, and $y_3$ (blue dots) and original (“true”) variables $x_1$, $x_2$, and $x_3$ (orange dashed line; not directly observed) used to generate $\eta$. (c), (e), and (g) Residues [see Eq. (5)] $u_i$, $i=1,2,3$ representing deviations from the model ODE's. In addition to the model variables, the model parameters are estimated at $\mathbf{p}=(2.99,4.98,0.00397,3.19)$ whereas the parameters used to integrate (22) and to generate the time series $\eta \left(t_n\right)$ are $(3,5,0.004,3.19)$. The parameter of the measurement function is estimated at $q=1.791$ [while $1.8$ was used in Eq. (23) to generate $\eta \left(t_n\right)$].

Figure 4

Adaption of the extended Rössler system (25) to the time series $\eta \left(t_n\right)$ given by Eq. (27). (a) The model output given by the measurement function $h$ [Eq. (28); dark blue dots] was adapted to the time series $\eta \left(t_n\right)$ (bright green circles). (b), (d), (f), (h), and (j) The estimated variables $y_1$, $y_2$, $y_5$, and $y_9$ (blue dots) and corresponding original (“true”) variables $x_1$, $x_2$, $x_5$, and $x_9$ (orange dashed line; not directly observed) used to generate $\eta$ (variables $y_3$, $y_4$, $y_6$, and $y_7$ not shown here exhibit a similar dynamics and coincide with the corresponding $x$ variables as well. (c), (e), (g), (i), and (k) Residues [see Eq. (5)] $u_i$, $i=1,2,5,8,9$ representing deviations from the model ODE's. In addition to the model variables, the model parameters are estimated to be $\mathbf{p}=(0.293,0.0801,4.08,1.96)$ whereas the parameters used to integrate (26) and to generate the time series $\eta \left(t_n\right)$ are $(0.3,0.1,4,2)$. The parameter of the measurement function is estimated at $q=1.69$ [while $1.7$ was used in Eq. (27) to generate $\eta \left(t_n\right)$]. For better visibility, the results are only shown from $t=0$ to 30 while the model time series (28) was adapted to the given data (27) from $t=0$ to 100.

Figure 5

The Colpitts oscillator model (17) is adapted to the time series (29). $p_2$ and $p_3$ are fixed to $p_2=0.08$ and $p_3=0.7$, $p_1$ and $p_4$ are estimated in addition to $\mathbf{y}\left(t_m\right)$ to ${\widehat{p}}_1=6.33$ and ${\widehat{p}}_4=4.97$ [$p_1=5$ and $p_3=6.3$ are used to integrate (18)]. The meaning of the line styles and colors is the same as in Fig. 2.

Figure 6

Contour lines of the minimized cost function (12) of the Colpitts oscillator model (17) when adapted to an $x_1$ time series (29) after transients decayed [(a), (b)] and including transient [(c), (d)]. In (a) and (c), the parameters $p_2$ and $p_3$ are fixed to $p_2=0.08$ and $p_3=0.7$ and the cost function is optimized for different combinations of $p_1$ and $p_4$ values resulting in a function $C(\widehat{\mathbf{y}}\left(t_m\right),p_1,p_4)$ (visualized by means of contour lines). The eigenvectors of the Hessian at the optimum $\widehat{\mathbf{p}}=({\widehat{p}}_1,{\widehat{p}}_4)=(6.29,5.00)$ are ${\mathbf{h}}_1=(-0.780,0.63){}^{\mathrm{T}}$ and ${\mathbf{h}}_2=(0.63,0.780){}^{\mathrm{T}}$ with eigenvalues ${\lambda }_1=4.5\times {10}^{-3}$ and ${\lambda }_2=0.547$, respectively. In (b) and (d), parameters $p_1$ and $p_4$ are fixed to $p_1=5$ and $p_4=6.3$, the cost function is minimized on a grid of $p_2$-$p_3$ values, and the contour plot shows $C(\widehat{\mathbf{y}}\left(t_m\right),p_2,p_3)$. The minimum is located at $\widehat{\mathbf{p}}=({\widehat{p}}_2,{\widehat{p}}_3)=(0.080,0.70)$ and the eigenvectors of the Hessian at $\widehat{\mathbf{p}}$ are ${\mathbf{h}}_1=(-0.0534,-0.9986){}^{\mathrm{T}}$ and ${\mathbf{h}}_2=(-0.9986,0.0534){}^{\mathrm{T}}$ with eigenvalues ${\lambda }_1=10.4$ and ${\lambda }_2=2.56\times {10}^3$, respectively. Note that ${\mathbf{h}}_1$ and ${\mathbf{h}}_2$ do not appear orthogonal in the plots due to the different scaling of the axes.

Figure 7

Parameter estimation results for the Colpitts oscillator (17) based on noisy observed $x_1$ time series (29) consisting of $15.150=2250$ samples. Six noise levels ranging from the noiseless case ($\infty$db) to 5.3 dB are considered as indicated in the legend in (b) that holds for all subfigures. The given time series was split into 15 parts, and for each segment parameters have been individually estimated. (a), (c), (e) Variance of the 15 estimated parameters logarithmically plotted vs. ${\log }_{10}(\alpha /(1-\alpha ))$, where $\alpha$ denotes the continuation parameter. (b), (d), (f) Deviations of the mean values of the estimated parameters $p_1$, $p_2$, and $p_3$ from the values used for generating the data.

Figure 8

Consistency measures vs fixed (detuned) parameter $p_1$ of the Colpitts oscillator (17). State variables and parameters $p_2$, $p_3$, and $p_4$ are estimated. A high consistency of model and data is indicated by $R$ values close to 1. The indicator $\overline{R}$ given by Eq. (35) uniquely selects the value $p_1=5$ that was also used for generating the (observed) time series.