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

全文:

1. ВВЕДЕНИЕ

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

Среди методов построения оптимального управления можно условно выделить два класса. В первый класс входят методы, основанные на принципе максимума Понтрягина, теории устойчивости Ляпунова, теории Флоке и теории игр [1–3]. Эти методы можно назвать локальными или тактическими, так как строится некоторое номинальное управление и управление, стабилизирующее движение в окрестности номинальной траектории. Во второй класс входят методы, основанные на принципе оптимальности Беллмана и динамическом программировании [4–6]. Эти методы можно назвать глобальными или стратегическими, так как здесь управление ищется в виде отображения из состояния в управляющие воздействия, причем управление должно «вести себя хорошо» в широком диапазоне значений аргумента.

В последнее время активно развивается и привлекает внимание исследователей раздел приближенного динамического программирования, называемый машинным обучением с подкреплением [7–9]. В обучении с подкреплением управление динамическими системами интерпретируется как взаимодействие агента с внешней средой, от которой агент за свои действия получает вознаграждения и стремится максимизировать суммарные вознаграждения. В механике космического полета средой может являться космический аппарат, агентом — управляющее программное обеспечение на аппарате, вознаграждением могут быть точность прилета в заданную область пространства и экономия топлива. Функция управления (отображение из состояния аппарата в управляющие воздействия) параметризуется, параметры ищутся так, чтобы удовлетворялось уравнение оптимальности Беллмана и / или вознаграждение, полученное агентом за весь полет, было максимальным в среднем по возможным начальным условиям старта. Результатом обучения является функция управления, которая способна направлять аппарат в заданную точку пространства. Эта функция может быть загружена на борт космического аппарата и может управлять им во время реального полета на основе состояния аппарата или оценок состояния аппарата.

Математически строгая и основанная на ляпуновском подходе теория обучения с подкреплением для построения управления изложена в монографии [9]. В частности, там описываются методы построения функции оптимального управления для классических вариантов функционалов — интегралов от квадратичной функции управления и невязок по состоянию. Формулируются теоремы о сходимости параметрически заданной функции управления к оптимальной. Методы излагаются для случая непрерывных динамических систем. Изложенная там теория имеет свои недостатки:

1. Делаются существенные предположения о динамике, функции управления и функции Беллмана, например условие богатства входного сигнала (постоянства возбуждения, persistence of excitation). Эти условия зачастую не удается доказать или проверить, а их невыполнение сказывается на сходимости к оптимальной функции управления.
2. Рассматриваются только функции вознаграждения, квадратичные относительно управления и положительно определенные относительно вектора невязки. Требуется обратимость матрицы квадратичной формы относительно управления.
3. В формулах коррекции параметров приближенных функции управления и функции Беллмана используются правые части уравнений движения.
4. Теория, построенная в монографии, касается непрерывных динамических систем. Для дискретных систем (например, в механике космического полета — случай движения с импульсами) необходимо проводить адаптацию этих методов.

Следует отметить, что существуют положительные примеры применения ляпуновского подхода к обучению нейросетевых моделей управления космическими аппаратами [10–13]. Эти примеры не относятся к обучению с подкреплением, но методы, изложенные в них, схожи с теми, что описываются в работе [9]. Ляпуновский подход к обучению нейросетевых моделей также известен в литературе под названиями детерминированное обучение [14] и нейродинамическое программирование [15].

В последние несколько лет область обучения с подкреплением пополнилась эффективными алгоритмами, зарекомендовавшими себя в разных областях, в том числе и в механике космического полета (см. обзор литературы по теме в публикации [16] в разделе Reinforcement learning). Эти численные методы основываются на алгоритмах приближенного динамического программирования, методах оптимизации функций с большим числом параметров и теории частично наблюдаемых марковских процессов принятия решений. Преимуществом этих методов является существенное сокращение математических предположений и значительный охват возможных решаемых задач. Примеры их применения показывают, что стратегии управления, создаваемые этими методами, естественным образом способны адаптироваться к неизвестным параметрам аппарата и внешней среды [17–20].

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

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

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

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

2. МАТЕМАТИЧЕСКИЕ ОСНОВЫ ОБУЧЕНИЯ С ПОДКРЕПЛЕНИЕМ

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

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

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

Будем считать, что агент взаимодействует со средой в дискретные моменты времени $k\in \{\mathrm{0,1,2,}\dots \}$. Пусть S — множество всех возможных состояний среды, A — множество всех возможных действий агента. Эти множества могут быть дискретными или непрерывными, конечными или бесконечными. Обозначим за $S_k$ состояние среды в момент k, $A_k$ — действие, которое производит агент в момент k.

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

$p:\mathcal{S}\times \mathcal{A}\times \mathcal{S}\to \left[\mathrm{0,1}\right],$

