On one problem of quadcopter control with given intermediate values of different parts of coordinates

Full Text

Введение

Все большую популярность получают беспилотные летательные аппараты в виде многовинтовых устройств, а чаще — квадрокоптеров. Квадрокоптеры обладают рядом преимуществ, таких как надежность и простота конструкции, большая стабильность, компактность и маневренность, малая взлетная масса. Область применения квадрокоптеров достаточно широка и не ограничена военной отраслью (см. [1–3]). Например, квадрокоптеры могут быть использованы как недорогое и эффективное средство для получения фото- и видеоизображений с воздуха. Квадрокоптер хорошо подходит для наблюдения и контроля объектов, территорий и зон, доступ к которым затруднен, или в условиях, непригодных для человека (см. [2–4]). Исследованиям различных задач, связанных с полетом (эксплуатацией) квадрокоптеров, посвящены, в частности, работы [5–7]. В [8] предложен способ управления квадрокоптером, оснащенным манипулятором. В [9–11] изложены основы моделирования и алгоритмы управления квадрокоптером.

В прикладных задачах при моделировании и проектировании движений различных механических систем, в частности, движения квадрокоптеров, возникают задачи управления движением с многоточечными промежуточными условиями (см. [12–14]). В подобных задачах управления наряду с классическими краевыми (начальное и конечное) условиями необходимо также учитывать многоточечные промежуточные условия (см. [15–17]). Эти задачи управления имеют важное прикладное и теоретическое значение.

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

1. Математическая модель динамики квадракоптера и постановка задачи

В качестве динамического объекта управления рассматривается модель геометрического симметричного беспилотного летательного аппарата (БПЛА) — квадрокоптера, построенного по классической четырехвинтовой схеме, структура которого представлена на рис. 1.

 

Рис. 1. Структурная схема БПЛА

 

Положение квадрокоптера с массой m в пространстве характеризуется координатами $x\mathrm{,}y\mathrm{,}z$ центра масс аппарата в неподвижной декартовой системе координат и тремя углами (углы Эйлера) $\phi \mathrm{,}\vartheta \mathrm{,}\psi$ поворота вокруг осей системы координат, жестко связанной с аппаратом, причем начало координат совпадает с центром масс аппарата. Здесь $\phi$ — угол тангажа (т. е. угол наклона квадрокоптера относительно оси $Ox$),  — угол крена (т. е. угол наклона БПЛА относительно оси $Oy$), $\psi$ — угол рыскания (т. е. угол наклона квадрокоптера относительно оси $Oz$). Движение квадрокоптера в задачах динамики управления полетом рассматривается как сложное — движение центра масс и движение вокруг центра масс [9–11]. Для составления математической модели квадрокоптера приняты следующие допущения [1–3]:

• квадрокоптер симметричен относительно осей x и y;
• рама квадрокоптера и его винты абсолютно жесткие;
• каждый двигатель располагается на конце стержня;
• тяга, создаваемая каждым винтом, перпендикулярна плоскости xy.

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

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

Угловые скорости винтов обозначим через ${\omega }_1\mathrm{,}{\omega }_2\mathrm{,}{\omega }_3\mathrm{,}{\omega }_4\mathrm{,}$ силы тяги, создаваемые вращением винтов — через $f_1\mathrm{,}{f}_2\mathrm{,}{f}_3\mathrm{,}{f}_4\mathrm{,}$ а моменты, которые образуются в результате вращения винтов — через $M_1\mathrm{,}{M}_2\mathrm{,}{M}_3\mathrm{,}{M}_4\mathrm{.}$ Для получения математической модели динамики квадрокоптера введем неподвижную систему координат $Oxyz\mathrm{,}$ а в центре масс $O_B$ квадрокоптера закрепим неподвижную систему координат $O_B{x}_B{y}_B{z}_B\mathrm{.}$ Обозначим через $x\mathrm{,}y\mathrm{,}z$ координаты центра масс квадрокоптера относительно системы координат $Oxyz\mathrm{.}$ Опишем вращательное движение квадрокоптера в инерциальной системе координат $Oxyz\mathrm{,}$ используя углы Эйлера $\phi \mathrm{,}\vartheta \mathrm{,}\psi \mathrm{.}$ Тогда пространственное положение квадрокоптера описывается двумя векторами

$\xi ={(x,y,z)}^T,\eta ={(\phi ,\vartheta ,\psi )}^T.$

Здесь и далее буква «T» в верхнем индексе означает операцию транспонирования. В системе координат $O_B{x}_B{y}_B{z}_B$ запишем векторы линейной скорости $V_B$ и угловой скорости v.  

$V_B={(V_{Bx},V_{By},V_{Bz})}^T,v={(p,q,r)}^T.$

А матрица поворота системы координат $O_B{x}_B{y}_B{z}_B$ от.осительно неподвижной инерциальной системы координат $Oxyz$ будет:

$R=\left(\begin{array}{ccc}{C}_{\psi }{C}_{\vartheta }& C_{\psi }{S}_{\vartheta }{S}_{\phi }-S_{\psi }{C}_{\phi }& C_{\psi }{S}_{\vartheta }{C}_{\phi }+S_{\psi }{S}_{\phi }\\ S_{\psi }{C}_{\vartheta }& S_{\psi }{S}_{\vartheta }{S}_{\phi }+C_{\psi }{C}_{\phi }& S_{\psi }{S}_{\vartheta }{C}_{\phi }-C_{\psi }{S}_{\phi }\\ -S_{\vartheta }& C_{\vartheta }{S}_{\phi }& C_{\vartheta }{C}_{\phi }\end{array}\right),$

где $C_{\alpha }=\cos \alpha ,$ $S_{\alpha }=\sin \alpha ,$ $\alpha =\phi ,\vartheta ,\psi .$   

Легко проверить, что матрица R ортогональна, следовательно, $R^{-1}=R^T,$ и является матрицей перехода из инерциальной системы координат $Oxyz$ в систему $O_B{x}_B{y}_B{z}_B\mathrm{.}$

Матрицы R и $R^{-1}$ используются для получения проекций вектора $V_B$ в инерциальной системе координат.

Матрица W преобразования из инерциальной системы координат $Oxyz$ в систему $O_B{x}_B{y}_B{z}_B$ и матрица $W^{-1}$ обратного преобразования для вектора угловой скорости v в данном случае имеют вид

$W=\left(\begin{array}{ccc}1& 0& -S_{\vartheta }\\ 0& C_{\phi }& C_{\vartheta }{S}_{\phi }\\ 0& -S_{\phi }& C_{\vartheta }{C}_{\phi }\end{array}\right)$

и

$W^{-1}=\left(\begin{array}{ccc}1& S_{\phi }{T}_{\vartheta }& C_{\phi }{T}_{\vartheta }\\ 0& C_{\phi }& -S_{\phi }\\ 0& \frac{S_{\phi }}{C_{\vartheta }}& \frac{C_{\phi }}{C_{\vartheta }}\end{array}\right).$

Здесь $T_{\alpha }=tg\alpha ,$ $\alpha =\phi ,\vartheta ,\psi .$ Заметим, что матрица W обратима тогда и только тогда, когда $\vartheta \ne \left(2k-1\right)\frac{\pi }{2},$ $k\in \mathbb{Z}\mathrm{.}$ Теперь, используя матрицу преобразования W и ее обратную $W^{-1}\mathrm{,}$ получаем

$\dot{\eta }=W^{-1}v,или\left(\begin{array}{c}\dot{\phi }\\ \dot{\vartheta }\\ \dot{\psi }\end{array}\right)=\left(\begin{array}{ccc}1& S_{\phi }{T}_{\vartheta }& C_{\phi }{T}_{\vartheta }\\ 0& C_{\phi }& -S_{\phi }\\ 0& \frac{S_{\phi }}{C_{\vartheta }}& \frac{C_{\phi }}{C_{\vartheta }}\end{array}\right)\left(\begin{array}{c}p\\ q\\ r\end{array}\right)$ (1.1)

и

$v=W_{\eta }\dot{\eta },или\left(\begin{array}{c}p\\ q\\ r\end{array}\right)=\left(\begin{array}{ccc}1& 0& -S_{\vartheta }\\ 0& C_{\phi }& C_{\vartheta }{S}_{\phi }\\ 0& -S_{\phi }& C_{\vartheta }{C}_{\phi }\end{array}\right)\left(\begin{array}{c}\dot{\phi }\\ \dot{\vartheta }\\ \dot{\psi }\end{array}\right).$

Благодаря указанному на рис. 1 выбору системы координат $O_B{x}_B{y}_B{z}_B\mathrm{,}$ квадрокоптер расположен симметрично относительно осей $O_B{x}_B$ и $O_B{y}_B\mathrm{.}$ Следовательно, инерционная матрица I будет диагональной

$I=\left(\begin{array}{ccc}{I}_{xx}& 0& 0\\ 0& I_{yy}& 0\\ 0& 0& I_{zz}\end{array}\right),$

где $I_{xx}\mathrm{,}$ $I_{yy}$ и $I_{zz}$ — моменты инерции квадрокоптера относительно осей $Ox\mathrm{,}$ $Oy$ и $Oz\mathrm{,}$ соответственно, и $I_{xx}=I_{yy}.$

Угловая скорость ${\omega }_i$ винта i $\left(i=\overline{\mathrm{1,4}}\right)$ создает силу тяги $f_i$ $\left(i=\overline{\mathrm{1,4}}\right),$ направленную по оси винта. Угловая скорость ${\omega }_i\mathrm{,}$ $\left(i=\overline{\mathrm{1,4}}\right)$ и угловое ускорение ${\dot{\omega }}_i\mathrm{,}$ $\left(i=\overline{\mathrm{1,4}}\right)$ винта i $\left(i=\overline{\mathrm{1,4}}\right)$ также создают крутящий момент $M_i$ $\left(i=\overline{\mathrm{1,4}}\right)$ вокруг оси винта. Сила тяги и крутящий момент определяются по формулам

