Синтез оптимального управления осесимметричным вращающимся летательным аппаратом с использованием модели спирального прогноза

Полный текст

Введение. Атмосфера оказывает существенное влияние на сельскохозяйственную деятельность человека. Одно из наиболее опасных явлений погоды – град, который уничтожает овощные и садовые культуры, виноград, побеги хлопчатника и др. Нередко от него погибают домашние птицы и мелкий рогатый скот. Для борьбы с градообразующими кучево-дождевыми облаками применяются специальные летательные аппараты (ЛА), распыляющие реагент (PbJ₂, ArJ), который вызывает кристаллизацию переохлажденных водяных капель. Занос реагента осуществляется в крупнокапельную зону, где происходит интенсивный рост градин [1]. Исходные данные для пусков ЛА готовятся по показаниям метеорологического локатора, например МРЛ-5 [2], и используются для выбора точек выброса реагента.

Со второй половины 50-х гг. XX в. для борьбы с градом в нашей стране, наряду с другими комплексами, применяются противоградовые ЛА серии “Эльбрус” [3], представляющие собой осесимметричные неуправляемые быстровращающиеся ЛА с активным стартом.

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

Для управления быстровращающимся ЛА можно использовать управляющую насадку на головную часть, оснащенную вращающимся блоком рулей с постоянным углом установки [4]. Общая схема применения ЛА с управляющей насадкой подразумевает перед пуском ввод данных о точке выброса реагента при помощи специального установщика. На траектории ЛА получает информацию о своем местоположении от космической радионавигационной системы (КРНС) и производит прогноз траектории. В зависимости от величины промаха спрогнозированной траектории относительно точки выброса реагента вырабатывается управляющая команда. Управление заключается в установке необходимого для устранения промаха углового положения блока рулей и сообщения последнему угловой скорости, равной по модулю и противоположной по направлению угловой скорости вращения ЛА. Таким образом, положение рулевых лопастей в полусвязанной с ЛА системе координат будет оставаться неизменным и последний будет совершать маневр. Вопросам движения ЛА с подобными насадками посвящен ряд публикаций в зарубежной печати [5–9], однако в них не затрагивается вопрос оптимального управления. В работе [10] рассмотрено приложение аналитического решения задачи о вращательном движении оси ЛА к проектированию систем коррекции, а в [11] исследован вопрос о движении ЛА с системой коррекции с одной степенью свободы и получены оценки областей достижимости.

В этой связи синтез оптимального управления вращающимися ЛА является актуальной и востребованной практикой научной задачей. Теория оптимального управления редко применяется на практике. И это связано не столько с объемом вычислений, сколько с необходимостью выбора начального приближения, обеспечивающего сходимость итерационных процедур нахождения управления. Отказ от оптимизации приводит к понижению роли выбора управления и использованию субъективных решений. Как правило, требуется найти управление сложной системой, да еще в реальном времени. Непосредственное применение теории и ее алгоритмов в реальном времени затруднительно. Однако к настоящему времени имеются обнадеживающие решения, позволяющие приблизить к повседневной практике идеи оптимизации. Это введение в рассмотрение академиком А. А. Красовским полуопределенного функционала обобщенной работы и использование для вычисления управления прогнозирующих моделей [12]. В результате двухточечная краевая задача сводится к двум задачам Коши (в прямом и обратном времени соответственно). Алгоритм реализуем в реальном времени. Простота реализации сопровождается, правда, потерей в терминальной точности, поскольку она определяется значениями весовых коэффициентов в целевом функционале. В ряде случаев удается повысить точность путем рассмотрения иерархии критериев А. А. Красовского и на основе упрощения полученного решения применить алгоритм последовательной оптимизации [13, 14]. Другим подходом является рассмотрение вспомогательной задачи оптимизации, облегчающей решение основной исходной задачи. В этом случае из решения основной задачи по принципу максимума формируется структура управления, коррекция параметров которой может производиться, в частности, из условия минимума как функционала с дополнительной интегральной частью, включающей затраты на управление, с применением принципа максимума, так и критерия А. А. Красовского с использованием алгоритма с прогнозирующей моделью. “Коль скоро структура оптимального управления тем или иным способом угадана, не так трудно установить, как правило, что она действительно такова”, – писал Р. Беллман [15]. Причем структура может быть выбрана с избытком, тогда в процессе вычислений происходит ее упрощение [16]. Этот алгоритм с успехом применялся для сложных задач управления: мостовым краном, автомобилем [17], автоматическим подводным аппаратом при наличии участка особого управления [18], для раскрытия трансформируемой антенны космического базирования, включающего разведение спиц, их выдвижение и натяжение фронтальной поверхности радиоотражающего сетеполотна с учетом колебательных процессов [19] и др. Была доказана теорема об эквивалентности основной и вспомогательной задач оптимизации [20].

