Задачи дифференциальной и топологической диагностики. Часть 8. Движение летательного аппарата и алгоритмы его диагностирования

Полный текст

Введение

Как уже отмечалось ранее [1; 8; 9], задача дифференциальной диагностики функционального состояния объектов управления может быть сведена к двум самостоятельным последовательно решаемым задачам: задаче контроля, т. е. установлению критерия наличия неисправности в системе, и задаче диагностирования, т. е. поиску происшедшей неисправности. Критерием наличия неисправности в системе может быть выход траектории объекта на некоторую заранее выбранную поверхность. Неисправность может произойти в любой заранее неизвестный момент времени движения объекта, в любой точке внутри данной поверхности контроля [10–12].

Исходной информацией при решении задачи контроля является математическая модель движения рассматриваемого объекта, ограниченная область ее начальных условий и априорный список математических моделей движения объекта с той или иной возможной неисправностью. По этой информации может быть выбрана поверхность контроля. Задача же диагностирования может быть решена путем последующего слежения за траекторией объекта после ее выхода на поверхность контроля. При этом необходимо, чтобы процесс диагностики совершался во время движения объекта, был осуществлен в течение весьма краткого интервала времени. Эти обстоятельства иногда не позволяют использовать довольно громоздкие алгоритмы теории идентификации и приводят к необходимости построения алгоритмов непрерывной экспресс-диагностики [13–15].

Рассматривается применение развиваемой методики диагностирования к теории летательных аппаратов.

В качестве численного эксперимента было рассмотрено движение летательного аппарата (ЛА), описываемого уравнениями, приведенными в [16]. ЛА находится в режиме планирующего спуска с высот, близких к орбитальным ($\cong$100 км), с начальной скоростью, близкой к первой космической скорости. Воспроизведем полнее и уточним уравнения движения ЛА и покажем, что они приводятся в некотором смысле к каноническому виду.

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

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

1. Cистемы, описывающие движение летательного аппарата

Эти уравнения имеют следующую структуру.

Динамические уравнения центра масс (штрихом в системах нормального вида обозначается производная по времени)