$f_i=k{\omega }_i^2,M_i=b{\omega }_i^2+I_M{\dot{\omega }}_i,i=\overline{\mathrm{1,4}},$

где k — постоянная подъемной силы, b — постоянная сопротивления (т. е. k и b — коэффициенты пропорциональности, характеризующие особенности роторов), $I_M$ — момент инерции винта (одинаковый для всех винтов). Влияние углового ускорения ${\dot{\omega }}_i$ $\left(i=\overline{\mathrm{1,4}}\right)$ настолько мало, что им можно пренебречь. Суммарная сила тяги f сил $f_1\mathrm{,}{f}_2\mathrm{,}{f}_3\mathrm{,}{f}_4$ направлена вдоль оси $O_B{z}_B\mathrm{.}$ Обозначим ${\tau }_B$ — вектор крутящего момента, который включает крутящие моменты ${\tau }_{\phi }\mathrm{,}{\tau }_{\vartheta }\mathrm{,}{\tau }_{\psi }\mathrm{,}$ соответствующие углам Эйлера $\phi \mathrm{,}\vartheta \mathrm{,}\psi \mathrm{.}$ Сила тяги и крутящий момент определяются по формулам

$f=\sum _{i=1}^4{f}_i=k\sum _{i=1}^4{\omega }_i^2,F_B=\left(\begin{array}{l}0\\ 0\\ f\end{array}\right),{\tau }_B=\left(\begin{array}{l}{\tau }_{\phi }\\ {\tau }_{\vartheta }\\ {\tau }_{\psi }\end{array}\right)=\left(\begin{array}{c}lk\left(-{\omega }_2^2+{\omega }_4^2\right)\\ lk\left(-{\omega }_1^2+{\omega }_2^2\right)\\ \sum _{i=1}^4{M}_i\end{array}\right),$ (1.2)

где I — расстояние от осей винта (т. е. центра роторов) до центра масс квадрокоптера.

Физический смысл формул (1.2) очевиден. Можно изменить угловые скорости четвертого и второго винтов, тем самым изменить угол крена. Таким же образом можно изменить угол наклона (тангажа), изменяя скорости вращения первого и третьего винтов. Для изменения угла рыскания можно уменьшить скорости вращения двух диагональных винтов (второго и четвертого) и увеличить скорости вращения двух других (первого и третьего) или наоборот, увеличить скорости вращения двух диагональных винтов (второго и четвертого) и уменьшить скорости вращения двух других (первого и третьего).

Для получения дифференциальных уравнений движения квадрокоптера воспользуемся вторым законом Ньютона и динамическими уравнениями Эйлера. Таким образом, относительно системы координат $O_B{x}_B{y}_B{z}_B$ получим (подробнее см. [9–11]) уравнения

$m{\dot{V}}_B+v\times m{V}_B=R^{T}G+F_B,$

где $G={(\mathrm{0,0,}-g)}^T,$ g — ускорение свободного падения. А относительно системы координат $Oxyz$ будем иметь

$m\ddot{\xi }=G+R{F}_B,$

или

$\left(\begin{array}{c}\ddot{x}\\ \ddot{y}\\ \ddot{z}\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ -g\end{array}\right)+\frac{f}{m}\left(\begin{array}{c}{C}_{\psi }{S}_{\vartheta }{C}_{\phi }+S_{\psi }{S}_{\phi }\\ S_{\psi }{S}_{\vartheta }{C}_{\phi }-C_{\psi }{S}_{\phi }\\ C_{\vartheta }{C}_{\phi }\end{array}\right).$ (1.3)

В системе координат $O_B{x}_B{y}_B{z}_B$ сумма углового ускорения инерции $I\dot{v}\mathrm{,}$ центростремительных сил $v\times Iv$ и гироскопических сил $\Gamma$ равна внешнему крутящему моменту

$I\dot{v}+v\times \left(Iv\right)+\Gamma ={\tau }_B,$ (1.4)

где

$\Gamma =I_r\left(\begin{array}{l}p\\ q\\ r\end{array}\right)\times \left(\begin{array}{l}0\\ 0\\ 1\end{array}\right){\omega }_{\Gamma },{\omega }_{\Gamma }={\omega }_1-{\omega }_2+{\omega }_3-{\omega }_4.$

Здесь $I_r$ — момент инерции квадрокоптера относительно его мгновенной оси вращения. Таким образом, из (1.4) получаем