которая для всех $s\in \mathcal{S}$, $a\in \mathcal{A}$, $s^{\text{'}}\in \mathcal{S}$ и $k\ge 0$ есть ¹

$p(s,a,s^{\text{'}})=\mathbb{P}(S_{k+1}=s^{\text{'}}|S_k=s,A_k=a).$

Будем считать, что вероятность перехода между состояниями не зависит от времени (такую среду называют стационарной).

Вознаграждение $R_k$ в момент времени k определяется распределением $d_R$, зависящем только от текущего состояния $S_k$, действия $A_k$ и, возможно, будущего состояния $S_{k+1}$ и не зависящего от времени k. Будем также считать, что вознаграждения равномерно по исходам ограничены сверху: $\left|R_k\right|$ ≤ R_max.

Введем начальное распределение состояний $d_0:\mathcal{S}\to \left[\mathrm{0,1}\right],$ то есть функцию, такую, что для всех s

$d_0\left(s\right)=\mathbb{P}(S_0=s).$

Опишем теперь правило, согласно которому в среде действует агент. Стратегией называется функция $\pi :\mathcal{S}\times \mathcal{A}\to \left[\mathrm{0,}1\right],$ которая для всех $s\in \mathcal{S}$, $a\in \mathcal{A}$ и $k\ge 0$ равна

$\pi (s,a)=\mathbb{P}(A_k=a|S_k=s).$

Таким образом, стратегия для каждого s определяет распределение на множестве действий, действие выбирается случайно в соответствии с этим распределением. Технически, действие выбирается в результате применения генератора случайных чисел, отвечающего данному распределению. Стратегия может быть детерминированной, в этом случае для всех s и a функция p(s, a) принимает значения из множества $\left\{\mathrm{0,1}\right\}$ и под стратегией можно понимать просто функцию a = p(s). Будем считать, что стратегия не зависит от времени t, а действие определяется лишь состоянием s, в котором находится среда.

Итак, рассмотрим процесс взаимодействия агента со средой. Сначала в соответствии с начальным распределением $d_0$ инициализируется начальное состояние $S_0~d_0$ (знак  означает, что случайная величина слева в выражении генерируется из распределения справа в выражении). Затем для этого состояния в соответствии со стратегией агент производит действие $A_0~\pi (S_0,\cdot )$. Далее, в соответствии с функцией перехода  среда переходит в новое состояние $S_1~p(S_0,A_0,\cdot )$. После этого агент получает вознаграждение $R_0~d_R(S_0,A_0,S_1)$. Далее агент производит новое действие $A_1~\pi (S_1,\cdot )$, среда переходит в новое состояние $S_2~p(S_1,A_1,\cdot )$, агент получает вознаграждение $R_1~d_R(S_1,A_1,S_2)$, и так далее. Взаимодействие со средой заканчивается либо в определенный момент времени $k=K$, либо при достижении средой некоторого особого состояния. Участок времени взаимодействия агента со средой от начального до финального состояния называется эпизодом.

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

$J\left(\pi \right)=\mathbb{E}\left(\left.\sum _{k=0}^{\infty }{R}_k\right|\pi \right),$

где вертикальная черта означает, что все действия, которые производит агент, производятся им в рамках стратегии p. Так как любой эпизод конечен, то и ряд представляет собой конечную сумму (формально можно считать, что $R_k=0$, начиная с некоторого момента времени, зависящего от исхода). Стратегия p* называется оптимальной, если

$\pi *\in \arg \underset{\pi \in \Pi }{\max }J\left(\pi \right),$

где $\Pi$ — множество рассматриваемых стратегий. В общем случае оптимальная стратегия может не существовать, а если существует, то может быть не единственной.

В определении целевой функции участвует сумма случайных величин. Чтобы сделать эту сумму ограниченной для любого исхода, вводят величину $\gamma \in \left[\mathrm{0,1}\right]$ и рассматривают целевую функцию $J\left(\pi \right)=\mathbb{E}\left(\left.\sum _{k=0}^{\infty }{\gamma }^k{R}_k\right|\pi \right),$ сумма в которой конечна с вероятностью единица, если y < 1 и $\left|R_k\right|$ ≤ R_max.

Функцией ценности состояния называется функция $v^{\pi }:\mathcal{S}\to ℝ$, для каждого $s\in \mathcal{S}$ равная

$v^{\pi }\left(s\right)=\mathbb{E}\left(\left.\sum _{k=0}^{\infty }{\gamma }^k{R}_{k+m}\right|S_m=s,\pi \right).$ (1)

Эта функция выражает среднее суммарное вознаграждение за эпизод, получаемое агентом, действуя в рамках стратегии p, стартуя из состояния s в момент m. Обратим внимание, что эта функция не зависит от m, и в ее определении допустимо брать $m=0$. Независимость функции ценности от m является следствием стационарности среды, стационарности вознаграждения и марковского свойства процесса.

Функцией ценности действия называется функция $q^{\pi }:\mathcal{S}\times \mathcal{A}\to ℝ$, для каждого s и a равная

$q^{\pi }(s,a)=\mathbb{E}\left(\left.\sum _{k=0}^{\infty }{\gamma }^k{R}_{k+m}\right|S_m=s,A_m=a,\pi \right).$

Эта функция выражает среднее суммарное вознаграждение за эпизод, получаемое агентом, стартуя из состояния s с действием a в момент m и далее действуя в рамках стратегии p. В этом выражении можно брать m произвольным.

В частично наблюдаемых марковских процессах принятия решений агент принимает решение о действиях не на основе состояний, а на основе наблюдений, которые являются функциями состояний (возможно случайными). Если наблюдение $o=\phi \left(s\right)$ является детерминированной и взаимно однозначной функцией состояния, то в качестве состояния можно выбрать наблюдение и ввести марковский процесс принятия решений, в котором стратегия может быть определена как

$\pi (o,a)=\mathbb{P}(A_k=a|O_k=o),$ (2)

где $O_k$ — вектор наблюдений в момент , а функция ценности —

$v^{\pi }\left(o\right)=\mathbb{E}\left(\left.\sum _{k=0}^{\infty }{\gamma }^k{R}_{k+m}\right|O_m=o,\pi \right).$ (3)

Если $o=\phi \left(s\right)+\xi$, где $\phi$ — взаимно однозначная детерминированная функция состояния, а $\xi$ — случайный вектор, то наблюдению  может соответствовать множество возможных состояний s. При достаточно малой дисперсии компонент $\xi$, допустимо рассматривать стратегию и функцию ценности как функции наблюдения по формулам (2)–(3), если это слабо влияет на результаты методов оптимизации стратегии. Если функция ценности используется только во время оптимизации стратегии, то есть играет только вспомогательную роль, то допустимо определять ее по формуле (1).

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

$o,h\to \phi \left(s\right)+\xi ,$

где $\phi$(s) — взаимно однозначная функция состояния, а $\xi$ — случайный вектор с малыми дисперсиями компонент. История наблюдений h может состоять из одного, двух или многих наблюдений, предшествующих наблюдению o. Этот вектор может быть равен и равносильной функции предыдущих наблюдений. В любом случае, стратегия вводится как отображение

$\pi (o,h,a)=\mathbb{P}(A_k=a|O_k=o,H_k=h),$

где $H_k$ обозначает вектор истории наблюдений к моменту k не включительно, а h — реализация этого вектора. Аналогично определяется функция ценности состояния

$v^{\pi }(o,h)=\mathbb{E}\left(\left.\sum _{k=0}^{\infty }{\gamma }^k{R}_{k+m}\right|O_m=o,H_m=h,\pi \right).$

Наблюдению и истории здесь все еще неоднозначно соответствует состояние, но оценка этого состояния тем точнее, чем меньше дисперсия компонент ошибки $\xi$, которая обычно связана с частотой наблюдений и объемом истории наблюдений.

3. МЕТОДЫ ОПТИМИЗАЦИИ СТРАТЕГИЙ

Существует множество алгоритмов поиска оптимальных стратегий. К классическим методам относятся методы итерации по стратегиям и ценности с известными условиями сходимости [6, 7]. Однако область применения этих методов ограничена средами с дискретным и конечным множеством состояний и агентами с дискретным и конечным множеством действий, в то время как для механики космического полета характерны задачи с непрерывными множествами состояний и действий. Для таких случаев разрабатываются методы приближенного динамического программирования, общую теорию о которых можно найти в монографиях [8, 15].

Обзор методов показывает, что среди них можно выделить два класса методов: основанные на использовании агентом функции перехода между состояниями (задачи планирования, model-based methods), и методы, в которых агент не использует эту функцию для расчета оптимальных действий (безмодельные методы, model-free methods). Характерными чертами первого класса методов являются эффективность и точность. Примером применения таких методов может служить работа по разработке гарантирующего синтеза управления для управления космическим аппаратом в окрестности неустойчивой точки либрации [21]. Методы же второго класса не используют специфические свойства динамики и потому более универсальны. Из наиболее известных — метод градиента глубокой детерминированной стратегии (Deep Deterministic Policy Gradient, DDPG) [22], метод асинхронного исполнителя–критика (Asynchronous Advantage Actor Critic, A3C) [23], метод оптимизации ближайшей стратегии (Proximal Policy Optimization, PPO) [24]. Перечисленные методы являются градиентными, среди безградиентных методов для оптимизации стратегий можно применять эволюционные алгоритмы [25, 26]. Обзор применения безмодельных методов к задачам механики космического полета можно найти в работе [16].

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

В современных методах обучения с подкреплением стратегия представляется в виде параметрически заданной функции с конечным числом параметров:

$\pi =\pi (a,s,\theta )$ или $a=\pi (s,\theta )$,

где $\theta$ — вектор оптимизируемых параметров. Задача поиска оптимальной стратегии сводится к оптимизации функционала

$J\left(\theta \right)=\mathbb{E}\left(\left.\sum _{k=0}^{\infty }{\gamma }^k{R}_k\right|\pi \left(\theta \right)\right).$

Рассмотрим один из вариантов оптимизации этого функционала, применяемый, например, методом PPO. Математическое ожидание в выражении для J(q) заменяется на выборочное среднее (среднее арифметическое). Производится серия испытаний Монте-Карло, в каждой серии среда инициализируется в начальном состоянии, агент при фиксированных значениях параметров q производит действия, получает за них вознаграждения R_t, и действует до конца эпизода. Так в серии испытаний получаются реализации суммарных вознаграждений за эпизод, их среднее дает оценку J(q). Далее значения параметров q с использованием конкретного метода оптимизации (например, PPO) корректируются в сторону повышения значения функционала J и процесс сбора данных повторяется снова. Процесс оптимизации останавливается, когда значение функционала перестает увеличиваться.

Интересно заметить, что для коррекции значений параметров в сторону повышения значения функционала J не требуется рассчитывать производные вознаграждений, действий, управляющих воздействий или состояний по параметрам q. Согласно теореме о градиенте стратегии (policy gradient theorem) [27] в силу стационарности среды и марковского свойства градиент J содержит только производную по стратегии p:

${\nabla }_{\theta }J\left(\theta \right)\propto {\mathbb{E}}_{S~d_0,A~\pi \left(\theta \right)}\left[q^{\pi }\right(S,A\left){\nabla }_{\theta }\ln \pi \right(\theta \left)\right],$

где $\propto$ означает пропорциональность. Это выражение используется в градиентных методах оптимизации стратегий.

4. МЕТОДИКА ПОСТРОЕНИЯ УПРАВЛЕНИЯ С ИСПОЛЬЗОВАНИЕМ ОБУЧЕНИЯ С ПОДКРЕПЛЕНИЕМ

Перейдем теперь к описанию методики построения управления. Пусть динамика механической системы описывается системой обыкновенных дифференциальных уравнений

$\dot{x}=f(t,x,u),t\in [t_0,t_f],$ (4)

где t — время, $x\in {ℝ}^n$ — вектор состояния, $u\in {ℝ}^m$ — вектор управления.

Рассматриваются функции управления вида $u=u(t,x)$ или, в более общем случае, вида $u=u(t,o,h)$, где $o=o(t,x)\in {ℝ}^k$ — вектор наблюдения, $h\in {ℝ}^l$ — вектор, характеризующий историю наблюдений. Назовем управление допустимым, если оно:

1) переводит каждое решение уравнений (4) с начальным условием из заданного множества $(t_0,x(t_0\left)\right)\in {\Omega }_0$ в множество заданных краевых условий $(t_f,x(t_f\left)\right)\in {\Omega }_f$, причем эти решения удовлетворяют промежуточным ограничениям $\left\{\right(t,x\left(t\right)),t_0<t<t_f\}\in {\Omega }_{\mathrm{int}}$,

2) удовлетворяет ограничениям