Теперь о сложности реальных процессов по сравнению с их математическим описанием. Для управления сложными процессами целесообразно использовать самоорганизующиеся оптимальные регуляторы с экстраполяцией (СОРЭ) А. А. Красовского [21, 22]. СОРЭ не требуют знания математической модели объекта, по результатам наблюдений сами формируют информационную модель. Такие алгоритмы апробированы для ряда технических объектов, в том числе для ЛА [13, 23], ускорителя электронов на основе высоковольтного тлеющего разряда [24], автомобиля [17] и др., а также для социально-экономических систем, гипотетически описываемых производственной функцией, моделью развития экономики В. В. Леонтьева [25], расчета стоимости опциона на покупку (продажу) акций [17], поддержки управленческих решений вузовской администрации [25]. Исследования показали, что СОРЭ можно использовать как непосредственно для управления процессом, так и для стабилизации его динамики на оптимальной траектории его математической модели, а также как компьютерные тренажеры при подготовке администраторов и в качестве компьютерного помощника в процессе принятия решения по управлению.

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

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

Далее исследуется управляемое пространственное движение осесимметричного вращающегося ЛА постоянной массы. Рассматривается терминальная задача перевода ЛА из произвольно заданного начального состояния в фиксированное конечное состояние.

Уравнения поступательного движения запишем в проекциях на оси траекторной системы координат (СК), а уравнения вращательного движения, чтобы избежать переменности моментов инерции, на оси полусвязанной СК [27–29]:

$\dot{V}=g\left(n_x-sin\theta \right)$,

$\dot{\theta }=\frac{g}{V}\left(n_y-cos\theta \right)$,

$\dot{\Psi }=-\frac{g}{Vcos\theta }{n}_z$,

${\dot{\omega }}_x=\frac{1}{I_x}\left[M_{x\mathrm{пс}}-{\omega }_{y\mathrm{пс}}{\omega }_{z\mathrm{пс}}\left(I_z-I_y\right)\right]$,

${\dot{\omega }}_{y\mathrm{пс}}=\frac{1}{I_y}\left[M_{y\mathrm{пс}}-{\omega }_{z\mathrm{пс}}{\omega }_x\left(I_x-I_z\right)-I_z{\omega }_{z\mathrm{пс}}\left({\omega }_x-{\omega }_{y\mathrm{пс}}tg\vartheta \right)\right]$,

${\dot{\omega }}_{z\mathrm{пс}}=\frac{1}{I_z}\left[M_{z\mathrm{пс}}-{\omega }_{y\mathrm{пс}}{\omega }_x\left(I_y-I_x\right)+I_y{\omega }_{y\mathrm{пс}}\left({\omega }_x-{\omega }_{y\mathrm{пс}}tg\vartheta \right)\right]$, (1.1)

$\dot{\gamma }={\omega }_x-{\omega }_{y\mathrm{пс}}tg\vartheta$,

$\dot{\psi }={\omega }_{y\mathrm{пс}}/cos\vartheta$,

$\dot{\vartheta }={\omega }_{z\mathrm{пс}}$,

$\dot{x}=Vcos\theta cos\Psi$,

$\dot{y}=Vsin\theta$,

$\dot{z}=-Vcos\theta sin\Psi$,

$n_x=A_{x\text{\hspace{0.17em}т}}/\left(mg\right)$, $n_y=A_{y\text{\hspace{0.17em}т}}/\left(mg\right)$, $n_z=A_{z\mathrm{т}}/\left(mg\right)$, $A_{x\text{\hspace{0.17em}т}}=\left(-c_x\left(M,\delta \right)-c_{x\text{\hspace{0.17em}р}}\left(M\right)S_{р}/S\right)qS$, $A_{y\mathrm{т}}=\left(c_y^{\delta }\left(M,\delta \right)\alpha +c_n^{\delta }\left(M\right){\delta }_{\mathrm{р}}sin{\gamma }_{\mathrm{р}}{S}_{\mathrm{р}}/S\right)qS$, $A_{z\mathrm{т}}=\left(-c_y^{\delta }\left(M,\delta \right)\beta -c_n^{\delta }\left(M\right){\delta }_{\mathrm{р}}cos{\gamma }_{\mathrm{р}}{S}_{р}/S\right)qS$,

$M_{xпс}=\left(m_x^{{\omega }_x}\left(M\right)\frac{\left({\omega }_x-{\omega }_{y\text{\hspace{0.17em}пс}}tg\vartheta \right)L}{V}+m_x\left(M\right)S_{р}/S\right)qSL$, $M_{y\text{\hspace{0.17em}пс}}=\left(m_z^{\delta }\left(M,\delta \right)\beta +m_z^{{\omega }_z}\left(M,\delta \right)\frac{{\omega }_{y\text{\hspace{0.17em}пс}}L}{V}+c_n^{\delta }\left(M\right){\delta }_{р}\frac{l_{р}}{L}\frac{S_{р}}{S}cos{\gamma }_{р}\right)qSL$, $M_{z\text{\hspace{0.17em}пс}}=\left(-m_z^{\delta }\left(M,\delta \right)\alpha +m_z^{{\omega }_z}\left(M,\delta \right)\frac{{\omega }_{z\text{\hspace{0.17em}пс}}L}{V}+c_n^{\delta }\left(M\right){\delta }_{р}\frac{l_{р}}{L}\frac{S_{р}}{S}sin{\gamma }_{р}\right)qSL$.

