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

Полный текст

Введение

Рассматривается задача оценки множества решений систем обыкновенных дифференциальных уравнений (ОДУ) с возмущениями. Такие задачи возникают при исследовании устойчивости с постоянно действующими возмущениями на конечном интервале времени [1–3], множеств достижимости (МД) управляемых систем [4–11], оценке выживаемости систем [12]. Оценка множеств решений для линейных систем ОДУ с возмущениями актуальна в наше время во многих современных практических сферах, к которым относятся различные системы управления реальными объектами, оценка характеристик объектов в баллистике при массовых вычислениях, управление беспилотными летательными аппаратами и автомобилями [13–14]. Постановка задачи оценки множества решений системы ОДУ с возмущающими воздействиями обусловлена тем, что для рассматриваемой модели отсутствует достаточное описание всех параметров, включающих действующие возмущения (управления) [1, 4, 9–10], а известны только границы изменения возмущений (управлений). В подобных случаях необходимо оценивать все множество решений, соответствующих возмущениям.

Универсального метода получения оценок не существует ввиду многообразия и сложности задач [9]. Существующие методы основаны на различных подходах (например, метод оптимальных двусторонних оценок МД на основе операций над эллипсоидами [4–9], подход работы [10]), в основе которых лежит многократное решение систем ОДУ с измененными значениями параметров, что позволяет учитывать все возможные значения коэффициентов этих систем под воздействием возмущений. В последние годы развиваются исследования в области разработки и применения методов адаптивной интерполяции в задачах с интервальными неопределенностями [15–17], в которых для каждого момента времени строятся кусочно-полиномиальные функции, интерполирующие зависимость решения задачи от точечных значений интервальных параметров с заданной точностью. Авторы метода адаптивной интерполяции применяют $k$-d-деревья для хранения информации о найденных значениях численного решения. С помощью этого метода решены несколько практически важных задач [15–17], в которых требуется определить верхние и нижние границы компонент множества решений.

Следует выделить также работы, описывающие подход, основанный на символьной формуле оператора сдвига вдоль траектории и вычислении множеств значений этой формулы по всем множествам управления и управляющим воздействиям [18] и его применение в ряде задач [19–23].

Иногда возникают трудности, связанные с применением многих методов оценки множеств решений с возмущающими и управляющими параметрами [13–14], которые заключаются в расширении этих оценок независимо от поведения самих решений, а также в возможном увеличении ошибок численных расчетов. Вычисление границы МД может иметь большую вычислительную сложность для систем высокой размерности как в линейном, так и нелинейном случаях. Некоторые методы не позволяют установить параметры и управление, ведущие на границу [13]. В настоящее время разработка надежных вычислительных алгоритмов для оценки МД остается актуальной [14].

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

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

1. Постановка задачи

Рассматривается следующая начальная задача для системы ОДУ:
\[ \begin{equation}
\frac{dy}{dt}=A(t)y(t)+u(t), \quad t \in [t_0, T];
\end{equation} \tag{1} \]
\[ \begin{equation}
y(t_0)=y_0 \in Y_0,
\end{equation} \tag{2} \]
где $y_0 \in Y_0 \subset \mathbb R^n$, $u: [t_0,T] \to [-1,1]$, $y: [t_0,T] \to \mathbb R^n$, $A(t): [t_0,T] \to \mathbb R^n \times \mathbb R^n$, $T>t_0$. Пусть для любого вектора начальных данных существует единственное решение этой задачи $y(t)$ на интервале $[t_0, T]$, $y(t)$ принадлежит соответствующему функциональному пространству (является непрерывной, непрерывно-дифференцируемой или абсолютно непрерывной функцией, в общем случае можно сказать — суммируемой функцией). Требуется найти включение (построить границы) множества решений $Y(t)\subseteq \mathbb R^n$ такое, что
\[ \begin{equation}
y(t) \in Y(t)
\end{equation} \tag{3} \]
на интервале $ [t_0, T] $ для всех решений $y(t)=y(t, t_0, y_0)$ задачи (1) с начальными данными (2) и возмущающим воздействием $u(t)$.