$\left\{\right(t,u(t,x(t\left)\right)),t_0\le t\le t_f\}\in {\Omega }_U$

или

$\left\{\right(t,u(t,o(t),h(t\left)\right)),t_0\le t\le t_f\}\in {\Omega }_U$.

Пусть определен функционал

$\mathcal{J}=\mathcal{J}(u,t_0,x(t_0),t_f,x(t_f\left)\right),$ (5)

который допустимому управлению u, переводящему решение с начальным условием $(t_0,x(t_0\left)\right)\in {\Omega }_0$ в краевое условие $(t_f,x(t_f\left)\right)\in {\Omega }_f$, ставит в соответствие число.

Рассматривается задача поиска функции управления $u=u(t,x)$ или, в более общем случае, $u=u(t,o,h)$, для каждых $(t_0,x(t_0\left)\right)\in {\Omega }_0$ и $(t_f,x(t_f\left)\right)\in {\Omega }_f$ оптимизирующей функционал (5).

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

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

1. Связать понятие состояния из теории обучения с подкреплением с понятием состояния механической системы. В общем случае состоянием из теории обучения с подкреплением можно считать $s=(t,x)$ или любую взаимно однозначную функцию от $(t,x)$. Если система (4) автономная, то состоянием можно считать x или любую взаимно однозначную функцию от x.
2. На области начальных условий $(t_0,x(t_0\left)\right)\in {\Omega }_0$ определить распределение вероятностей ${\mathcal{D}}_0$, в соответствии с которым в серии испытаний Монте-Карло будут генерироваться начальные условия для обучения стратегии. Это может быть равномерное распределение, нормальное распределение или какое-либо другое распределение. Это распределение должно отражать ожидания разработчиков и сопровождающих миссию касательно того, в каких областях пространства и с какой вероятностью будет находиться аппарат в момент начала действия управления.
3. Так как алгоритмы обучения с подкреплением, рассматриваемые в настоящей работе, оперируют средами с дискретным временем, следует перейти от динамики с непрерывным временем (4) к динамике с дискретным временем. Это значит, что необходимо определить отображение $(t_k,x_k,u_k)\to (t_{k+1},x_{k+1})$ для каждого дискретного шага k. Например, это отображение можно определить с использованием схемы интегрирования Эйлера

$x_{k+1}=x_k+hf(t_k,x_k,u_k),t_{k+1}=t_k+h,$

задав некоторый малый шаг h. В общем случае шаг на участке может вычисляться с использованием любого численного метода интегрирования. Если управление представляет собой импульсы скорости, под дискретным шагом можно понимать целую траекторию между импульсами, рассчитываемую методом интегрирования. Алгоритмам обучения с подкреплением обычно требуется информация о том, является ли новое состояние $k+1$ финальным. Например, это состояние может быть финальным, если $t_{k+1}=t_f$ или если в момент времени $t_{k+1}<t_f$ происходит событие, означающее конец траектории движения (столкновение с небесным телом, выход за допустимые границы движения и т. п.). Поэтому в общем случае следует задавать отображения вида $(t_k,x_k,u_k)\to (t_{k+1},x_{k+1},d_{k+1})$, где $d_{k+1}=0$, если состояние $k+1$ не является финальным, и $d_{k+1}=1$, если это состояние является финальным. Интервал времени от начального состояния до финального является эпизодом в терминах обучения с подкреплением.

4. Определить функцию вознаграждения, которая будет представлять собой отображение $(t_k,x_k,u_k,t_{k+1},x_{k+1},d_{k+1})\to r_k.$ Здесь при расчете сигнала вознаграждения учитывается то, в каком состоянии находилась система до управления, вектор управления, состояние, в которое перешла система, и информация о том, является ли новое состояние финальным. В качестве функции вознаграждения можно рассматривать функции вида

$r_k=-\alpha \left|u_k\right|-\rho (t_{k+1},x_{k+1},{\Omega }_f),$ (6)

где r — это определенное исследователем расстояние от точки $(t_{k+1},x_{k+1})$ до множества W_f краевых условий, а a — задаваемая постоянная. Например, если ищется управление, переводящее космический аппарат в состояние, характеризуемое положением $r_f$ и скоростью $v_f$, и вектор состояния представляет собой положение и скорость аппарата $x=[r,v]$, то в качестве функции p можно взять $\rho (x,{\Omega }_f)=|r-r_f|+|v-v_f|$. Еще один вариант функции вознаграждения:

$r_k=-\alpha \left|u_k\right|+\rho (t_k,x_k,{\Omega }_f)-\rho (t_{k+1},x_{k+1},{\Omega }_f).$

В этом случае на каждом шаге сигнал вознаграждения содержит информацию о приближении к краевым условиям. Заметим, что в этом случае при a = 0 суммарные вознаграждения за эпизод совпадают с величиной, на которую улучшается расстояние до краевых условий:

$R=\sum _{k=1}^K{r}_k=\rho (t_0,x_0,{\Omega }_f)-\rho (t_{K+1},x_{K+1},{\Omega }_f),$

где K — число шагов в эпизоде. В некоторых работах функция вознаграждения состоит из трех слагаемых:

$r_k=r_{k,\text{int}}+r_{k,f}+r_{k,\text{good}},$

где $r_{k,\text{int}}$ состоит из значения функционала $\mathcal{J}$ и вознаграждения за удовлетворение ограничениям на траекторию, $r_{k,f}$ — вознаграждение за удовлетворение краевым условиям (подсчитывается при $d_{k+1}=1$), и $r_{k,\text{good}}$ — вознаграждение за следование в окрестности «хорошей траектории» (например, оптимальной траектории, полученной в упрощенной постановке).

5. Определить модель восприятия, то есть отображение $(t_k,x_k)\to o_k$ из состояния в вектор наблюдений или измерений. В простейшем случае наблюдение может совпадать с состоянием: $o_k=x_k$. В общем случае это отображение моделирует работу датчиков и результат применения навигационных процедур, и потому именно наблюдение и история наблюдений используются для расчета действия и управляющих воздействий. Так наблюдением может являться оценка состояния, которая может моделироваться как $o_k=x_k+{\xi }_k$, где x_k — случайный вектор. Наблюдением может быть изображение, в этом случае стратегия будет отображать наблюдение-изображение непосредственно в управляющие воздействия, минуя стадию навигации (оценки состояния по изображениям).
6. Наконец, определить модель управления — параметрически заданное отображение $o_k\to u_k$ или $(o_k,h_k)\to (u_{k+1},h_{k+1})$. Модель управления состоит из композиции двух отображений — из наблюдения в действие $o_k\to a_k$ или $(o_k,h_k)\to (a_{k+1},h_{k+1})$ и из действия в управление $a_k\to u_k$. Отображение из наблюдения в действие обычно строят на основе нейросетевых моделей. Так, в случае отображений вида $o_k\to a_k$, чаще всего используют многослойные нейронные сети прямого распространения, например сети с одним скрытым слоем

$a=A_2{\phi }_1(A_{1}o+b_1)+b_2$

или двумя скрытыми слоями

$a=A_3{\phi }_2\left(A_2{\phi }_1\right(A_{1}o+b_1)+b_2)+b_3,$

где $A_1$, $A_2$, $A_3$, $b_1$, $b_2$, $b_3$ — матрицы и векторы оптимизируемых параметров, j₁, j₂ — активационные функции. Выбор размеров матриц и векторов параметров, а также активационных функций остается за исследователем. Универсальные теоремы аппроксимации [28–31] утверждают, что выбором достаточно большого числа параметров и произвольных активационных функцией из широкого множества нелинейных функций можно добиться аппроксимации любой гладкой функции. Однако эти теоремы не говорят о том, сколько параметров следует брать, и каковы их значения. Выбор числа параметров и активационных функций может значительно влиять на точность аппроксимации и качество получаемой стратегии. В случаях $(o_k,h_k)\to (u_{k+1},h_{k+1})$ можно использовать рекуррентные нейронные сети, а $h_k$ считать скрытым состояние рекуррентного слоя, несущим в себе информацию об истории наблюдений. Вместо рекуррентных слоев можно использовать слои прямого распространения, но на вход сети подавать конкатенацию векторов наблюдений. В общем случае модель не обязана быть нейросетевой. Модель управления может быть построена на основе ляпуновского управления, выведенного в рамках упрощенной модели движения системы. В таком управлении чаще всего есть параметры, которые можно сделать обучаемыми (оптимизируемыми). Некоторые алгоритмы обучения с подкреплением для коррекции параметров стратегии используют также параметрическую модель функции ценности. Эту модель как правило выбирают нейросетевой, причем архитектура (число слоев, активационные функции, число нейронов в слоях) обычно совпадает с архитектурой стратегии. Что касается отображения из действия в управление, вводят функцию $u=\psi \left(a\right)$. В механике космического полета функция y часто вводится, чтобы ограничивать значения, которые может принимать управление, а также в целях нормировки и масштабирования значений a.

Процесс обучения представлен на рис. 1. Вначале инициализируются параметры модели управления $\pi (o,h,\theta )$, где q содержит все оптимизируемые параметры (разбирается наиболее общий случай управления по наблюдениям и истории наблюдений). Далее начинается оценивание этой стратегии в серии испытаний Монте-Карло. Для этого в соответствии с заданным распределением ${\mathcal{D}}_0$ в области начальных условий W₀ инициализируется состояние системы $(t_0,x_0)$. Затем модель восприятия сопоставляет этому состоянию наблюдение $o_0$. Наблюдение подается на вход стратегии, получается вектор управления $u_0$ и история наблюдений (или скрытое состояние рекуррентного слоя) $h_1$. Состояние и управление подаются на вход отображения дискретного шага системы и функции вознаграждения, получаются новое состояние, вознаграждение и флаг, сигнализирующий о конце или продолжении эпизода. Если эпизод не закончен, цикл повторяется заново для нового состояния. Если эпизод закончен, инициализируется новое состояние. Этот цикл (серии испытаний Монте-Карло) продолжается много раз, в результате накапливается история взаимодействия со средой: состояния s, действия a, вознаграждения r и флаги завершения эпизода d. Эта история подается на вход алгоритму оптимизации, который корректирует параметры стратегии (этим алгоритмом может быть PPO, DDPG, A3C, генетический алгоритм и любой другой алгоритм оптимизации для обучения с подкреплением). Параметры стратегии корректируются, и серия испытаний Монте-Карло повторяется для исправленной стратегии.

Рис. 1. Процесс обучения стратегии.

5. ОЦЕНКА КАЧЕСТВА РАБОТЫ ФУНКЦИИ УПРАВЛЕНИЯ

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

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

Теорема. Пусть $X_1$, …, $X_n$ — независимые случайные величины, для которых выполнено a_i ≤ X_i ≤ ≤ b_i с вероятностью единица. Тогда для среднего выборочного $\overline{X}=(1/n)\sum _{i=1}^n{X}_i$ справедлива оценка

$\mathbb{P}\left(\right|\overline{X}-\mathbb{E}\overline{X}|\ge \epsilon )\le 2\exp \left(-\frac{2{\epsilon }^2{n}^2}{{\displaystyle \sum _{i=1}^n{(b_i-a_i)}^2}}\right).$

Следствие. Если $X_i$ одинаково распределены, a_i ≤ X_i ≤ b_i для всех i и $\mathbb{E}{X}_1=$m, то

$\mathbb{P}\left(\right|\overline{X}-\mu |\ge \epsilon )\le 2\exp \left(-\frac{2{\epsilon }^{2}n}{{(b-a)}^2}\right).$