Здесь $n={\left(n_x{n}_y{n}_z\right)}^{Т}$ – вектор перегрузок с проекциями на оси траекторной СК, m – масса ЛА, I_x, I_y I_z – осевой и экваториальные моменты инерции ЛА; $V=\sqrt{V_x^2+V_y^2+V_z^2}$ – модуль вектора скорости с проекциями на оси связанной СК, $\theta$ – угол наклона траектории; $\gamma$– угол пути, $A_{\mathrm{т}}={\left(A_{x\mathrm{т}}{A}_{y\mathrm{т}}{A}_{z\mathrm{т}}\right)}^{Т}$– вектор аэродинамических сил в траекторной системе координат; $M_{пс}={\left(M_{xпс}{M}_{yпс}{M}_{zпс}\right)}^{Т}$– вектор аэродинамических моментов в полусвязанной системе координат; g – ускорение силы тяжести; ${\omega }_{пс}={\left({\omega }_{yпс}tg\vartheta {\omega }_{yпс}{\omega }_{zпс}\right)}^{Т}$ – вектор угловой скорости вращения в полусвязанной системе координат; $\psi$, $\vartheta$, $\gamma$ – углы рыскания, тангажа и крена ЛА соответственно; $\delta$ – угол нутации ЛА, x, y, z – координаты центра масс ЛА в стартовой системе координат, радиус Земли принимается равным радиусу сферы, аппроксимирующему геоид, R_з = 6356766 м; $q=\rho V^2/2$ – скоростной напор, M – число Маха, $\rho$ – плотность воздуха, S – площадь миделя, $S_{р}$– площадь рулевой поверхности, L – длина ЛА; аэродинамические коэффициенты: $c_{\text{x}}$– лобового сопротивления, $c_y^{\delta }$и $c_n^{\delta }$ – производных коэффициентов подъемной силы и нормальной силы рулевой поверхности по углу атаки руля; $m_z^{\delta }$и $m_z^{{\omega }_z}$ – производная коэффициента опрокидывающего момента по углу нутации и коэффициент демпфирующего момента. Из главных векторов аэродинамических сил и моментов исключены сила и момент Магнуса, поскольку их вклад в движение ЛА не превышает ошибок, вызванных приближенным определением основных аэродинамических коэффициентов. Углы атаки $\alpha$ и скольжения $\beta$, в силу их малости для устойчивых быстровращающихся ЛА, а также малости углов пути и рыскания, определяются выражениями: $\alpha =\vartheta -\theta$, $\beta =\left(\psi -\Psi \right)/\cos \theta$.

Блок рулей характеризуется площадью рулевой поверхности S_p, продольной координатой, отсчитываемой от центра масс ЛА l_p, и углом установки рулей ${\delta }_p$. За управление принимается угол наклона блока рулей ${\gamma }_p$ от вертикальной плоскости (рис. 1), принимающий значения от 0 до 2p.

 

Рис. 1. Схема ЛА с системой коррекции вращающимся блоком рулей. ЦМ – центр масс.

 

Система уравнений (1.1) может быть получена из уравнений (4.56)–(4.64) из работы [28], если положить в них массу ЛА постоянной, Землю в пределах траектории плоской и не вращающейся, а угол крена принять равным нулю.

Блок рулей кинематически развязан по вращению с корпусом ЛА по осевому вращению и является ротором бесколлекторного электродвигателя постоянного тока. На участке коррекции, имея заданный системой управления угол ${\gamma }_p$, по данным датчиков угловой скорости и положения блока рулей, управляющий контроллер электродвигателя изменяет напряжение питания таким образом, чтобы угловая скорость вращения блока рулей равнялась угловой скорости вращения ЛА, а угол наклона блока рулей соответствовал бы заданному. Равенство угловых скоростей вращения ЛА и блока рулей приводит к практически неизменному положению последнего в пространстве.

Таким образом, состояние ЛА описывается 12 переменными, из которых система управления может оказывать целенаправленное влияние на линейные и угловые ускорения, причем ускорения вдоль оси Х появляются на участке управления, зависят только от факта управления и не зависят от других параметров, а линейные и угловые ускорения относительно полусвязанных осей Y и Z оказываются функционально связанными величинами, т. е. некоторому линейному ускорению однозначно соответствует угловое.