$\dot{v}=I^{-1}(-\left(\begin{array}{l}p\\ q\\ r\end{array}\right)\times \left(\begin{array}{l}{I}_{xx}p\\ I_{yy}q\\ I_{zz}r\end{array}\right)-I_r\left(\begin{array}{l}p\\ q\\ r\end{array}\right)\times \left(\begin{array}{l}0\\ 0\\ 1\end{array}\right){\omega }_{\Gamma }+{\tau }_B),$ (1.5)

Имея угловые ускорения в системе координат $O_B{x}_B{y}_B{z}_B\mathrm{,}$ можно получить угловые ускорения в системе координат $Oxyz\mathrm{,}$ используя матрицу преобразования $W^{-1}$ и ее производную по времени.

$\ddot{\eta }=\frac{d}{dt}\left(W^{-1}v\right)=\frac{d}{dt}\left(W^{-1}\right)v+W^{-1}\dot{v}$

$=\left(\begin{array}{ccc}0& \dot{\phi }{C}_{\phi }{T}_{\vartheta }+\frac{\dot{\vartheta }{S}_{\phi }}{C_{\vartheta }^2}& -\dot{\phi }{S}_{\phi }{C}_{\vartheta }+\frac{\dot{\vartheta }{C}_{\phi }}{C_{\vartheta }^2}\\ 0& -\dot{\phi }{S}_{\phi }& -\dot{\phi }{C}_{\phi }\\ 0& \frac{\dot{\phi }{C}_{\phi }}{C_{\vartheta }}+\frac{\dot{\phi }{S}_{\phi }{T}_{\vartheta }}{C_{\vartheta }}& -\frac{\dot{\phi }{S}_{\phi }}{C_{\vartheta }}+\frac{\dot{\vartheta }{C}_{\phi }{T}_{\vartheta }}{C_{\vartheta }}\end{array}\right)v+W^{-1}\dot{v}.$

Из (1.1), (1.3), (1.5) получим следующую систему уравнений относительно компонент векторов:

$\begin{array}{r}\ddot{x}=\frac{f}{m}{C}_{\psi }{S}_{\vartheta }{C}_{\phi }+\frac{f}{m}{S}_{\psi }{S}_{\phi },\ddot{y}=\frac{f}{m}{S}_{\psi }{S}_{\vartheta }{C}_{\phi }-\frac{f}{m}{C}_{\psi }{S}_{\phi },\ddot{z}=-g+\frac{f}{m}{C}_{\vartheta }{C}_{\phi },\\ \dot{\phi }=p+\frac{S_{\phi }{S}_{\vartheta }}{C_{\vartheta }}q+\frac{C_{\phi }{S}_{\vartheta }}{C_{\vartheta }}r,\dot{\vartheta }=C_{\phi }q-S_{\phi }r,\dot{\psi }=\frac{S_{\phi }}{C_{\vartheta }}q+\frac{C_{\phi }}{C_{\vartheta }}r,\\ \dot{p}=\frac{\left(I_{yy}-I_{zz}\right)qr}{I_{xx}}-I_r\frac{q}{I_{xx}}{\omega }_{\Gamma }+\frac{{\tau }_{\phi }}{I_{xx}},\dot{q}=\frac{\left(I_{zz}-I_{xx}\right)pr}{I_{yy}}-I_r\frac{p}{I_{yy}}{\omega }_{\Gamma }+\frac{{\tau }_{\vartheta }}{I_{yy}},\\ \dot{r}=\frac{\left(I_{xx}-I_{yy}\right)pq}{I_{zz}}-I_r\frac{q}{I_{zz}}{\omega }_{\Gamma }+\frac{{\tau }_{\psi }}{I_{zz}}.\end{array}$ (1.6)

Линеаризуем математическую модель квадрокоптера (1.6). Линеаризацию произведем в окрестности начала координат, считая, что углы $\phi \mathrm{,}\vartheta \mathrm{,}\psi$ — малы, соответствующие синусы равны нулю и косинусы равны единице. Далее введя обозначения

$\begin{array}{l}{x}_1=x,x_2=\dot{x},x_3=y,x_4=\dot{y},x_5=z,x_6=\dot{z},\\ x_7=\phi ,x_8=\vartheta ,x_9=\psi ,x_{10}=p,x_{11}=q,x_{12}=r,\end{array}$

окончательно получим систему линеаризованных уравнений, описывающих динамику линейной модели квадрокоптера:

$\begin{array}{l}{\dot{x}}_1=x_2,{\dot{x}}_2=g{x}_8,{\dot{x}}_3=x_4,{\dot{x}}_4=-g{x}_7,{\dot{x}}_5=x_6,{\dot{x}}_6=u_1,\\ {\dot{x}}_7=x_{10},{\dot{x}}_8=x_{11},{\dot{x}}_9=x_{12},{\dot{x}}_{10}=\frac{u_2}{I_{xx}},{\dot{x}}_{11}=\frac{u_3}{I_{yy}},{\dot{x}}_{12}=\frac{u_4}{I_{zz}},\end{array}$ (1.7)