Это неравенство означает, что проведя n независимых измерений случайной величины X, ограниченной промежутком $[a,b]$, мы получим, что вероятность отклонения среднего выборочного $\overline{X}$ от истинного математического ожидания m более чем на e не превосходит p = = $2\exp (-2{\epsilon }^{2}n/{(b-a)}^2)$. Неравенство Хефдинга удобно записывать в виде доверительного интервала:

$\mathbb{E}X=\overline{X}\pm \epsilon$ с вероятностью не менее $1-p$.

Величину e можно назвать точностью определения математического ожидания случайной величины, а $1-p$ называется уровнем доверия.

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

Ниже в табл. 1 для случая $a=0$, $b=1$ приводятся формулы для расчета величин e, p, n, когда даны любые две из этих величин, а также значения объема измерений n для различных значений e, p:

$\begin{array}{l}n(\epsilon ,p)=\frac{1}{2{\epsilon }^2}\ln \left(\frac{2}{p}\right),\\ p(\epsilon ,n)=2\exp (-2{\epsilon }^{2}n),\\ \epsilon (p,n)=\sqrt{\frac{1}{2n}\ln \left(\frac{2}{p}\right)}.\end{array}$

Если вместо промежутка $\left[\mathrm{0,1}\right]$ рассматривается промежуток $[a,b]$, то объем выборки следует увеличить в ${(b-a)}^2$ раз, чтобы получить те же значения p, e оценки величины.

Таблица 1. Таблица значений n, e, p в случае a = 0, b = 1

Теперь применим теорему Хефдинга к оценке вероятности интересуемого исследователя события A. Примером такого события может быть то, что невязка по краевым условиям в конце эпизода управления увеличится. Другой пример — движение космического аппарата вышло в фазовом пространстве за допустимые границы, и миссия потеряна.

Пусть $\mathcal{I}\left(A\right)$ — индикатор события A, то есть $\mathcal{I}\left(A\right)=1$, если A произошло, и $\mathcal{I}\left(A\right)=0$, если A не произошло. Пусть $p_A$ — неизвестная вероятность того, что A происходит. Введем независимые одинаково распределенные случайные величины $X_i$, имеющие распределение Бернулли $\text{Be}\left(p_A\right)$. Тогда $\mu =\mathbb{E}{X}_i=p_A$, и так как 0 ≤ X_i ≤ 1, то согласно неравенству Хефдинга

$\mathbb{P}\left(\right|\overline{X}-p_A|\ge \epsilon )\le 2\exp \left(-2{\epsilon }^{2}n\right).$

Вероятность того, что ошибка оценки $p_A$ превысит e, равна $p=2\exp (-2{\epsilon }^{2}n)$. Например, для e = 0.01, $p=0.01$ получается $n=26492$. Это значит, что проведя $n=26492$ измерений случайной величины X, мы получим, что вероятность отклонения истинной вероятности события A от оценки $\overline{X}$ более чем на $1\%$ не превосходит $1\%$. Таким образом, неравенство Хефдинга можно рассматривать как инструмент оценки требуемого числа испытаний для расчета вероятности события.

В заключение отметим, что вероятность событий и распределение исследуемых случайных величин зависит от распределения ${\mathcal{D}}_0$ на множестве начальных условий. Это распределение может быть разным во время обучения стратегии и во время ее тестирования.

6. АВТОРСКАЯ ПРОГРАММНАЯ БИБЛИОТЕКА KIAM_RL

Описанная выше методика построения управления механическими системами была воплощена автором в виде программной библиотеки kiam_rl. Библиотека написана на языке Python и на момент написания статьи состоит из двух модулей: routines.py и ppo_hyperparameters_tuning.py. Далее следует описание их возможностей.

Модуль routines.py содержит базовые классы и функции для создания сред и моделей стратегий. Модуль содержит абстрактный класс RLProblem, расширяющий возможности пакета gymnasium [33], позволяющего создавать стандартные среды, с которыми оперируют популярные программные библиотеки алгоритмов обучения с подкреплением, например stable-baselines3 [34], который содержит реализации алгоритмов PPO, DDPG и др. Класс RLProblem содержит несколько методов:

1) observation_space, возвращающий gymnasium.spaces.Box-объект, определяющий множество состояний s;

2) action_space, возвращающий gymnasium.spaces.Box-объект, определяющий множество действий a;

3) initialize, возвращающий начальные время $t_0$ и фазовое состояние $x_0$;

4) equations_of_motion, принимающий на вход время t, фазовое состояние x, вектор управления u и возвращающий правые части уравнений движения, то есть вектор $dx/dt$;

5) step, принимающий на вход время t, фазовое состояние x, вектор управления u и возвращающий результаты дискретного шага: новое время $t^{\text{'}}$, новое состояние $x^{\text{'}}$ и флаги конца эпизода, первый из которых сигнализирует о нормальном завершении эпизода, второй — о вынужденном завершении эпизода;

6) reward, принимающий на вход момент времени t, состояние x, новый момент времени $t^{\text{'}}$, новое состояние $x^{\text{'}}$, вектор управления u, флаги конца эпизода и возвращающий вознаграждение;

7) action2u, принимающий на вход действие a и возвращающий управление u;

8) perception_model, принимающий на вход время t и состояние x и возвращающее вектор наблюдения o.

Пользователь создает класс, описывающий его задачу, наследуя класс RLProblem, и самостоятельно наполняет указанные методы отображениями согласно постановке его задачи. На основе этого класса с использованием Environment модуля routines.py создается стандартный gymnasium.Env-объект среды, который можно далее использовать в сочетании с популярными библиотеками алгоритмов оптимизации стратегий. Таким образом, класс RLProblem представляет собой интерфейс между механическими аспектами задачи и обучением с подкреплением и помогает создавать и использовать стандартизированные среды.

Программный вид модели управления не стандартизирован, его пользователь создает самостоятельно. Если предполагается использовать алгоритмы библиотеки stable-baselines3, то модели должны создаваться с помощью пакета pytorch [35]. В этом может помочь класс ActorCriticNetworks модуля routines.py. Пользователь создает класс, наследуя его из класса ActorCriticNetworks, и самостоятельно определяет модель стратегии в конструкторе.

При создании и оптимизации моделей стратегии пользователю необходимо задать опции алгоритмов: архитектуру, глубину и ширину нейросетевых моделей, скорость обучения, число эпизодов для оценивания функционала  и многие другие параметры. Эти параметры называются гиперпараметрами, чтобы отличать их от параметров модели. Результаты оптимизации стратегии могут существенно зависеть от выбранных значений гиперпараметров, поэтому зачастую встает задача оптимизации гиперпараметров. Поиск разумных значений гиперпараметров пользователь может производить вручную. Существуют и автоматические процедуры оптимизации гиперпараметров, основанные на методах оптимизации вычислительно затратных функций, например байесовской оптимизации [36], древесно-структурированной оценки Парзена [37]. Существуют программные библиотеки, позволяющие оптимизировать гиперпараметры, например Optuna [38], Ray [39], BoTorch [40], Hyperopt [41]. Модуль ppo_hyperparameters_tuning.py использует процедуры оптимизации гиперпараметров на основе библиотеки Optuna и предполагает, что оптимизация параметров модели осуществляется алгоритмом PPO.