При полете малого осесимметричного вращающегося ЛА, стремящегося достичь зоны градообразования внутри облака, надо использовать алгоритмы, формирующие управление в реальном времени. Поскольку условия применения таких ЛА изменчивы в зависимости от явлений погоды, то реализация программных траекторий не продуктивна. Задачу управления можно решать по-разному. В работе [30] с помощью концепции обратных задач динамики по заданной достигающей цели траектории полета ЛА (как в вертикальной плоскости, так и в пространстве) в виде степенного полинома от продольной дальности, коэффициенты которого определяются из решения краевой задачи, вычисляется значение перегрузки в текущий момент времени. Это решение затруднительно использовать при необходимости ситуационного облета запретной области. Уравнения (1.1) примем в качестве имитационной модели.

В данной статье предложено формировать оптимальное управление с использованием алгоритма спирального прогноза [26], который допускает реализацию в реальном времени.

2. Основной результат. Для целей управления земную нормальную систему координат будем считать инерциальной и рассмотрим уравнения динамики ЛА в связанных осях [26]:

$\dot{v}=\Omega v+g\left(n-{\epsilon }_2\right)$, (2.1)

${\dot{\epsilon }}^T=\Omega {\epsilon }^T$, (2.2)

${\left(\dot{\phi }\dot{h}\dot{\lambda }\right)}^T=d^{*}\epsilon v$,

где

$\Omega =\left[\begin{array}{ccc}0& {\omega }_z& -{\omega }_y\\ {\omega }_z& 0& {\omega }_x\\ {\omega }_y& -{\omega }_x& 0\end{array}\right]$, ${\epsilon }^T=\left[\begin{array}{ccc}{\epsilon }_{11}& {\epsilon }_{21}& {\epsilon }_{31}\\ {\epsilon }_{12}& {\epsilon }_{22}& {\epsilon }_{32}\\ {\epsilon }_{13}& {\epsilon }_{23}& {\epsilon }_{33}\end{array}\right]$,

$d^{\ast }=diag\left(\frac{1}{R_3+h},\mathrm{1,}\frac{1}{\left(R_3+h\right)cos\varphi }\right)$, ${\epsilon }_2={\left({\epsilon }_{21}{\epsilon }_{22}{\epsilon }_{23}\right)}^T$,

${\epsilon }_{11}=cos\psi cos\vartheta$, ${\epsilon }_{21}=sin\upsilon$, ${\epsilon }_{31}=-sin\psi cos\vartheta$,

${\epsilon }_{12}=sin\psi sin\gamma -cos\psi sin\vartheta cos\gamma$, ${\epsilon }_{22}=cos\vartheta cos\gamma$,

${\epsilon }_{32}=cos\psi sin\gamma +sin\psi sin\vartheta cos\gamma$, ${\epsilon }_{13}=sin\psi \cos \gamma +\cos \psi sin\vartheta sin\gamma$,

${\epsilon }_{23}=-cos\vartheta sin\gamma$, ${\epsilon }_{33}=cos\psi \cos \gamma -sin\psi sin\vartheta sin\gamma$,

$\epsilon v=v_g$, $v={\left(v_x{v}_y{v}_z\right)}^T$– скорость ЛА с проекциями на оси связанной СК, $v_g$ – скорость ЛА в земной СК, $\varphi$, $\lambda$ – широта и долгота ЛА, h – высота полета, R_з – радиус Земли.

Среди достоинств такой формы описания кинематики углового движения можно отметить отсутствие особых точек, т. е. моделируемые угловые движения не имеют ограничений. Избыточность вычисляемых параметров (девять вместо трех) частично компенсируется, если одну из строк матрицы $\epsilon =D_{н}^{св}$определять через алгебраические дополнения.

Управление зададим в виде ${\dot{y}}_i=u_i$, $i=\overline{\mathrm{1,5}}$, $y={\left(n_y{n}_z{\omega }_y{\omega }_z{t}_f\right)}^{Т}$, ${\dot{n}}_y=u_1$, ${\dot{n}}_z=u_2$, ${\dot{\omega }}_y=u_3$, ${\dot{\omega }}_z=u_4$, ${\dot{t}}_f=u_5$. Здесь $n_x$, $n_y$, $n_z$ – компоненты вектора перегрузки $n$в связанных осях. Уравнения (2.1), (2.2) имеют аналитическое решение [26, 31].

Задача оптимизации системы (2.1), (2.2) по классическому критерию по принципу максимума с решением двухточечной краевой задачи методом Ньютона продемонстрировала вычислительные затруднения [32], связанные с заданием начальных условий для сопряженных переменных. И хотя удалось выявить локальные зоны сходимости, рекомендовать его для реализации в реальном времени преждевременно, даже с учетом концепции “гибких кинематических траекторий” [33]. Поэтому для разработки алгоритма, реализуемого в реальном времени, рассмотрим целевой функционал в виде критерия Красовского:

$I=V_f\left(x,t_f\right)+\underset{t_0}{\overset{t_f}{\int }}{f}_0\left(x,u,u_o,t\right)dt$, (2.3)

где x – n-мерный вектор состояния, u – m-мерный вектор управления, $u_o$ – оптимальное значение вектора u, $f_0$, $V_f$ – заданные положительно-определенные функции своих аргументов, имеющие непрерывные частные производные по x, t, а функция $f_0$ – еще и по u. Примем $V_f\left(x,t_f\right)=\mathrm{0,5}{\left(\Delta l_f\Delta h_f\Delta z_f\right)}^T\rho \left(\Delta l_f\Delta h_f\Delta z_f\right)$, $x_0=x\left(t_0\right)$, $\Delta l_f=l\left(t_f\right)-l_f$, $\Delta h_f=h\left(t_f\right)-h_f$, $\Delta z_f=z\left(t_f\right)-z_f$, $f_0\left(x,u,t\right)=Q\left(x,t\right)+0.5\left(u^{\text{T}}{K}^{-1}u+u_{\text{o}}^{\text{T}}{K}^{-1}{u}_{\text{o}}\right)$, $K=diag\left(k_1^2,k_2^2,\dots ,k_5^2\right)$, $\rho =diag\left({\rho }_1,{\rho }_2,{\rho }_3\right)$; $l_f$, $h_f$, $z_f$ – заданные конечные значения координат точки доставки полезной нагрузки.

В соответствии с алгоритмом с прогнозирующей моделью (на прогнозе $u_i=0$, $i=\overline{\mathrm{1,5}}$) управления формируются в виде $u_i=-k_i^2\left(\partial V/\partial y_i\right)$. При $Q\left(x,t\right)=0$ уравнение Ляпунова $\dot{V}=-Q\left(x,t\right)$ принимает вид $\dot{V}=0$ и $\partial V/\partial y_i=\partial V_f/\partial y_i$, отсюда $u_i=-k_i^2\left(\partial V_f/\partial y_i\right)$. Учет условий полета на всем интервале оптимизации (ограничения по высоте, перегрузкам, аэродинамическим углам и др.) может осуществляться посредством соответствующего выбора подынтегральной функции $Q\left(x,t\right)$. При $Q\left(x,t\right)\ne 0$путем введения дополнительной переменной ${\dot{x}}_{n+1}=Q\left(x,t\right)$имеем $V_{fQ}=V_f+x_{n+1}\left(t_f\right)$, и задача сводится к предыдущему случаю. Если ограничения действуют на ограниченных временных интервалах (например, включением в $Q\left(x,t\right)$функции штрафа за превышения некоторыми компонентами вектора состояния предельных значений), то приходим к алгоритму с прогнозирующей моделью при нулевом управлении. При этом в прямом времени прогноз производится по аналитическим формулам спирального прогноза, а при численном обратном прогнозировании вдоль модели (2.1), (2.2) добавляются уравнения для сопряженных переменных [34].

Алгоритмы отработки оптимальных траекторных задающих воздействий могут строиться на основе различных систем стабилизации [27]. Если маневрирование ЛА не является интенсивным, то могут применяться оптимальные регуляторы. Для случая быстрых интенсивных маневров можно использовать алгоритмы обратных задач динамики. Продуктивным, особенно при изменении аэродинамических характеристик ЛА в полете, является применение СОРЭ.

Если выбор порядка информационного полинома, формируемого по дискретно получаемым наблюдениям (самоорганизация информационной модели), происходит периодически на заданном интервале времени [0, t_c] (цикле), то необходимая для вычисления управления оценка вектора x состояния внутри этого цикла (в данном случае это вектор выходных параметров и их производных) производится с привлечением текущих измерений с помощью фильтра Калмана. В работе [35] рекомендован выбор соответствующих начальных значений элементов матрицы ковариаций фильтра Калмана, где, в частности, приведены числовые значения этой матрицы при ее порядке от 3 по 7. На борту хранятся в памяти значения элементов ковариационных матриц для разных порядков полиномиальной модели. Это симметричная матрица квазидиагонального вида, диагональные и соседние с ними элементы которой компактно можно представить следующим образом [17]:

$P_{l-\mathrm{1,}l-1}={\left[\left(l-i-1\right)!\right]}^2{C}_{2(l-1)}^{2(l-1-i)},l\ge \mathrm{1,}i=\overline{\mathrm{0,}l-\mathrm{1,}}$

$P_{l-\mathrm{1,}l-i+1}=-\frac{1}{2}\left(l-i-1\right)!\left(l-i\right)!C_{2(l-1)}^{2(l-1-i)+1},l\ge \mathrm{2,}i=\overline{\mathrm{1,}l-\mathrm{1,}}$

а все остальные элементы равны нулю.