Для построения включения (3) множества решений системы ОДУ с начальными данными воспользуемся оператором Коши $Y(t_0,t) : \mathbb R^n \to \mathbb R^n$, сопоставляющим со значением каждого решения $y(t)$ системы ОДУ в точке $t=t_0$ значение этого же решения в точке $t$:
\[ \begin{equation}
Y(t_0,t) y(t_0)=y(t).
\end{equation} \tag{4} \]

 Выполняются следующие свойства оператора Коши:
\[ \begin{equation*}
Y(t_0,t_0)=0, \quad Y(t_0,t)=Y^{-1}(t,t_0), \quad Y(t_0,\tau)Y(\tau,t)=Y(t_0,t).
\end{equation*} \]

Если система ОДУ (1) является линейной системой с суммируемой правой частью, то формула (4) записывается в виде следующих форм:
\[ \begin{equation}
y(t)=y(t,t_0,y_0)=\Phi(t)\biggl(\Phi^{-1}(t_0)y_0 +\int_{t_0}^t \Phi^{-1}(\tau)u(\tau)d\tau \biggr),
\end{equation} \tag{5} \]
\[ \begin{equation}
y(t)=y(t,t_0,y_0)=\exp(tA)\biggl(\exp^{-1}(t_0A)y_0 +\int_{t_0}^t \exp^{-1}(\tau A)u(\tau)d\tau \biggr).
\end{equation} \tag{6} \]

Формула (5) справедлива для систем линейных ОДУ с переменными коэффициентами, формула (6) выполняется для систем с постоянными матрицами коэффициентов. Здесь $u(\tau)=(u_1(\tau), \dots ,u_n(\tau))$ — вектор внешних воздействий (возмущений); $y_0=(y_{1,0}, \dots, y_{n,0})$ — вектор начальных данных; $\Phi(t)$ — матрица фундаментальных решений линейной однородной системы, полученной из системы (1); $\Phi^{-1}(t)$ — обратная матрица к матрице фундаментальных решений; $\exp(tA)$ — матрица фундаментальных решений линейной однородной системы с постоянными коэффициентами; $\exp^{-1}(tA)$ — обратная матрица к матрице фундаментальных решений линейной однородной системы с постоянными коэффициентами.

На основе формул (5) или (6) нельзя сразу же получить числовые значения. Необходимо выбрать алгоритм и организовать вычисления; величины $\exp(tA)$, $\exp^{-1}(tA)$, $\Phi(t)$, $\displaystyle \int\Phi(t)\Phi^{-1}(t_0)$ — символьные величины, для определения значений которых нужно реализовать соответствующие алгоритмы.

Представим формулы (5) и (6) как соотношения, связывающие решение и неопределенные воздействия. На их основе необходимо выбрать алгоритм и организовать вычисления, позволяющие получить верхние и нижние оценки множества всех решений. При этом величины, подобные $\exp(tA)$, $\exp^{-1}(tA)$, $\Phi(t)$, $\displaystyle \int\Phi(t)\Phi^{-1}(t_0)$, воспринимаются как символьная формула — выражение.

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

2. Описание численно-символьного алгоритма

Первая часть алгоритма описывает построение символьных формул и организацию вычислений по этим формулам включенных в процесс вычисления оценок областей значений. Следует отметить, что оценивать области значений в этой части формулы Коши не требуется, так как полагаем, что начальные данные и все коэффициенты заданы точно. Далее в описании алгоритма мы не будем ссылаться на номера формул (5) или (6), так как эти действия на всех шагах выполняются как для формул (5), так и для формул (6).

Шаг 1. Определяется символьная формула решений фундаментальной матрицы решений матричной системы ОДУ 
\[ \begin{equation}
\frac{d\Phi}{dt}=A(t)\Phi.
\end{equation} \tag{7} \]
Для нахождения матричного решения системы (7) используется метод нахождения полиномиального решения линейной функциональной системы либо метод нахождения рационального решения линейной функциональной системы [25], либо метод построения матричной экспоненты, если $A$ — постоянная матрица. В итоге получаем символьные формулы решений. Эти формулы входят как компоненты в первый член формулы Коши $\Phi(t)\Phi^{-1}(t_0)Y_0$.

Шаг 2. На этом шаге символьная формула $\Phi(t)$ преобразуется к символьным формулам $\Phi(t_0)$, $\Phi^{-1}(t_0)y_0$, выполняется подстановка символьных переменных и их преобразование. Первый член символьной формулы Коши $\Phi(t)\Phi^{-1}(t_0)y_0$ конструируется на основе выражений, полученных после этой подстановки. Размерность произведения матричных и векторных функций равна размерности решения. Под интеграл во втором члене суммы входит также вектор возмущающих воздействий.

