Вывод уравнений движения ЛА типа "самолет"
Допущения:
1) Конструкция ЛА абсолютно жесткая (ЛА - твердое тело);
2) Атмосфера спокойная (W=0 - скорость ветра);
3) Масса ЛА постоянная т.к. заданное время переходных процессов меньше 60 с;
4) Земля плоская и не вращается;
5) Угол тангажа не равен +/-90 градусов;
6) Спарведлива гипотеза стационарности (согласно этой гипотезе, аэродинамическая сила и аэродинамический момент в нестационарном режиме полета полностью определяются величинами параметров в данный момент времени).
Системы координат и матрицы перехода

MatrixForm[
      M1 = {{Cos[Bteta[t]], -Sin[Bteta[t]], 0}, {Sin[Bteta[t]], 
            Cos[Bteta[t]], 0}, {0, 0, 
            1}}]

MatrixForm[
      M2 = {{Cos[BPsi[t]], 0, Sin[BPsi[t]]}, {0, 1, 0}, {-Sin[BPsi[t]], 
            0, Cos[BPsi[
                t]]}}]

MatrixForm[
      M = M2 . M1]

MatrixForm[
      ReplaceAll[
        M, {BPsi[t] -> 0, 
          Bteta[t] -> 
            0}]]

\(\(\(MatrixForm[
      MM1 = {{Cos[LPsi[t]], 0, Sin[LPsi[t]]}, {0, 1, 0}, {\(-Sin[LPsi[t]]\), 
            0, Cos[LPsi[
                t]]}}]\)\(\ \)\( (*Матрица\ поворота\ на\ угол\ рыскания*) \)\)\)

MatrixForm[
      MM2 = {{Cos[LTeta[t]], -Sin[LTeta[t]], 0}, {Sin[LTeta[t]], 
            Cos[LTeta[t]], 0}, {0, 0, 
            1}}]

MatrixForm[
      MM3 = {{1, 0, 0}, {0, Cos[LGamma[t]], -Sin[LGamma[t]]}, {0, 
            Sin[LGamma[t]], 
            Cos[LGamma[
                t]]}}]

MatrixForm[
      MM = MM1 . MM2 . 
          MM3]

MatrixForm[
      SM1 = {{Cos[Alpha[t]], -Sin[Alpha[t]], 0}, {Sin[Alpha[t]], 
            Cos[Alpha[t]], 0}, {0, 0, 
            1}}]

MatrixForm[
      SM2 = {{Cos[Betta[t]], 0, Sin[Betta[t]]}, {0, 1, 
            0}, {-Sin[Betta[t]], 0, 
            Cos[Betta[
                t]]}}]

MatrixForm[
      SM = Transpose[
          SM2 . SM1]]

Кинематические уравнения поступательного движения

R={L[t], H[t], 
        dZ[t]}

vg={vgx[t], vdy[t], 
        vgz[t]}

\(\(\(MatrixForm[D[R, t]]\)\(==\)\(MatrixForm[
        vg]\)\(\ \)\( (*Кинематическое\ уравнение\ движения\ центра\ масс\ самолета\ в\ векторной\ форме*) \)\)\)

MatrixForm[
      vzem = {vk[t], 0, 
          0}]

\(\(\(MatrixForm[
      vg = M . vzem]\)\(\ \)\( (*Проецируем\ вектор\ земной\ скорости\ на\ нормальные\ оси\ координат*) \)\)\)

\(\(\(MatrixForm[D[R, t]]\)\(==\)\(MatrixForm[
        vg]\)\(\ \)\( (*Кинематическое\ уравнение\ движения\ центра\ масс\ самолета\ в\ векторной\ форме*) \)\)\)

ukpd1 = 
      D[R[[1]], t] == vg[[1]]

ukpd2 = 
      D[R[[2]], t] == vg[[2]]

ukpd3 = 
      D[R[[3]], t] == vg[[3]]

Динамические уравнения поступательного движения

v={vx[t], vy[t], 
        vz[t]}

Clear[fx, fy, fz]

f={fx, fy, 
        fz}

w={wx[t], wy[t], 
        wz[t]}

MatrixForm[
        massa*(D[v, t] + Cross[w, v])]==MatrixForm[
        f]

udpd1 = (massa*(D[v, t] + 
                  Cross[w, v]))[[1]] == 
        f[[1]]

udpd2 = (massa*(D[v, t] + 
                  Cross[w, v]))[[2]] == 
        f[[2]]