где $u_1=\frac{f}{m}-g,$ $u_2={\tau }_{\phi },$ $u_3={\tau }_{\vartheta },$ $u_4={\tau }_{\psi }$  — управляющие воздействия. Далее будем обозначать

Непосредственной проверкой можно убедиться, что система (1.7) является вполне управляемой (см. [18, гл. 6]).

Пусть заданы начальное (при $t=t_0$) и конечное (при $t=T$) состояния системы (1.7)

$x\left(t_0\right)={\left(x_1\left(t_0\right),\dots ,x_{12}\left(t_0\right)\right)}^T,x\left(T\right)={\left(x_1\left(T\right),\dots ,x_{12}\left(T\right)\right)}^T,$ (1.8)

а также значения

$x_5\left(t_1\right),x_1\left(t_2\right),x_3\left(t_2\right)$ (1.9)

части координат фазового вектора в заданные промежуточные моменты времени $t=t_1$ и $t=t_2$ ($t_0<t_1<t_2<T$).

Рассмотрим следующую задачу: требуется найти условия, при которых существует программное управляющее воздействие $U\left(t\right),$ $t\in \left[t_0,T\right]$ и программное движение $x\left(t\right),$ удовлетворяющие системе (1.7) и условиям (1.8) и (1.9), а также построить их.

2. Решение задачи

Решение уравнения (1.7) с помощью формулы Коши запишем следующим образом (см. [18, §5], [19, с. 13]):

$x\left[t\right]=X\left[t,t_0\right]x\left(t_0\right)+{\int }_{t_0}^{t}H\left[t,\tau \right]U\left(\tau \right)d\tau ,$ (2.1)

где $X\left[t,\tau \right]$ — нормированная фундаментальная матрица решения однородной части уравнения (1.7), a $H\left[t,\tau \right]$ — импульсно-переходная матрица, которая имеет вид:

$H[t,\tau ]={\left(\begin{array}{cccccccccccc}0& 0& 0& 0& t-\tau & 1& 0& 0& 0& 0& 0& 0\\ 0& 0& -\frac{g{\left(t-\tau \right)}^3}{6I_{xx}}& -\frac{g{\left(t-\tau \right)}^2}{2I_{xx}}& 0& 0& \frac{t-\tau }{I_{xx}}& 0& 0& \frac{1}{I_{xx}}& 0& 0\\ \frac{g{\left(t-\tau \right)}^3}{6I_{yy}}& \frac{g{\left(t-\tau \right)}^2}{2I_{yy}}& 0& 0& 0& 0& 0& \frac{t-\tau }{I_{yy}}& 0& 0& \frac{1}{I_{yy}}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& \frac{t-\tau }{I_{zz}}& 0& 0& \frac{1}{I_{zz}}\end{array}\right)}^{\text{}T}.$

Применяя подходы, приведенные в работах [15], [16, гл. 3], [17], для определения законов управления с учетом условий (1.8) и (1.9) получим следующие интегральные соотношения

${\int }_{t_0}^T\frac{g{\left(T-t\right)}^3}{6I_{yy}}{u}_3\left(t\right)dt=C_1,{\int }_{t_0}^T\frac{g{\left(T-t\right)}^2}{2I_{yy}}{u}_3\left(t\right)dt=C_2,{\int }_{t_0}^T-\frac{g{\left(T-t\right)}^3}{6I_{xx}}{u}_2\left(t\right)dt=C_3,$

${\int }_{t_0}^T-\frac{g{\left(T-t\right)}^2}{2I_{xx}}{u}_2\left(t\right)dt=C_4,{\int }_{t_0}^T\left(T-t\right)u_1\left(t\right)dt=C_5,{\int }_{t_0}^T{u}_1\left(t\right)dt=C_6,$

${\int }_{t_0}^T\frac{T-t}{I_{xx}}{u}_2\left(t\right)dt=C_7,{\int }_{t_0}^T\frac{T-t}{I_{yy}}{u}_3\left(t\right)dt=C_8,{\int }_{t_0}^T\frac{T-t}{I_{zz}}{u}_4\left(t\right)dt=C_9,$

${\int }_{t_0}^T\frac{1}{I_{xx}}{u}_2\left(t\right)dt=C_{10},{\int }_{t_0}^T\frac{1}{I_{yy}}{u}_3\left(t\right)dt=C_{11},{\int }_{t_0}^T\frac{1}{I_{zz}}{u}_4\left(t\right)dt=C_{12},$

${\int }_{t_0}^{t_2}\frac{g{\left(t_2-t\right)}^3}{6I_{yy}}{u}_3\left(t\right)dt=C_{13},{\int }_{t_0}^{t_2}-\frac{g{\left(t_2-t\right)}^3}{6I_{xx}}{u}_2\left(t\right)dt=C_{14},{\int }_{t_0}^{t_1}\left(t_1-t\right)u_1\left(t\right)dt=C_{15},$