В данной работе для стабилизации в целях упрощения вместо СОРЭ будет применен пропорционально-интегрально-дифференциальный (ПИД) регулятор

$u=r_{1}e+r_2\underset{0}{\overset{t}{\int }}ed\tau +r_3\dot{e}$.

Если решать обратную задачу оптимизации с целью получения целевого функционала для управления в виде ПИД-регулятора, то оказывается, что надо использовать СОРЭ. В работе [36] показано, что ПИД-регулятор является частным случаем СОРЭ. При этом измеряется интеграл

$z\left(t\right)=\underset{0}{\overset{t}{\int }}ed\tau$

от рассогласования е в трех точках, а информационный полином выбирается без самоорганизации аппроксимирующим в виде параболы. Оптимальное по критерию А. А. Красовского (2.3) ($V_f=0.5{\beta }_1{x}_1^2\left(t_c\right)$, $f_o=0.5{\beta }_2{x}_1^2$, $t_0=0$, $t_f=t_c$, $x_1=z$) управление формируется регулятором в виде $u=-K{p}_{l+1}$, где ${\beta }_1$, ${\beta }_2$, $t_c$заданы, а сопряженные переменные определяются в виде

$p_j=\frac{t_c^{j-1}}{\left(j-1\right)!}\sum _{i=0}^l{a}_i{t}_c^i\left({\beta }_1+\frac{{\beta }_2{t}_c}{l+i}\right)$,

полученном методом математической индукции [13]. Здесь u – управление, где l – порядок полинома (l=2, так как используется парабола); ${\beta }_{1,} {\beta }_2, K$ – весовые коэффициенты в целевом функционале, – интервал оптимизации, $a_i, i=\overline{0,2}$ – коэффициенты полинома $z\left(t\right)=a_2{t}^2+a_{1}t+a_0$. В этом частном случае структура управления совпадает со структурой ПИД-регулятора, а коэффициенты ПИД-регулятора вычисляются однозначно через коэффициенты функционала и величину интервала оптимизации:

$r_1=\frac{1}{24}K{\beta }_2{t}_c^4$, $r_2=-K\frac{t_c^2}{2}\left({\beta }_1+\frac{{\beta }_2{t}_c}{3}\right)$, $r_3=-\frac{1}{120}K{\beta }_2{t}_c^5$.

3. Моделирование. При моделировании необходимые для определения управления частные производные $\partial V_f/\partial y_i$ в алгоритме со спиральным прогнозом определяются численно, и управления вычисляются в виде $u_i=-k_i^2\left(\Delta V_f/\Delta y_i\right)$, $i=\overline{\mathrm{1,5}}$. Приращения $\Delta y$ принимались такими: $\Delta y={\left(\Delta n_x\Delta n_y\Delta n_z\Delta {\omega }_y\Delta {\omega }_z\Delta t_f\right)}^T=$${\left(0.20.20.010.0050.1\right)}^T$.

В расчетах для упрощения на вход имитационной модели через инерционное звено подается управление в виде ПД-регулятора (пропорционально-дифференциальный):

${\gamma }_p=k_{p\text{ 1}}\Delta n_y+k_{\text{p 2}}\Delta n_z+k_{\text{p 3}}\Delta {\omega }_y+k_{\text{p 4}}\Delta {\omega }_z+k_{\text{p 5}}\Delta {\dot{n}}_y+k_{\text{p 6}}\Delta {\dot{n}}_z+k_{\text{p 7}}\Delta {\dot{\omega }}_y+k_{\text{p 8}}\Delta {\dot{\omega }}_z$.

Здесь $\Delta n_y=n_y-n_{ys}$, $\Delta n_z=n_z-n_{zs}$, $\Delta {\omega }_y={\omega }_y-{\omega }_{ys}$, $\Delta {\omega }_z={\omega }_z-{\omega }_{zs}$, $\Delta {\dot{n}}_y=u_1-{\dot{n}}_{ys}$, $\Delta {\dot{n}}_z=u_2-{\dot{n}}_{zs}$, $\Delta {\dot{\omega }}_y=u_3-{\dot{\omega }}_{ys}$, $\Delta {\dot{\omega }}_z=u_4-{\dot{\omega }}_{zs}$, где ${\dot{n}}_{ys}=\left[n_{ys}\left(t+\Delta t\right)-n_{ys}\left(t\right)\right]/\Delta t$, ${\dot{n}}_{zs}=\left[n_{zs}\left(t+\Delta t\right)-n_{zs}\left(t\right)\right]/\Delta t$, ${\dot{\omega }}_{ys}=\left[{\omega }_{ys}\left(t+\Delta t\right)-{\omega }_{ys}\left(t\right)\right]/\Delta t$, ${\dot{\omega }}_{zs}=\left[{\omega }_{zs}\left(t+\Delta t\right)-{\omega }_{zs}\left(t\right)\right]/\Delta t$. Проекции перегрузки в алгоритме со спиральным прогнозом рассматриваются в связанных осях, а в имитационной модели – в траекторной СК. Поэтому n (т. е. $n_T$) надо перевести в $n_{св}$ с помощью матрицы перехода $n_{св}=D_{св}^{Т}{n}_{Т}$, $D_{св}^{Т}=D_{св}^{Н}{D}_{Н}^{Т}$, $D_{св}^{Н}={\left(D_{Н}^{св}\right)}^{Т}$,

