Регрессионная идентификация линейных непрерывных систем управления

Теоретическая часть

Описание непрерывных стационарных систем управления в пространстве состояний имеет вид

$$\frac{dX(t)}{dt}=AX(t)+BU(t),\\ Y(t)=CX(t)+DU(t),$$

где: $$X(t)$$ — $$n$$ -мерный вектор состояния;

$$U(t)$$ — $$r$$ -мерный вектор управления (входные управляющие воздействия);

$$Y(t)$$ — $$m$$ -мерный вектор выхода системы, $$А$$ – матрица состояния размера $$n\times n$$ ;

$$В$$ — матрица входа размера $$n\times r$$ ;

$$С$$ — матрица выхода размера $$m\times n$$ ;

$$D$$ — матрица обхода размера $$m\times r$$.

Для строго реализуемых систем матрица $$D = 0$$.

В задачу регрессионной оценки (идентификации) линейной стационарной системы (12.1) может входить определение ее структуры и параметров по наблюдаемым данным — входным и выходным сигналам функционирующей системы. Если структура системы задана, т. е. известны дифференциальные уравнения, описывающие систему, то в задачу входит определение ее параметров — коэффициентов дифференциальных уравнений — матриц действительных чисел $$A,\mbox{ }B,\mbox{ }C,\mbox{ }D$$.

Для применения регрессионного анализа систему с непрерывным временем (12.1) следует представить в дискретной форме:

$$X[(k+1)T]=A_dX(kT)+B_dU(kT),\\ Y(kT)=CX(kT)+DU(kT),$$

где:

$$Т$$ — шаг квантования (период дискретизации) по времени;

$$k$$ — целые числа, $$A_{d}$$, $$B_{d}$$ ;

$$С,\mbox{ }D$$ — матрицы дискретной системы тех же размеров, что и для исходной непрерывной системы [7].

Матрицы $$A_{d}$$, $$B_{d}$$ имеют следующий вид:

$$A_d=e^{AT},$$ $$B_d=\int\limits_0^T e^{A\tau}Bd\tau.$$

Если матрица $$А$$ непрерывной системы не вырожденная, то матрицу можно представить в виде

$$B_d=A^{-1}(e^{AT}-E)B,$$

где $$Е$$ — единичная матрица $$n$$ -го порядка [3].

Поскольку в уравнениях дискретной системы шаг квантования по времени $$Т$$ входит слева и справа, то его часто опускают и систему записывают в виде

$$X(k+1)=A_dX(k)+B_dU(k),\\ Y(k)=CX(k)+DU(k).$$

По известной матрице $$A_{d}$$ дискретной системы можно определить матрицу непрерывной системы, логарифмируя обе части уравнения (12.3). При этом следует иметь в виду, что в (12.3) следует применять матричный экспоненциал, а для обратного преобразования — матричный логарифм. В системе MATLAB имеются функции матричного экспоненциала и матричного логарифма (см. $$help\mbox{ }expm$$ и $$help\mbox{ }logm$$ ).

В случае неособенной матрицы $$А$$ непрерывной системы уравнение (12.5) можно разрешить относительно матрицы $$В$$ в виде

$$B=(A^{-1}(e^{AT}-E))^{-1}B_d=(A^{-1}(A_d-E))^{-1}B_d.$$

Для достаточно малых значений шага квантования $$Т$$ можно воспользоваться следующей приближенной формулой:

$$B=\frac{B_d}{T}.$$

Таким образом, оценка матриц $$А$$ и $$В$$ непрерывной системы (12.1) осуществляется через оценку соответствующих матриц дискретной системы (12.2). При этом существуют более универсальные методы преобразования дискретной системы к непрерывной [11, 14].

Систему разностных уравнений запишем в виде скалярных уравнений:

$$x_1(k+1)=a_{d11}x_1(k)+a_{d12}x_2(k)+...+a_{d1n}x_n(k)+b_{d11}u_1(k)+...+b_{d1r}u_r(k),\\ x_2(k+1)=a_{d21}x_1(k)+a_{d22}x_2(k)+...+a_{d2n}x_n(k)+b_{d21}u_1(k)+...+b_{d2r}u_r(k),\\ ................................................\\ x_i(k+1)=a_{di1}x_1(k)+a_{di2}x_2(k)+...+a_{din}x_n(k)+b_{di1}u_1(k)+...+b_{dir}u_r(k),\\ ................................................\\ x_n(k+1)=a_{dn1}x_1(k)+a_{dn2}x_2(k)+...+a_{dnn}x_n(k)+b_{dn1}u_1(k)+...+b_{dnr}u_r(k),\\$$

В задачу входит оценка (идентификация) параметров $$a_{dij}$$, $$b_{diq}$$, $$i,j=\overrightarrow{1,n}$$, $$q=\overrightarrow{1,r}$$.

Выпишем из (12.8) промежуточное уравнение:

$$x_i(k+1)=a_{di1}x_1(k)+a_{di2}x_2(k)+...+a_{din}x_n(k)+b_{di1}u_1(k)+...+b_{dir}u_r(k).$$

В ходе эксперимента нужно запомнить $$g$$ решений $$u(k)$$, $$x(k)$$, $$x(k+1)$$ для идентификации параметров дискретной системы, причем $$g \ge n^2 + (r n) + 1 $$ [1]. Решения системы (12.8) определяются в результате подачи на нее некоторых входных сигналов (управляющие воздействия). Чтобы проверить, насколько построенная модель точно имитирует или предсказывает данные наблюдений, необходимо сравнить их при одинаковых воздействиях. Эта процедура называется верификацией модели [7].

Составим матрицу $$W_{ik}$$ из элементов $$x(k),\mbox{ }u(k)$$ размерностью $$(n+r)\times 1$$. Эта матрица будет представлять собой совокупность входных воздействий и выхода системы на момент дискретного времени $$k$$:

$$W_{ik}=[x_1(k),x_2(k),...,x_n(k),u_1(k),u_2(k),...,u_r(k)]^T.$$

Составим матрицу $$Ф_{i}$$ из искомых коэффициентов уравнения (12.9):

$$Ф_i^T=[a_{di1},a_{di2},...,a_{din},b_{di1},b_{di2},...,b_{dir}]_{1\times (n+r)},\mbox{ }i=\overrightarrow{1,n}.$$

С учетом (12.10), (12.11) запишем уравнение (12.9) в матричном виде:

$$x_i(k+1)=W_{ik}^TФ_i.$$

Сформируем вектор $$\chi_i$$ из элементов правой части уравнения (12.12) после $$g$$ испытаний (значений, решений, наблюдений):

$$\chi_i= \left[\begin{array}{c} x_{(1)i}(k+1)\\ x_{(2)i}(k+1)\\ \vdots\\ x_{(g)i}(k+1)\\ \end{array}\right].$$

Выражая $$\chi_i$$ через $$W_{ik}$$, $$Ф_i$$, получим $$\chi_i=W_kФ_i$$, где

$$\chi_i= \left[\begin{array}{cccccccc} x_{(1)1}(k)x_{(1)2}(k)...x_{(1)n}(k)u_{(1)1}(k)u_{(1)2}(k)...u_{(1),r}(k)\\ x_{(2)1}(k)x_{(2)2}(k)...x_{(2)n}(k)u_{(2)1}(k)u_{(2)2}(k)...u_{(2),r}(k)\\ ........................\\ x_{(g)1}(k)x_{(g)2}(k)...x_{(g)n}(k)u_{(g)1}(k)u_{(g)2}(k)...u_{(g),r}(k)\\ \end{array}\right].$$

Размерность матрицы $$W_{k}$$ равна $$g\times (n + r)$$.

Считая, что вектор $$\chi_i$$ и матрица $$W_{k}$$ известны (входные и выходные сигналы), можно применить метод наименьших квадратов, в соответствии с которым получим следующее нормальное уравнение относительно искомых параметров уравнения (12.9):

$$(W_k^TW_k)\hat Ф_i=W_k^T\chi_i.$$

Если матрица $$W_k^TW_k$$ невырожденная, то оптимальная оценка параметров дискретной системы определяется в виде

$$\hat Ф_i^{*}=(W_k^TW_k)^{-1}W_k^T\chi_i.$$

