Параметрична ідентифікація лінійної динамічної системи

Введення

Шановні читачі. В даний час процесам ідентифікації динамічних систем приділяється багато уваги. На цю тему написано багато дисертацій, дипломів і наукових публікацій. У різній літературі написано багато чого про ідентифікацію, наведено моделі та методи. Але все це для обивателя стає ясним не відразу. Я спробую в цій статті пояснити як вирішувати задачу параметричної ідентифікації, коли технічна система (об'єкт) описується системою диференціальних рівнянь за допомогою методу МНК.

Трохи з теорії

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

{{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_2y_m— це фазові координати x_1і (або) x_2обчислені за математичної моделі об'єкта. Але ж у моделі об'єкта, представленої у вигляді системи диференціальних рівнянь, x_1x_2не виражені в явному вигляді. Для того, щоб позбутися від цього безумства необхідно вирішити дану систему диференціальних рівнянь з заданими початковими умовами.

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

Sa_i^*=y_i
Потім підставляючи замість yзначення «виміряних» фазових координат, знаходимо оцінку параметра a_i^*для кожного моменту часу t_i.

Де взяти ці «виміряні» значення фазових координат?

Взагалі ці значення беруться з експерименту. Але так як ми ніякої експеримент не проводили, то візьмемо ці значення чисельного рішення системи диференціальних рівнянь методом Рунге-Кутта 4-5 порядку. Виберемо параметр a=-0,7

Рішення знайдемо вбудованими функціями пакета MatLab. Детальніше тут. Рішення даним методом показано нижче.

код MatLabfunction 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^*та «істинного» значення викликано помилками при розв'язанні системи диференціальних рівнянь методом Рунге-Кутта.

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.