Вторая часть шага 2 алгоритма описывает преобразование подынтегрального выражения $\Phi^{-1}(\tau)u(\tau)$, нахождение первообразной преобразованного выражения, затем подстановку первообразной в символьную формулу Коши и оценивание границ области значений этой формулы.

Это достаточно сложная для реализации часть всего алгоритма. Для внешнего воздействия $u(t)$ известны только границы его значений. Поэтому интегрирование неизвестной функции, зависящей от времени, невыполнимо. Можно преобразовать подынтегральное символьное выражение к виду, для которого определяется первообразная. Затем с использованием этой первообразной, оценивается область значений второй части формулы (1), а значит оценивается область значений формулы (1).

Шаг 3. На этом шаге оценивается область значений формулы 
\[ \begin{equation*}
\Phi(t)\int_{t_{0}}^t \Phi^{-1}(\tau)u(\tau)d\tau.
\end{equation*} \]
В ней изменяется один из сомножителей подынтегрального выражения $u(\tau),$ $ | u(\tau)| \leqslant 1$ для любого $\tau$, формулы фундаментальной матрицы решений и обратной матрицы не изменяются.

Чтобы вычислить границы множества значений 
\[ \begin{equation*}
\bigcup_{u(\tau) \in U}\int_{0}^{1}\Phi^{-1}(\tau)u(\tau)d\tau
\end{equation*} \]
при $u(\tau) \in [u_{\min}, u_{\max}]$, в алгоритме определяются экстремальные значения подынтегрального выражения $\displaystyle \Phi(t)\int_0^t\Phi^{-1} (\tau)u(\tau)d\tau$ в фиксированных точках времени $t_k$, $k=1,\dots,n$. Выбор значений возмущающих (управляющих) воздействий как решений задач 
\[ \begin{equation*}
\max_{u_{\min}\leqslant u(\tau) \leqslant u_{\max}}\Phi^{-1}(\tau)u(\tau)
\quad \text{и}
\quad 
\min_{u_{\min}\leqslant u(\tau) \leqslant u_{\max}}\Phi^{-1}(\tau)u(\tau)
\end{equation*} \]
в фиксированные моменты времени эквивалентен тому, что возмущающее воздействие выбирается из условия приближения подынтегрального выражения сверху, затем — из условия приближения снизу. При этом используются методы, в которых применяются только значения целевой функции и не вычисляются производные [26].

Для графического представления границы множества значений и ее оценки строится интерполяционный полином по точкам максимума подынтегрального выражения 
\[ \begin{equation*}
\Phi(t)\Phi^{-1}(\tau)u(\tau)=\Phi(t_{\rm const})\Phi^{-1}(\tau)u(\tau),
\end{equation*} \]
вычисленным в точках $t_{k}$; значение $t=t_{\rm const}$ фиксировано, значения переменной $\tau$ изменяются в точках $t_k$:
\[ \begin{equation*}
 \max_{u_{\min}\leqslant u(t_1) \leqslant u_{\max}}\Phi^{-1}(t_1)u(t_1),\quad 
 \dots, 
 \quad \max_{u_{\min}\leqslant u(t_n) \leqslant u_{\max}}\Phi^{-1}(t_n)u(t_n).
\end{equation*} \]
Для интерполирования максимальных значений применяется кубический интерполяционный сплайн [27], обозначенный здесь как $\operatorname{MaxICS}(t)$, $t_0\leqslant t \leqslant t_n$. Кубический интерполяционный сплайн $\operatorname{MinICS}(t)$, $t_0\leqslant t \leqslant t_n$, 
строится аналогично по точкам
\[ \begin{equation*}
 \min_{u_{\min}\leqslant u(t_1) \leqslant u_{\max}}\Phi^{-1}(t_1)u(t_1),\quad 
 \dots, 
\quad 
\min_{u_{\min}\leqslant u(t_n) \leqslant u_{\max}}\Phi^{-1}(t_n)u(t_n).
\end{equation*} \]
Так как все входящие в подынтегральное выражение функции векторные, данная процедура повторяется для каждой компоненты построенных векторных функций.