$V{\text{'}}_{y_i}=-W_{e_{y_i}}+g_{y_i}+\frac{F_{y_i}}{m}, i=1,2,3,$ (1.1)

где $V_{y_i}$ — проекции вектора путевой скорости ЛА $V$ на оси системы координат $M_{y_i}$, связанной с географической вертикалью и ориентированной в азимуте в ортодромической координатной сетке, $W_{e_{y_i}}$ — проекции составляющей ускорения точки $M$, обусловленной кривизной и вращением Земли на те же оси, $g_{y_i}$ — проекции ускорения силы тяжести на рассматриваемую подвижную систему координат, $F_{y_i}$ — проекции силы $F$, действующей на ЛА, где $F=A+T$, $A$ — аэродинамическая сила, $T$ — сила тяги двигателя, $m$ — масса ЛА.

Кинематические уравнения движения центра масс имеют вид:

$\begin{array}{c}\sigma {\text{'}}_1=\frac{V_{y_1}}{r\cos {\sigma }_2}\mathrm{,}\\ \sigma {\text{'}}_2=\frac{V_{y_2}}{r}\mathrm{,}\\ r\text{'}=V_{y_3}\mathrm{.}\end{array}$ (1.2)

Здесь ${\sigma }_1$ и ${\sigma }_2$ — ортодромические долгота и широта центра масс (точки $M$) летательного аппарата, $r$ — радиус-вектор этой точки в системе $O_{y_i}$.

Уравнения движения ЛА вокруг центра масс в проекциях на оси $M_s=M_{s_1}\mathrm{,}{M}_{s_2}\mathrm{,}{M}_{s_3}$ системы, жестко связанной с ЛА, таковы:

$\begin{array}{c}{I}_{s_1}\frac{d{\omega }_{s_1}}{dt}+(I_{s_3}-I_{s_2}){\omega }_{s_2}{\omega }_{s_3}-I_{s_2{s}_3}({\omega }_{s_2}^2-{\omega }_{s_3}^2)=M_{s_1}\end{array}$

$\begin{array}{c}{I}_{s_2}\frac{d{\omega }_{s_2}}{dt}+(I_{s_1}-I_{s_3}){\omega }_{s_1}{\omega }_{s_3}-I_{s_2{s}_3}(\frac{d{\omega }_{s_3}}{dt}+{\omega }_{s_2}{\omega }_{s_3})=M_{s_2}\\ \end{array}$ (1.3)

$I_{s_3}\frac{d{\omega }_{s_3}}{dt}+(I_{s_2}-I_{s_1}){\omega }_{s_1}{\omega }_{s_2}-I_{s_2{s}_3}(\frac{d{\omega }_{s_2}}{dt}-{\omega }_{s_1}{\omega }_{s_2})=M_{s_3}$

Здесь $I_{s_1}\mathrm{,}{I}_{s_2}\mathrm{,}{I}_{s_3}$ — главные, а $I_{s_2{s}_3}$ — центробежный моменты инерции, ${\omega }_{s_i}\mathrm{,}i=1,2,3$ — угловые скорости ЛА; все — в проекциях на оси $M_s$.

Так как

${\omega }_s={\omega }_c+{\alpha }^{\text{'}}+{\beta }^{\text{'}}+{\gamma }_{c^{\text{'}}}\mathrm{,}$ (1.4)

где ${\omega }_c$ — угловая скорость траекторной системы координат, $\alpha$ — угол атаки, $\beta$ — угол скольжения, ${\gamma }_c$ — угол скоростного крена, то имеем еще одну группу кинематических уравнений:

${\omega }_s{D}_{sc}{\omega }_c+{\gamma }_{c^{\text{'}}}{D}_{sn}\left(\begin{array}{c}0\\ 1\\ 0\\ \end{array}\right)+{\beta }^{\text{'}}\left(\begin{array}{c}0\\ 0\\ 1\\ \end{array}\right)+{\alpha }^{\text{'}}\left(\begin{array}{c}1\\ 0\\ 0\\ \end{array}\right)\mathrm{,}$ (1.5)

в которой сами уравнения можно разрешить относительно $\alpha \text{'}\mathrm{,}\beta \text{'}\mathrm{,}\gamma {\text{'}}_c$. Здесь $D_{sc}\mathrm{,}{D}_{sn}$ — матрицы перехода [17; 18].

Кроме того, имеем группу уравнений, выражающих величину ${\omega }_c$ через $U\mathrm{,}{\sigma }_1\mathrm{,}{\sigma }_2\mathrm{,}{\psi }_c$ и $\theta$, где $U$ — угловая скорость вращения Земли, ${\psi }_c$ — угол скоростного курса, $\theta$ — угол наклона траектории:

${\omega }_c=D_{cy}\left[U_y+\psi {\text{'}}_c\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)-\sigma {\text{'}}_2\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)\right]+\sigma {\text{'}}_1{D}_{c\zeta }$ (1.6)

где $D_{cy}\mathrm{,}{D}_{c\zeta }$ — соответствующие матрицы перехода [19; 20].

Из определения углов ${\psi }_c$ и $\theta$ следуют соотношения

$\begin{array}{c}\psi {\text{'}}_c=-\frac{V{\text{'}}_{y_1}\cos {\psi }_c+V{\text{'}}_{y_2}\sin {\psi }_c}{V\cos \theta }\mathrm{,}\\ {\theta }^{\text{'}}=\frac{V{\text{'}}_{y_3}\cos \theta +\sin \theta (V{\text{'}}_{y_1}\sin {\psi }_c-V{\text{'}}_{y_2}\cos {\psi }_c)}{V}\mathrm{,}\end{array}$ (1.7)

где $V$ — абсолютная величина вектора путевой скорости ЛА.

Уравнения (1.1)–(1.7) могут быть представлены в виде одного уравнения нормального вида на своем 14-мерном фазовом многообразии как

$x^{\text{'}}=K\left(x\right)$, (1.8)

где x — 14-мерный вектор [21]:

$x=\left(V_{y_1}\mathrm{,}{V}_{y_2}\mathrm{,}{V}_{y_3}\mathrm{,}{\psi }_c\mathrm{,}\theta \mathrm{,}r\mathrm{,}\lambda \mathrm{,}\phi \mathrm{,}{\omega }_{s_1}\mathrm{,}{\omega }_{s_2}\mathrm{,}{\omega }_{s_3}\mathrm{,}\alpha \mathrm{,}\beta \mathrm{,}{\gamma }_c\right),$ (1.9)

Здесь (для простоты) рассмотрен случай, когда полюс ортодромии лежит на оси вращения Земли, т. е. ${\sigma }_1=\lambda$, ${\sigma }_2=\phi$, где $\lambda$ и $\phi$ — геоцентрические долгота и широта центра масс ЛА [22; 23].

2. Система управления летательного аппарата

Аэродинамические силы и моменты, действующие на ЛА, определяются следующими выражениями:

$\begin{array}{c}X=c_x\frac{\rho V^2}{2}S, Y=c_y\frac{\rho V^2}{2}S, Z=c_z\frac{\rho V^2}{2}S\mathrm{,}\\ M_{s_1}=\frac{\rho V^2}{2}S{b}_a{m}_z, M_{s_2}=\frac{\rho V^2}{2}Sl{m}_x, M_{s_3}=\frac{\rho V^2}{2}Sl{m}_y.\end{array}$ (2.1)

 

где $\rho$ — плотность воздуха на высоте полета, $S$ — характерная площадь ЛА, $b_a$ — средняя аэродинамическая хорда, $l$ — размах крыльев, $c_x\mathrm{,}{c}_y\mathrm{,}{c}_z$ — аэродинамические коэффициенты сил, $m_x\mathrm{,}{m}_y\mathrm{,}{m}_z$ — аэродинамические коэффициенты моментов.

Будем рассматривать аэродинамические коэффициенты в виде

$\begin{array}{c}{c}_x=c_{x\alpha }\left(M\mathrm{,}\alpha \right)+c_{xTP}\left(M\mathrm{,}h\right),\\ c_y=c_{y\alpha }\left(M\mathrm{,}\alpha \right)+c_y^{{\delta }_b}\left(M\right){\delta }_b\mathrm{,}\\ c_z=c_z^{\beta }\left(M\mathrm{,}\alpha \right)\beta +c_z^{{\delta }_H}\left(M\mathrm{,}\alpha \right){\delta }_H\mathrm{,}\end{array}$ (2.2)

$\begin{array}{l}{m}_x=m_{s_1}=m_{s_1\alpha }\left(M\mathrm{,}\alpha \right)+m_{s_1}^{{\omega }_{s_1}}\left(M\right)\frac{b_a}{V}{\omega }_{s_1}+m_{s_1}^{{\delta }_b}\left(M\right){\delta }_b\mathrm{,}\\ m_y=m_{s_2}=m_{s_2}^{\beta }\left(M\mathrm{,}\alpha \right)\beta +m_{s_2}^{{\omega }_{s_2}}\left(M\right)\frac{l}{2V}{\omega }_{s_2}+m_{s_2}^{{\omega }_{s_3}}\left(M\right)\frac{l}{2V}{\omega }_{s_3}+m_{s_2}^{{\delta }_H}\left(M\mathrm{,}\alpha \right){\delta }_H+m_{s_2}^{{\delta }_E}\left(M\right){\delta }_E\mathrm{,}\\ \\ m_z=m_{s_3}=m_{s_3}^{\beta }\left(M\mathrm{,}\alpha \right)\beta +m_{s_3}^{{\omega }_{s_2}}\left(M\right)\frac{l}{2V}{\omega }_{s_2}+m_{s_3}^{{\omega }_{s_3}}\left(M\right)\frac{l}{2V}{\omega }_{s_3}+m_{s_3}^{{\delta }_H}\left(M\mathrm{,}\alpha \right){\delta }_H\mathrm{.}\\ \end{array}$ (2.3)

Здесь ${\delta }_b\mathrm{,}{\delta }_E\mathrm{,}{\delta }_H$ имеют вполне ясный геометрический смысл, а именно отклонение рулей высоты, элеронов и рулей направления, соответственно.

Структура системы управления отклонением рулей высоты, направления и элеронов зависит от выбранной программы движения. В данной работе рассмотрено управления следующего вида:

${\delta }_u=f\left({\delta }_{u\max }\mathrm{,}{\delta }_u^0\right), u=b\mathrm{,}E\mathrm{,}H\mathrm{,}$ (2.4)

где ${\delta }_{u\max }$ — максимальная величина отклонения величин ${\delta }_b\mathrm{,}{\delta }_E\mathrm{,}{\delta }_H$, соответственно, а $f$ — некоторая (не обязательно гладкая) функция.

Переменные ${\delta }_u^0$ определяются по следующим формулам:

$\begin{array}{l}{\delta }_b^0=r_1{\omega }_{s_1}-r_2\Phi \left[1\right]+f\left({\alpha }_m\mathrm{,}{\sigma }_b\right)\\ {\delta }_E^0=k_1{\omega }_{s_2}-k_2\Phi \left[2\right]+f\left({\gamma }_m\mathrm{,}{\sigma }_E\right)+k_7\Phi \left[3\right]\\ {\delta }_H^0=l_1{\omega }_{s_3}-l_2\Phi \left[3\right]+f\left({\beta }_m\mathrm{,}{\sigma }_H\right)+l_7\Phi \left[2\right]\end{array}$ (2.5)

Здесь $\Phi \left[1\right]={\theta }^{\text{'}}$, $\Phi \left[2\right]={\gamma }_n-{\gamma }_{s^{\text{'}}}$, $\Phi \left[3\right]=-{\beta }^{\text{'}}$ — сигналы, поступающие с гироплатформы [24]: значения углов тангажа (${\theta }^{\text{'}}$), скольжения ($-{\beta }^{\text{'}}$) и разности между программным (${\gamma }_n$) и вычисленным (${\gamma }_{s^{\text{'}}}$) значениями угла крена.

Сигналы с постоянными коэффициентами $r_i\mathrm{,}{k}_j\mathrm{,}{l}_s$ формируются в зависимости от углового движения ЛА. Сигналы $f\left({\sigma }_u\right)$ формируются в зависимости от характеристик траекторного движения ЛА. Здесь $f\left({\sigma }_u\right)$ также некоторая (не обязательно гладкая) функция. Сигналы формируются в следующем виде:

$\begin{array}{l}{\sigma }_b=r_3(\alpha \text{'}-{\alpha }_n)+r_4(\widetilde{\alpha \text{'}}-{\alpha }_n)+r_5{\int }_{t_0}^t(\alpha \text{'}-{\alpha }_n)dt\mathrm{,}\\ {\sigma }_E=k_3(\sigma {\text{'}}_{2n}^0-{\sigma }_{2n})+k_4(\widetilde{\sigma {\text{'}}_2^0}-{\tilde{\sigma }}_{2n})+k_5{\int }_{t_0}^t({\sigma }_2^0-{\sigma }_2)dt+k_6(\psi {\text{'}}_s-{\psi }_{sn}),\\ {\sigma }_H=l_3(\sigma {{\text{'}}_{2n}}^0-{\sigma }_{2n})+l_4(\tilde{\sigma }{{\text{'}}_2}^0-{\tilde{\sigma }}_{2n})+l_5{\int }_{t_0}^t(\sigma {{\text{'}}_2}^0-{\sigma }_2)dt+l_6(\psi {\text{'}}_s-{\psi }_{sn}).\end{array}$ (2.6)

Здесь величины $K^{\text{'}}$ (со штрихом) — суть вычисленные значения переменных, $K_n$ (с индексом “n”) — программные значения этих переменных, $r_i\mathrm{,}{k}_j\mathrm{,}{l}_s$ — постоянные коэффициенты, ${\alpha }_m\mathrm{,}{\beta }_m\mathrm{,}{\gamma }_m$ — постоянные, ограничивающие значения рассматриваемых сигналов траекторного слежения.

Значения коэффициентов $r_i\mathrm{,}{k}_j\mathrm{,}{l}_s$ приведены в табл. 2.1.

 

Таблица 2.1

Table 2.1

	1	2	3	4	5	6	7
r	20	0	10	0	0	0	0
k	0,4	0,3	0	0	0	0	1
l	7	7	0	0	0	0	0,3

 

Уравнения (1.8) совместно с (1.9)–(2.6) и с учетом табл. 2.1 могут быть представлены в следующем виде:

$\begin{array}{l}x\text{'}\text{'}=X\left(x\right)+A\left(x\right)\xi \mathrm{,}\\ \xi =\Phi \left(\delta \right)\\ \delta =Bx*+f\left(\sigma \right)\\ \sigma =Cx**\mathrm{.}\end{array}$ (2.7)

Здесь

$x*=({\omega }_{s_1}\mathrm{,}{\omega }_{s_2}\mathrm{,}{\omega }_{s_3}\mathrm{,}\theta *\mathrm{,}{\gamma }_n-{\gamma }_s^{*}\mathrm{,}-\beta )$ (2.8)

и

$x**=(\alpha *-{\alpha }_n\mathrm{,}{\psi }_s^{*}-{\psi }_n\mathrm{,}{\sigma }_s^{*}-{\sigma }_n)$ (2.9)

— шестимерный и трехмерный векторы, составляющие которых представляют значения некоторых координат фазового вектора состояния , вычисленные на ЛА в процессе полета или сигналов, поступающих с гироплатформы;

 

$\delta =\left({\delta }_b^0\mathrm{,}{\delta }_E^0\mathrm{,}{\delta }_H^0\right)$ (2.10)

 

— трехмерный вектор переменных вида (2.5);

 

$\Phi \left(\delta \right)=\left(f\right({\delta }_{b\max }\mathrm{,}{\delta }_b^0),f({\delta }_{E\max }\mathrm{,}{\delta }_E^0),f({\delta }_{H\max }\mathrm{,}{\delta }_H^0\left)\right)$ (2.11)

 

и

$f\left(\sigma \right)=\left(f\right({\alpha }_m\mathrm{,}{\sigma }_b),f({\gamma }_m\mathrm{,}{\sigma }_E),f({\beta }_m\mathrm{,}{\sigma }_H\left)\right)$ (2.12)

— трехмерные векторы допустимых нелинейных функций, имеющих вид (2.5)–(2.6), где

$\sigma =\left({\sigma }_b\mathrm{,}{\sigma }_E\mathrm{,}{\sigma }_H\right)$

— трехмерный вектор сигналов (2.6);

$B=\left(\begin{array}{cccccc}{r}_1& 0& 0& -r_2& 0& 0\\ 0& k_1& 0& 0& -k_2& k_7\\ 0& 0& l_1& 0& l_7& -l_2\end{array}\right)$ (2.14)

 и

$C=\left(\begin{array}{ccc}{r}_3& 0& 0\\ 0& k_6& 0\\ 0& 0& l_6\end{array}\right)$ (2.15)

— постоянные матрицы;

в настоящей работе в численном эксперименте коэффициенты $k_6$ и $l_6$ принимались равными нулю.

Уравнение

$x\text{'}\text{'}=X\left(x\right)+A\left(x\right)\xi \mathrm{,}$ (2.16)

где x — 14-мерный фазовый вектор (1.9) системы (2.7), $\xi =\left({\delta }_b\mathrm{,}{\delta }_E\mathrm{,}{\delta }_H\right)$ — 3-мерный вектор управления, X(x) и A(x) — матрицы-функции с довольно громоздкими коэффициентами — здесь не приводятся [25].

Система (2.7) представляет в настоящей работе систему с прямым перекрестным управлением по отклонениям рулей высоты, элеронов и направления.

3. Численный эксперимент и анализ пространства неисправностей

Мы моделировали движение ЛА, представляющее собой планирующий спуск с высот, близких к орбитальным ( $\cong$100 км), с начальной скоростью, близкой к первой космической скорости. Программное движение определялось заданием программного угла атаки и угла крена.

Численный эксперимент заключался в моделировании различных неисправностей в системе управления (2.5) и (2.6). Был составлен список возможных неисправностей, в соответствии с классификацией неисправностей, приведенной в [26].

Первый список возможных неисправностей включал следующие неисправности, происходящие в первом канале управления (2.5), т. е. в канале управления рулями высоты ${\delta }_b$.

3.1. Априорный список неисправностей №1

Список состоял из 5 следующих неисправностей.

1) Отказ датчика угловой скорости ${\omega }_{s_1}$. Неисправность моделировалась путем обнуления коэффициента $r_1$ в матрице $B$ в момент времени $t=t_0$, $t_0$ — момент возникновения неисправности.

Таким образом, при неисправности 1 из априорного списка неисправностей 1 выполняется условие

$r_1\left(t\right)=0, t\ge t_0\mathrm{.}$ (3.1)

2) Отказ при формировании сигнала ${\sigma }_b$. Моделируется путем обнуления коэффициента $r_3$ в матрице $C$:

$r_3\left(t\right)=0, t\ge t_0\mathrm{.}$ (3.2)

3) Нарушение симметрии функции $f\left({\delta }_{b\max }\mathrm{,}{\delta }_b^0\right)$. Моделируется путем замены $f\left({\delta }_{b\max }\mathrm{,}{\delta }_b^0\right)$ на $f({\delta }_{b\max }\mathrm{,}{\delta }_b^0+{\delta }^{\text{'}})$, т. е. сдвига графика функции $f$ по оси $x$.

4) Заклинивание управляющего органа (руля высоты). Моделируется как ${\delta }_b\left(t\right)={\delta }_b\left(t_0\right)$ при $t\ge t_0$, где $t_0$ — момент возникновения неисправности.

5) Активный отказ управляющего органа (руля высоты). Значение ${\delta }_b$ в момент возникновения неисправности $t_0$ скачком меняется на ${\delta }_{b\max }$ — максимально возможное значение ${\delta }_b$.

Таким образом, при неисправности 5 из априорного списка неисправностей №1 выполняется условие

 ${\delta }_b\left(t\right)={\delta }_{b\max }, t\ge t_0\mathrm{.}$ (3.3)

В процессе полета тяжелого ЛА, уравнения движения которого рассмотрены выше, характерно наличие двух существенно отличных по временным характеристикам движений. Это движения вокруг центра масс с постоянными времени порядка минут.

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

Численное интегрирование исправной и соответствующих неисправных систем (2.7) проводилось с шагом $h=0,8$ с, характерным для движения относительно центра масс тяжелого ЛА с данными уравнениями движения.

3.2. О поверхности контроля

Вектор контроля $y\left(t\right)$ для данной системы состоял из одной компоненты — угла атаки $\alpha$. Множество начальных условий представляло собой сферу $S_{0,1}^{x_0}$ радиуса 0,1 (рад) в пространстве фазовых переменных с центром в точке

$x_0=\left(V_{y_1}^0\mathrm{,}{V}_{y_2}^0\mathrm{,}{V}_{y_3}^0\mathrm{,}{\psi }_c^0\mathrm{,}{\theta }^0\mathrm{,}{r}^0\mathrm{,}{\lambda }^0\mathrm{,}{\phi }^0\mathrm{,}{\omega }_{s_1}^0\mathrm{,}{\omega }_{s_2}^0\mathrm{,}{\omega }_{s_3}^0\mathrm{,}{\alpha }^0\mathrm{,}{\beta }^0\mathrm{,}{\gamma }_c^0\right)$ (3.4)

Здесь используются следующие постоянные:

$\begin{array}{l}{V}_{y_1}^0=7350 км, V_{y_2}^0=0, V_{y_3}^0=0,\\ {\lambda }^0=0, {\phi }^0=45°, {\omega }_{s_1}^0=0, {\omega }_{s_2}^0=0, {\omega }_{s_3}^0=0,\\ {\alpha }^0=0,519 рад, {\gamma }_c^0=0, {\beta }^0=0,\\ r^0=646572 м (h=100 км)\end{array}$ (3.5)

Промежуток времени процесса построения был выбран в следующих пределах:

$[0; 2000 с]$ (3.6)

Для такого множества начальных условий $X^0$, вектора контроля $y\left(0\right)$ и априорного списка неисправностей методом статистических испытаний [27] была получена поверхность контроля ${\pi }_k$ — отрезок $[0,499; 0,539]$ (в радианах). При исправном движении ЛА значение угла атаки находится внутри ${\pi }_k$. Поверхность контроля построена с доверительной вероятностью, не меньшей 0,95 (доказанное утверждение).

3.3. Об обнаружении неисправностей

Путем численного интегрирования системы (7) моделировалось исправное функционирование объекта (полет ЛА), затем возникновение некоторой (j-й) неисправности из априорного списка в момент времени $t_0$ и дальнейшее функционирование вплоть до выхода траектории на поверхность ${\pi }_k$. Времена выхода на ${\pi }_k$ для разных неисправностей из списка приведены в табл. 3.1.

 

Таблица 3.1

Table 3.1

Номер неисправности	1	2	3	4	5
Время выхода на поверхность контроля ${\pi }_k$ (с)	28	118	16	36	2,4

 

По выходе траектории на ${\pi }_k$ включается алгоритм диагностирования с вектором диагностирования $z=\alpha$, т. е. содержащим ту же компоненту фазового вектора, что и вектор контроля.

Характеристиками алгоритма являются время диагностирования $\tau$ и число измерений $N$ (так как численное интегрирование осуществлялось с шагом $h=0,8$ с, то $N$ и $\tau$ связаны соотношением $\tau =Nh$).

Численный эксперимент показал, что для всех номеров $i$ из априорного списка неисправностей алгоритм диагностирования правильно определяет априорные неисправности при $N=3$ ($\tau =2,4$ с). Моделировалось также обнаружение неисправностей с вектором диагностирования

$\overline{z}=\left({\omega }_{s_1}^0\mathrm{,}{\omega }_{s_2}^0\mathrm{,}{\omega }_{s_3}^0\right)$ (3.7)

Численный эксперимент также показал, что для $\overline{z}$ все неисправности из априорного списка определялись однозначно за число измерений $N=3$.

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

Определение 3.1 Диагностическим пространством назовем совокупность датчиков, которой становится в соответствие набор возможных опорных неисправностей $H_j$ с их окрестностями $O_j$, подвергающихся диагностированию посредством опорных невырожденных неисправностей из класса возможных.

Уточним математическую структуру всего диагностического пространства:

$(M;O_1\mathrm{,}{O}_2\mathrm{,}\dots \mathrm{,}{O}_l;A_1\mathrm{,}{A}_2\mathrm{,}{A}_3),$ (3.8)

где M — множество рассмотренных неисправностей $H_1\mathrm{,}\dots \mathrm{,}{H}_l$ вместе с их окрестностями $O_1\mathrm{,}\dots \mathrm{,}{O}_l$. Аксиомы $A_1\mathrm{,}{A}_2\mathrm{,}{A}_3$ определяются следующим образом:

$\begin{array}{l}{A}_1:\forall H_j\in M\exists O_j(H_j\subset O_j);\\ A_2:\forall O_j\exists H_j(H_j\subset O_j);\\ A_3:H_j\subset O_j\cap O_k\Rightarrow \exists O_{\mu }(H_j\subset O_{\mu }\subset O_j\cap O_k\vee O_{\mu }\subset O_j\cup O_k).\end{array}$ (3.9)

Аксиома $A_1$ утверждает, что окрестности $O_j$, являющиеся подмножествами множества $M$, покрывают все $M$, а из аксиомы $A_2$ следует, что эти окрестности не пусты.

Аксиома $A_3$ позволяет обеспечить непрерывный процесс приближения к элементу $H_j$:

$H_j=\underset{\mu \to \infty }{lim}{O}_{\mu }\iff \forall O_{\mu } \exists \overline{M}(\mu >\overline{M}\Rightarrow H_j\subset O_{\mu }),$ (6)

каждая окрестность $O_{\mu }$ которого содержит $H_j$ и “близкие” к $H_j$ непредвиденные и не содержащиеся в рассматриваемом наборе неисправности, которые, если они возникнут в окрестностях $O_j$, надо уметь диагностировать посредством априорного набора неисправностей.

Если же под элементом x понимать не только событие $H_j$, но и непредвиденное событие (не включенное в список $H$ возможных неисправностей), которое может произойти в любой точке $M$ и которое требуется диагностировать посредством $H_j$, то аксиомы математической структуры диагностического пространства (3.8) могут быть записаны в следующем виде:

$A_1:\forall x\in M\exists O_i(x\in O_i),$ (3.11)

то есть окрестности покрывают все $M$;

$A_2:\forall O_i\exists x(x\in O_i)$ (3.12)

то есть окрестности не пусты;

$A_3:x\in O_i\cap O_j\Rightarrow \exists O_k(x\in O_k\subset O_i\cap O_j\vee O_k\subset O_i\cup O_j)$ (3.13)

то есть окрестности можно измельчать и обеспечивать процесс приближения к элементу x.

Предел последовательности $\left\{O_k\right\}$ можно определить как элемент

$x=\underset{k\to \infty }{lim}{O}_k\mathrm{,}$ (3.14)

каждая окрестность $O_j\left(x\right)$ которого содержит $H_j$ и “близкие” к $H_j$ непредвиденные и не содержащиеся в наборе неисправности, которые, если они возникнут в окрестностях $O_j\left(H_j\right)$, надо уметь диагностировать посредством априорного набора неисправностей:

$x=\underset{k\to \infty }{lim}{O}_k\iff \forall O_k\left(x\right)\exists \overline{K}(k>\overline{K}\Rightarrow x\in O_k\mathrm{,}{H}_j\subset O_k\mathrm{,}j\dots \mathrm{,}l)$ (3.15)

Математическая структура диагностического пространства (3.8) позволяет обеспечить ситуацию, при которой решения рассматриваемой динамической системы и системы с неисправностями в системе управления с одинаковыми начальными условиями $x^0$ из некоторого ограниченного пространства будут отличаться друг от друга, а решения рассматриваемой системы с опорной неисправностью $H_j$ и с не предусмотренной априорным списком (непредвиденной) неисправностью из окрестности $O_j$ этой опорной неисправности с одинаковыми начальными условиями $x^0$ будут “близкими”, то есть “мало” отличаться друг от друга, и они могут быть диагностированы как опорные неисправности.

Это во всяком случае будет справедливо для непересекающихся окрестностей $O_j$, $j=1,\dots \mathrm{,}l$, или окрестностей, пересекающихся только вдоль прямых по оси времени.

А вот теперь переходим к неисправностям не из априорного списка.

3.4. Определение неисправностей не из априорного списка

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

Были рассмотрены две неисправности не из списка.

1. Постепенное ухудшение качества показаний датчика угловой скорости ${\omega }_{s_1}^0$ вплоть до полного исчезновения поступающего с него сигнала.

Эта неисправность моделировалась путем уменьшения $r_1$ до нуля по линейному закону:

$r_1\left(t\right)=r_{1N}{c}_1(t-t_0), t>t_0\mathrm{,}$ (3.16)

где $r_{1N}$ — номинальное значение коэффициента $r_1$, $c_1$ — отрицательная константа. Достигнув нуля, значение $r_1$ больше не меняется. Значение $r_1=0$ достигается за время $\tau \cong 30÷40$ с, т. е. после выхода неисправной системы на поверхность контроля. Таким образом, эта неисправность близка к неисправности из априорного списка неисправностей 1, но не совпадает с ней, так как в момент выхода на ${\pi }_k$ коэффициент $r_1$ еще не равен нулю.

2. Постепенный отказ в формировании сигнала ${\sigma }_b$. Моделируется как

$r_3\left(t\right)=r_{3N}{c}_3(t-t_0), t>t_0\mathrm{.}$ (3.17)

Здесь $r_{3N}$ — номинальное значение коэффициента $r_3$, $c_3$ — некоторая отрицательная константа.

Достигнув нуля, величина $r_3$ далее не меняется. Значение $r_3=0$ достигается за время $\tau \cong 140÷150$ с, т. е. после выхода системы на ${\pi }_k$. Эта неисправность близка к неисправности из априорного списка неисправностей №1.

Линейный закон в последних двух случаях 1 и 2 из раздела 3.4 взят в принципе на предварительном этапе. Ну а в дальнейших исследованиях предполагается также и нелинейная аппроксимация рассматриваемых величин.

Обнаружение неисправностей 1 и 2 из раздела 3.4 (т. е. не из априорного списка неисправностей №1) моделировалось следующим образом: возникновение неисправности 1 (неисправности 2) в момент $t_0$, включение алгоритма на выходе на ${\pi }_k$ (время выхода на ${\pi }_k$ для неисправности 1 — 78 с, для неисправности 2 — 224 с), включение алгоритма с одним из векторов диагностирования $z$, $\overline{z}$ и выбор минимума из $S_j^N$, $j=1,\dots \mathrm{,5}$.

Обнаружением рассматриваемой неисправности 1 (неисправности 2) из раздела 3.4 в данном случае является определение случившейся неисправности как неисправности 1 (неисправности 2) из априорного списка неисправностей 1.

Численное моделирование показало, что для $N=5$ ($\tau =4$ с) алгоритм с векторами диагностирования $z$ и $\overline{z}$ правильно обнаруживает рассматриваемые неисправности 1 и 2.

3.5 Априорный список неисправностей №2

Каждая неисправность из этого списка характеризовалась наличием функций $f_j$, $j=1,\dots \mathrm{,17}$, где $f_j$ содержала j-й набор значений коэффициентов $r_1\mathrm{,}{r}_3\mathrm{,}{k}_1\mathrm{,}{k}_2\mathrm{,}{k}_3\mathrm{,}{l}_1\mathrm{,}{l}_2\mathrm{,}{l}_3$ в цепях формирования сигналов ${\delta }_u^0$ из (2.5). Различные комбинации приведены в табл. 3.2. В ней индекс $Nom$ означает номинальное значение коэффициента. Из таблицы видно, что неисправности с номерами 1–3 относятся к первому каналу управления из (2.5), формирующему величину ${\delta }_b$, с номерами 4–10 — ко второму (формирующему величину ${\delta }_E$), а с номерами 11–17 — к третьему (формирующему величину ${\delta }_H$).

По этим опорным неисправностям проводилось затем определение несписочных неисправностей в каналах управления. В этом случае распознавался только номер канала: при номере i минимального из чисел $S_j^N$, $j=1,\dots \mathrm{,17}$, равном 1, 2 или 3, неисправность определялась как случившаяся в первом канале управления (${\delta }_b$), при номере $i=4,\dots \mathrm{,10}$, — во втором (${\delta }_E$), при номере $i=11,\dots \mathrm{,17}$, — в третьем (${\delta }_H$).

 

Таблица 3.2

Table 3.2

N	$r_1$	$r_2$	$r_3$	$r_4$	$r_5$	$r_6$	$r_7$	$r_8$
1	0	$Nom$
2	$Nom$	0
3	0	0
4			0	$Nom$	$Nom$
5			$Nom$	0	$Nom$
6			$Nom$	$Nom$	0
7			0	0	$Nom$
8			$Nom$	0	0
9			0	$Nom$	0
10			0	0	0
11						0	$Nom$	$Nom$
12						$Nom$	0	$Nom$
13						$Nom$	$Nom$	0
14						0	0	$Nom$
15						$Nom$	0	0
16						0	$Nom$	0
17						0	0	0

 

Поверхность контроля ${\pi }_k$ в данном случае при векторе контроля $y=\alpha$ и множестве начальных условий была получена такой же, как и для априорного списка неисправностей №1 (множество начальных условий такое же, как и для априорного списка неисправностей №1):

${\pi }_k:[0,499; 0,539].$ (3.18)

Численное моделирование показало, что диагностирование неисправностей с векторами диагностирования $z=\alpha$ и $\overline{z}=\left({\omega }_{s_1}^0\mathrm{,}{\omega }_{s_2}^0\mathrm{,}{\omega }_{s_3}^0\right)$ за число измерений $N=8$ ($\tau =6,4$ с) привело к правильному определению номера канала, в котором произошла неисправность. Набор несписочных неисправностей, для которых проводилось определение номера неисправного канала, приведен в табл. 3.3.

Номинальные значения коэффициентов в цепях управления (2.7) приведены в табл. 2.1.

При этом были выбраны следующие значения постоянных:

$r_1=20, r_3=10, k_1=7, k_2=1,$ (3.19)

$k_3=0,3, l_1=7, l_2=1, l_3=0,7,$ (3.20)

 

Таблица 3.3

Table 3.3

Номер канала	Коэффициенты	в цепи	управления
1	$r_1=5$	$r_3=10$
1	$r_1=1$	$r_3=5$
1	$r_1=0$	$r_3=2$
2	$k_1=0,5$	$k_2=1$	$k_3=0,3$
2	$k_1=1$	$k_2=0$	$k_3=0$
2	$k_1=3$	$k_2=0$	$k_3=1$
3	$l_1=3$	$l_2=1$	$l_3=0,7$
3	$l_1=1$	$l_2=0$	$l_3=0,3$
3	$l_1=1$	$l_2=0$	$l_3=0$

 

3.6. Обнаружение неисправностей без применения метода поверхности контроля

Возможно обнаружение неисправностей по алгоритму диагностирования без построения самой поверхности контроля ${\pi }_k$. При этом:

1. Под номером 0 в априорный список неисправностей вносится исправная система.
2. Алгоритм диагностирования в принципе включается циклически, с некоторым интервалом времени $\Delta t$.
3. Если обнаружена так называемая “неисправность” под номером 0, то есть система исправна, продолжается функционирование объекта до момента нового включения алгоритма диагностирования.
4. Если обнаружена неисправность с номером $i\ne 0$, выдается сообщение о наличии этой неисправности.

Численное моделирование диагностики с циклическим включением алгоритма показало, что при интервалах включения алгоритма $\Delta t=10,\mathrm{20,30}$ с все вышеперечисленные неисправности из априорного списка неисправностей 1 и априорного списка неисправностей 2 были правильно определены.

Алгоритм диагностирования работал с векторами диагностирования $z=\alpha$ и $\overline{z}=\left({\omega }_{s_1}\mathrm{,}{\omega }_{s_2}\mathrm{,}{\omega }_{s_3}\right)$ и числом измерений $N=5$ ($\tau =4$ с).

Таким образом, с помощью предлагаемого алгоритма диагностирования численным экспериментом показана возможность диагностирования неисправностей датчиков управляющих сигналов, формирующих СУ движением ЛА и, в частности, датчиков управляющих сигналов с гиростабилизированной платформы [28–30].

Таблица 3.2 представляет собой в некотором роде опорные неисправности. В этой таблице неисправности 5, 7, 8, 10, 12, 14, 15, 17 формируются, в частности, отказами датчиков сигналов управления с гиростабилизированной платформы. Эти неисправности правильно диагностируются [31; 32].

По опорным неисправностям (табл. 3.2) проводилось диагностирование не списочных неисправностей (табл. 3.3) в каналах управления движением ЛА (неисправности 5, 6, 8, 9 формируются в табл. 3.3 и с помощью отказов датчиков сигналов управления с гиростабилизированной платформы). В этом случае распознавался только номер канала управления движением летательного аппарата. Численное моделирование показало, что диагностирование неисправностей за вполне приемлемое время привело к правильному определению номера сигнала управления, в котором произошла неисправность.

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

4. Диагностика в условиях измерения части фазового вектора

Для практического применения алгоритма диагностирования неисправностей, как показано в предыдущем разделе о численном эксперименте, требуется знание начальных условий для всего 14-мерного фазового вектора состояния x. Это несколько затрудняет возможность практического применения алгоритма диагностирования неисправностей.

Из системы уравнений (7) может быть выделена динамическая подсистема из 3-х уравнений (рассматриваемая на алгебре Ли $so\left(3\right)$) относительно угловых скоростей ${\omega }_{s_i}\mathrm{,}i=1,2,3$. В нормальном виде она может быть представлена как

$\begin{array}{l}{\omega }_{s_1}=\frac{\rho V_e^2}{2}\frac{S{b}_a}{I_{s_1}}\left(m_{s_1}^{\alpha }+m_{s_1}^{{\delta }_e}{\delta }_e+m_{s_1}^{{\omega }_{s_1}}\frac{b_a}{V_e}{\omega }_{s_1}\right)+\frac{I_{s_2}-I_{s_3}}{I_{s_1}}{\omega }_{s_2}{\omega }_{s_3}+\frac{M_{s_1}^{\partial }}{I_{s_1}}\mathrm{,}\\ \\ {\omega }_{s_2}=\frac{\rho V_e^2}{2}\frac{SL}{I_{s_2}}\left(m_{s_2}^{\beta }(\beta -{\beta }_e)+m_{s_2}^{{\delta }_E}{\delta }_E+m_{s_2}^{{\delta }_H}{\delta }_H+m_{s_2}^{{\omega }_{s_2}}\frac{L}{2V_e}{\omega }_{s_2}\right)+\\ +\frac{I_{s_3}-I_{s_1}}{I_{s_2}}{\omega }_{s_1}{\omega }_{s_3}+\frac{M_{s_2}^{\partial }}{I_{s_2}}\mathrm{,}\\ {\omega }_{s_3}=\frac{\rho V_e^2}{2}\frac{SL}{I_{s_3}}\left(m_{s_3}^{\beta }(\beta -{\beta }_e)+m_{s_3}^{{\delta }_H}{\delta }_H+m_{s_3}^{{\omega }_{s_2}}\frac{L}{2V_e}{\omega }_{s_2}+m_{s_3}^{{\omega }_{s_3}}\frac{L}{2V_e}{\omega }_{s_3}\right)+\\ +\frac{I_{s_1}-I_{s_2}}{I_{s_3}}{\omega }_{s_1}{\omega }_{s_2}+\frac{M_{s_3}^{\partial }}{I_{s_3}}\mathrm{.}\end{array}$ (4.1)

Эти уравнения имеют место при совпадении главных осей инерции с так называемыми “строительными осями”, т. е. при тензоре инерции следующего вида:

$I=\left(\begin{array}{ccc}{I}_{s_1}& 0& 0\\ 0& I_{s_2}& 0\\ 0& 0& I_{s_3}\end{array}\right)\mathrm{.}$ (4.2)

Здесь $\rho \mathrm{,}L\mathrm{,}S\mathrm{,}{I}_{s_i}$ — вполне определенные физические константы, $m_{s_i}$ — медленно меняющиеся коэффициенты аэродинамических моментов, которые можно считать постоянными на времени работы алгоритма диагностирования. Величины ${\omega }_{s_i}$ наблюдаются и, таким образом, начальные условия для уравнения (1) известны. Величины ${\delta }_e\mathrm{,}{\delta }_E\mathrm{,}{\delta }_H$ известны в любой момент времени, так как это — формируемое управление.

Таким образом, измеряя величины $V_e$, $\beta -{\beta }_e$, присутствующие в правой части, можно замкнуть систему уравнений (1) и численно ее интегрировать на некотором промежутке времени (диагностирования) $\left[t_0\mathrm{,}t\right]$ с начальными условиями ${\omega }_{s_i}{t}_0$.

Численный эксперимент представлял собой диагностирование неисправностей из априорного списка неисправностей 1 в условиях неточных измерений величин $V_e$, $\beta -{\beta }_e$ (измеренные значения отличались (при каждом измерении) на $5÷10\%$ от действительных). Измеряемый вектор $z\left(t\right)=\left({\omega }_{s_1}\mathrm{,}{\omega }_{s_2}\mathrm{,}{\omega }_{s_3}\right)$, компоненты расчетного вектора $z_j$ получались путем численного интегрирования системы (4.1). Алгоритм диагностирования работал циклически (под 0 в него была включена исправная система) с интервалом включения $15÷20$ с и правильно определял происшедшее в системе управления ЛА неисправности за $15÷20$ измерений ($8÷14$ с).

Заключение

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

Об авторах

Максим Владимирович Шамолин

Институт механики Московского государственного университета имени М.В. Ломоносова

Автор, ответственный за переписку.
Email: shamolin@rambler.ru
ORCID iD: 0000-0002-9534-0213

доктор физико-математических наук, профессор, ведущий научный сотрудник Института механики, эксперт РАН

Россия, г. Москва

Список литературы

Дополнительные файлы