7. ПРИМЕРЫ ПРИМЕНЕНИЯ МЕТОДИКИ

В данном разделе демонстрируется применение методики построения управления. Рассматриваются два примера — задача стабилизации движения в простой динамической системе и астродинамическая задача поддержания движения аппарата в окрестности неустойчивой гало-орбиты вокруг точки либрации. Обе задачи были с использованием описанной выше библиотеки kiam_rl.

Простая динамическая система

Ставится задача построения управления динамической системой

$\dot{x}=u$

на интервале времени $t\in \left[\mathrm{0,1}\right]$, с начальными условиями $x_0\in [-\mathrm{1,1}]$, краевым условием $x\left(1\right)=0$ и управлением $u=u\left(x\right)\in [-\mathrm{1,1}]$. Оптимизируемый функционал:

$\mathcal{J}={\int }_0^1\left|x\right|dt\to \min ,$

его оптимизация равносильна минимизации времени достижения системой состояния $x=0$. Заметим, что оптимальным управлением в этой задаче является функция

$u\left(x\right)=-\text{sign}\left(x\right),$

где $\text{sign}$ — функция знака. Легко показать, что соответствующей функцией ценности (функцией, сопоставляющей начальному условию значение функционала $\mathcal{J}$) является функция

$v\left(x\right)=-x^2/2.$

Для поиска управления воспользуемся описанной выше методикой. Будем считать, что состоянием из теории обучения с подкреплением является переменная x, то есть $s=x$. На области начальных условий ${\Omega }_0=[-\mathrm{1,1}]$ определим равномерное распределение вероятностей. Дискретный шаг определим так:

$x_{k+1}=x_k+h{u}_k,t_{k+1}=t_k+h,$

и будем считать состояние $x_{k+1}$ финальным, если $t_{k+1}=1$. Величину шага по времени выберем $h=0.01$. В качестве функции вознаграждения выберем

$r_k=-\left|x_k\right|h.$

В этом случае будет искаться управление, максимизирующее среднее значение величины

$R=-\sum _{k=1}^K\left|x_k\right|h,$

где $K=1/h=100$ — число шагов. Будем считать наблюдение состоянием, то есть $o=x$. В качестве модели управления рассмотрим

$u=\max \left(\min \right(a\mathrm{,1}),-1),$

$a(x,\theta )={\theta }_3\text{th}({\theta }_{1}x+{\theta }_2)+{\theta }_4,$

где действие a моделируется как полносвязная нейронная сеть прямого распространения с одним скрытым слоем и гиперболическим тангенсом в роли функции активации. Для обучения этой модели воспользуемся алгоритмом обучения PPO. Этот алгоритм требует также задания модели для функции ценности, определим ее следующим образом:

$v=v(x,w)=A_{2}th(A_{1}x+b_1)+b_2,$

где матрицы и векторы $A_1\in {ℝ}^{5\times 1}$, $b_1\in {ℝ}^5$, $A_2\in {ℝ}^{1\times 5}$, $b_2\in ℝ$. Это тоже полносвязная нейронная сеть с одним скрытым слоем и гиперболическим тангенсом в роли активационной функции, но с пятью нейронами на скрытом слое.

Обучение моделей произведем с помощью реализации метода PPO из библиотеки stable-baselines3. Объем выборки для аппроксимации среднего значения функционала r (опция n_steps) выберем равным 1000. Число итераций градиентного метода для коррекции весов нейросетевых моделей (опция n_epochs) выберем равным 100, а скорость обучения (опция learning_rate) – 0.01. В общем случае объем выборки для аппроксимации среднего значения функционала влияет на устойчивость процесса сходимости, а коэффициент скорости сходимости и число итераций градиентного метода оптимизации влияют на устойчивость и скорость сходимости к оптимальному решению. Обучение будем производить на центральном процессоре и завершим, когда число дискретных шагов достигнет 100000 (опция total_timesteps).

В результате оптимизации получается модель

$\begin{array}{c}u(x,{\theta }^{*})=\max \left(\min \right(1.5035\cdot \text{th}(-21.8299x+\\ +0.0612)-\mathrm{0.0725,1}),-1),\end{array}$

которая приближает оптимальную функцию

$-\text{sign}\left(x\right)=1\cdot \text{th}(-\infty \cdot x+0)+0.$

На рис. 2 и 3 показаны графики теоретически оптимальной и приближенной функций управления и ценности. Графики показывают близость приближенных функций к теоретическим. Эта близость регулируется шагом дискретизации h и богатством и удачностью выбора параметрических моделей управления и функции ценности.

Рис. 2. Оптимальная (синий цвет) и приближенная (красный цвет) функции управления.

Рис. 3. Оптимальная (синий цвет) и приближенная (красный цвет) функции ценности.

Поддержание движения в окрестности гало-орбиты

Рассмотрим теперь задачу поддержания движения космического аппарата в окрестности неустойчивой гало-орбиты вокруг точки либрации $L_1$ системы Земля – Луна. В качестве модели движения выберем круговую ограниченную задачу трех тел, а уравнения движения запишем во вращающейся системе координат: начало правой системы координат поместим в центр Земли, ось x направим вдоль направления Земля – Луна, ось z направим вдоль угловой скорости орбитального движения Луны вокруг Земли. В качестве единицы расстояния выберем расстояние между Луной и Землей, а единицы частоты — орбитальную частоту движения Луны вокруг Земли. Уравнения движения в этом случае запишутся следующим образом:

$\dot{x}=v_x,\dot{y}=v_y,\dot{z}=v_z,$

${\dot{v}}_x=2v_y+U_x,{\dot{v}}_y=-2v_x+U_y,{\dot{v}}_z=U_z,$

где

$U(x,y,z)=\frac{{(x-\mu )}^2+y^2}{2}+\frac{1-\mu }{r_1}+\frac{\mu }{r_2}+\frac{\mu (1-\mu )}{2},$

$r_1=\sqrt{x^2+y^2+z^2},r_2=\sqrt{{(x-1)}^2+y^2+z^2},$

а $U_x$, $U_y$, $U_z$ означают частные производные функции $U=U(x,y,z)$ по x, y, z соответственно. Здесь $\mu =m_{\text{M}}/(m_{\text{E}}+m_{\text{M}})$ — массовый параметр, $m_{\text{E}}$ — масса Земли, $m_{\text{M}}$ — масса Луны. Используются следующие значения массового параметра, единицы расстояния $\text{DU}$, единицы скорости  и единицы времени $\text{TU}$:

$\mu =1.215058446035100\cdot {10}^{-2}$, $\text{DU}=384405$ км,

$\text{VU}=1.024540192302405$ км/с,

$\text{TU}=4.342564574695797$ дней².

Рассматривается движение космического аппарата вблизи гало-орбиты вокруг точки либрации $L_1$ с максимальной z-координатой равной $z_{\max }=34981$ км (рис. 4).

Рис. 4. Гало-орбита (синий цвет), в окрестности которой рассматривается движение аппарата. Серым цветом показаны проекции орбиты на плоскости xy, xz, yz.

Начальное условие, отвечающее этой орбите, есть

$x_{ref\mathrm{,0}}=[x_{ref\mathrm{,0}},y_{ref\mathrm{,0}},z_{ref\mathrm{,0}},v_{x,ref\mathrm{,0}},v_{y,ref\mathrm{,0}},v_{z,ref\mathrm{,0}}],$