Зная $$\hat Ф_i^{*}$$, для всех индексов $$i=\overline{1,n}$$ можно найти коэффициенты системы уравнений (12.8), т. е. матрицы $$A_{d}$$, $$B_{d}$$, а затем, например, по формулам (12.3), (12.7) найти матрицы $$А,\mbox{ }В$$ непрерывной системы (12.1).

В случае, когда матрица $$W_k^TW_k$$ вырожденная или плохо обусловленная, при решении нормального уравнения (12.6) прибегают к псевдообращению, например, используют псевдообратную матрицу Мура–Пенроуза.

Расчет вектора состояния дискретной системы можно произвести по соотношению

$$X(kT)=A_d^kX(0)+\sum\limits_{j=0}^{k-1}A_d^{k-1-j}B_dU(j),$$

где $$X(0)$$ — начальный вектор состояния системы в момент времени, равный нулю [3].

Для регрессионной оценки матриц уравнения выхода $$C$$ и $$D$$ используются те же способы и приемы, которые были описаны для оценки матриц $$А$$ и $$В$$. В случае невырожденной матрицы $$W_k^TW_k$$ одновременная оценка матриц $$C$$ и $$D$$ может быть получена с помощью матричного уравнения следующего вида:

$$\overline{CD}=((W_k^TW_k)^{-1}W_k^TY_k)^T,$$

где

$$\overline{CD}= \left[\begin{array}{cccccccc} c_{11}c_{12}...c_{1n}d_{11}d_{12}...d_{1r}\\ c_{21}c_{22}...c_{2n}d_{21}d_{22}...d_{2r}\\ ........................\\ c_{m1}c_{m2}...c_{mn}d_{m1}d_{m2}...d_{mr}\\ \end{array}\right],\\ \\ \\ Y(k)= \left[\begin{array}{cccc} y_{(1)1}(k)y_{(1)2}(k)...y_{(1)m}(k)\\ y_{(2)1}(k)y_{(2)2}(k)...y_{(2)m}(k)\\ ............\\ y_{(g)1}(k)y_{(g)2}(k)...y_{(g)m}(k)\\ \end{array}\right].$$

Примечание. При верификации модели на вход модели и на реальную систему подают сигналы, которые не использовались при идентификации параметров и сравнении выходных сигналов модели и реальной системы.

Практическая часть

В практической части лабораторной работы будут использовать модели стационарных систем управления. Сигналы входа и выхода будут определяться для известной системы, при этом будем считать, что матрицы $$А$$, $$В$$ уравнения (12.1) и матрицы $$A_{d}$$, $$B_{d}$$ уравнения (12.2) известны. Но для оценки этих матриц будут использоваться только входные воздействия и значения векторов состояния.

Пример 1. Произведите регрессионную оценку матриц $$A$$, $$B$$ системы управления, состоящую из последовательного соединения трех инерционных звеньев с параметрами $$k_{1} = 1, T_{1} = 1, k_{2} = 2, T_{2} = 2, k_{3} = 3, T_{3} = 0.5$$. Входное воздействие на систему примите в виде $$u(t)=e^{-0.5t}\cos(2t)$$. Начальное состояние системы примите нулевым по всем переменным состояния.

Изобразим схему объекта управления, представленного через передаточные функции. Схема показана на рис. 12.1.

(рис 12.1) Схема объекта управления

Связь между входом и выходом, например, первого звена определяется соотношением, выраженным через изображение по Лапласу:

$$X_1(s)=W_1(s)U(s)=\frac{k_1U(s)}{T_1s+1}.$$

Передаточные функции определены для нулевых начальных условий. Поэтому комплексную переменную s формально можно заменить на производную, т. е.

$$(T_1s+1)X_1(s)=k_1U(s);\mbox{ }T_1sX_1(s)+X_1(s)=k_1U(s)\mbox{ }\{s=d/dt\}\\ T_1\frac{dx_1(t)}{dt}+x_1(t)=k_1u(t); \to \frac{dx_1(t)}{dt}=-\frac{1}{T_1}x_1(t)+\frac{k_1}{T_1}u(t).$$

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

