Линейное приближение сложных систем по идентификации
Этот пример показывает, как получить линейные приближения комплексной, нелинейной системы с помощью линейной идентификации модели. Подход основан на выборе входного сигнала, который возбуждает систему. Линейное приближение получается путем подбора линейной модели к симулированному отклику нелинейной модели для выбранного входного сигнала.
Этот пример использует Simulink ®, Control System Toolbox™ и Simulink Control Design™.
Введение
Во многих ситуациях линейная модель получается путем упрощения более сложной нелинейной системы при определенных локальных условиях. Для примера модель динамики самолета с высокой точностью может быть описана подробной моделью Simulink. Общим подходом, принятым для ускорения симуляции таких систем, изучения их локального поведения относительно рабочей точки или проектирования компенсаторов, является их линеаризация. Если мы выполним линеаризацию исходной модели аналитически относительно рабочей точки, это, в целом, даст модель порядка, такого высокого или близкого к количеству состояний в исходной модели. Этот порядок может быть излишне высоким для класса входов, который предполагается использовать для анализа или разработки системы управления. Следовательно, мы можем рассмотреть альтернативный подход, ориентированный на сбор входно-выходных данных из симуляции системы и использование его, чтобы вывести линейную модель справедливого порядка.
Аналитическая линеаризация F14 модели
Рассмотрим модель 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 улучшить его подгонку данным. Для этой оценки мы выбираем метод поиска «Levenberg-Marquardt» и изменяем максимальное количество допустимых итераций на 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')