[Из песочницы] Параметрическая идентификация линейной динамической системы11.03.2016 13:33

Введение

Уважаемые читатели. В настоящее время процессам идентификации динамических систем уделяется много внимания. На эту тему написано много диссертаций, дипломов и научных публикаций. В различной литературе написано много чего про идентификацию, приведены различные модели и методы. Но всё это для обывателя становится ясным не сразу. Я попытаюсь в этой статье объяснить как решать задачу параметрической идентификации, когда техническая система (объект) описывается системой дифференциальных уравнений, с помощью метода МНК.

Немного из теории

Для начала нужно понять, что такое динамическая система. Если говорить как можно проще, то это система, параметры которой изменяются во времени. Подробнее здесь. Практически любую динамическую систему можно описать дифференциальным уравнением какого-либо порядка, например:

${{dx_1\over dt}=ax_1+bx_2} \\ {dx_2\over dt}=cx_1+dx_2} .$

Данная система дифференциальных уравнений характеризуется своими параметрами. В нашем случае это a, b, c и d. Они могут быть как статическими так и динамическими.

Что означают эти коэффициенты?

Применительно к реальным физическим динамическим системам эти коэффициенты дифференциального уравнения имеют конкретную физическую привязку. Например в системе ориентации и стабилизации космического аппарата данные коэффициенты могут играть различную роль: коэффициент статической устойчивости КА, коэффициент эффективности бортового управления, коэффициент способности изменять траекторию и т.п. Подробнее здесь.

Так вот задача параметрической идентификации это определение этих самых коэффициентов параметров a, b, c и d.

Задача наблюдения и измерения

Стоит отметить, что для решения задачи параметрической идентификации необходимо получить «измерения» одной (или всех) фазовой координаты (в нашем случае это x1 и (или) x2).

Для того чтобы система была идентифицируема, она должна быть наблюдаема. То есть ранг матрицы наблюдаемости должен быть равен порядку системы. Подробнее про наблюдаемость здесь.

Наблюдение процессов, происходящих в объекте, происходит следующим образом:

$y=Hx+\xi(t).$

где:

у — вектор наблюдаемых параметров;
H — матрица связи параметров состояния и наблюдаемых параметров;

$\xi(t)$ — помеховая составляющая (в ней спрятаны все погрешности наблюдения);

Подробнее про вектора и матрицы

Динамическую систему которую мы описывали выше, можно представить в векторно-матричной форме:

${dx\over dt}=Ax+Bu+\varepsilon(t)$

где:

x — вектор состояния динамической системы. В общем случае он имеет бесконечную размерность. Но когда мы имеем дело с конкретной моделью, то уже говорим не о «векторе состояния», а о «векторе фазовых координат». В нашем примере он имеет две составляющие
$x=\begin{pmatrix}x_1\\x_2\end{pmatrix}$
A — переходная матрица. Содержит те самые коэффициенты которые нам хотелось бы найти.
$A=\begin{pmatrix}a &b\\c&d\end{pmatrix}$
u — управляющее воздействие. В нашем примере оно равно 0. То есть наш объект не управляем. Прошу обратить внимания что это просто пример, не имеющий ничего общего с какой-то реально существующей системой.
B — матрица управления.

$\varepsilon(t)$ — помеховая составляющая.

Измерение процессов, происходящих в объекте, описывается следующим образом:

${{z}=y+n(t)} \\ {z}=n_1(t)y+n_2(t)} .$

Как мы видим погрешность измерения может быть как аддитивной (в первом случае), так и мультипликативной (во втором)

Задача идентификации

Рассмотрим решение задачи параметрической идентификации в случае когда не известен один коэффициент. Перейдем к конкретному примеру. Пусть дана следующая система:

${{dx_1\over dt}=ax_1+x_2} \\ {dx_2\over dt}=0.0225x_1-0.3x_2} \\ {x_1(0)=20 \ x_2(0)=20}$

Видно, что параметры равны b = 1, c = 0.0225 и d = -0.3. Параметр a нам неизвестен. Попробуем дать его оценку с помощью метода наименьших квадратов.

Задача состоит в следующем: по имеющимся выборочным данным наблюдений за выходным сигналами с интервалом дискретизации Δt требуется оценить значения параметра, обеспечивающего минимум величины функционала невязки между модельными и фактическими данными.

$J=(y-y_m)^T(y-y_m)=e^Te=\sum\limits_{j=1}^{N}e_j^2\longrightarrow\min$

где $e_j=y_i-y_m_i$ — невязка, определённая как разность между выходом исследуемого объекта и реакцией, вычисленной по математической модели объекта.