$$\frac{dx_1(t)}{dt}=-\frac{1}{T_1}x_1(t)+\frac{k_1}{T_1}u(t);\\ \frac{dx_2(t)}{dt}=\frac{k_2}{T_2}x_1(t)-\frac{1}{T_2}x_2(t);\\ \frac{dx_3(t)}{dt}=\frac{k_3}{T_3}x_2(t)-\frac{1}{T_3}x_3(t);$$

Систему (12.21) можно представить в матричном виде. Тогда матрицы системы управления будут иметь следующий вид:

$$A= \left[\begin{array}{cccc} -\frac{1}{T_1}00\\ \frac{k_2}{T_2}-\frac{1}{T_2}0\\ 0\frac{k_3}{T_3}-\frac{1}{T_3} \end{array}\right], B= \left[\begin{array}{cccc} \frac{k_1}{T_1}\\ 0\\ 0 \end{array}\right], C=[0,0,1],D=0.\\$$

С учетом числовых значений параметров инерционных звеньев будем иметь

$$A= \left[\begin{array}{cccc} -100\\ 1-0.50\\ 062 \end{array}\right], B= \left[\begin{array}{cccc} 1\\ 0\\ 0 \end{array}\right], C=[0,0,1],D=0.$$

Программный код решения примера:

function LAB12;
clc,close all
%% Параметры инерционных звеньев
k1 = 1; T1 = 1;
k2 = 2; T2 = 2;
k3 = 3; T3 = 0.5;
%%%------------------------------
%% Матрицы непрерывной системы
global A B 
A = [-1/T1, 0, 0;
k2/T2, -1/T2, 0;
0, k3/T3, -1/T3]
B = [k1/T1; 0; 0]
C = [0, 0, 1]
D = 0
%%%------------------------------
%% Размерность системы управления
n = length(A);
%% Размерность управления
r = size(B, 2);
 %% Преобразование к дискретной системе
Ts = 0.01; % шаг квантования
%% Расчет матрицы Ad
Ad = expm(A*Ts)
%% Расчет матрицы Bd
if abs(det(A)) < 1e-10
Bd = inv(A)*(Ad - eye(size(A)))*B
 
else
syms tau

Bd2 = int(expm(-A*tau)*B, tau, 0, Ts);
Bd = double(Bd2)
end
%%%------------------------------
 %% Время наблюдений дискретной системы
Kend = 100; %% число шагов квантования
%% Вид входного воздействия
% % U(t) = exp(-0.5*t).*cos(2*t); 
%% Решение разностного уравнения
X0 = zeros(length(Ad),1);
Xk = zeros(length(Ad),Kend);
for k = 1 : Kend
sm = zeros(length(Ad), 1);
   for J = 0 : k-1
sm = sm + Ad^(k-1-J)*Bd*exp(-0.5*J*Ts)*cos(2*J*Ts);
   end
Xk(:,k) = (Ad^k)*X0 + sm;

Uk(1:r,k) = exp(-0.5*k*Ts)*cos(2*k*Ts);
end
 
%% Хранение массива KSI
KSI = [Xk(:, 2 : end)]';
 
%% Хранение массива Wk
Wk = [(Xk(:, 1 : end-1))',(Uk(1 : end-1))'];
 