udpd3 = (massa*(D[v, t] + 
                  Cross[w, v]))[[3]] == 
        f[[3]]

Кинематические уравнения вращательного движения

\(proekc = \(\({D[LGamma[t], t], 0, 
          0}\)\(+\)\(Transpose[MM2 . MM3] . {0, D[LPsi[t], t], 
            0}\)\(+\)\(Transpose[MM3] . {0, 0, 
            D[LTeta[t], 
              t]}\)\(\ \)\( (*Проецируем\ угловые\ скорости\ рыскания, \ тангажа\ и\ крена\ на\ оси\ связанной\ системы\ координат*) \)\)\)

ukvd1 = w[[1]] == proekc[[1]]

ukvd2 = w[[2]] == proekc[[2]]

ukvd3 = w[[3]] == proekc[[3]]

Динамические уравнения вращательного движения

MatrixForm[
      TI = {{Ixx, -Ixy, -Ixz}, {-Iyx, 
            Iyy, -Iyz}, {-Izx, -Izy, 
            Izz}}]

Ixz = Izx = Iyz = Izy=0

Iyx = Ixy

MatrixForm[TI]

MatrixForm[
      K = TI . w]

Clear[Mx, My, Mz]

\(\(\(M\)\(=\)\({Mx, My, 
        Mz}\)\(\ \)\( (*Результирующий\ момент\ аэродинамических\ сил\ и\ тяги\ в\ связанной\ системе\ координат \((Момент\ силы\ тяжести, \ центробежных\ и\ кориолисовых\ сил\ инерции\ вокруг\ центра\ масс\ ЛА\ равен\ нулю)\)*) \)\)\)

lchudvd = D[K, t]+Cross[w, 
          K]

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

udvd1 = M[[1]] == lchudvd[[1]]

udvd2 = M[[2]] == lchudvd[[2]]

udvd3 = M[[3]] == lchudvd[[3]]

MatrixForm[
      polnsyst = {ukpd1, ukpd2, ukpd3, udpd1, udpd2, udpd3, ukvd1, ukvd2, 
          ukvd3, udvd1, udvd2, 
          udvd3}]

Уточним зависимости для сил и моментов

Clear[Cy]
Cy[Alpha[t], Mah[t], ugolstrel, ugolVobr, dzak[t], 
        drv[t]] = D[
          Cy[Alpha[t], Mah[t], ugolstrel, ugolVobr, dzak[t], drv[t]], 
          Alpha[t]]*(Alpha[t] - 
            Alpha0)

lyambda=(l^2)/
        S

Cx = Cx0+(Cy[Alpha[t], Mah[t], ugolstrel, ugolVobr, 
                dzak[t], drv[t]]^2)/(lyambda*
              Pi)

\(\(\( (*q = \((ro[H[t]]*\((vk[t])\)^2)\)/
            2\ Скоростной\ напор*) \)\(
\)\(q = \((ρ*\((vk[t])\)^2)\)/2\)\)\)

Maer={kmx*q*s*l, 
        kmy*q*s*l, 
        kmz*q*s*ba}

\(\(\( (*faer = {\(-Cx\)*q*s, 
          Cy[Alpha[t], Mah[t], ugolstrel, ugolVobr, dzak[t], drv[t]]*q*s, 
          Cz[Betta[t], drn[t]]*q*
            s}*) \)\(\ \)\( (*Аэродинамическая\ сила\ в\ скоростной\ системе\ координат*) \)\(
\)\(\(faer\)\(=\)\({\(-Cx\)*q*s, 
        Cy[Alpha[t]]*q*s, 
        Cz[Betta[t], drn[t]]*q*
          s}\)\( (*Аэродинамическая\ сила\ в\ скоростной\ системе\ координат*) \)\)\)\)

MatrixForm[
      faersv = SM . 
          faer]

MatrixForm[
      Gsv = Transpose[MM] . {0, -G, 
            0}]

Mx = Maer[[1]]

My = Maer[[2]]

Mz = 
      Maer[[3]] - P[dsg[t], vk[t], H[t]]*yp*Cos[fip]

fx = 
      faersv[[1]] + Gsv[[1]] + 
        P[dsg[t], vk[t], H[t]]*Cos[fip]

fy = 
      faersv[[2]] + Gsv[[2]] + 
        P[dsg[t], vk[t], H[t]]*Sin[fip]

fz = faersv[[3]] + Gsv[[3]]