$n_{\text{св \hspace{0.17em}x}}=d_{11}{n}_x+d_{12}{n}_y+d_{13}{n}_z$, $n_{\text{св \hspace{0.17em}}y}=d_{21}{n}_x+d_{22}{n}_y+d_{23}{n}_z$, $n_{\text{св \hspace{0.17em}}z}=d_{31}{n}_x+d_{32}{n}_y+d_{33}{n}_z$;

$d_{11}=sin\vartheta sin\theta +cos\left(\psi -\Psi \right)cos\vartheta cos\theta$, $d_{12}=sin\vartheta cos\theta -cos\left(\psi -\Psi \right)cos\vartheta sin\theta$,

$d_{13}=-cos\vartheta sin\left(\psi -\Psi \right)$,

$d_{21}=cos\gamma cos\vartheta sin\theta -cos\left(\psi -\Psi \right)cos\gamma sin\vartheta cos\theta +sin\left(\psi -\Psi \right)sin\gamma cos\theta$,

$d_{22}=cos\gamma cos\vartheta cos\theta +cos\left(\psi -\Psi \right)cos\gamma sin\vartheta sin\theta -sin\left(\psi -\Psi \right)sin\gamma sin\theta$,

$d_{23}=cos\left(\psi -\Psi \right)sin\gamma -sin\left(\psi -\Psi \right)cos\gamma sin\upsilon$,

$d_{31}=-sin\gamma cos\vartheta sin\theta +sin\left(\psi -\Psi \right)cos\gamma cos\theta +cos\left(\psi -\Psi \right)sin\gamma sin\vartheta cos\theta$,

$d_{32}=-sin\gamma cos\vartheta cos\theta -sin\left(\psi -\Psi \right)cos\gamma sin\theta -cos\left(\psi -\Psi \right)sin\gamma sin\vartheta sin\theta$,

$d_{33}=cos\left(\psi -\Psi \right)cos\gamma -sin\left(\psi -\Psi \right)sin\gamma sin\vartheta$

($d_{ij}$, $i,j=\overline{\mathrm{1,3}}$, – элементы матрицы $D_{св}^{Т}$). Эти проекции перегрузки отмечены в алгоритме управления подстрочным индексом s.

Проекции угловой скорости в алгоритме со спиральным прогнозом рассматриваются в связанных осях, а в имитационной модели – в полусвязанной СК. Поэтому ${\omega }_{пс}$ надо перевести в ${\omega }_{св}$ (т. е. $\omega$) с помощью матрицы перехода $\omega =D_{св}^{пс}{\omega }_{пс}$, ${\omega }_{\text{x}}={\omega }_{пс x}$, ${\omega }_y=cos\gamma {\omega }_{псy}-sin\gamma {\omega }_{пс z}$, ${\omega }_z=sin\gamma {\omega }_{пс y}+cos\gamma {\omega }_{пс z}$. Здесь $o_{11}=1$, $o_{12}=0$, $o_{13}=0$, $o_{21}=0$, $o_{22}=cos\gamma$, $o_{23}=-sin\gamma$, $o_{31}=0$, $o_{23}=sin\gamma$, $o_{33}=cos\gamma$($o_{ij}$, $i,j=\overline{\mathrm{1,3}}$, – элементы матрицы $D_{св}^{пс}$). Эти проекции угловой скорости отмечены в алгоритме управления подстрочным индексом s.

Успешность решения задачи сильно зависит от малых маневренных возможностей стабилизированного вращением ЛА. Поэтому на носителе перед пуском производится оценка области достижимости путем моделирования динамики имитационной модели при различных углах поворота блока управления (${\gamma }_p\in \left[\mathrm{0,}\text{\hspace{0.17em}}2\pi \right]$). При нахождении ЛА в области управляемости можно осуществлять пуск, и алгоритм управления со спиральным прогнозом сформирует оптимальные значения перегрузки и угловой скорости, реализация которых посредством ПИД-регулятора (СОРЭ) обеспечит доставку реагента в досягаемую точку с допустимой точностью.