Для многих задач возмущающее воздействие [5] может не быть непрерывной, а также дифференцируемой функцией. Возможные варианты выбора возмущающих (управляющих) воздействий: кусочно-гладкая функция, производная которой интегрируема по Лебегу, либо абсолютно непрерывная функция. Поэтому необходимо оценить точность построения интерполяционной или сглаживающей сплайн-функции $\operatorname{MaxICS}(t)$, $t_0 \leqslant t \leqslant t_n$, вычисленной по точкам максимальных значений подынтегральных выражений, а также сплайн-функции $\operatorname{MinICS}(t)$, $t_0 \leqslant t \leqslant t_n$, вычисленной по точкам минимальных значений подынтегральных выражений. Компоненты подынтегрального выражения — непрерывные функции, зависящие от $t$.

Определение. Для любой функции $f(t)$, заданной на отрезке $T$, верхняя грань колебаний функции [27] по всем подотрезкам из $T$ длиной меньше $\delta$ является модулем непрерывности этой функции:
\[ \begin{equation*}
\omega(f,\delta) = \sup \{|f(t_1)-f(t_2)|, \; t_1, t_2 \in T, \; |t_1-t_2|<\delta \}.
\end{equation*} \]

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

Утверждение. Пусть на отрезке $[t_0,T]$ заданы непрерывная функция $y(t)$ и такая последовательность сеток $\Delta_k$, что $\| \Delta_k \| \to 0$ при $k \to \infty$. Если кубический сплайн дефекта 2 интерполирует функцию $y(t)$ в узлах сетки $\Delta_k$ и его производная в этих узлах равна нулю, то последовательность $\{S_{\Delta_k(t)} \}$ сходится к $y(t )$ равномерно на $[t_0, T]$. Выполняется неравенство
\[ \begin{equation}
\| y(t)- S_{\Delta_k(t)} \| \leqslant 2\omega \Bigl( \frac{\| \Delta \}}{2}, y \Bigr).
\end{equation} \tag{8} \]

Доказательство утверждения и оценки (8) проводится методами из [27]. Дополнительным ограничением интерполяционного сплайна является требование равенства производной сплайна нулю в узлах сетки. Поскольку по сути задачи нас более всего интересуют оценки значений подынтегрального выражения в узлах сетки, равенство производной интерполяционного сплайна нулю в добавленных узлах сетки не оказывает существенного воздействия и не усложняет реализацию алгоритма. В противном случае возможно построение сглаживающего сплайна [27].

В выбранных узлах $t_k$, $k=1,\dots, n$, функции $\operatorname{MaxICS}(t)$ и $\operatorname{MinICS}(t)$ принимают экстремальные значения, совпадающие со значениями какой-либо траектории задачи (1), (2) в этих узлах. В промежуточных точках интерполяции к значениям интерполирующих функций нужно добавить или отнять оценку ошибки интерполяции (вычисляемая величина) согласно (8). Ниже на шаге 4 описан способ контроля за оценкой отклонения решений в пучке траекторий. Поскольку используется символьная формула точного решения, ошибок приближенного интегрирования системы при реализации этого алгоритма нет. Далее в алгоритме строятся формулы первообразных 
\[ \begin{equation*}
\operatorname{IMaxICS}(t)=\int_{t_{0}}^t \operatorname{MaxICS}(\tau) d\tau,
\quad 
\operatorname{IMinICS}(t)=\int_{t_{0}}^t \operatorname{MinICS}(\tau) d\tau.
\end{equation*} \]

Шаг 4. Для контроля точности алгоритма оценим максимальные отклонения любых двух решений, принадлежащих пучку траекторий
\[ \begin{equation*}
y(t)=\Phi(t)\biggl(\Phi^{-1}(t_0)y_0+\int_{t_0}^t \Phi^{-1}(\tau)u(\tau)d\tau \biggr).
\end{equation*} \]
Изменяющаяся величина в формуле решения — возмущающее воздействие $u(t)$, остальные параметры решений совпадают. Цель шага 4 заключается в демонстрации эквивалентности максимальных отклонений любых двух решений и максимальных отклонений решений, вычисляемых в данном методе по формуле Коши. При этом максимальные отклонения в первом из рассматриваемых примеров не будут превосходить вычисляемые в методе отклонения произвольных значений двух решений.

