Имитационное моделирование
Цель: определить вероятность того, что система самолет-автопилот тангажа ЖОС с астатическим законом управления с рассчитанными фиксированными передаточными числами на данном режиме полета не будет удовлетворять требуемым характеристикам переходного процесса из-за неточности определения коэффициентов линейной модели (±10%) и рассчитать при этом наиболее вероятное время регулирования и перерегулирование.

Коэффициенты линейной модели самолета:

c1=0.88;
c2=0.404;
c3=7.028;
c4=0.987;
c5=0.251;
c6=2.640;

Рассчет передаточных чисел автопилота:

OmegaA=Sqrt[c1*c4+c2] (*Собственная частота короткопериодического движения самолета*)

TA=1/OmegaA (*Постоянная времени короткопериодического движения самолета*)

Tv=1/c4 (*Постоянная времени форсирующего звена*)

Kc=c3*c4/(c1*c4+c2)

ash=2*(c1+c5-c4)/c3

bsh=((c1+c4+c5)^2-4*(c1*c4+c2))/(c3^2)

mu=4*(-ash/2+Sqrt[(ash^2)/4-bsh] )(*Передаточное число автопилота по сигналу угловой скорости тангажа*)

Kcsh=Kc/(1+mu*Kc)

TAsh=TA/Sqrt[1+mu*Kc]

iv=1.2/(Kcsh*Tv) (*Передаточное число автопилота по сигналу угла тангажа*)

nu=0.09/(Kcsh*Tv^2) (*Передаточное число автопилота по сигналу интеграла угла тангажа*)

Tnu=iv/nu

Получение передаточной и переходной функции. Построение переходного процесса.

MatrixForm[M1={{s+c1,c5*s+c2,0,c3},{-1,s+c4,0,0},{-1,0,s,0},{mu,0,nu*(Tnu*s+1)/s,-1}}] (*Матрица состояния свободного самолета*)

MatrixForm[M2={{s+c1,c5*s+c2,0,c3},{-1,s+c4,0,0},{-1,0,0,0},{mu,0,nu*(Tnu*s+1)/s,-1}}] (*Матрица, в которой столбец соответствующий углу тангажа заменен на столбец соответствующий рулю высоты*)

W=Simplify[Det[M2]/Det[M1]] (*Расчет передаточной функции, как отношения определителей матриц М2 и М1*)

h=Simplify[InverseLaplaceTransform[W/s,s,t]] (*Получение переходной функции - обратное преобразование Лапласа*)

hf[t_]=Abs[h]

Plot[h,{t,0,20},PlotRange->{0,1.2}]

Аналитическое определение времени регулирования и перерегулирования:

kol=10000; (*Количество промежутков разбиения  диапазона, на котором ищется решение*)
dt=20/kol; (*Шаг разбиения*)

Ust=hf[50000]; (*Определение установившегося значения*)
hmax=hf[0]; (*Начальное значение максимального значения*)
For[i=1,i<kol+1,{If[Sign[Abs[hf[dt*(i-1)]-Ust]-0.05]≠Sign[Abs[hf[dt*(i)]-Ust]-0.05],Treg=dt*(i-0.5)],If[hmax<hf[dt*i],hmax=hf[dt*i]]};i++] (*Численное решение уравнений h(t)=0.95 и h(t)=1.05, объединенных в одно |h(t)-1|-0.05=0 и определение максимума*)
sigma=Round[100*(hmax-Ust)/Ust,0.01]; (*Расчет перерегулирования*)
Print["Treg=",Treg,"  с."] 
Print["Sigma=",sigma,"%."]

Soch=4; (*Количество уровней в каждом факторе равное Soch-1*)
koefstep=(1.1-0.9)/(Soch-2) (*Шаг изменения коэффициентов линейной модели*)

Полный факторный эксперимент:
в последнем вложенном цикле (m) происходит: расчет коэффициентов линейной модели самолета для текущего сочетания уровней факторов, расчет передаточной и переходной функции системы, аналитическое определение времени регулирования и перерегулирования и их оценка.