Невязка складывается из неточностей структуры модели, погрешностей измерений и неучтённых взаимодействий среды и объекта. Однако, независимо от природы возникающих ошибок, метод наименьших квадратов минимизирует сумму квадратичной невязки для дискретных значений. В принципе, МНК не требует никакой априорной информации о помехе. Но для того, чтобы полученные оценки обладали желательными свойствами, будем предполагать, что помеха является случайным процессом типа белого шума.

Оценка $a^*$ по методу наименьших квадратов, минимизирующая критерий J, находится из условия существования минимума функционала:

$J=\min\limits_{a} J=J |_{a=a^*}$

Важным свойством оценок по МНК является существование только одного локального минимума, совпадающего с глобальным. Поэтому оценка $a^*$ является единственной. Ее значение определяется из условия экстремума функционала J:

${{\partial J\over \partial a}|_{a=a^*}=2(y-y_m)(y-y_m)^{\prime} =0$

То есть необходимо от функционала $J=(y-y_m)^2$ взять производную по a и приравнять ее к нулю.

Обращаю внимание, что $y$ — это «измеренные» значения фазовых координат $x_1$ и (или) $x_2$ , а $y_m$ — это фазовые координаты $x_1$ и (или) $x_2$ вычисленные по математической модели объекта. Но ведь в модели объекта, представленной в виде системы дифференциальных уравнений, $x_1$ и $x_2$ не выражены в явном виде. Для того, чтобы избавиться от этого безумия необходимо решить данную систему дифференциальных уравнений с заданными начальными условиями.

Решать можно как «вручную», так и используя какое-либо программное обеспечение. Ниже будет показано решение в MatLab. В итоге должна получится система алгебраических уравнений для каждого момента времени $t_i$ :

$Sa_i^*=y_i$

Затем подставляя вместо $y$ значения «измеренных» фазовых координат, находим оценку параметра $a_i^*$ для каждого момента времени $t_i$ .

Где взять эти «измеренные» значения фазовых координат?

Вообще эти значения берутся из эксперимента. Но так как мы никакой эксперимент не проводили, то возьмем эти значения из численного решения нашей системы дифференциальных уравнений методом Рунге-Кутта 4–5 порядка. Выберем параметр $a=-0,7$

Решение найдем встроенными функциями пакета MatLab. Подробнее здесь. Решение данным методом показано ниже.

код MatLab

function F=diffurafun (t, x)
F=[-0.7*x (1)+x (2);0.0225*x (1)-0.3*x (2)];
end
%===============================================================%
% формирование вектора начальных условий
X0=[20;20];
% формирование интервала интегрирования
interval=[0 50];
% обращение к функции оde 45
[T, X]=ode45(@diffurafun, interval, X0);
% вывод графика решения
figure; plot (T, X (:,1),'r-', T, X (:,2),'b--'); % вывод графика
legend ('Параметр х1',' Параметр х2');
grid on;

График

На данном графике красной сплошной линией обозначена фазовая координата $x_1$ , а синей пунктирной линией обозначена фазовая координата $x_2$

Ниже показан код на MatLab и график.

Код MatLab

% dx/dt = A*x — линейная динамическая система без управления и помех
% для начала нам необходимо решить аналитически данную СДУ
% обозначим тип переменных
syms x (t) y (t) a
% решим систему при заданных начальных условиях
S = dsolve (diff (x) == a*x + 1*y,'x (0)=20', diff (y) == 0.0225*x — 0.3*y,'y (0)=20');
% выберем решение первой фазовой координаты, так как именно в его уравнении
% содержится искомый параметр а
x (t) = S.x;
% найдем частную производную первого уравнения по параметру, а (в
% соответствии с методом МНК)
f=diff (x (t),'a');
% теперь немного упростим получившееся выражение
S1=simplify (f);
% зададим переменной t массив значений T
t=T;
% найдем выражения, содержашие параметр, а для каждого момента времени
SS=eval (S1);
% теперь в цикле, подставляя в каждое выражение значение «измеренной»
% первой фазовой координаты, определим параметр, а для каждого момента
% времени T. Значения «измеренной» фазовой координаты берем из решения СДУ
% методом Рунге-Кутта 4-ого порядка
for i=2:81
SSS (i)=solve (SS (i)==X (i,1), a);
end
ist=zeros (length (T),1);
ist (1: length (T))=-0.7;
figure; plot (T, SSS,'b--', T, ist,'r-');
legend ('оценка параметра а','истинное значение');
grid on;

На графике синей пунктирной линией обозначена оценка параметра $a^*$ , а красной сплошной линией обозначено непосредственно «истинное» значение параметра модели $a=-0,7$ . Мы видим, что примерно на 3,5 секунде процесс стабилизируется. Небольшое расхождение оценки параметра $a^*$ и «истинного» значения вызвано ошибками при решении системы дифференциальных уравнений методом Рунге-Кутта.