В этом примере показано, как спроектировать устойчивый контроллер для автопилота 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