Этот пример показывает, как спроектировать робастный контроллер для двухконтурного автопилота, который управляет скоростью тангажа и вертикальным ускорением планера. Контроллер является устойчивым к изменениям усиления и фазы в многоканальном цикле обратной связи.
Динамика планера и автопилота моделируются в Simulink. См. раздел Настройка автопилота с двумя циклами для получения дополнительной информации об этой модели, а также для получения дополнительной информации о проекте и опциях.
open_system('rct_airframe2')
Как и в примере Настройка автопилота с двумя контурами, обрезайте планер для
и. 
Условие обрезки соответствует нулю нормального ускорения и крутящего момента (
и
устойчивого). Использование 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:
---------------------------------
Operating point search report for the Model rct_airframe2.
(Time-Varying Components Evaluated at time t=0)
Operating point specifications were successfully met.
States:
----------
(1.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/Position
x: 0 dx: 984
x: -3.05e+03 dx: 0
(2.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/Theta
x: 0 dx: -0.00972
(3.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/U,w
x: 984 dx: 22.7
x: 0 dx: 2.46e-11 (0)
(4.) rct_airframe2/Airframe Model/Aerodynamics & Equations of Motion/ Equations of Motion (Body Axes)/q
x: -0.00972 dx: -1.72e-16 (0)
(5.) rct_airframe2/MIMO Controller
x: 0.000654 dx: -0.009
x: 4.13e-19 dx: 0.0303
Inputs:
----------
(1.) rct_airframe2/delta trim
u: 0.000436 [-Inf Inf]
Outputs: None
----------
Разработайте контроллер H-бесконечности, который реагирует на изменение шага вертикального ускорения в течение 1 секунды. Используйте композицию со смешанной чувствительностью с низким весом wS что штрафует низкочастотную ошибку отслеживания и высокий вес wT что обеспечивает надлежащее развертывание. Создайте дополненный объект с azref,delta как входы и отфильтрованные wS*e,wT*az,e,q как выходы. (Для получения информации о составе со смешанной чувствительностью смотрите Mixed-Sensitivity Loop Shaping.)
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]
Для полей на выходах объекта используйте multiloop margin для расчета одновременных независимых изменений в обоих выходных каналах.
[~,MM] = diskmargin(loopsgn*G*Knom)
MM =
struct with fields:
GainMargin: [0.4994 2.0025]
PhaseMargin: [-36.9261 36.9261]
DiskMargin: 0.6678
LowerBound: 0.6678
UpperBound: 0.6691
Frequency: 23.6855
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 (6dB) и 35 степени на входе объекта и выходах объекта. Более того, MMIO указывает, что этот цикл обратной связи может выдерживать изменения коэффициента усиления в 1,45 раза или изменения фазы в 21 степень, влияющие на все входы и выходы объекта сразу.
Затем исследуйте влияние неопределенности усиления и фазы на эффективность. Используйте umargin управляйте проектирующим блоком, чтобы смоделировать коэффициент изменения усиления 1,4 (3dB) и изменение фазы 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.029 5.24 5.304 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.0167 24.0167]
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
7.82e-14 + 2.37e+01i -3.30e-15 2.37e+01 -1.28e+13
7.82e-14 - 2.37e+01i -3.30e-15 2.37e+01 -1.28e+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.6855
Обратите внимание, что изменения усиления и фазы, вызванные этим возмущением, лежат на контуре области значений безопасного усиления/фазы, вычисленных 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;
Комментируя блок наихудшего возмущения, можно симулировать ответ с этим возмущением или без него. Результаты показаны ниже и подтверждают дестабилизирующий эффект изменения усиления и фазы, вычисленные diskmargin.
Фигура 1: Номинальная характеристика.
Фигура 2: Реакция с дестабилизирующим возмущением усиления/фазы.
diskmargin | diskmarginplot | musyn | umargin