Линеаризированная модель планера

Динамика планера и автопилота моделируются в 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: Реакция с дестабилизирующим возмущением усиления/фазы.