где

$C_1=x_1\left(T\right)-x_1\left(t_0\right)-x_2\left(t_0\right)\left(T-t_0\right)-\frac{g}{2}{\left(T-t_0\right)}^2{x}_8\left(t_0\right)-\frac{g}{6}{\left(T-t_0\right)}^3{x}_{11}\left(t_0\right),$

$C_2=x_2\left(T\right)-x_2\left(t_0\right)-x_8\left(t_0\right)g\left(T-t_0\right)-\frac{g}{2}{\left(T-t_0\right)}^2{x}_{11}\left(t_0\right),$

$C_3=x_3\left(T\right)-x_3\left(t_0\right)-x_4\left(t_0\right)\left(T-t_0\right)+\frac{g}{2}{\left(T-t_0\right)}^2{x}_7\left(t_0\right)+\frac{g}{6}{\left(T-t_0\right)}^3{x}_{10}\left(t_0\right),$

$C_4=x_4\left(T\right)-x_4\left(t_0\right)+x_7\left(t_0\right)g\left(T-t_0\right)+\frac{g}{2}{\left(T-t_0\right)}^2{x}_{10}\left(t_0\right),$

$C_5=x_5\left(T\right)-x_5\left(t_0\right)-x_6\left(t_0\right)\left(T-t_0\right),$

$C_j=x_j\left(T\right)-x_j\left(t_0\right),j=\mathrm{6,10,11,12,}$

$C_7=x_7\left(T\right)-x_7\left(t_0\right)-x_{10}\left(t_0\right)\left(T-t_0\right),C_8=x_8\left(T\right)-x_8\left(t_0\right)-x_{11}\left(t_0\right)\left(T-t_0\right),$

$C_9=x_9\left(T\right)-x_9\left(t_0\right)-x_{12}\left(t_0\right)\left(T-t_0\right),$

$C_{13}=x_1\left(t_2\right)-x_1\left(t_0\right)-x_2\left(t_0\right)\left(t_2-t_0\right)-\frac{g}{2}{\left(t_2-t_0\right)}^2{x}_8\left(t_0\right)-\frac{g}{6}{\left(t_2-t_0\right)}^3{x}_{11}\left(t_0\right),$

$C_{14}=x_3\left(t_2\right)-x_3\left(t_0\right)-x_4\left(t_0\right)\left(t_2-t_0\right)+\frac{g}{2}{\left(t_2-t_0\right)}^2{x}_7\left(t_0\right)+\frac{g}{6}{\left(t_2-t_0\right)}^3{x}_{10}\left(t_0\right),$

$C_{15}=x_5\left(t_1\right)-x_5\left(t_0\right)-\left(t_1-t_0\right)x_6\left(t_0\right).$

В векторно-матричной форме полученные интегральные соотношения будут иметь вид

${\int }_{t_0}^{T}H\left(t\right)U\left(t\right)dt=\sigma \left(t_0,t_1,t_2,T\right),$ (2.2)

где $\sigma \left(t_0,t_1,t_2,T\right)={\left(C_1,\dots ,C_{15}\right)}^T;$ элементы $15\times 4$ матрицы $H\left(t\right)=H[T,t]$ определяются формулами

$h_1^{\left(3\right)}[T,t]=\frac{g{\left(T-t\right)}^3}{6I_{yy}},h_2^{\left(3\right)}[T,t]=\frac{g{\left(T-t\right)}^2}{2I_{yy}},h_3^{\left(2\right)}[T,t]=-\frac{g{\left(T-t\right)}^3}{6I_{xx}},$

$h_4^{\left(2\right)}[T,t]=-\frac{g{\left(T-t\right)}^2}{2I_{xx}},h_5^{\left(1\right)}[T,t]=\left(T-t\right),h_6^{\left(1\right)}[T,t]=\mathrm{1,}{h}_7^{\left(2\right)}[T,t]=\frac{T-t}{I_{xx}},$

$h_8^{\left(3\right)}[T,t]=\frac{T-t}{I_{yy}},h_9^{\left(4\right)}[T,t]=\frac{T-t}{I_{zz}},h_{10}^{\left(2\right)}[T,t]=\frac{1}{I_{xx}},h_{11}^{\left(3\right)}[T,t]=\frac{1}{I_{yy}},$