ocen=0; (*Переменная, соответсвующая количеству режимов, на которых получаются неудовлетворительные характеристики переходных процессов*)
VTsr={}; (*Массив для записи времени регулирования на неудовлетворительных режимах*)
VSsr={}; (*Массив для записи перерегулирования на неудовлетворительных режимах*)
For[i=1,i<Soch, (*Цикл для изменения коэффициента с1*)
{For[j=1,j<Soch, (*Цикл для изменения коэффициента с2*)
{For[k=1,k<Soch, (*Цикл для изменения коэффициента с3*)
{For[l=1,l<Soch, (*Цикл для изменения коэффициента с4*)
{For[m=1,m<Soch, (*Цикл для изменения коэффициента с5*)
{c1c=c1*(0.9+koefstep*(i-1));
c2c=c2*(0.9+koefstep*(j-1));
c3c=c3*(0.9+koefstep*(k-1));
c4c=c4*(0.9+koefstep*(l-1));
c5c=c5*(0.9+koefstep*(m-1));
M1={{s+c1c,c5c*s+c2c,0,c3c},{-1,s+c4c,0,0},{-1,0,s,0},{mu,0,nu*(Tnu*s+1)/s,-1}};(*Матрица состояния свободного самолета*)
M2={{s+c1c,c5c*s+c2c,0,c3c},{-1,s+c4c,0,0},{-1,0,0,0},{mu,0,nu*(Tnu*s+1)/s,-1}};(*Матрица, в которой столбец соответствующий углу тангажа заменен на столбец соответствующий рулю высоты*)
W=Det[M2]/Det[M1];(*Расчет передаточной функции, как отношения определителей матриц М2 и М1*)
h=InverseLaplaceTransform[W/s,s,t];(*Получение переходной функции - обратное преобразование Лапласа*)
hf[t_]=Abs[h];
Ust=hf[50000];(*Определение установившегося значения*)
hmax=hf[0];(*Начальное значение максимального значения*)
For[q=1,q<kol+1,{If[Sign[Abs[hf[dt*(q-1)]-Ust]-0.05]≠Sign[Abs[hf[dt*(q)]-Ust]-0.05],Treg=dt*(q-0.5)],If[hmax<hf[dt*q],hmax=hf[dt*q]]};q++](*Численное решение уравнений h(t)=0.95 и h(t)=1.05, объединенных в одно |h(t)-1|-0.05=0 и определение максимума*)
sigma=Round[100*(hmax-Ust)/Ust,0.01];(*Расчет перерегулирования*)If[Treg>5|| sigma>5,{Print["Режим с коэффициентами с1=",c1c,", c2=",c2c,", c3=",c3c,", c4=",c4c,", c5=",c5c," не удовлетворяет требованиям к переходному процессу, т.к. Treg=",Treg," c, Sigma=",sigma,"%."],ocen++,VTsr=Append[VTsr,Treg],VSsr=Append[VSsr,sigma]}]}; (*Вывод текущего результата моделирования*)
m++]};
l++]};
k++]};
j++]};
i++]
Ver=(ocen/((Soch-1)^5))*100.; (*Расчет вероятности того, что система не удовлетворяет заданным характеристикам переходного процесса*)
Tregsr=Mean[VTsr]; (*Расчет наиболее вероятного времени регулирования*)
sigmasr=Mean[VSsr];(*Расчет наиболее вероятного перерегулирования*)
Print["Вероятность того, что система самолет – рассчитанный автопилот на данном режиме полета не будет удовлетворять заданным характеристикам переходного процесса равна ",Ver,"%."]
If[Ver>0,Print["При этом наиболее вероятное время регулирования – ",Tregsr,"с, а наиболее вероятное перерегулирование - ",sigmasr,"%. Для решения этой проблемы следует уменьшить число nu(передаточное число АП по сигналу интеграла угла тангажа), либо увеличить число iv(передаточное число АП по сигналу угла тангажа), либо уменьшать nu и увеличивать iv одновременно, но на меньшие величины, чем в первых двух случаях, т.е. таким образом сдвнуть Tnu ещё дальше от Tv - увеличить Tnu."]]

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

Mver={{3,7.41},{5,5.12},{7,4.14}}

ListLinePlot[Mver,PlotRange->{0,10}]