Algorithm of parameter identification of second-order dynamical system with small perturbations according to experimental data

Abstract

The problem of identifying the parameters of the second-order linear dynamical systems with small perturbations according to experimental data is solved by reducing it to the problem of minimizing a quadratic functional with constraints in the form of nonlinear algebraic equations. The calculations were carried out on three experimental data sets corresponding to the three types of singular points of a linear system: the saddle, the node, and the center.

Full Text

Рассмотрим линейную систему обыкновенных дифференциальных уравнений c малыми возмущениями вида

$\left\{\begin{array}{c}\frac{d{x}_1}{dt}={\theta }_1{x}_1+{\theta }_2{x}_2+{\epsilon }_1{x}_1\\ \frac{d{x}_2}{dt}={\theta }_3{x}_1+{\theta }_4{x}_2+{\epsilon }_2{x}_2\end{array}\right.\text{ (1)}$

где 𝑥_𝑖∈ ℝ,𝑖=1,2 – зависимые переменные, 𝑡∈[0,𝑏] – независимая переменная, 𝑏>0, 𝜃_𝑘∈ℝ, 𝑘=1,…,4 – неизвестные параметры, 𝜀_𝑗 (𝑗=1,2) – достаточно малые вещественные параметры.

Обозначим через 𝑥_𝑗(𝑡,𝜃) – 𝑗-компоненту решения системы (1), 𝑗=1,2, зависящую от векторного параметра 𝜃=𝑐𝑜𝑙𝑢𝑚𝑛(𝜃₁,𝜃₂,𝜃₃,𝜃₄).

Пусть при некоторых фиксированных значениях 𝜃_𝑘, 𝑘=1,…,4, решение системы (1) удовлетворяет задаче Коши с начальными условиями:

$x_1^{\left(1\right)}=x_1\left(0\right),x_2^{\left(1\right)}=x_2\left(0\right).\text{ (2)}$

Пусть так же по переменной 𝑥 с шагом $\tau \text{=}\frac{b}{N}$ на равномерной сетке

𝑡₁=0,…,𝑡_𝑖+1=𝑡_𝑖+𝜏,…,𝑡_𝑁=𝑏, (3)

для экспериментальных данных справедливы соотношения

𝑥̃^(𝑖)=𝑥^(𝑖)+𝜀^(𝑖), 𝑖=1,…,𝑁, (4)

где ${\tilde{x}}^{\left(i\right)}=column\left({\tilde{x}}_1^{\left(i\right)},{\tilde{x}}_2^{\left(i\right)}\right),x^{\left(i\right)}=column\left(x_1^{\left(i\right)},x_2^{\left(i\right)}\right),x_j^{\left(i\right)}=x_j\left(t_i,\theta \right)$ - значение компоненты решения системы (1) в точке $t_i$ при фиксированном значении векторного параметра $\theta ,{\epsilon }^{\left(i\right)}=column\left({\epsilon }_{i\mathrm{1,}}\mathrm{...}{\epsilon }_{iN}\right)$ – вектор, элементы которого являются случайными величинами, имеющими стандартное нормальное распределение, то есть ${\epsilon }_{ij}$∈ N (0, 1), 𝑖 =$\overline{\mathrm{1,}N}$, 𝑗=1,2.

Ставится задача идентификации параметров системы вида (1), заключающаяся в нахождении таких оценок ${\widehat{\theta }}_k$ параметров ${\theta }_k$, k=1,…4, при которых решение задачи (1) приближается экспериментальными данными $\left\{{\tilde{x}}^{\left(i\right)},i=\overline{\mathrm{1,}N}\right\}$ в смысле метода наименьших квадратов [1].

Заменим уравнение (1) симметричной разностной схемой [2] на сетке (3)

$\begin{array}{l}\frac{x_1^{\left(i+1\right)}-x_1^{\left(i\right)}}{\tau }=\frac{1}{2}\left(f_{i\mathrm{,1}}+f_{i+\mathrm{1,1}}\right),\\ \frac{x_2^{\left(i+1\right)}-x_2^{\left(i\right)}}{\tau }=\frac{1}{2}\left(f_{i\mathrm{,2}}+f_{i+\mathrm{1,2}}\right),\text{ (5)}\\ i=\overline{\mathrm{1,}N-1},\end{array}$

где

$\begin{array}{l}{f}_{i\mathrm{,1}}=f_1\left(x^{\left(i\right)},\theta \right)\equiv {\theta }_1{x}_1^{\left(i\right)}+{\theta }_2{x}_2^{\left(i\right)}+{\epsilon }_1{x}_1^{\left(i\right)},\\ f_{i\mathrm{,2}}=f_2\left(x^{\left(i\right)},\theta \right)\equiv {\theta }_3{x}_1^{\left(i\right)}+{\theta }_4{x}_2^{\left(i\right)}+{\epsilon }_2{x}_2^{\left(i\right)}.\end{array}$

Вводя обозначения [3]

$\begin{array}{l}z=column\left(x^{\left(N\right)},\theta \right),x^{\left(N\right)}=column\left(x^{\left(1\right)}\mathrm{,...,}{x}^{\left(N\right)}\right),\\ \tilde{z}=column\left({\tilde{x}}^{\left(N\right)},o_4\right),{\tilde{x}}^{\left(N\right)}=column\left({\tilde{x}}^{\left(1\right)}\mathrm{,...,}{\tilde{x}}^{\left(N\right)}\right),\end{array}$

$o_4$– нулевой вектор размерности 4,

получим

$x^{\left(N\right)}=H_{1}z,{\tilde{x}}^{\left(N\right)}={\tilde{H}}_{1}z,$

где

$H_1=\left[I_{2N}⋮O_{2N\times 4}\right],{\text{ H}}_2=\frac{1}{N}{H}_1^T{H}_1,$

здесь 𝐼_2𝑁 – единичная (2𝑁×2𝑁)-матрица, 𝑂_2𝑁×4 – нулевая (2𝑁×4)-матрица.

Тогда, согласно [3–5], задача идентификации параметров может быть сформулирована как задача минимизации квадратичного функционала с ограничениями в виде нелинейных алгебраических уравнений

$\left\{\begin{array}{ll}\underset{z}{\min m\left(z\right),}& m\left(z\right)=\frac{1}{2}(H_{2}z,z)-(H_2\tilde{z},z)\\ & T\left(z\right)=h\\ & g\left(z\right)=0\end{array}\right.\text{ (6)}$

Здесь введены следующие обозначения: 𝑇=[𝑇₁ 𝑇₂,…,𝑇_𝑁,𝑇_𝜃] – (2×(2𝑁+4)) – матрица, 𝑇_𝑖, 𝑖=$\overline{\mathrm{2,}N}$– нулевые (2×2) – матрицы, 𝑇_𝜃 – (2×4) – матрица;

$g\left(z\right)=column\left(g_1\left(z\right)\mathrm{,...,}{g}_N\left(z\right)\right),g_i\left(z\right)column\left(g_{i1}\left(z\right),g_{i2}\left(z\right)\right),i=\overline{\mathrm{1,}N},$

$\begin{array}{l}{g}_{i1}\left(z\right)=g_{i1}\left(x^{\left(N\right)},\theta \right)\equiv \\ \equiv \left(1+\frac{\tau }{2}\left({\theta }_1+{\epsilon }_1\right)\right)x_1^{\left(i\right)}+\frac{\tau }{2}{\theta }_2{x}_2^{\left(i\right)}-\left(1-\frac{\tau }{2}\left({\theta }_1+{\epsilon }_1\right)\right)x_1^{\left(i+1\right)}+\frac{\tau }{2}{\theta }_2{x}_2^{\left(i+1\right)},\text{ (7)}\end{array}$

$\begin{array}{l}{g}_{i2}\left(z\right)=g_{i2}\left(x^{\left(N\right)},\theta \right)\equiv \\ \equiv \frac{\tau }{2}{\theta }_3{x}_1^{\left(i\right)}+\left(1+\frac{\tau }{2}\left({\theta }_4+{\epsilon }_2\right)\right)x_2^{\left(i\right)}+\frac{\tau }{2}{\theta }_3{x}_1^{\left(i+1\right)}-\left(1-\frac{\tau }{2}\left({\theta }_4+{\epsilon }_2\right)\right)x_2^{\left(i+1\right)}\text{. (8)}\end{array}$

Заметим, что разностная схема (5) с учетом обозначений (7)-(8) может быть записана в виде

$g\left(z\right)=0$

Для решения задачи (6) воспользуемся алгоритмом [3], основанным на аппроксимации исходной задачи последовательностью квадратичных задач минимизации с линейными ограничениями. На каждом шаге разреженная система линейных алгебраических уравнений большой размерности решалась с использованием метода сопряженных градиентов [4].

Вычисления проводились на трех наборах экспериментальных данных 𝑥̃^(𝑖), $i=\overline{\mathrm{1,}N}$, отличающихся от приближенных решений с начальными данными 𝑥₁(0)=1, 𝑥₂(0)=1 системы (1) на случайные величины 𝜀_𝑖𝑗∈𝑁(0,1) (𝑗=1,2), имеющие стандартное нормальное распределение. Каждое приближенное решение системы (1) получено при фиксированном наборе параметров 𝜃_𝑘, 𝑘=1,…,4, соответствующих трем типам особой точки системы (1): седло, узел, центр. Точность вычислений 𝛿 полагалась равной 0.001. В результате вычислений для трех различных наборов экспериментальных данных получены оценки ${\widehat{\theta }}_k$ для параметров 𝜃_𝑘, 𝑘=1,…,4, а также компоненты 𝑥₁(𝑡, $\widehat{\theta }$), 𝑥₂(𝑡, $\widehat{\theta }$) приближенных решений системы (1) .

Приведем графики компонент 𝑥₁(𝑡, $\widehat{\theta }$), 𝑥₂(𝑡, $\widehat{\theta }$) приближенных решений системы (1) и соответствующие им экспериментальные данные.

Случай 1. Нулевое положение равновесия системы (1) – седло.

Получены следующие оценки параметров 𝜃₁,𝜃₂,𝜃₃,𝜃₄:

${\widehat{\theta }}_1$=2.02298,  ${\widehat{\theta }}_2$=0.0001,  ${\widehat{\theta }}_3$=0.00002, ${\widehat{\theta }}_1$=−2.90002

Случай 2. Нулевое положение равновесия системы (1) – узел.

Получены следующие оценки параметров 𝜃1,𝜃2,𝜃3,𝜃4:

${\widehat{\theta }}_1$=−3.14132, ${\widehat{\theta }}_2$=0.01895, ${\widehat{\theta }}_3$=1.03559, ${\widehat{\theta }}_4$=0.02129.

Случай 3. Нулевое положение равновесия системы (1) – центр.

Получены следующие оценки параметров 𝜃1,𝜃2,𝜃3,𝜃4:

${\widehat{\theta }}_1$=0.07274,  ${\widehat{\theta }}_2$=4.01077,  ${\widehat{\theta }}_3$=−3.98702,  ${\widehat{\theta }}_4$=−0.02925.

Рис. 1. Экспериментальные данные и график компоненты ${\mathit{x}}_{\mathbf{1}}\left(t,\widehat{\theta }\right)$ приближенного решения системы (1).

Рис. 2. Экспериментальные данные и график компоненты ${\mathit{x}}_{\mathbf{2}}\left(t,\widehat{\theta }\right)$ приближенного решения системы (1).

Рис. 3. Экспериментальные данные и приближение фазовой траектории системы (1) с начальными данными 𝑥₁(0)=1,𝑥₂(0)=1.

Рис. 4. Экспериментальные данные и график компоненты ${\mathit{x}}_{\mathbf{1}}\left(t,\widehat{\theta }\right)$ приближенного решения системы (1).

Рис. 5. Экспериментальные данные и график компоненты ${\mathit{x}}_{\mathbf{2}}\left(t,\widehat{\theta }\right)$ приближенного решения системы (1).

Рис. 6. Экспериментальные данные и приближение фазовой траектория системы (1) с начальными данными 𝑥₁(0)=1, 𝑥₂(0)=1.

Рис. 7. Экспериментальные данные и график компоненты ${\mathit{x}}_{\mathbf{1}}\left(t,\widehat{\theta }\right)$ приближенного решения системы (1).

Рис. 8. Экспериментальные данные и график компоненты ${\mathit{x}}_{\mathbf{2}}\left(t,\widehat{\theta }\right)$ приближенного решения системы (1).

Рис. 9. Экспериментальные данные и приближение фазовой траектории системы (1).