В этом примере показано, как спроектировать устойчивый контроллер для автопилота 2D цикла, который управляет уровнем тангажа и вертикальным ускорением корпуса. Контроллер устойчив против усиления и изменений фазы многоканальной обратной связи.
Динамика корпуса и автопилот моделируются в Simulink. Смотрите Настройку Автопилота 2D Цикла для получения дополнительной информации об этой модели и для дополнительного проекта и настройки опций.
open_system('rct_airframe2')
Как в Настройке в качестве примера Автопилота 2D Цикла, обрежьте корпус для
и
. Условие для обрезки соответствует нулевому нормальному ускорению и подаче момента (
и
устойчивый). Используйте findop вычислить соответствующие условия работы. Затем линеаризуйте модель корпуса при этом условии для обрезки.
opspec = operspec('rct_airframe2'); % Specify trim condition % Xe,Ze: known, not steady opspec.States(1).Known = [1;1]; opspec.States(1).SteadyState = [0;0]; % u,w: known, w steady opspec.States(3).Known = [1 1]; opspec.States(3).SteadyState = [0 1]; % theta: known, not steady opspec.States(2).Known = 1; opspec.States(2).SteadyState = 0; % q: unknown, steady opspec.States(4).Known = 0; opspec.States(4).SteadyState = 1; % controller states unknown, not steady opspec.States(5).SteadyState = [0;0]; op = findop('rct_airframe2',opspec); G = linearize('rct_airframe2','rct_airframe2/Airframe Model',op); G.InputName = 'delta'; G.OutputName = {'az','q'};
Operating point search report:
---------------------------------
opreport =
Operating point search report for the Model rct_airframe2.
(Time-Varying Components Evaluated at time t=0)
Operating point specifications were successfully met.
States:
----------
<strong>Min</strong> <strong>x</strong> <strong>Max</strong> <strong>dxMin</strong> <strong>dx</strong> <strong>dxMax</strong>
<strong>_____</strong> <strong>__________</strong> <strong>_____</strong> <strong>_____</strong> <strong>___________</strong> <strong>_____</strong>
(1.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/Position
0 0 0 -Inf 984 Inf
-3048 -3048 -3048 -Inf 0 Inf
(2.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/Theta
0 0 0 -Inf -0.0097235 Inf
(3.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/U,w
984 984 984 -Inf 22.69 Inf
0 0 0 0 2.4587e-11 0
(4.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/q
-Inf -0.0097235 Inf 0 -1.7215e-16 0
(5.) rct_airframe2/MIMO Controller
-Inf 0.00065361 Inf -Inf -0.0089973 Inf
-Inf 3.7592e-20 Inf -Inf 0.030259 Inf
Inputs:
----------
<strong>Min</strong> <strong>u</strong> <strong>Max</strong>
<strong>____</strong> <strong>__________</strong> <strong>___</strong>
(1.) rct_airframe2/delta trim
-Inf 0.00043574 Inf
Outputs: None
----------
Спроектируйте контроллер H-бесконечности, который отвечает на ступенчатое изменение на вертикальном ускорении менее чем 1 секунда. Используйте формулировку смешанной чувствительности с весом lowpass wS это штрафует низкочастотную ошибку отслеживания и highpass вес wT это осуществляет соответствующий спад. Создайте увеличенный объект с azref,delta как входные параметры и отфильтрованный wS*e,wT*az,e,q как выходные параметры. (Для получения информации о формулировке смешанной чувствительности смотрите, что Цикл Смешанной Чувствительности Формирует.)
sb = sumblk('e = azref-az'); wS = makeweight(1e2,5,0.1); wS.u = 'e'; wS.y = 'we'; wT = makeweight(0.1,50,1e2); wT.u = 'az'; wT.y = 'waz'; Paug = connect(G,wS,wT,sb,{'azref','delta'},{'we','waz','e','q'});
Вычислите оптимальный контроллер H-бесконечности, использующий hinfsyn. В этой настройке существует два измерения e,q и одно управление delta.
[Knom,~,gam] = hinfsyn(Paug,2,1); gam
gam =
0.5928
Проверьте коэффициент усиления разомкнутого контура по wS,wT.
f = figure(); sigma(Knom*G,wS,'r--',1/wT,'g--'), grid legend('open-loop gain','> wS at low freq','< 1/wT at high freq')
Постройте ответ с обратной связью.
CL = feedback(G*Knom,diag([1 -1])); step(CL(:,1)), grid
Вычислите дисковые поля во входе объекта и выходных параметрах. То, что az цикл использует отрицательную обратную связь в то время как q цикл использует положительную обратную связь, поэтому измените знак усиления цикла q цикл перед использованием diskmargin.
loopsgn = diag([1 -1]);
Исследуйте поля во входе объекта.
DM = diskmargin(Knom*loopsgn*G)
DM =
struct with fields:
GainMargin: [0.3918 2.5524]
PhaseMargin: [-47.2105 47.2105]
DiskMargin: 0.8740
LowerBound: 0.8740
UpperBound: 0.8740
Frequency: 28.8842
WorstPerturbation: [1x1 ss]
Для полей на объекте выходные параметры используйте многоконтурное поле с учетом одновременных, независимых изменений обоих выходных каналов.
[~,MM] = diskmargin(loopsgn*G*Knom)
MM =
struct with fields:
GainMargin: [0.4994 2.0025]
PhaseMargin: [-36.9262 36.9262]
DiskMargin: 0.6678
LowerBound: 0.6678
UpperBound: 0.6691
Frequency: 23.6845
WorstPerturbation: [2x2 ss]
Наконец, вычислите поля против одновременных изменений во входе объекта и выходных параметров.
MMIO = diskmargin(loopsgn*G,Knom)
MMIO =
struct with fields:
GainMargin: [0.6866 1.4565]
PhaseMargin: [-21.0565 21.0565]
DiskMargin: 0.3717
LowerBound: 0.3717
UpperBound: 0.3725
Frequency: 24.8030
WorstPerturbation: [1x1 struct]
Находящиеся на диске запасы по амплитуде и фазе превышают 2 (6 дБ) и 35 градусов во входе объекта и объекте выходные параметры. Кроме того, MMIO указывает, что эта обратная связь может противостоять изменениям усиления фактором 1.45 или изменениям фазы 21 градуса, влияющего на все вводы и выводы объекта целиком.
Затем исследуйте удар усиления и неопределенности фазы на эффективности. Используйте umargin блок системы управления, чтобы смоделировать усиление изменяет фактор 1,4 (3 дБ) и фазовый переход 20 градусов в области каждого канала обратной связи. Используйте getDGM соответствовать диску неопределенности к этим суммам усиления и фазового перехода.
GM = 1.4; PM = 20; DGM = getDGM(GM,PM,'balanced'); ue = umargin('e',DGM); uq = umargin('q',DGM); Gunc = blkdiag(ue,uq)*G; Gunc.OutputName = {'az','q'};
Восстановите модель с обратной связью и произведите усиление и неопределенность фазы, чтобы измерить удар на переходной процесс.
CLunc = feedback(Gunc*Knom,loopsgn); step(CLunc(:,1),3) grid
График показывает значительную изменчивость в эффективности с большим перерегулированием и установившейся ошибкой для некоторых комбинаций изменений фазы и усиления.
Восстановите синтез контроллера с учетом усиления и изменения фазы с помощью musyn. Этот синтез оптимизирует эффективность однородно для заданной области неопределенности фазы и усиления.
Punc = connect(Gunc,wS,wT,sb,{'azref','delta'},{'we','waz','e','q'});
[Kr,gam] = musyn(Punc,2,1);
gam
D-K ITERATION SUMMARY:
-----------------------------------------------------------------
Robust performance Fit order
-----------------------------------------------------------------
Iter K Step Peak MU D Fit D
1 51.06 26.33 26.62 4
2 7.028 5.24 5.303 8
3 1.681 1.156 1.157 4
4 0.9705 0.9702 0.9822 10
5 0.962 0.9619 0.9622 10
6 0.9625 0.9623 0.9629 10
Best achieved robust performance: 0.962
gam =
0.9619
Сравните произведенные переходные процессы для номинальных и устойчивых контроллеров. Устойчивый проект уменьшает и перерегулирование и установившиеся ошибки и дает более сопоставимую эффективность через смоделированную область значений неопределенности фазы и усиления.
CLr = feedback(Gunc*Kr,loopsgn);
rng(0) % for reproducibility
step(CLunc(:,1),3)
hold
rng(0)
step(CLr(:,1),3)
grid
Current plot held
Устойчивый контроллер достигает этой эффективности путем увеличения (номинальных) дисковых полей на объекте выход и, до меньшей степени, дискового поля ввода-вывода. Например, сравните находящиеся на диске поля на объекте выходные параметры для номинальных и устойчивых проектов. Используйте diskmarginplot видеть изменения полей с частотой.
Lnom = loopsgn*G*Knom; Lrob = loopsgn*G*Kr; clf diskmarginplot(Lnom,Lrob) title('Disk margins at plant outputs') grid legend('Nominal','Robust')
Исследуйте поля против изменений одновременные изменения при вводах и выводах объекта с устойчивым контроллером.
MMIO = diskmargin(loopsgn*G,Kr)
MMIO =
struct with fields:
GainMargin: [0.6492 1.5404]
PhaseMargin: [-24.0166 24.0166]
DiskMargin: 0.4254
LowerBound: 0.4254
UpperBound: 0.4263
Frequency: 19.6041
WorstPerturbation: [1x1 struct]
Вспомните, что с номинальным контроллером, обратная связь могла противостоять изменениям усиления фактора 1,45 или изменениям фазы 21 градуса, влияющего на все вводы и выводы объекта целиком. С устойчивым контроллером те поля увеличиваются приблизительно до 1,54 и 24 градусов.
Просмотрите область значений одновременного усиления и изменений фазы, которые устойчивый проект может выдержать при всех вводах и выводах объекта.
diskmarginplot(MMIO.GainMargin)
diskmargin возвращает самое маленькое изменение в усилении и фазе, которая дестабилизирует обратную связь. Это может быть проницательно, чтобы ввести это возмущение в нелинейной симуляции, чтобы далее анализировать последствия для действительной системы. Например, вычислите возмущение дестабилизации на объекте выходные параметры для номинального контроллера.
[~,MM] = diskmargin(Lnom);
WP = MM.WorstPerturbation;
bode(WP)
title('Smallest destabilizing perturbation')
Худшее возмущение является диагональным, динамическим возмущением, которое умножает ответ разомкнутого контура на объекте выходные параметры. С этим возмущением система с обратной связью становится нестабильной с незатухающим полюсом на частоте w0 = MM.Frequency.
damp(feedback(WP*G*Knom,loopsgn))
Pole Damping Frequency Time Constant
(rad/seconds) (seconds)
-1.88e-03 1.00e+00 1.88e-03 5.33e+02
-2.55e-02 1.00e+00 2.55e-02 3.92e+01
-3.20e-02 1.00e+00 3.20e-02 3.12e+01
-5.16e-02 1.00e+00 5.16e-02 1.94e+01
-1.25e-01 1.00e+00 1.25e-01 7.98e+00
-3.85e+00 + 5.04e+00i 6.07e-01 6.34e+00 2.60e-01
-3.85e+00 - 5.04e+00i 6.07e-01 6.34e+00 2.60e-01
-8.38e+00 + 1.17e+01i 5.81e-01 1.44e+01 1.19e-01
-8.38e+00 - 1.17e+01i 5.81e-01 1.44e+01 1.19e-01
6.48e-14 + 2.37e+01i -2.74e-15 2.37e+01 -1.54e+13
6.48e-14 - 2.37e+01i -2.74e-15 2.37e+01 -1.54e+13
-2.95e+01 1.00e+00 2.95e+01 3.39e-02
-3.33e+01 1.00e+00 3.33e+01 3.00e-02
-6.04e+01 + 2.28e+01i 9.36e-01 6.46e+01 1.66e-02
-6.04e+01 - 2.28e+01i 9.36e-01 6.46e+01 1.66e-02
-3.86e+01 + 5.36e+01i 5.85e-01 6.61e+01 2.59e-02
-3.86e+01 - 5.36e+01i 5.85e-01 6.61e+01 2.59e-02
-1.10e+03 1.00e+00 1.10e+03 9.05e-04
w0 = MM.Frequency
w0 = 23.6845
Обратите внимание на то, что усиление и изменения фазы, вызванные этим возмущением, лежат на контуре области значений безопасных изменений усиления/фазы, вычисленных diskmargin. Чтобы подтвердить этот результат, вычислите частотную характеристику худшего возмущения на частоте w0, преобразуйте его в усиление и изменение фазы, и визуализируйте его наряду с областью значений безопасного усиления и изменений фазы, представленных многоконтурным дисковым полем.
DELTA = freqresp(WP,w0); clf diskmarginplot(MM.GainMargin) title('Range of stable gain and phase variations') hold on plot(20*log10(abs(DELTA(1,1))),abs( angle(DELTA(1,1))*180/pi),'ro') plot(20*log10(abs(DELTA(2,2))),abs( angle(DELTA(2,2))*180/pi),'ro')
Чтобы симулировать эффект этого худшего возмущения на полной нелинейной модели в Simulink, вставьте его как блок перед блоком контроллера, как сделано в модифицированной модели rct_airframeWP.
close(f)
open_system('rct_airframeWP')
Здесь блок MIMO Controller установлен в номинальный контроллер Knom. Чтобы симулировать нелинейный ответ с этим контроллером, вычислите отклонение для обрезки и q начальное значение от условий работы op.
delta_trim = op.Inputs.u + [1.5 0]*op.States(5).x; q_ini = op.States(4).x;
Путем комментария блока WorstPerturbation в и, можно симулировать ответ с или без этого возмущения. Результаты показывают ниже и подтверждают эффект дестабилизации усиления и изменения фазы, вычисленного diskmargin.
Рисунок 1: Номинальный ответ.
Рисунок 2: Ответ с дестабилизацией возмущения усиления/фазы.
umargin | diskmargin | diskmarginplot | musyn