В этом примере показано, как сконструировать надежный контроллер для двухконтурного автопилота, который управляет скоростью тангажа и вертикальным ускорением планера. Контроллер надежен против усиления и фазовых изменений в многоканальном контуре обратной связи.
Динамика планера и автопилот смоделированы в 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 в качестве выходных данных. (Сведения о формулировке смешанной чувствительности см. в разделе Формирование петель смешанной чувствительности.)
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.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 устанавливается в номинальный контроллер Knom. Чтобы смоделировать нелинейный отклик с помощью этого контроллера, вычислите отклонение обрезки и q начальное значение из рабочего состояния op.
delta_trim = op.Inputs.u + [1.5 0]*op.States(5).x; q_ini = op.States(4).x;
Комментируя входящий и исходящий блок Purturbation, можно смоделировать ответ с этим возмущением или без него. Результаты показаны ниже и подтверждают дестабилизирующий эффект усиления и изменения фазы, вычисленный diskmargin.
Рис. 1: Номинальная реакция.
Рис. 2: Реакция с дестабилизирующим усилением/фазовым возмущением.
diskmargin | diskmarginplot | musyn | umargin