Выполните управление и статический анализ устойчивости на линеаризированных самолетах

Этот пример использует:

В этом примере показано, как преобразовать самолет в линейную модель пространства состояний инварианта времени (LTI) для линейного анализа.

Этот пример описывает:

Импорт и заполнение данных из файла DATCOM.
Построение самолета из данных DATCOM.
Вычисление статической устойчивости самолета.
Линеаризация самолета вокруг начального состояния.
Валидация статического анализа устойчивости с динамической характеристикой.
Изоляция передаточной функции «лифт на тангаж» и разработка контроллера обратной связи для лифта.

Определение самолета и состояния

Этот пример использует файл DATCOM, созданный для самолета Sky Hogg.

Во-первых, импортируйте выходной файл DATCOM с помощью datcomimport.

allData = datcomimport('astSkyHoggDatcom.out', false, 0);
skyHoggData = allData{1}

skyHoggData = struct with fields:
          case: 'SKYHOGG BODY-WING-HORIZONTAL TAIL-VERTICAL TAIL CONFIG'
          mach: [0.1000 0.2000 0.3000 0.3500]
           alt: [1000 3000 5000 7000 9000 11000 13000 15000]
         alpha: [-16 -12 -8 -4 -2 0 2 4 8 12]
         nmach: 4
          nalt: 8
        nalpha: 10
         rnnub: []
        hypers: 0
          loop: 2
          sref: 225.8000
          cbar: 5.7500
         blref: 41.1500
           dim: 'ft'
         deriv: 'deg'
        stmach: 0.6000
        tsmach: 1.4000
          save: 0
         stype: []
          trim: 0
          damp: 1
         build: 1
          part: 0
       highsym: 1
       highasy: 0
       highcon: 0
          tjet: 0
        hypeff: 0
            lb: 0
           pwr: 0
          grnd: 0
         wsspn: 18.7000
         hsspn: 5.7000
        ndelta: 5
         delta: [-20 -10 0 10 20]
        deltal: []
        deltar: []
           ngh: 0
        grndht: []
        config: [1x1 struct]
       version: 1976
            cd: [10x4x8 double]
            cl: [10x4x8 double]
            cm: [10x4x8 double]
            cn: [10x4x8 double]
            ca: [10x4x8 double]
           xcp: [10x4x8 double]
           cma: [10x4x8 double]
           cyb: [10x4x8 double]
           cnb: [10x4x8 double]
           clb: [10x4x8 double]
           cla: [10x4x8 double]
         qqinf: [10x4x8 double]
           eps: [10x4x8 double]
      depsdalp: [10x4x8 double]
           clq: [10x4x8 double]
           cmq: [10x4x8 double]
          clad: [10x4x8 double]
          cmad: [10x4x8 double]
           clp: [10x4x8 double]
           cyp: [10x4x8 double]
           cnp: [10x4x8 double]
           cnr: [10x4x8 double]
           clr: [10x4x8 double]
       dcl_sym: [5x4x8 double]
       dcm_sym: [5x4x8 double]
    dclmax_sym: [5x4x8 double]
    dcdmin_sym: [5x4x8 double]
      clad_sym: [5x4x8 double]
       cha_sym: [5x4x8 double]
       chd_sym: [5x4x8 double]
      dcdi_sym: [10x5x4x8 double]

Затем подготовьте интерполяционные таблицы DATCOM.

Интерполяционные таблицы DATCOM могут иметь отсутствующие значения из-за того, что таблицы заполняют только одно значение для всего столбца.

Эти отсутствующие данные представлены в интерполяционных таблицах как значение 99999 и могут быть заполнены с помощью «предыдущего» метода заполнения.

В этом примере, $C_{y β}$ , $C_{n β}$ , $C_{Lq}$ , и $C_{mq}$ иметь отсутствующие данные.

skyHoggData.cyb = fillmissing(skyHoggData.cyb, "previous", "MissingLocations", skyHoggData.cyb == 99999);
skyHoggData.cnb = fillmissing(skyHoggData.cnb, "previous", "MissingLocations", skyHoggData.cnb == 99999);
skyHoggData.clq = fillmissing(skyHoggData.clq, "previous", "MissingLocations", skyHoggData.clq == 99999);
skyHoggData.cmq = fillmissing(skyHoggData.cmq, "previous", "MissingLocations", skyHoggData.cmq == 99999);

При заполненных недостающих данных могут быть построены самолеты.

Во-первых, самолет готовят с желаемым именем самолета.

Опционально имя самолета может быть извлечено из поля «case» на struct DATCOM путем передачи пустого объекта самолета.

skyHogg = Aero.FixedWing();
skyHogg.Properties.Name = "Sky_Hogg";
skyHogg.DegreesOfFreedom = "3DOF";
[skyHogg, cruiseState] = datcomToFixedWing(skyHogg, skyHoggData);

datcomToFixedWing преобразует все совместимые данные из struct datcom в объект fixwing и его состояние. Однако возвращенное состояние все еще нуждается в обработке, чтобы получить желаемые начальные условия самолета.

В этом примере окружение, масса, инерция, воздушная скорость и центр давления нуждаются в регулировании.

h = 2000;
[T, a, p, rho] = atmoscoesa(convlength(h, "ft", "m"));
cruiseState.AltitudeMSL = h;
cruiseState.Environment.Temperature = convtemp(T,"K", "R");
cruiseState.Environment.SpeedOfSound = convvel(a, "m/s", "ft/s");
cruiseState.Environment.Pressure = convpres(p, "pa", "psf");
cruiseState.Environment.Density = convdensity(rho, "kg/m^3", "slug/ft^3");