\(\(\(MatrixForm[
      polnsystut = {ukpd1, ukpd2, ukpd3, udpd1, udpd2, udpd3, ukvd1, ukvd2, 
          ukvd3, udvd1, udvd2, 
          udvd3}]\)\(\ \)\( (*Уточненная\ полная\ система\ уравнений\ движения\ ЛА*) \)\)\)

vektorsost=Flatten[{R, v, 
          w, {LPsi[t], LTeta[t], 
            LGamma[t]}}]

MatrixForm[
      SystFK = Flatten[
          Simplify[
            Solve[polnsystut, D[vektorsost, t], 
              Sort -> 
                False]]]]

MatrixForm[
      systema = 
        Table[SystFK[[i, 1]] == 
            SystFK[[i, 2]], {i, 1, 
            12}]]

Представим нелинейную нестационарную модель движения ЛА в матричном виде (Overscript[X, ·])[t]=Ax[t]+Bu[t]

MatrixForm[
      MatrA = Table[
          D[systema[[i, 2]], vektorsost[[j]]], {i, 1, 
            12}, {j, 1, 
            12}]]

MatrixForm[
      Perest = {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 
            0, 0, 0, 0, 0}, {0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 
            0, 1, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 
            0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}, {0, 0, 1, 0, 0, 0, 0, 
            0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 
            0, 0, 0, 1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
            0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 
            0, 0, 0, 0, 1}}]

U={drv[t], deler[t], drn[t], 
        dsg[t]}

MatrixForm[
      MatrB = Table[
          D[systema[[i, 2]], U[[j]]], {i, 1, 12}, {j, 
            1, 4}]]

MatrixForm[MatrARazd = Perest . MatrA]

Матрица состояния "свободного" ЛА разделяется на 4 матрицы. Две из них характеризуют динамику продольного движения - Apr и бокового движения - Abok, соответственно. Две другие характеризуют влияние боковго движения на продольное и продольного на бокове. Т.к. ЛА типа "самолет" имеет одну ось симметрии, то разделение на продольное и боковое движение возможно и коэффициенты матриц взаимосвязи движений пренебрежимо малы по сравнению с коэффициентами матриц Apr и  Abok.

MatrixForm[
      Apr = Table[
          MatrARazd[[i, j]], {i, 1, 6}, {j, 1, 6}]]

MatrixForm[
      Abok = Table[
          MatrARazd[[i, j]], {i, 7, 12}, {j, 7, 12}]]

a1 = 0.635; 
a2 = 5.47; 
a3 = 2.72; 
a4 = 0.269; 
a5 = 3.26; 
a6 = 0.709; 
a7 = 0.043;

b1 = 3.1; 
b2 = 20.2; 
b3 = 17.6; 
b4 = 0.072; 
b5 = -0.518; 
b6 = 0.057; 
b7 = 0.065;

c1 = 0.82; 
c2 = 11.0; 
c3 = 8.8; 
c4 = 1.2; 
c5 = 0.3; 
c6 = 7.3;

Допустимость изолированного исследования короткопериодического и длиннопериодического движения самолета:

If[c2 > 0,   Print[Dopustimo], 
      Print[Nedopustimo]]

Допустимость изолированного исследования "быстрого" бокового и спирального движения самолета:

If[
      b1*(a1*a4 + a2) - b2*b6 + a1*b2*b7 - 
          a6*(a4*b6 + a2*b7) ≥ 9*b2*b4,   Print[
        Dopustimo], Print[Nedopustimo]]

Допустимость изолированного исследования составляющей "рыскание-скольжение" и креновой составляющей "быстрого" бокового движения самолета:

If[
      b1*(a1*a4 + a2) >= 
        0.9*(a1*b2*b7 - b2*b6 - a6*(a4*b6 + a2*b7)),   Print[
        Dopustimo], Print[Nedopustimo]]

MatrixForm[
      Apr = {{-c1 - c5, -c2 + c4*c5, -e3 + c5*e2, 
            0, -c5*c10}, {1, -c4, -e2, 
            0, -c10}, {0, -c8, -e1, 0, -c7}, {0, -c6, 
            c11, 0, c6}, {1, 0, 0, 0, 0}}]

MatrixForm[
      Xpr = {wz[t], Alpha[t], dV[t], dH[t], 
          LTeta[t]}]

MatrixForm[
      Abok = {{-b1, -a6, -b2, 0, 0}, {-b6, -a1, -a2, 
            0, 0}, {b7, a8, -a4, b4, 0}, {1, -a9, 0, b8, 0}, {0, a10, 
            0, -b9, 0}}]