$h_{12}^{\left(4\right)}[T,t]=\frac{1}{I_{zz}},h_{13}^{\left(1\right)}[t_1,t]=\left\{\begin{array}{ll}{t}_1-t& приt_0\le t\le t_1\\ 0& приt_1<t\le T\end{array}\right.,$

$h_{14}^{\left(3\right)}[t_2,t]=\left\{\begin{array}{ll}-\frac{g{\left(t_2-t\right)}^3}{6I_{yy}}& приt_0\le t\le t_2\\ 0& приt_2<t\le T\end{array}\right.\text{},h_{15}^{\left(2\right)}[t_2,t]=\left\{\begin{array}{ll}-\frac{g{\left(t_2-t\right)}^3}{6I_{xx}}& приt_0\le t\le t_2\\ 0& приt_2<t\le T\end{array}\right.$

(нижний индекс j функции $h_j^i$ обозначает номер строки, а верхний индекс i — номер столбца матрицы), все остальные не приведенные элементы матрицы $H\left(t\right)$ равны нулю.

Из (2.2) следует справедливость следующего утверждения.

Утверждение 2.1. Система (1.7) с условиями (1.8) и (1.9) вполне управляема тогда и только тогда, когда для любого вектора $\sigma \left(t_0,t_1,t_2,T\right)$ можно найти управление $U\left(t\right),$ $t\in \left[t_0,T\right],$ удовлетворяющее условию (2.2).

Управляющее воздействие $U\left(t\right),$ удовлетворяющее интегральному соотношению (2.2), представим в виде (см. [16, c. 116--118], [19, §2])

$U\left(t\right)=H^T\left(t\right)Q^{-1}\sigma \left(t_0,t_1,t_2,T\right)+e\left(t\right),$ (2.3)

где

$Q={\int }_{t_0}^{T}H\left(t\right)H^T\left(t\right)dt,$ (2.4)

а вектор-функция $e\left(t\right)$ удовлетворяет соотношению ${\int }_{t_0}^{T}H\left(t\right)e\left(t\right)dt=0.$

Таким образом, решение задачи можно сформулировать в виде следующей теоремы, аналогичной теореме, доказанной в [19, гл. I].

Теорема 2.1. Для того чтобы существовало программное управление (2.3) и соответствующее ему решение системы (1.7), удовлетворяющее условиям (2.1) (или (2.2)), необходимо и достаточно, чтобы определяемая соотношением (2.4) матрица $Q$ была не особой или чтобы ранги матрицы Q и расширенной матрицы $\left\{Q,\sigma \right\}$ совпадали.  Здесь матрица Q имеет размерность $\left(15\times 15\right)$ и

$\det Q=\frac{g^{12}{T}^{23}{\left(T-t_1\right)}^3{t}_1^3{\left(T-t_2\right)}^{14}{t}_2^{14}}{20808145575327301632000000I_{xx}^{10}{I}_{yy}^{10}{I}_{zz}^4}.$

Отсюда видно, что матрица Q не особая.

Подставляя выражение для функции управления $U\left(t\right)$ в формулу Коши (2.1), получим соответствующее программное движение. Затем, если найденные выражения функции управления $U\left(t\right)$ подставить в (1.7) и проинтегрировать эти уравнения при заданных начальном и промежуточных значениях фазового вектора на каждом промежутке времени, получим законы движения.

Иллюстрацией вышеизложенного служит следующий пример.

Пример 2.1. Пусть $t_0=\mathrm{0,}$ $t_1=\mathrm{2,}$ $t_2=\mathrm{6,}$ $T=\mathrm{12,}$

$x\left(0\right)={\left(\mathrm{0,2,0,5,0,4,6,2,1,1,3,1,0,0,0}\right)}^T,$

$x\left(T\right)={\left(\mathrm{0,5,0,5,0,0,0,0,0,0,0,4,50,80,100}\right)}^T,$

$x_5\left(t_1\right)=\mathrm{50,}{x}_1\left(t_2\right)=\mathrm{80,}{x}_3\left(t_2\right)=100.$    

По вышеприведенным формулам вычислим значения элементов матриц $H^T\left(t\right),$ $Q^{-1}$ и вектора $\sigma$ и подставим их значения в формулы (2.3). Для простоты предполагая, что $e\left(t\right)=\mathrm{0,}$ получим явные выражения для управляющих воздействий в следующем виде:

$u_1\left(t\right)=\left\{\begin{array}{ll}\frac{116}{3}-\frac{53t}{2}& при\le t\le 2\\ -\frac{284}{15}+\frac{23t}{10}& при<t\le 12\end{array}\right.,$

$u_2\left(t\right)=\left\{\begin{array}{ll}-\frac{22021}{1512}+\frac{2935t}{168}-\frac{172535t^2}{36288}+\frac{33445t^3}{93312}& при\le t\le 6\\ \frac{8719}{72}-\frac{76145t}{1512}+\frac{79235t^2}{12096}-\frac{176125t^3}{653184}& при<t\le 12\end{array}\right.,$

$u_3\left(t\right)=\left\{\begin{array}{ll}-\frac{149225}{10584}+\frac{41539t}{3528}-\frac{646183t^2}{254016}+\frac{15107t^3}{93312}& при\le t\le 6\\ \frac{29327}{1176}-\frac{81967t}{10584}+\frac{20017t^2}{28224}-\frac{1757t^3}{93312}& при<t\le 12\end{array}\right.,$

$u_4\left(t\right)=-\frac{25}{8}+\frac{31t}{48}.$

Графики управляющих воздействий приведены ниже (см. рис. 2–5).

 

Рис. 2. График функции управления ${\mathit{u}}_{\mathbf{1}}\mathbf{\left(}\mathit{t}\mathbf{\right)}$

 

Рис. 3. График функции управления ${\mathit{u}}_{\mathbf{2}}\mathbf{\left(}\mathit{t}\mathbf{\right)}$

 

Рис. 4. График функции управления ${\mathit{u}}_{\mathbf{3}}\mathbf{\left(}\mathit{t}\mathbf{\right)}$

 

Рис. 5. График функции управления ${\mathit{u}}_{\mathbf{4}}\mathbf{\left(}\mathit{t}\mathbf{\right)}$

 

Если подставить найденные выражения функции управления $U\left(t\right)$ в уравнения (1.7) и проинтегрировать эти уравнения при заданных начальном и промежуточных значениях фазового вектора на каждом промежутке времени, получим законы фазового вектора. Считаем целесообразным привести только законы движения для геометрических координат, которые имеют следующие явные виды:

$x_1\left(t\right)=\left\{\begin{array}{l}2t+\frac{49t^2}{5}+\frac{49t^3}{10}-\frac{29845t^4}{10368}+\frac{41539t^5}{86400}-\frac{646183t^6}{18662400}+\frac{105749t^7}{111974400}при0\le t\le 6\\ \frac{7378}{25}-\frac{25673t}{75}+\frac{27293t^2}{150}-\frac{23177t^3}{540}+\frac{29327t^4}{5760}-\frac{81967t^5}{259200}+\\ +\frac{20017t^6}{2073600}-\frac{12299t^7}{111974400}при6<t\le 12\end{array}\right.,$

$x_3\left(t\right)=\left\{\begin{array}{l}5t-\frac{147t^2}{5}-\frac{49t^3}{30}+\frac{154147t^4}{25920}-\frac{4109t^5}{2880}+\frac{241549t^6}{1866240}-\frac{46823t^7}{11197440}при0\le t\le 6\\ -\frac{10256}{5}+\frac{35971t}{15}-\frac{18389t^2}{15}+\frac{89299t^3}{270}-\frac{427231t^4}{8640}+\\ +\frac{106603t^5}{25920}-\frac{110929t^6}{622080}+\frac{7045t^7}{2239488}при6<t\le 12\end{array}\right.,$

$x_5\left(t\right)=\left\{\begin{array}{ll}\frac{1}{12}\left(48t+232t^2-53t^3\right)& при0\le t\le 2\\ \frac{1}{60}\left(-2304+3696t-568t^2+23t^3\right)& при2<t\le 12\end{array}\right..$

Графики фазовой вектор-функции $x\left(t\right)$ по геометрическим координатам $x_1\left(t\right),$ $x_3\left(t\right)$ и $x_5\left(t\right)$ при $t\in \left[\mathrm{0,12}\right]$ представлены на рис. 6–8. Траектория движения изображена на рис. 9.

 

Рис. 6. График функции ${\mathit{x}}_{\mathbf{1}}\mathbf{\left(}\mathit{t}\mathbf{\right)}$

 

Рис. 7. График функции ${\mathit{x}}_{\mathbf{3}}\mathbf{\left(}\mathit{t}\mathbf{\right)}$

 

Рис. 8. График функции ${\mathit{x}}_{\mathbf{5}}\mathbf{\left(}\mathit{t}\mathbf{\right)}$

 

Рис. 9. Траектория движения

About the authors

Vanya R. Barseghyan

Institute of Mechanics of National Academy of Science of the Republic of Armenia; Yerevan State University

Author for correspondence.
Email: barseghyan@sci.am
ORCID iD: 0000-0001-6518-3694

Doctor of Physical and Mathematical Sciences, Professor. Leading Scientific Researcher of Institute of Mechanics of NAS of RA; Professor of Mathematics and Mechanics Department, Yerevan State University

Armenia, 24B Marshal Baghramyan Ave., Yerevan 0019; 1 Alec Manukyan St., Yerevan 0025

Tamara A. Simonyan

Yerevan State University

Email: simtom09@gmail.com
ORCID iD: 0000-0003-0434-6183

Candidate of Physical and Mathematical Sciences, Associate Professor, Associate Professor of Mathematics and Mechanics Department

Armenia, 1 Alec Manukyan St., Yerevan 0025

Aram G. Matevosyan

Yerevan State University

Email: matevosaram@gmail.com
ORCID iD: 0009-0006-3432-2253

Candidate of Physical and Mathematical Sciences, Associate Professor, Associate Professor of Mathematics and Mechanics Department

Armenia, 1 Alec Manukyan St., Yerevan 0025

References

Supplementary files