cruiseState.U = 169.42;
cruiseState.Mass = 1299.214;
cruiseState.Inertia.Variables = [5787.969 0 117.64;0 6928.93 0;-117.64 0 11578.329];
cruiseState.CenterOfPressure = [0.183, 0, 0];

Вычисление статической устойчивости

Выполнение статического анализа устойчивости помогает определить динамическую устойчивость системы без вычисления динамической системы отклика.

stability = staticStability(skyHogg, cruiseState)

stability=6×8 table
             U            V           W         Alpha        Beta           P           Q            R    
          ________    _________    ________    ________    _________    _________    ________    _________

    FX    "Stable"    ""           ""          ""          ""           ""           ""          ""       
    FY    ""          "Neutral"    ""          ""          ""           ""           ""          ""       
    FZ    ""          ""           "Stable"    ""          ""           ""           ""          ""       
    L     ""          ""           ""          ""          "Neutral"    "Neutral"    ""          ""       
    M     "Stable"    ""           ""          "Stable"    ""           ""           "Stable"    ""       
    N     ""          ""           ""          ""          "Neutral"    ""           ""          "Neutral"

При статически устойчивых силах и моментах, с возмущениями к прямым и вертикальным скоростям и возмущениями к углу атаки, динамическая устойчивость системы имеет тенденцию к колебанию установившегося состояния при возмущении передней и вертикальной скоростей.

Чтобы проверить это поведение, используйте Control System Toolbox™.

Линеаризация самолета

Чтобы использовать инструменты в Control System Toolbox™, линеаризируйте самолет вокруг состояния.

Это делается с помощью метода линеаризации с таким же круизным состоянием, как и раньше.

linSys = linearize(skyHogg, cruiseState)

linSys =
 
  A = 
                  XN          XD           U           W           Q
   XN              0           0           1           0           0
   XD              0           0           0           1           0
   U               0   1.319e-06   -0.001714  -0.0007608           0
   W               0           0   -0.002705     -0.3319       2.557
   Q               0           0     0.02379      -2.495      -26.41
   Theta           0           0           0           0           1
 
               Theta
   XN     -2.586e-07
   XD         -2.957
   U         -0.5617
   W      -4.867e-08
   Q               0
   Theta           0
 
  B = 
               Delta
   XN              0
   XD              0
   U      -0.0004314
   W        -0.02084
   Q          -2.321
   Theta           0
 
  C = 
             XN     XD      U      W      Q  Theta
   XN         1      0      0      0      0      0
   XD         0      1      0      0      0      0
   U          0      0      1      0      0      0
   W          0      0      0      1      0      0
   Q          0      0      0      0      1      0
   Theta      0      0      0      0      0      1
 
  D = 
          Delta
   XN         0
   XD         0
   U          0
   W          0
   Q          0
   Theta      0
 
Continuous-time state-space model.

Валидация статической устойчивости с динамической характеристикой

С помощью построенной линейной модели пространства состояний можно построить график динамического поведения системы.

Чтобы проверить результаты статической устойчивости с динамическим поведением системы, постройте график модели пространства состояний с учетом начальных условий.

Чтобы вызвать возмущение в системе, к сигналу лифта добавляется 5 степень шаг в течение 1 секунды.

x0 = getState(cruiseState, linSys.OutputName);
t = linspace(0, 50, 500);
u = zeros(size(t));
u(t > 1 & t < 2) = 5;
lsim(linSys,u,t,x0)

Как ожидалось из статического анализа устойчивости, скорость воздуха и скорость тангажа стабильны при реагировании на небольшие возмущения в лифте.

Изоляция отклика лифт-тангаж

В дополнение к статической верификации устойчивости, изоляция управляющих поверхностей к их предполагаемой динамической реакции может помочь спроектировать контроллеры, характерные для отдельных поверхностей.

В этом случае имеется только одна управляющая поверхность - лифт.

Лифт управляет тангажом ответом самолета. Чтобы показать ответ тангажа, изолируйте вход лифта от угла тангажа до передаточной функции лифта.

linSysElevatorTF = tf(linSys(6,1))

linSysElevatorTF =
 
  From input "Delta" to output "Theta":
           -2.321 s^3 - 0.7222 s^2 - 0.001232 s - 8.935e-09
  ------------------------------------------------------------------
  s^5 + 26.75 s^4 + 15.19 s^3 + 0.03932 s^2 + 0.008227 s + 5.712e-08
 
Continuous-time transfer function.

step(linSysElevatorTF)

Figure contains an axes. The axes with title From: Delta To: Theta contains an object of type line. This object represents linSysElevatorTF.

Как видно на графике шага, ответ тангажа на вход лифта имеет нежелательную колебательную природу и большую установившуюся ошибку.

Путем добавления ПИД обратной связи контроллера к входу лифта может быть достигнута гораздо более желательная характеристика тангажа.

C = pidtune(linSysElevatorTF, "PID")

C =
 
        1 
  Ki * ---
        s 

  with Ki = -0.000524
 
Continuous-time I-only controller.

elevatorFeedback = feedback(linSysElevatorTF * C, 1)

elevatorFeedback =
 
  From input to output "Theta":
                                                                 
           0.001217 s^3 + 0.0003786 s^2 + 6.459e-07 s + 4.684e-12
                                                                 
  ------------------------------------------------------------------------
                                                                          
  s^6 + 26.75 s^5 + 15.19 s^4 + 0.04053 s^3 + 0.008605 s^2 + 7.03e-07 s   
                                                                          
                                                               + 4.684e-12
                                                                          
 
Continuous-time transfer function.

step(elevatorFeedback)

Figure contains an axes. The axes with title From: In(1) To: Theta contains an object of type line. This object represents elevatorFeedback.

Документация