Решение задачи управления полетом ЛА выполнялось с использованием предложенного алгоритма при следующих начальных условиях: $x={\left(270-0.350125000000040000\right)}^T$, конечных условиях $x_f=1900$ м, $y_f=3000$м, $z_f=50$ м, при $d=0.1м$, $S=0.00785{м}^2$, $S_{р}=0.0016{м}^2$, $m=15.5кг$, $I_x=0.0208кг {\text{\hspace{0.17em}м}}^2$, $I_y=I_z=0.2219кг\cdot {\text{\hspace{0.17em}м}}^2$, $L=0.49\text{\hspace{0.17em}\hspace{0.17em}м}$, $c_x=0.1486$, $c_{x\text{\hspace{0.17em}р}}=0.0219$, $c_n^{\delta }=1.8683$, $c_y^{\delta }=1.8146$, $m_z^{\delta }=0.2647$, $m_z^{{\omega }_z}=0.0386$, $m_x=0.0001$, $m_x^{{\omega }_x}=-0.0007$, $l_{р}=0.25\text{\hspace{0.17em}\hspace{0.17em}м}$, ${\delta }_{р}=0.087$. Численное интегрирование уравнений имитационной модели выполнялось методом Эйлера с шагом 0.0001 с при обращении к алгоритму спирального прогноза для получения оптимальных значений перегрузок и угловых скоростей через каждые 100 шагов. Время счета составило 3 с.

Вначале алгоритм спирального прогноза формировал свою траекторию независимо от динамики модели с выдачей значений проекций перегрузки и угловой скорости для ПД-регулятора, с выхода которого поступали значения, усредненные на интервалах в 20 шагов ${\gamma }_{р}\left(t\right)$, в имитационную модель. Принимались следующие значения параметров критерия оптимальности: $t_f-t_0=7$с, ${\rho }_1=0.02$, ${\rho }_2=0.04$, ${\rho }_3=0.04$, $k_1^2={10}^{-3}$, $k_2^2={10}^{-2}$, $k_3^2={10}^{-4}$, $k_4^2={10}^{-8}$, $k_5^2={10}^{-3}$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}1}=0.02$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}2}=0.2$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}3}=0.2$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}4}=0.02$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}5}=0.02$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}6}=0.2$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}7}=0.2$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}8}=0.02$. На рис. 2 представлен график , характеризующий изменение величины интервала оптимизации для траектории, формируемой с помощью алгоритма со спиральным прогнозом. При этом терминальные значения требуемых координат получились следующими: $x\left(t_f\right)=1959.4$м, $y\left(t_f\right)=3052.7$м, $z\left(t_f\right)=49.1$м. На рис. 3 помещена зависимость ${\gamma }_{р}\left(t\right)$.

 

Рис. 2. Зависимость $t_f\left(t\right)$

 

Рис. 3. Зависимость ${\gamma }_{р}\left(t\right)$

 

Этот результат уже можно использовать для доставки распылителя в окрестность заданной точки. Однако с увеличением длины интервала оптимизации траектория модели все больше отличается от траектории оптимального эталона.

Далее, алгоритм со спиральным прогнозом применялся из текущего положения ЛА, выдавая проекции перегрузки и угловой скорости в ПД-регулятор для введения ${\gamma }_{р}\left(t\right)$ в имитационную модель. Параметры равнялись: ${\rho }_1=1$, ${\rho }_2=1$, ${\rho }_3=40$, $k_1^2=5\cdot {10}^{-4}$, $k_2^2=5\cdot {10}^{-3}$, $k_3^2={10}^{-4}$, $k_4^2={10}^{-8}$, $k_5^2={10}^{-5}$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}1}=0.02$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}2}=0.1$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}3}=0.2$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}4}=0.02$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}5}=0.02$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}6}=0.1$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}7}=0.2$, $k_{\text{p\hspace{0.17em}\hspace{0.17em}}8}=0.02$. Имитационная модель посредством ПД-регулятора приводит ЛА в точку с координатами $x\left(t_f\right)=1919.7$м, $y\left(t_f\right)=3035.8$м, $z\left(t_f\right)=48.8$м. На рис. 4, 5 помещены графики $t_f\left(t\right)$ и ${\gamma }_{р}\left(t\right)$соответственно.

 

Рис. 4. Зависимость $t_f\left(t\right)$

 

Рис. 5. Зависимость ${\gamma }_{р}\left(t\right)$

 

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

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

Об авторах

Е. А. Знаменский

БГТУ “ВОЕНМЕХ” им. Д.Ф. Устинова

Автор, ответственный за переписку.
Email: eugenznam@mail.ru
Россия, Санкт-Петербург

Д. С. Кабанов

БГТУ “ВОЕНМЕХ” им. Д.Ф. Устинова

Email: kabanovds@mail.ru
Россия, Санкт-Петербург

С. А. Кабанов

БГТУ “ВОЕНМЕХ” им. Д.Ф. Устинова

Email: kaba-sa@mail.ru
Россия, Санкт-Петербург

Е. Н. Никулин

БГТУ “ВОЕНМЕХ” им. Д.Ф. Устинова

Email: yevgeniy.nikulin.47@mail.ru
Россия, Санкт-Петербург

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

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