\(\(\(MatrixForm[
      Xbok = {wx[t], wy[t], Betta[t], LGamma[t], 
          LPsi[t]}]\)\(\ \)\( (*Вектор\ состояния\ бокового\ движения*) \)\)\)

Рассмотрим свободный самолет:

dxpr = D[Xpr, t]

lxpr = LaplaceTransform[dxpr, t, s]

LaplaceTransform[Cosh[t], t, s]

lxpr = Xpr . DiagonalMatrix[{s, s, s, s, s}]

UrPr = MatrixForm[lxpr - Apr . Xpr]

MatrixForm[
      MatrUrPr = 
        Table[D[UrPr[[1, i]], Xpr[[j]]], {i, 1, 
            5}, {j, 1, 5}]]

MatrixForm[MatrUrPr = ReplaceAll[MatrUrPr, c10 -> 0]]

MatrUrPr[[1]] = 
      Simplify[MatrUrPr[[1]] + 
          MatrUrPr[[2]]*c5]

MatrixForm[
      MatrUrPr]

dxbok = D[Xbok, t]

lxbok = LaplaceTransform[dxbok, t, s]

lxbok = Xbok . DiagonalMatrix[{s, s, s, s, s}]

UrBok = MatrixForm[lxbok - Abok . Xbok]

Abok . Xbok

MatrixForm[
      MatrUrBok = 
        Table[D[UrBok[[1, i]], Xbok[[j]]], {i, 1, 
            5}, {j, 1, 
            5}]]

MatrixForm[
      MatrUrKPD = 
        Delete[MatrUrPr, {{1, 3}, {1, 4}, {1, 5}, {2, 3}, {2, 4}, {2, 
              5}, {3}, {4}, {5}}]]

моментное возмущение

Walmz = Det[{{MatrUrKPD[[1, 1]], 
              57.3}, {MatrUrKPD[[2, 1]], 0}}]/
        Det[MatrUrKPD]

halmz = InverseLaplaceTransform[Walmz/s, s, t]

Plot[halmz, {t, 0, 10}, PlotRange -> {0, 7}]

Рассмотрим моментное возмущение для кренового изолированного движения свободного самолета:

MatrUrKren = 
      Delete[MatrUrBok, {{2}, {3}, {5}, {1, 2}, {1, 3}, {1, 5}, {4, 2}, {4, 
            3}, {4, 5}}]; 
MatrixForm[
      MatrUrKren = ReplaceAll[MatrUrKren, {b8 -> 0}]]

\(Wkrenmx = \(\(Det[{{MatrUrKren[\([\)\(1, 1\)\(]\)], 
              57.3}, {MatrUrKren[\([\)\(2, 1\)\(]\)], 0}}]/
        Det[MatrUrKren]\)\(\ \)\( (*Передаточная\ функция\ крен/
            Mx*) \)\)\)

hkrenmx = InverseLaplaceTransform[Wkrenmx/s, s, t]

Plot[hkrenmx, {t, 0, 10}, PlotRange -> {0, 100}]

MatrUrRSk = 
      Delete[MatrUrBok, {{1}, {4}, {5}, {2, 1}, {2, 4}, {2, 5}, {3, 1}, {3, 
            4}, {3, 5}}]; 
MatrixForm[
      MatrUrRSk = ReplaceAll[MatrUrRSk, {a8 -> 1}]]

Wrskmy = Det[{{MatrUrRSk[[1, 1]], 
              57.3}, {MatrUrRSk[[2, 1]], 0}}]/
        Det[MatrUrRSk]

hrskmy = InverseLaplaceTransform[Wrskmy/s, s, t]

Plot[hrskmy, {t, 0, 15}, PlotRange -> {0, 16}]

\(\(\(
\)\(
\)\(
\)\(АП\ крена\ с\ ЖОС\)\)\)

asd

treg = 2;

mu = (9.48 - b1*treg)/b3/treg

ie = 22.5/b3/treg/treg

MatrixForm[
      Udv  =  {{(s + b1)*s,  b3},  {-(mu*s + ie),  
            1}}]

MatrixForm[
      Matgam2  =  {{0,  b3},  {-ie,  
            1}}]

Wgamma = Det[Matgam2]/Det[
          Udv]

hgamma = InverseLaplaceTransform[Wgamma/s, s, t]

Plot[hgamma, {t, 0, 10}, PlotRange -> {0, 1.2}];

\(Задача : \ определить\ оптимальные\ значения\ передаточных\ чисел\ автопилота\ при\ изменении\ коэф . \ в\ лин . \ модели\ в\ диапаз\)