Домашнее задание. Расчет порога срабатывания.

Коэффициенты линейной модели самолета для различных режимов:

c1 = {0.560, 0.560, 0.496, 0.496, 0.670, 0.670, 0.590, 0.590, 1.100, 1.100, 0.940, 0.940, 0.760, 0.760, 0.880, 0.880, 0.668, 0.668, 0.793, 0,793};
c2 = {1.635, 0.922, 1.273, 0.718, 1.864, 0.655, 1.200, 0.245, 1.864, 0.569, 2.259, 0.0726, 2.640, 0.357, 3.337, 0.404, 1.742, 0.056, 2.980, 0.681};
c3 = {1.468, 1.468, 1.142, 1.142, 2.570, 2.570, 2.071, 2.071, 7.827, 7.827, 5.975, 5.975, 5.617, 5.617, 7.028, 7.028, 4.496, 4.496, 6.340, 6.340};
c4 = {0.765, 0.765, 0.868, 0.868, 1.053, 1.053, 0.944, 0.944, 1.122, 1.122, 0.962, 0.962, 0.864, 0.864, 0.987, 0.987, 0.671, 0.671, 0.875, 0.875};
c5 = {0.160, 0.160, 0.168, 0.168, 0.195, 0.195, 0.107, 0.107, 0.314, 0.314, 0.269, 0.269, 0.226, 0.226, 0.251, 0.251, 0.191, 0.191, 0.226, 0.226};

MPS={};
MTreg={};
Mperereg={};
MT={};

For[rezim=1,rezim≤20,rezim++,
{OmegaA=Sqrt[c1[[rezim]]*c4[[rezim]]+c2[[rezim]]]; (*Собственная частота короткопериодического движения самолета*)
TA=1/OmegaA; (*Постоянная времени короткопериодического движения самолета*)
Tv=1/c4[[rezim]]; (*Постоянная времени форсирующего звена*)
Kc=c3[[rezim]]*c4[[rezim]]/(c1[[rezim]]*c4[[rezim]]+c2[[rezim]]);
ash=2*(c1[[rezim]]+c5[[rezim]]-c4[[rezim]])/c3[[rezim]];
bsh=((c1[[rezim]]+c4[[rezim]]+c5[[rezim]])^2-4*(c1[[rezim]]*c4[[rezim]]+c2[[rezim]]))/((c3[[rezim]])^2);
mu=Abs[-ash/2+Sqrt[(ash^2)/4-bsh]]; (*Передаточное число автопилота по сигналу угловой скорости тангажа*)
Kcsh=Kc/(1+mu*Kc);
TAsh=TA/Sqrt[1+mu*Kc];
(*Передаточное число автопилота по сигналу угла тангажа - iv, gередаточное число автопилота по сигналу интеграла угла тангажа - nu:*)
If[Tv<(10*TAsh),{iv=1/(Kcsh*Tv), nu=0.09/(Kcsh*Tv^2)},{iv=0.5/(Kcsh*Tv), nu=5/(Kcsh*Tv )}];
Tnu=iv/nu;
M1={{s+c1[[rezim]],c5[[rezim]]*s+c2[[rezim]],0,c3[[rezim]]},{-1,s+c4[[rezim]],0,0},{-1,0,s,0},{mu,0,nu*(Tnu*s+1)/s,-1}};(*Матрица состояния свободного самолета*)
M2={{s+c1[[rezim]],c5[[rezim]]*s+c2[[rezim]],0,c3[[rezim]]},{-1,s+c4[[rezim]],0,0},{-1,0,0,0},{mu,0,nu*(Tnu*s+1)/s,-1}};(*Матрица, в которой столбец соответствующий углу тангажа заменен на столбец соответствующий заданному углу тангажа*)
W=Simplify[Det[M2]/Det[M1]]; (*Расчет передаточной функции, как отношения определителей матриц М2 и М1*)
Wch=W/.s->I*om;
F1=(Abs[Wch])^2;
Fin=1;
Fout=F1*Fin;
pl=NIntegrate[Fout,{om,0,Infinity}];
tochnost=6;
om0=0;
step=0.1;
While[step≥10^-tochnost,{While[(NIntegrate[Fout,{om,0,om0}]-0.95*pl)<0,om0=om0+step],om0=om0-step,step=step/10}];
om0=om0-step*5;
kpkv=0.1;(*Поправочный коэффициент*)
T=kpkv*Pi/om0; (*Период квантования*)
h=InverseLaplaceTransform[W/s,s,t];
hf[t_]=Abs[h];
(*Нахождение времени регулирования и перерегулирования*)
kol=50000; (*Количество промежутков разбиения  диапазона, на котором ищется решение*)
dt=50/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 и определение максимума*)
sigmaper=Round[100*(hmax-Ust)/Ust,0.01]; (*Расчет перерегулирования*)
T1=Round[T*100000]/100000; (*Период квантования, округленный до 5 знаков после запятой*)
Nrasch=Ceiling[Treg/T1]; (*Количество периодов квантования за время регулирования Treg*)
Nekstr=4; (*Степень полинома экстраполяции*)
MOtkl={}; (*Массив отклонений*)
hizm=Table[0,{Nekstr}]; (*Массив последних Nekstr измеренных (реальных) значений используемый при экстраполяции*)
For[i=Nekstr-1,i<Nrasch,i++,
{Ttek=T1*i;
hTtek=hf[Ttek];
hizm[[1]]=hf[Ttek-T1*3];
hizm[[2]]=hf[Ttek-T1*2];
hizm[[3]]=hf[Ttek-T1*1];
hizm[[4]]=hf[Ttek-T1*0];
hekstr=hizm[[4]]+T1*(55*hizm[[4]]-59*hizm[[3]]+37*hizm[[2]]-9*hizm[[1]])/24;(*Экстраполяция*)
MOtkl=Append[MOtkl,Abs[hf[Ttek+T1]-hekstr]];
}];
MO=Mean[MOtkl]; (*Математическое ожидание отклонения истинного значения угла тангажа, полученного во время кратное периоду квантования Т, от экстраполированного значения*)
sigma=StandardDeviation[MOtkl]; (*Среднее квадратическое отклонение отклонения*)
b=3; (*1...5 - доверительная граница*)
PS=MO+b*sigma; (*Порог срабатывания*)
MT=Append[MT,T1*1.];
MPS=Append[MPS,PS];
MTreg=Append[MTreg,Treg];
Mperereg=Append[Mperereg,sigmaper];
}]

For[i=1,i≤20,i++,Print["Режим №",i,": время регулирования Tрег = ",MTreg[[i]],"с, перерегулирование σ = ",Mperereg[[i]],"%, период квантования T = ",MT[[i]],"с, порог срабатывания равен ",MPS[[i]],"."]]
Print[]
Print["Наибольшее значение порога срабатывания равно ",Max[MPS],", среднее значение порога срабатывания равно ",Mean[MPS],"."]
Print[]
Print["Δ = ",Mean[MPS]+3*StandardDeviation[MPS]]