В этом примере показано, как получить линейные аппроксимации комплексной нелинейной системы посредством идентификации линейной модели. Подход основан на выборе входного сигнала, возбуждающего систему. Линейную аппроксимацию получают подгонкой линейной модели к моделируемому отклику нелинейной модели для выбранного входного сигнала.
В этом примере используются Simulink ®, Control System Toolbox™ и Simulink Control Design™.
Во многих ситуациях линейная модель получается упрощением более сложной нелинейной системы при определенных локальных условиях. Например, модель с высокой точностью динамики самолета может быть описана детальной моделью Simulink. Общим подходом, принятым для ускорения моделирования таких систем, изучения их локального поведения относительно рабочей точки или проектирования компенсаторов, является линеаризация их. Если выполнить линеаризацию исходной модели аналитически относительно рабочей точки, это, в общем случае, даст модель порядка, равного или близкого к числу состояний в исходной модели. Этот порядок может быть излишне высоким для класса вводов, который предполагается использовать для анализа или проектирования системы управления. Поэтому мы можем рассмотреть альтернативный подход, ориентированный на сбор данных ввода-вывода из моделирования системы и использование его для получения линейной модели только правильного порядка.
Рассмотрим модель F14. Это уже линейная модель, но содержит блоки производных и источники возмущений, которые могут повлиять на характер ее выхода. Мы можем «линеаризовать» его между одним входным портом и двумя выходными портами следующим образом:
open_system('idF14Model') IO = getlinio('idF14Model') syslin = linearize('idF14Model',IO)
3x1 vector of Linearization IOs:
--------------------------
1. Linearization input perturbation located at the following signal:
- Block: idF14Model/Pilot
- Port: 1
- Signal Name: Stick Input
2. Linearization output measurement located at the following signal:
- Block: idF14Model/Gain5
- Port: 1
- Signal Name: Angle of Attack
3. Linearization output measurement located at the following signal:
- Block: idF14Model/Pilot G force
- Port: 1
- Signal Name: Pilot G force
syslin =
A =
Transfer Fcn Derivative Transfer Fcn Derivative1
Transfer Fcn -0.6385 0 689.4 0
Derivative 1 -1e+05 0 0
Transfer Fcn -0.00592 0 -0.6571 0
Derivative1 0 0 1 -1e+05
Actuator Mod 0 0 1.424 0
Alpha-sensor 0.001451 0 0 0
Stick Prefil 0 0 0 0
Pitch Rate L 0 0 1 0
Proportional 0 0 -0.8156 0
Actuator Mod Alpha-sensor Stick Prefil Pitch Rate L
Transfer Fcn -1280 0 0 0
Derivative 0 0 0 0
Transfer Fcn -137.7 0 0 0
Derivative1 0 0 0 0
Actuator Mod -20 2.986 -39.32 -1.67
Alpha-sensor 0 -2.526 0 0
Stick Prefil 0 0 -22.52 0
Pitch Rate L 0 0 0 -4.144
Proportional 0 -1.71 22.52 0.9567
Proportional
Transfer Fcn 0
Derivative 0
Transfer Fcn 0
Derivative1 0
Actuator Mod -3.864
Alpha-sensor 0
Stick Prefil 0
Pitch Rate L 0
Proportional 0
B =
Stick Input
Transfer Fcn 0
Derivative 0
Transfer Fcn 0
Derivative1 0
Actuator Mod 0
Alpha-sensor 0
Stick Prefil 1
Pitch Rate L 0
Proportional 0
C =
Transfer Fcn Derivative Transfer Fcn Derivative1
Angle of Att 0.001451 0 0 0
Pilot G forc -3106 3.106e+08 7.083e+04 -7.081e+09
Actuator Mod Alpha-sensor Stick Prefil Pitch Rate L
Angle of Att 0 0 0 0
Pilot G forc 0 0 0 0
Proportional
Angle of Att 0
Pilot G forc 0
D =
Stick Input
Angle of Att 0
Pilot G forc 0
Continuous-time state-space model.
syslin - модель с 2 выходами, одним входом и 9 состояниями. Это происходит потому, что путь линеаризации от входа «Stick Input» к двум выходам имеет 9 состояний в исходной системе. Мы можем проверить это, используя operpoint:
operpoint('idF14Model')
Operating point for the Model idF14Model.
(Time-Varying Components Evaluated at time t=0)
States:
----------
(1.) idF14Model/Actuator Model
x: 0
(2.) idF14Model/Aircraft Dynamics Model/Transfer Fcn.1
x: 0
(3.) idF14Model/Aircraft Dynamics Model/Transfer Fcn.2
x: 0
(4.) idF14Model/Controller/Alpha-sensor Low-pass Filter
x: 0
(5.) idF14Model/Controller/Pitch Rate Lead Filter
x: 0
(6.) idF14Model/Controller/Proportional plus integral compensator
x: 0
(7.) idF14Model/Controller/Stick Prefilter
x: 0
(8.) idF14Model/Dryden Wind Gust Models/Q-gust model
x: 0
(9.) idF14Model/Dryden Wind Gust Models/W-gust model
x: 0
x: 0
Inputs: None
----------
Может ли этот порядок быть уменьшен при сохранении точности ответов для выбранного входного сигнала прямоугольной формы («Stick input»)?
Смоделировать модель и записать входную квадратную волну (u) и выходы «Угол атаки» (y1) и «Сила пилота G» (y2) в течение 0:30 секунд. Эти данные после интерполяции на равномерно разнесенный вектор времени (время выборки 0,0444 секунды) сохраняются в файле «idF14SimData.mat».
load idF14SimData Z = iddata([y1, y2],u,Ts,'Tstart',0); Z.InputName = 'Stick input'; Z.InputUnit = 'rad/s'; Z.OutputName = {'Angle of attack', 'Pilot G force'}; Z.OutputUnit = {'rad', 'g'}; t = Z.SamplingInstants; subplot(311) plot(t,Z.y(:,1)), ylabel('Angle of attack (rad)') title('Logged Input-Output Data') subplot(312) plot(t,Z.y(:,2)), ylabel('Pilot G force (g)') subplot(313) plot(t,Z.u), ylabel('Stick input (rad/s)') axis([0 30 -1.2 1.2]) xlabel('Time (seconds)')
Оценка моделей состояния пространства порядков от 2 до 4 с помощью ssest команда. Мы конфигурируем оценку для использования фокуса «моделирования» и выбираем, чтобы не оценивать компонент возмущения модели.
opt = ssestOptions('Focus','simulation'); syslin2 = ssest(Z, 2, 'DisturbanceModel', 'none', opt); syslin3 = ssest(Z, 3, 'DisturbanceModel', 'none', opt); syslin4 = ssest(Z, 4, 'DisturbanceModel', 'none', opt);
Сравнение производительности линеаризованной модели syslin и три идентифицированные модели к данным. Обратите внимание, что syslin является моделью SS, пока syslin2, syslin3 и syslin4 являются моделями IDSS.
syslin.InputName = Z.InputName;
syslin.OutputName = Z.OutputName; % reconcile names to facilitate comparison
clf
compare(Z, syslin, syslin2, syslin3, syslin4)
На графике видно, что модель 3-го порядка (syslin3) достаточно хорошо работает как линейная аппроксимация динамики самолета о его стандартных (t = 0) условиях эксплуатации. Он подходит к данным несколько лучше, чем полученный аналитической линеаризацией (syslin). Если исходный idF14Model является линейным, почему его линеаризация не дает результата, syslin, произвести 100% подгонку к данным? Для этого есть две причины:
На измеряемые выходные сигналы влияют порывы ветра, что означает, что регистрируемые выходные сигналы не являются просто функцией входного сигнала Stick. Есть нарушения, влияющие на него.
Блок «Pilot G force» использует производные блоки, линеаризация которых зависит от значения постоянной времени «c». «c» должен быть маленьким (мы используем 1e-5), но не нулевым. Ненулевое значение для c вносит аппроксимационную ошибку в линеаризацию.
Рассмотрим параметры модели. syslin3 что, по-видимому, довольно хорошо отражает ответы:
syslin3
syslin3 =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2 x3
x1 -1.006 -2.029 0.5842
x2 8.284 -19.39 -5.611
x3 2.784 -12.63 -6.956
B =
Stick input
x1 0.2614
x2 5.512
x3 3.606
C =
x1 x2 x3
Angle of att -8.841 -0.5347 1.402
Pilot G forc -86.42 -15.85 66.12
D =
Stick input
Angle of att 0
Pilot G forc 0
K =
Angle of att Pilot G forc
x1 0 0
x2 0 0
x3 0 0
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: none
Disturbance component: none
Number of free coefficients: 18
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using SSEST on time domain data "Z".
Fit to estimation data: [98.4;97.02]%
FPE: 2.367e-05, MSE: 0.1103
Мы также можем использовать подход уменьшения порядка линеаризованной модели syslin и уточнение параметров уменьшенной модели, чтобы наилучшим образом соответствовать данным Z. Чтобы определить хорошее значение для уменьшенного порядка, мы используем hsvd:
[S, BalData] = hsvd(syslin); clf; bar(S)
Гистограмма показывает, что сингулярные значения довольно малы для состояний 4 и выше. Следовательно, 3-й порядок может быть оптимальным для уменьшения.
sysr = balred(syslin,3,BalData)
opt2 = bodeoptions; opt2.PhaseMatching = 'on';
clf; bodeplot(sysr,syslin,opt2)
sysr =
A =
x1 x2 x3
x1 -2.854 -7.61 -54.04
x2 -0.9714 2.341 9.123
x3 0.6979 -7.203 -24.08
B =
Stick input
x1 -137.7
x2 -869
x3 -506.7
C =
x1 x2 x3
Angle of att -0.0005063 -0.0008826 -0.001016
Pilot G forc -0.005926 -0.04692 -0.1646
D =
Stick input
Angle of att -0.03784
Pilot G forc -1.617
Continuous-time state-space model.
График бода показывает хорошую верность до 10 рад/с. sysr может эмулировать ответы так же хорошо, как и исходная модель с 9 состояниями, как показано compare график:
compare(Z, sysr, syslin)
Уточним параметры sysr улучшить его соответствие данным. Для этой оценки выбираем метод поиска «Левенберга-Марквардта» и изменяем максимальное количество допустимых итераций на 10. Этот выбор был сделан на основе проб и ошибок. Мы также включаем отображение прогресса оценки.
opt.Display = 'on'; opt.SearchMethod = 'lm'; opt.SearchOptions.MaxIterations = 10; sysr2 = ssest(Z, sysr, opt) compare(Z, sysr2)
sysr2 =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2 x3
x1 -4.048 -7.681 -54.01
x2 -0.4844 1.549 8.895
x3 -0.2398 -6.777 -25.78
B =
Stick input
x1 -137.7
x2 -869
x3 -506.7
C =
x1 x2 x3
Angle of att -0.0003361 -0.0004964 0.00215
Pilot G forc -0.01191 -0.03599 -0.1434
D =
Stick input
Angle of att 0.003022
Pilot G forc 0.6438
K =
Angle of att Pilot G forc
x1 0 0
x2 0 0
x3 0 0
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: yes
Disturbance component: none
Number of free coefficients: 20
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using SSEST on time domain data "Z".
Fit to estimation data: [98.78;97.03]%
FPE: 1.434e-05, MSE: 0.1097
Уточненная модель sysr2 вполне подходит под отклик модели F14 (около 99% для первого выпуска и около 97% для второго).
Мы показали альтернативный подход к аналитической линеаризации для получения линейных аппроксимаций сложных систем. Результаты получаются для конкретного входного сигнала и, строго говоря, применимы только к этому входу. Чтобы улучшить применимость результатов к различным профилям ввода, мы могли бы выполнить несколько моделирований с использованием различных типов входных данных. Затем мы могли бы объединить результирующие наборы данных в один набор данных из нескольких экспериментов (см. iddata/merge) и использовать его для оценок. В этом примере мы использовали сложную, но линейную систему для удобства. Реальная польза этого подхода будет видна для нелинейных систем.
Мы также показали подход к уменьшению порядка линейной системы при сохранении уменьшенной модели верной моделируемой реакции исходной модели Simulink. Роль сокращенной модели sysr должен был предоставить первоначальное предположение для расчетной модели sysr2. Подход также подчеркивает тот факт, что любая линейная система, включая системы другого класса, может использоваться в качестве начальной модели для оценок.
bdclose('idF14Model')