%% Регрессионная оценка параметров дискретной системы
if abs(det(Wk'*Wk)) < 1e-10
F = (pinv(Wk'*Wk)*Wk'*KSI)'; 
 
else
F = (inv(Wk'*Wk)*Wk'*KSI)' ;   
end
%% Идентифицированные матрицы
Adr = F(:, 1 : n)
Bdr = F(:, n+1 : end)
 
%% Обратное преобразование - оценка матриц А, В
global Areg Breg
Areg = logm(Adr)/Ts
if Ts <= 1e-1 
Breg = Bdr/Ts
 
elseif Ts > 1e-1  abs(det(Areg)) <= 1e-10 
Breg = inv(inv(Areg)*(Adr - eye(size(Adr))))*Bdr    
    
else
    sd2 = ss(Adr, Bdr, C, D, Ts);
        if abs(min(eig(Adr)) ) < 1e-3
    sreg = d2c(sd2,'tustin');
        else
    
sreg = d2c(sd2,'zox');        
    end
[Areg, Breg, C, D] = ssdata(sreg)    
end
 
%%% Верификация модели
%%% Анализ систем при ступенчатом воздействии
global Um
Um = 12;
Xreg0 = zeros(length(Areg), 1);
X0 = zeros(length(A), 1);
T = [0, 12];
[treg, Xreg] = ode23(@fun, T, Xreg0);
[t, X] = ode23(@fun0, T, X0);
 
h12 = figure(1);
set(h12, 'name','Верификация при ступенчатом воздействии');
line(treg, Xreg, 'linew', 2);
line(t, X,'marker', 'o');
grid on
N = 2*length(Areg);
MN = [ [1:N/2],[1:N/2]]';
str1 = '\bf\it\fontname{times} x\rm\bf_';
str2 = num2str(MN);
str3 = '(\itt\rm\bf)';
legend(strcat(str1, str2, str3), 'location','best');
 
title('\bf\fontsize{10}\fontname{times}Результат верификации двух моделей систем управления');
xlabel('\bf\fontsize{12}\fontname{times}\it -  -  -  -  -  -  -  t -  -  -  -  -  -  -');
ylabel('\bf\fontsize{12}\fontname{times}\it X\rm\bf(\itt\rm\bf) ');
 %%%------------------------------
function f = fun0(t,X)
global A B Um
f = A*X + B*Um;
 
function f = fun(tr, Xr)
global Areg Breg Um
f = Areg*Xr + Breg*Um;

В программе преобразование матриц дискретной системы к непрерывной системе может выполняться при ряде условий, в частности, с помощью специализированных функций $$ss$$ (создание класса объекта – дискретной модели), $$d2c$$ (преобразование дискретной модели к непрерывной) и $$ssdata$$ (извлечение матриц системы). В функции $$d2c$$ используются методы преобразования $$zox$$ – экстаполятор нулевого порядка, $$tustin$$ – экстраполятор на основе билинейной аппроксимации Тастина [11, 14]. Метод $$zox$$ применяется, когда собственные числа матрицы состояния дискретной системы не лежат вблизи нуля комплексной плоскости.

Результат выполнения программы

A =
   -1.0000         0         0
    1.0000   -0.5000         0
         0    6.0000   -2.0000

B =
     1
     0
     0
C =
     0     0     1

D =
     0

Ad =
    0.9900         0         0
    0.0099    0.9950         0
    0.0003    0.0593    0.9802

Bd =
    0.0101
   -0.0001
    0.0000

Adr =
    0.9900    0.0000   -0.0000
    0.0099    0.9950    0.0000
    0.0003    0.0593    0.9802

Bdr =
    0.0101
   -0.0001
    0.0000

Areg =
   -1.0000    0.0000   -0.0000
    1.0000   -0.5000    0.0000
    0.0000    6.0000   -2.0000

Breg =
    1.0050
   -0.0050
    0.0001

function LAB122;
clc,close all
%% Параметры инерционных звеньев
k1 = 1; T1 = 1;
k2 = 2; T2 = 2;
k3 = 3; T3 = 0.5;
%%%------------------------------
%% Матрицы непрерывной системы
global A B 
A = [-1/T1, 0, 0;
k2/T2, -1/T2, 0;
0, k3/T3, -1/T3];
B = [k1/T1; 0; 0];
C = [0, 0, 1]
D = 0;
%%%------------------------------
%% Размерность системы управления
n = length(A);
%% Размерность управления
r = size(B, 2);
%% Размерность выхода
m = size(C, 1);
 
%% Преобразование к дискретной системе
Ts = 0.01; % шаг квантования
%% Расчет матрицы Ad
Ad = expm(A*Ts);
%% Расчет матрицы Bd
if abs(det(A)) < 1e-10
Bd = inv(A)*(Ad - eye(size(A)))*B;
 
else
syms tau
 
Bd2 = int(expm(-A*tau)*B, tau, 0, Ts);
Bd = double(Bd2);
end
%%%------------------------------
 
%% Время наблюдений дискретной системы
Kend = 100; %% число шагов квантования
%% Вид входного воздействия
% % U(t) = exp(-0.5*t).*cos(2*t); 
 
%% Решение разностного уравнения
X0 = zeros(length(Ad),1);
Xk = zeros(length(Ad),Kend);

for k = 1 : Kend
sm = zeros(length(Ad), 1);
   for J = 0 : k-1
sm = sm + Ad^(k-1-J)*Bd*exp(-0.5*J*Ts)*cos(2*J*Ts);
   end
Xk(:,k) = (Ad^k)*X0 + sm;
Uk(1:r,k) = exp(-0.5*k*Ts)*cos(2*k*Ts);
end

%% Массив выходных значений системы 
Ykd = C*Xk; %% как результат измерения

%% Хранение массива KSI
KSI = [Xk(:, 2 : end)]';
 
%% Хранение массива Wk
Wk = [(Xk(:, 1 : end-1))',(Uk(1 : end-1))'];
%% Регрессионная оценка параметров дискретной системы
if abs(det(Wk'*Wk)) < 1e-10
F = (pinv(Wk'*Wk)*Wk'*KSI)'; 
 
else
F = (inv(Wk'*Wk)*Wk'*KSI)' ;   
end

%% Идентифицированные матрицы
Adr = F(:, 1 : n);
Bdr = F(:, n+1 : end);
 
%% Хранение массива Yk
Yk = Ykd(:, 2:end);
 
%% Хранение массива YWk
YWk = Xk(:, 1:end-1);
 
%% Регрессионная оценка матриц дискретной системы
if abs(det(YWk*YWk')) < 1e-10
FY = (pinv(YWk*YWk')*YWk*Yk'); 
 
else
FY = (inv(YWk*YWk')*YWk*Yk'); 
end
 
%% Идентифицированная матрица выхода
Cdr = FY' %% для дискретной системы
Creg = Cdr %% для непрерывной системы
 
%% Обратное преобразование - оценка матриц А, В
global Areg Breg
Areg = logm(Adr)/Ts;
 
if Ts <= 1e-1 
Breg = Bdr/Ts;
 
elseif Ts > 1e-1  abs(det(Areg)) <= 1e-10 
Breg = inv(inv(Areg)*(Adr - eye(size(Adr))))*Bdr;    
 
else
    sd2 = ss(Adr, Bdr, C, D, Ts);
        if abs(min(eig(Adr)) ) < 1e-3
    sreg = d2c(sd2,'tustin');
        else
sreg = d2c(sd2,'zox');        
        end
 [Areg, Breg, C, D] = ssdata(sreg);    
end
 
%%% Верификация модели
%%% Анализ систем при ступенчатом воздействии
global Um
Um = 12;
Xreg0 = zeros(length(Areg), 1);
X0 = zeros(length(A), 1);
T = [0, 12];
[treg, Xreg] = ode23(@fun, T, Xreg0);
Yreg = Xreg*Creg'; %% выход системы
 
[t, X] = ode23(@fun0, T, X0);
Y = X*C'; %% выход системы
 
h12 = figure(1);
set(h12, 'name','Верификация при ступенчатом воздействии');
line(treg, Yreg, 'linew', 2);
line(t, Y, 'marker', 'o','color','r');
str1 = ...
 '\bf\fontsize{11}\fontname{times}\itY_r_e_g\rm\bf(\itt\rm\bf)';
str2 ='\bf\fontsize{11}\fontname{times}\itY\rm\bf(\itt\rm\bf)';
legend(str1, str2, 'location','best');
 
title('\bf\fontsize{10}\fontname{times}Результат верификации двух моделей по выходу');
xlabel('\bf\fontsize{12}\fontname{times}\it -  -  -  -  -  -  -  t -  -  -  -  -  -  -');
ylabel('\bf\fontsize{12}\fontname{times}\it Y\rm\bf(\itt\rm\bf) '); grid on

 
%%%------------------------------
function f = fun0(t,X)
%% Для заданной системы
global A B Um
f = A*X + B*Um;
 
function f = fun(tr, Xr)
%% Для идентифицированной системы 
global Areg Breg Um
f = Areg*Xr + Breg*Um;

C =
     0     0     1
Cdr =
    0.0003    0.0592    0.9802
Creg =
    0.0003    0.0592    0.9802