Преобразуем разность двух решений этого пучка, выходящих из одной начальной точки, но имеющих разные управляющие воздействия:
\[ \begin{multline*}
y_{1}(t)-y_{2}(t)= \Phi(t)\Phi^{-1}(t_0)y_0+\Phi(t)\int^t_{t_0}\Phi^{-1}(\tau)u_1(\tau)d\tau - \Phi(t)\Phi^{-1}(t_0)y_0-\Phi(t)\int^t_{t_0}\Phi^{-1}(\tau)u_2(\tau)d\tau =\\\
=\Phi(t)\biggl(\int^t_{t_0}\Phi^{-1}(\tau)u_1(\tau)d\tau - \int^t_{t_0}\Phi^{-1}(\tau)u_2(\tau)d\tau\biggr)=\Phi(t)\int_{t_0}^{t} \Phi^{-1}(\tau)(u_{1}(\tau)-u_{2}(\tau))d\tau.
\end{multline*} \]
Используя свойства интеграла Лебега в этом выражении, оценим покомпонентно модуль разности двух решений:
\[ \begin{equation*}
|y_1(t)-y_2(t) \leqslant 2 | \Phi(t) |
\int_{t_0}^{t} |\Phi^{-1} (\tau)| |u_1(\tau)-u_{2}(\tau)|d\tau. 
\end{equation*} \]
Далее, 
\[ \begin{equation*}
| \Phi^{-1}(\tau)(u_{1}(\tau)-u_2(\tau)| \leqslant |\Phi^{-1}(\tau)| |{u_{1}(\tau)-u_{2}(\tau)}| \leqslant 2 | \Phi^{-1}(\tau)|.
\end{equation*} \]

В этом неравенстве учтено, что $u_{1}(t) \in [-1, 1]$, $u_{2}(t) \in [-1, 1]$. Отсюда для любого $t$ справедлива оценка $|u_1(t)-u_2(t)| \leqslant 2$.

Это означает, что отклонение любых двух решений пучка траекторий покомпонентно оценивается так:
\[ \begin{equation}
|y_1(t)-y_2(t)| \leqslant 2|\Phi(t)|\int_{t_0}^t |\Phi^{-1}(\tau)|d\tau.
\end{equation} \tag{9} \]

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

Метод, представленный в этой статье, оценивает формулу множества точных решений как формулу объединения всех выражений, полученных по формуле Коши:
\[ \begin{equation}
\bigcup_{u\in [-1,1]} y(t,t_0,y_0)=
\bigcup_{u\in [-1,1]} \Phi(t)\biggl(\Phi^{-1}(t_0)y_0 + \int^t_{t_0}\Phi^{-1}(\tau)u(\tau)d\tau \biggr).
\end{equation} \tag{10} \]
Каждый элемент объединения (10) является точным решением задачи (1), (2). Основная сложность реализации алгоритма заключается в оценке области значений на основе символьных формул этого выражения. Эта сложность преодолена — этап описан на шаге 3.

При этом значение оценок выражения вычисляется сразу в выбранный момент времени. Нет необходимости накапливать (суммировать) множества значений выражений, что возникает при реализации любого численного метода оценки множеств решений. Это свойство является важной характеристикой алгоритма, составленного на основе формулы Коши. Явное достоинство методов, основанных на символьных формулах решений (10), заключается также в отсутствии накопления ошибок со всех временных шагов. В (10) отсутствуют какие-либо промежуточные численные значения, появляющиеся в любом численном методе вида 
\[ \begin{equation*}
y_k+h\Psi(y_k,y_{k-1}, \dots, y_{k-m}),\quad 0 \leqslant k \leqslant N,\quad 0 \leqslant m \leqslant k.
\end{equation*} \]

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

Шаг 5. Печать таблиц значений и построение графиков решений.

Для анализа вычислительной сложности символьного метода на основе оператора Коши отметим следующее. Шаг 3 алгоритма требует пересчета векторной численно-символьной величины $\Phi(t)\Phi(\tau)u(\tau)$ для каждого момента времени $t$ и определения ее экстремальных значений для всех $\tau$, пробегающих узлы сетки от 0 до $t$. Следовательно, асимптотическая зависимость числа операций от количества узлов сетки будет оцениваться как $O( N(N-1)/2)$, где $N$ — число узлов сетки. Поскольку система является линейной и в методе не используется рекуррентная зависимость решения в каждый момент времени от неопределенных параметров (известна формула решения — формула Коши), в каждом узле сетки выполняется число операций с асимптотикой $O(n)$, где $n$ — размерность системы.

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

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

3. Примеры вычисления оценок множеств решений

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

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

Рассмотрим управляемую систему, описываемую системой ОДУ четвертого порядка
\[ \begin{equation}
\begin{array}{ll}
\dfrac{dy_1}{dt}=y_2, &
\dfrac{dy_2}{dt}=-y_1+y_2+u, \\
\dfrac{dy_3}{dt}=y_4, &
\dfrac{dy_4}{dt}=y_1-y_2-w.
\end{array}
\end{equation} \tag{11} \]
В системе (11) $y_1$, $y_3$ — компоненты вектора состояния, описывающие положение управляемого объекта; $y_2$, $y_4$ — компоненты вектора, описывающие величину скорости управляемого объекта; на возмущения $u(t)$ и $w(t)$ наложены ограничения 
\[ \begin{equation*}
u(t), w(t) \in U=\{|u| \leqslant 1, |w | \leqslant 1\}. 
\end{equation*} \]
Начальные данные для системы (11) следующие: $y_i(t_0)=1$, $i= 1, 2, 3, 4$. В уравнения движения справа входит вектор равнодействующей всех сил в системе.

Выбор возмущающих или управляющих воздействий системы (11) заранее неизвестен. Поэтому ставится задача определить включение $Y_h(t)$ для пучка траекторий точных решений системы (11):
\[ \begin{equation*}
Y_h(t) \supseteq \{y(t,t_0,u,w) \mid u \in U , \; w \in U\}.
\end{equation*} \]

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

На основе символьного метода проведены оценки множеств достижимости задачи (11).

На рис. 1, 2 изображены верхние и нижние границы включений области достижимости для первой и третьей компоненты вектора решений, то есть в проекциях на оси $y_1$, $y_3$. На этих рисунках включены только два объекта: верхняя и нижняя границы множества решений в проекции на координатные плоскости.

Рис. 1. Верхняя и нижняя границы множества достижимости первой компоненты вектора состояний на интервалах $[0{, \,}4]$ и $[0{, \,}8]$ для управляемой системы с возмущениями (11)
[Figure 1. The upper and lower bounds of the reachable set of the first component of the state vector are defined on the intervals $[0{, \,}4]$ and $[0{, \,}8]$ for the controlled system with perturbations (11)]

Рис. 2. Верхняя и нижняя границы множества достижимости третьей компоненты вектора состояний на интервалах $[0{, \,}4]$ и $[0{, \,}8]$ для управляемой системы с возмущениями (11)
[Figure 2. The upper and lower bounds of the reachable set of the third component of the state vector are defined on the intervals $[0{, \,}4]$ and $[0{, \,}8]$ for the controlled system with perturbations (11)]

Рис. 3. Верхняя и нижняя границы множества достижимости первой компоненты вектора состояний на интервале $[0{, \,}16]$ для управляемой системы с возмущениями (12)
[Figure 3. The upper and lower bounds of the reachable set of the first component of the state vector are defined on the intervals $[0{, \,}16]$ for the controlled system with perturbations (12)]

Рис. 4. Верхняя и нижняя границы множества достижимости первой компоненты вектора состояний в задаче вертикального взлета на интервале $[0{, \,}10]$
[Figure 4. The upper and lower bounds of the set of attainable values of the first component of the state vector in the vertical takeoff problem within the interval $[0{, \,}10]$]

Оценка множеств решений выполнена также для управляемой системы четвертого порядка [6]:
\[ \begin{equation}
\begin{array}{ll}
\dfrac{dy_1}{dt}=2y_2+y_4, & \dfrac{dy_2}{dt}=-2y_1+y_3+u_1, 
\\
\dfrac{dy_3}{dt}=y_2+10y_4, & \dfrac{dy_4}{dt}=y_1-2y_3+u_2
\end{array}
\end{equation} \tag{12} \]
с возмущающими силами $|u_1(t)| \leqslant 1$, $|u_2(t)| \leqslant 1$. В [6] построены графики роста объема эллипсоида траекторий системы, входящих в пучок точных решений задачи (12). Полученные численно-символьным методом результаты оценки множеств решений всех компонент вектора состояния задачи (12) совпадают с качественными оценками в [6]. Это подтверждает график на рис. 3.

Рассмотрим задачу оценки множества достижимости при плоском движении самолета, при котором траектории его центра масс расположена в некоторой фиксированной вертикальной плоскости, служащей плоскостью материальной симметрии самолета. Силами, действующим на самолет, являются сила тяги винта, направленная по оси винта и составляющая с хордой крыла постоянный угол, сила тяжести и аэродинамические силы. Совокупность последних может быть приведена к главному вектору. Движение летательного аппарата в вертикальной плоскости на некотором интервале времени можно описать следующей линейной системой обыкновенных дифференциальных уравнений с постоянными коэффициентами [28]:
\[ \begin{equation}
\begin{array}{ll}
\dfrac{dy}{dt}= a_{13}v+a_{14}\theta, & 
\dfrac{dx}{dt}= a_{23}v+a_{24}\theta,\\
\dfrac{dv}{dt}=a_{35} |u | +a_{30}, &
\dfrac{d\theta}{dt}=a_{45}u+a_{40},
\end{array}
\end{equation} \tag{13} \]
где $y$, $x$ — координаты местоположения летательного аппарата, $v$ — скорость, $\theta$ — угол наклона траектории, $a_{ij}$ — динамические коэффициенты. Система (13) получена в [28] из нелинейных дифференциальных уравнений движения летательного аппарата заменой правых частей функциями линейными относительно фазовых координат и управлений. Управляющее воздействие $u$ — значения угла атаки, удовлетворяющие ограничению
\[ \begin{equation*}
-a_{\max}\leqslant u \leqslant a_{\max},\quad a_{\max} = 0.25.
\end{equation*} \]
Динамические коэффициенты в работе [28] заданы как величины, имеющие размерности. Они принимают значения $a_{13}=0$, $a_{14}=1000$ м/с, $a_{23}=1$, $a_{24}=-1$ м/с, $a_{35}=-50$ м/с$^2$, $a_{30}=-20$ м/с$^2$, $a_{45}=0.25$ 1/с, $a_{40}=0$. Требуется построить оценку области достижимости летательного аппарата в вертикальной плоскости $Oxy$ для заданного момента времени.

При построении границ множества достижимости использовался символьный метод, описанный в этой статье. Для символьного приближения интеграла в формуле Коши применялось два алгоритма. Один алгоритм расчета основан на аппроксимации управляющего воздействия кусочно-постоянной функцией, второй — на приближении интеграла символьной квадратурной формулой. Динамические коэффициенты, например скорости, принимают большие значения, что означает возможность быстрого роста границ множества значений координат летательного аппарата. Численные значения границ множества решений этой задачи подтверждаются при сравнении со значениями формул точных решений (рис. 4).

Заключение

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

Разработан численно-символьный алгоритм построения формул, описывающих точные решения систем ОДУ, организующий упрощения этих формул и вычисления их экстремальных значений. В алгоритм также входят неравенства (9), связывающие максимальные отклонения решений и возмущающие воздействия на шаге 4 алгоритма. Результаты этого этапа нужны для контроля за точностью численно-символьного алгоритма. Построенный численно-символьный алгоритм применим на всем интервале времени, на котором существуют решения всех систем ОДУ, входящих в поставленную задачу оценки множеств решений. 

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

Конкурирующие интересы. Конкурирующих интересов не имею.
Авторский вклад и ответственность. Я несу полную ответственность за предоставление окончательной версии рукописи в печать. Окончательная версия рукописи мною одобрена.
Финансирование. Работа выполнена в рамках проекта, поддержанного за счет средств гранта № 20–41–240002, предоставленного РФФИ, Правительством Красноярского края, Краевым фондом науки.

Об авторах

Александр Алексеевич Рогалев

Институт космических и информационных технологий Сибирского федерального университета

Автор, ответственный за переписку.
Email: gogoba88@mail.ru
ORCID iD: 0000-0003-2176-9639
SPIN-код: 1313-8673
Scopus Author ID: 57212867169

старший преподаватель, каф. информационных систем

Россия, 660074, Красноярск, ул. Академика Киренского, 26 к/1

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

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