$x\text{ref},\text{0}=\mathrm{0,826890333820514}\text{\hspace{0.33em}y}\text{ref},\text{0}=0.z_{\text{ref,0}}=0.091$

$v_{x,\text{ref},\text{0}}=\mathrm{0,}{v}_{y,\text{ref},\text{0}}=\mathrm{0.205889408677437,}{v}_{z,\text{ref},\text{0}}=0.$

Период орбиты равен $P_{\text{ref}}=2.78227853520921$ безразмерных единиц времени, то есть приблизительно 12 дней. Параметризуем точки орбиты параметром $\tau \in \left[\mathrm{0,}{P}_{\text{ref}}\right]$ так, что $x_{ref}\left(\tau \right)$ — фазовое состояние на орбите в момент времени t = t и $x_{ref}\left(0\right)=x_{ref\mathrm{,0}}$. Под окрестностью орбиты будем понимать область фазового пространства

$\begin{array}{c}{\Omega }_{vic}=\{x=[r,v]\in {ℝ}^6:\exists \tau \in [\mathrm{0,}{P}_{ref}\left]\right|r-r_{ref}\left(\tau \right)|\le R_{vic},\\ |v-v_{ref}(\tau \left)\right|\le V_{vic}\},\end{array}$

где r_ref(t) и v_ref(t) — положение и скорость в векторе x_ref (t), а размеры окрестности по положению и скорости выбраны равными $R_{vic}=100$ км, $V_{vic}=0.1$ м/с. Наконец, будем считать, что состояние аппарата в моменты управления известно со среднеквадратичными ошибками ${\sigma }_r=1$ км по положению и ${\sigma }_v=0.01$ м/с по скорости.

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

Сформулируем и решим эту задачу в терминах теории обучения с подкреплением. Состояние из теории обучения с подкреплением определим как фазовый вектор механической системы, то есть $s=x=[x,y,z,v_x,v_y,v_z]$.

Инициализация состояния в окрестности орбиты происходит следующим образом. Орбита дискретизируется, то есть представляется в виде конечного числа точек x_ref (t_i), $i=\mathrm{1,}\dots \mathrm{,1000}$, с равномерным по t разбиением. Случайным образом выбирается одна из 1000 точек x_ref (t_i) и в соответствии с равномерным распределением в шаре радиуса 100 км вокруг нее генерируется отклонение по положению dr, а в шаре радиуса 0.1 м/с — скорость dv. Начальное состояние определяется как $s=x_{ref}\left({\tau }_i\right)+[\delta r,\delta v]$.

Траектория аппарата получается численным интегрированием уравнений движения. Для этого используется реализация метода Рунге – Кутты 8-го порядка с адаптивным шагом DOP853 [42]. Дискретный шаг представляет собой интегрирование уравнений движения на интервале времени $[\mathrm{0,}{P}_{ref}/4]$. Будем считать на этапе обучения, что эпизод состоит из одного шага. Из-за неустойчивости движения выбор большего числа шагов в эпизоде приводит к быстрому удалению траектории от орбиты и процесс обучения на получаемых данных затрудняется. Тестирование обученной стратегии можно осуществлять на эпизодах с большим числом шагов.

Функция вознаграждения определяется следующим образом:

$r_1=\max (-\mathrm{1,}-|\Delta x_1|\cdot {10}^3),$

где $\left|\Delta x_1\right|.$ — это минимальное расстояние в фазовом пространстве от проинтегрированного состояния $x_1$ до орбиты, то есть

$\left|\Delta x_1\right|=|x_1-x_{ref}({\tau }_{i^{*}}\left)\right|,i^{*}={\mathrm{argmin}}_i|x_1-x_{ref}({\tau }_i\left)\right|.$

Нормировочный коэффициент ${10}^3$ и ограничение снизу -1 выставляются для того, чтобы значения вознаграждения лежали в интервале $[-\mathrm{1,0}]$. Нормировка и ограничение вознаграждения положительно влияют на процесс обучения и аппроксимацию функции ценности. Кроме того, ограничение функции вознаграждения позволяет определять объем выборки для доверительного оценивания суммарных вознаграждений (см. раздел 5), и чем меньше дисперсия суммарного вознаграждения, тем более устойчивым является процесс обучения. Так как эпизод состоит из одного шага, то суммарное вознаграждение за эпизод $R=r_1$.

Наблюдением считается вектор

