Линеаризовавшая модель корпуса

Динамика корпуса и автопилот моделируются в Simulink. Смотрите Настройку Автопилота 2D Цикла (Control System Toolbox) для получения дополнительной информации об этой модели и для дополнительного проекта и настройки опций.

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:
---------------------------------

 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 секунда. Используйте формулировку смешанной чувствительности с весом 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.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 (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.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.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    
  7.02e-14 + 2.37e+01i    -2.96e-15       2.37e+01        -1.43e+13    
  7.02e-14 - 2.37e+01i    -2.96e-15       2.37e+01        -1.43e+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;

Путем комментария блока WorstPerturbation в и, можно симулировать ответ с или без этого возмущения. Результаты показывают ниже и подтверждают эффект дестабилизации усиления и изменения фазы, вычисленного diskmargin.

Рисунок 1: Номинальный ответ.

Рисунок 2: Ответ с дестабилизацией возмущения усиления/фазы.