A =
x1 x2 x3
x1 0.77993 0.54376 -0.013931
x2 -0.12996 0.073872 0.19929
x3 -0.1067 0.016064 -0.47392
B =
расход газа
x1 -0.022796
x2 -0.048006
x3 -0.034112
C =
x1 x2 x3
температура -4.6742 -0.54702 0.0027609
D =
расход газа
температура 0
K =
температура
x1 -0.20139
x2 0.075372
x3 0.02383
x(0) =
x1 0.10598
x2 0.033411
x3 -0.0020287
Estimated using N4SID from data set zdanv
Loss function 0.00753932 and FPE 0.00791011
Sampling interval: 0.08
Функция pem оценивает параметры обобщенной многомерной линейной модели:>> zpem=pem(zdanv)
State-space model: x(t+Ts) = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2 x3 x4 x5
x1 0.79771 0.52932 0.017441 -0.043818 -0.055481
x2 -0.089871 -0.056516 0.34751 0.4433 -0.14175
x3 0.12179 -0.31736 -0.29564 -0.43427 0.29463
x4 0.096387 -0.0086843 -0.063652 0.64711 -0.21098
x5 -0.012304 0.025405 -0.18701 0.69982 1.0514
B =
расход газа
x1 -0.021221
x2 -0.053799
x3 -0.040284
x4 -0.00059208
x5 -0.032983
C =
x1 x2 x3 x4 x5
температура -4.675 -0.54794 0.010925 0.068154 -0.1202
D =
расход газа
температура 0
K =
температура
x1 -0.22248
x2 0.037523
x3 0.024536
x4 0.11512
x5 -0.068061
x(0) =
x1 0.10985
x2 0.0039575
x3 0.071289
x4 -0.15615
x5 -0.15783
Estimated using PEM from data set zdanv
Loss function 0.00664935 and FPE 0.00720346
Sampling interval: 0.08
2. 9. Функции преобразования моделей
Для дальнейшего использования, полученных моделей при анализе и синтезе систем автоматизации технологических процессов в пакете System Identification Toolbox MATLAB имеются специальные функции, позволяющие выполнить преобразование этих моделей из тета-формата (внутреннего вида матричных моделей, являющегося дискретным) в другие виды и в частности из дискретной в непрерывную модель в виде передаточной функции.
Функция th2arx преобразует модель тета-формата в ARX модель:
Функция имеет синтаксис:
>> [A,B]=th2arx(darx)
где darx - модель тета-формата
A =
Columns 1 through 8
1.0000 -1.0105 0.3552 -0.0347 -0.1432 0.1302 -0.0128 -0.0858
Column 9
0.0630
B =
0 0.1367 0.0733
Функция th2ff вычисляет частотные характеристики и соответствующие стандартные отклонения по модели в тета-формате.
В качестве аргумента функции может выступать любая из рассмотренных ранее моделей, например darx:
>> [g,phiv]=th2ff(darx)
IDFRD model g.
Contains Frequency Response Data for 1 output and 1 input
and SpectrumData for disturbances at 1 output
at 129 frequency points, ranging from 0.1 rad/s to 39.27 rad/s.
Output Channels: температура
Input Channels: расход газа
Sampling time: 0.08
Estimated from data set zdanv using ARX.
IDFRD model phiv.
Contains SpectrumData for 1 signal
at 126 frequency points, ranging from 0.1 rad/s to 39.27 rad/s.
Output Channels: температура
Sampling time: 0.08
Estimated from data set zdanv using ARX.
Функция th2poly преобразует матрицу модели тета-формата в матрицы обобщенной (многомерной) линейной модели:
>> [A,B,C,D,K,lan,T]=th2poly(zpem)
A = 1.0000 -2.1441 1.5148 -0.3280 0.1302 -0.1081
B = 0 0.1322 -0.0698 -0.0946 0.0197 0.0668
C = 1.0000 -1.1083 -0.0108 0.1446 0.3438 -0.0371
D =1
K = 1
lan = 0.0069
T = 0.0800
Здесь параметр lan определяет интенсивность шума наблюдений.
Функция th2ss преобразует тета-модель в модель для переменных состояния. В качестве аргумента функции может выступать любая из рассмотренных ранее моделей, например darmax:
>> [A,B,C,D,K,x0]=th2ss(darmax)
A =
0.8733 1.0000
-0.1567 0
B =
0.1331
0.1028
C =
1 0
D =
0
K =1.0587
-0.1701
x0 =
0
0
Функция th2tf преобразует модель тета-формата многомерного объекта в вектор передаточных функций, связанных с выбранным входом:
>> [num,den]=th2tf(zn4s)
num = 0 0.1327 0.1566 0.0575
den = 1.0000 -0.3799 -0.2810 0.0749
Команда tf служит для представления передаточной функции в виде отношения:
>> zzn4s=tf(num,den,0.08)
Возвращает:
Transfer function:
0.1327 z^2 + 0.1566 z + 0.0575
------------------------------------
z^3 - 0.3799 z^2 - 0.281 z + 0.07493
Sampling time: 0.08
Функция thd2thc преобразует дискретную модель в непрерывную
Например: преобразуем дискретную модель тета-формата zn4s (модель переменных состояния в канонической форме при произвольном числе входов и выходов) в непрерывную модель и представим ее в виде передаточной функции. Для этого необходимо сначала выполнить функцию thd2thc преобразующую дискретную модель в непрерывную, затем выполнить функцию th2tf преобразующую модель тета-формата многомерного объекта в вектор передаточных функций, связанных с выбранным входом, а затем команду tf для представления передаточной функции в виде отношения:
>> sn4s=thd2thc(zn4s);
>> [num,den]=th2tf(sn4s);
>> sysn4s=tf(num,den)
Transfer function:
Transfer function:
-0.891 s^2 + 77.33 s + 746.9
---------------------------------
s^3 + 32.39 s^2 + 308.9 s + 891.7
Для обратного преобразования непрерывной модели в дискретную модель существует функция thc2thd.Функция th2zp рассчитывает нули, полюса и статические коэффициенты передачи (коэффициенты усиления) модели тета - формата zn4s многомерного объекта:
>> [zepo,k]=th2zp(zn4s)
zepo =1.0000 61.0000 21.0000 81.0000
-0.5898 + 0.2921i 0.2754 + 0.1282i -0.4946 0.6660
-0.5898 - 0.2921i 0.3531 0.6365 0.0651
Inf Inf 0.2380 0.2121
k =
1.0000
0.8376
0.0648
Информацию о нулях и полюсах модели zn4s можно получить с помощью команды
>> [zero,polus]=getzp(zepo)
zero =
-0.5898 + 0.2921i
-0.5898 - 0.2921i
polus =
-0.4946
0.6365
0.2380
С помощью команды zpplot можно построить график нулей и полюсов модели zn4s
>>zpplot(zpform(zepo))
На рис. 2. 10. представлен график нулей (обозначены кружком) и полюсов (обозначены крестиком) модели zn4s, который получен с помощью команды zpplot.
Рис. 2.10. Графики нулей и полюсов модели zn4s
Данные графика показывают, что модель является устойчивой: полюса модели находятся внутри окружности с радиусом, равным 1 и проходящей через точку с координатами (–1; j0).
2. 10. Проверка адекватности модели
Одним из важных этапов идентификации объектов автоматизации является проверка качества модели по выбранному критерию близости выхода модели и объекта, т.е проверка ее адекватности. В пакете System Identification Toolbox MATLAB в качестве такого критерия принята оценка адекватности модели fit, которая рассчитывается по формуле: fit = norm (yh – y)/
, где norm – норма вектора; yh и y – выходы модели и объекта соответственно; N – количество элементов массива данных.Для проверки адекватности полученных ранее моделей воспользуемся функцией:
>> compare(zdane,zn4s,zpem,zoe,zbj,darx,darmax).
где: zdane – выход объекта;
zn4s,zpem,zoe,zbj,darx,darmax – выходы моделей zn4s,zpem,zoe,zbj,darx,darmaxРис. 2. 11. Графики выходов объекта и моделей.
Результатом выполнения команды является вывод графика выходов объекта и построенных моделей (Рис. 2. 11). На графике цветными линиями представлены выходы полученных моделей и значения критерия адекватности, выраженного в процентах. Наилучшие показатели имеют модели darx, zn4s и zpem.
Для проверки адекватности модели zn4s воспользуемся функцией:
>>compare(zdane,zn4s)
Результат выполнения команды является вывод графика объекта на рис. 2. 12.
Рис. 2. 12. Графики выходов объекта и модели zn4s.В пакете System Identification Toolbox MATLAB имеется возможность прогнозировать ошибку моделирования при заданном входном воздействии u(t) и известной выходной координате объекта y(t). Оценивание производится методом прогноза ошибки Preictive Error Method, сокращенно PEM, который заключается в следующем. Пусть модель исследуемого объекта имеет вид так называемой обобщенной линейной модели
y(t) = W(z) u(t) + v(t),
где W(z) – дискретная передаточная функция любой из ранее рассмотренных моделей. При этом шум v(t) может быть представлен как
v(t) = H(z) e(t),
где e(z) – дискретный белый шум, который собственно и характеризует ошибку модели; H(z) – некоторый полином от z, приводящий дискретный белый шум к реальным помехам при измерении выходных параметров объекта.
Из данных выражений следует, что
e(t) = H-1(z) [y(t) – W(z) u(t)].
Функция resid вычисляет остаточную ошибку e для заданной модели, а также r – матрицу значений автокорреляционной функции процесса e(t) и значения взаимокорреляционой функции между остаточными ошибками e(t) и выходами объекта автоматизации y(t) вместе с соответствующими 99 %-ми доверительными коридорами. Кроме указанных значений выводятся графики данных функций. В качестве примера сравним остаточные ошибки и соответствующие корреляционные функции для полученных моделей darx и zbj, имеющих максимальную и минимальную оценки адекватности с помощью команд:
>> [e,r]=resid(zdan,darx)
>> [e1,r1]=resid(zdan,zbj)
Приведенные графики (рис. 2. 13 и 2 14) характеризуют равномерное распределение остаточных ошибок во всем диапазоне изменения интервалов времени τ. Причем значения остаточных ошибок
для модели darx практически в два раза больше, чем для модели zbj. Для вывода графиков необходимо выполнить команду resid(r).