$o=\left[\right(x^{\text{'}}-x_{ref}\left({\tau }_{i^{*}}\right))\cdot {10}^3,\cos {\phi }_{i^{*}},\sin {\phi }_{i^{*}}]\in {ℝ}^8,$

где оценка состояния $x^{\text{'}}=x+\xi$, $\xi ~\mathcal{N}\left(\mathrm{0,}{\Sigma }_{nav}\right)$ — вектор ошибок определения состояния, ${\Sigma }_{\text{nav}}=\text{diag}({\sigma }_r,{\sigma }_r,{\sigma }_r,{\sigma }_v,{\sigma }_v,{\sigma }_v)$ — ковариационная матрица ошибок навигации, $i^{*}={\mathrm{argmin}}_i|x^{\text{'}}-x_{ref}({\tau }_i\left)\right|$ — номер ближайшей к $x^{\text{'}}$ точки на орбите, ${\phi }_{i^{*}}=2\pi \cdot {\tau }_{i^{*}}/P_{ref}$, ${10}^3$ — нормировочный коэффициент. Первые шесть компонент этого вектора представляют собой отклонение в фазовом пространстве от ближайшей точки на орбите. Две последние компоненты определяют абсолютное положение этой точки на орбите в пространстве.

Выбор такого вектора наблюдения вызван необходимостью нормирования входного вектора в функцию стратегии для сходимости процесса обучения. Выбор $o=x+\xi$ в качестве вектора наблюдения приводит к тому, что значения этого вектора распределены в относительно узкой окрестности орбиты и стратегия во время обучения не способна различить близкие значения этого вектора. Поэтому здесь предлагается использовать масштабированную локальную информацию об отклонении от орбиты и информацию об абсолютном положении точки, в окрестности которой происходит наблюдение. Выбор компонент cos j и sin j обусловлено желанием обеспечить периодичность управления по этому углу.

Действие и импульсы скорости свяжем равенством

$\Delta v=a\cdot 6\cdot {10}^{-4}/VU.$

Эта нормировка делается для того, чтобы значения компонент действия a лежали в пределах $[-\mathrm{1,1}]$, что способствует сходимости процесса оптимизации. Компоненты импульсов скорости лежат в пределах от $-0.6$ до $0.6$ м/с. Нормировочный коэффициент подбирался автором вручную апостериорно. Модель для действий определим следующим образом:

$a=A_{2}th(A_{1}o+a_1)+a_2,$

где $o\in {ℝ}^n$ — вектор наблюдения, а матрицы и векторы $A_1\in {ℝ}^{16\times 8}$, $a_1\in {ℝ}^{16}$, $A_2\in {ℝ}^{3\times 16}$, $a_2\in {ℝ}^3$ — обучаемые параметры модели. Модель функции ценности определим с похожей архитектурой:

$V=b_{2}th(B_{1}o+b_1)+c_2,$

где $B_1\in {ℝ}^{16\times 8}$, $b_1\in {ℝ}^{16}$, $b_2\in {ℝ}^{16}$, $c_2\in ℝ$. Модель действий имеет таким образом 16 ⋅ 8 +16 + 3 ⋅ 16 + + 3 = 195 параметров, модель функции ценности имеет $16\cdot 8+16+1\cdot 16+1=161$ параметр.

Обучение моделей произведем с помощью реализации метода PPO из библиотеки stable-baselines3. Объем выборки для аппроксимации среднего значения функционала R (опция n_steps) выберем равным 10000. Число итераций градиентного метода для коррекции весов нейросетевых моделей (опция n_epochs) выберем равным 30, а скорость обучения (опция learning_rate) – 0.005. Обучение будем производить на центральном процессоре и завершим, когда число дискретных шагов достигнет 10 млн (опция total_timesteps).

На рис. 5 показано среднее суммарное вознаграждение за эпизод как функция шага. Вознаграждение в среднем растет по мере обучения моделей и в конце оптимизации колеблется на величину $\pm 0.01$ и в среднем равно $0.29$. На рис. 6 изображен график зависимости от номера шага среднеквадратичного отклонения модели функции ценности от средних суммарных вознаграждений. Отклонение падает с ростом числа шагов, финальное значение равно $7.7\cdot {10}^{-4}$. Указано также значение отклонения после первого плато значений: $5.0\cdot {10}^{-3}$.

Рис. 5. Среднее вознаграждение за эпизод в зависимости от шага.

Рис. 6. Среднеквадратичное отклонение модели функции ценности от средних суммарных вознаграждений.

Оценка качества обученной модели управления производилась в эпизодах с последовательными 4 шагами (полный виток в окрестности орбиты). Оценивалась вероятность выхода аппарата за окрестность орбиты в результате маневрирования, а также затраты характеристической скорости на каждый импульс. Так как значение вероятности принадлежит интервалу $\left[\mathrm{0,1}\right]$, для оценивания ее с точностью $0.1\%$ на уровне доверия $99.9\%$ согласно табл. 1 достаточно произвести 3800452 измерения (шага), что соответствует 950113 эпизодам. Максимальное значение импульса скорости равно $\sqrt{3\cdot {0.6}^2}\le 1.04$ м/с, поэтому для оценки затрат скорости с точностью до 0.001 м/с достаточно смоделировать $950113\cdot {1.04}^2=1027643$ эпизода.

В итоге было смоделировано 1 028 000 эпизодов. Результаты оценок промаха по положению мимо окрестности орбиты и затраты характеристической скорости приведены в табл. 2–5. Здесь $q_0$ — минимальное значение величины из встреченных; $q_{0.25}$ — величина, ниже которой находятся 25% встреченных величин; $q_{0.5}$ — медиана; $q_{0.75}$ — величина, выше которой находятся 25% встреченных величин; $q_1$ — максимальное значение величины из встреченных; m — среднее арифметическое всех величин. Средние затраты характеристической скорости для каждого импульса равны:

1) первый импульс: $0.229\pm 0.001$ м/с с вероятностью не менее 99.9%,

2) второй импульс: $0.265\pm 0.001$ м/с с вероятностью не менее 99.9%,

3) третий импульс: $0.192\pm 0.001$ м/с с вероятностью не менее 99.9%,

4) четвертый импульс: $0.143\pm 0.001$ м/с с вероятностью не менее 99.9%.

Таблица 2. Квантили и средние значения распределений промаха по положению мимо орбиты Dr и затраты характеристической скорости Dv после первого импульса

Таблица 3. Квантили и средние значения распределений промаха по положению мимо орбиты Dr и затраты характеристической скорости Dv после второго импульса

Таблица 4. Квантили и средние значения распределений промаха по положению мимо орбиты Dr и затраты характеристической скорости Dv после третьего импульса

Таблица 5. Квантили и средние значения распределений промаха по положению мимо орбиты Dr и затраты характеристической скорости Dv после четвертого импульса

Так как суммарные затраты топлива за виток ограничены величиной 4.16 м/с, то с учетом $n=1028000$ в соответствии с неравенством Хефдинга можно рассмотреть e = 0.008 и $p=2\exp (-2{\epsilon }^{2}n)=0.001$, то есть $0.829\pm 0.008$ м/с с вероятностью не менее $99.9\%$. Экстраполяция этих значений на один год дает $25.21\pm 0.25$ м/с. Для сравнения, поддержание близкой по размерам орбиты в миссии ARTEMIS в годовом выражении потребовало порядка 7.39 м/с [43]. Затраты на поддержание можно уменьшить, если добавить в вознаграждение затраты скорости на импульс и выбрать коэффициент перед соответствующим слагаемым. В настоящей работе такое исследование не проводится.

Результаты оценки вероятности выйти за пределы окрестности орбиты по положению (отклонение более 100 км) получились следующими:

1) после первого импульса: $0.4\%\pm 0.1\%$ с вероятностью не менее 99.9%,

2) после второго импульса: $0.2\%\pm 0.1\%$ с вероятностью не менее 99.9%,

3) после третьего импульса: $0.05\%\pm 0.1\%$ с вероятностью не менее 99.9%,

4) после четвертого импульса: $0.1\%\pm 0.1\%$ с вероятностью не менее 99.9%.

Вероятность выйти за пределы окрестности орбиты на произвольном шаге равна $0.2\%\pm 0.1\%$ с вероятностью не менее $99.9\%$.

ЗАКЛЮЧЕНИЕ

Сформулирована методика сведения общей задачи оптимального управления космическим аппаратом к задаче машинного обучения с подкреплением. Методика состоит из нескольких шагов: 1) определение состояния и действия в терминах переменных механической задачи; 2) определение распределения начальных фазовых состояний; 3) определение дискретного шага системы; 4) определение функции вознаграждения; 5) определение модели восприятия; 6) определение модели управления. Приведены строгие подходы к оценке математического ожидания случайных величин и вероятности наступления событий на основе неравенства Хефдинга. Рассмотрены два примера применения методики: 1) к простой динамической системе, 2) к задаче поддержания неустойчивой орбиты вокруг точки либрации. В первом примере показана близость получаемых решений к теоретическим значениям. Во втором примере приведены результаты обучения моделей управления, строго оценены средние затраты характеристической скорости и вероятность неудачи поддержания орбиты.

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

² Большое число знаков после запятой приводится для воспроизводимости результатов.

作者简介

М. Широбоков

编辑信件的主要联系方式.
Email: shirobokov@keldysh.ru
俄罗斯联邦, Москва

参考

补充文件

附件文件

动作

1. JATS XML

下载

2. Fig. 1. The strategy learning process.

下载 (29KB)

索引源数据

3. Fig. 2. Optimal (blue) and approximate (red) control functions.

下载 (10KB)

索引源数据

4. Fig. 3. Optimal (blue) and approximate (red) value functions.

下载 (11KB)

索引源数据

5. Fig. 4. Halo orbit (blue), in the vicinity of which the motion of the spacecraft is considered. The projections of the orbit on the xy, xz, and yz planes are shown in gray.

下载 (23KB)

索引源数据

6. Fig. 5. Average reward per episode depending on the step.

下载 (16KB)

索引源数据

7. Figure 6. The standard deviation of the value function model from the average total rewards.

下载 (17KB)

索引源数据