Журнал вычислительной математики и математической физики, 2021, T. 61, № 5, стр. 865-877

TT-QI: ускоренная итерация функции ценности в формате тензорного поезда для задач стохастического оптимального управления

А. И. Бойко^1, *, И. В. Оселедец^1, 2, **, Г. Феррер¹

¹ Сколковский институт науки и технологий
121205 Москва, Большой бульвар, 30, стр. 1, Россия

² ИВМ РАН
119333 Москва, ул. Губкина, 8, Россия

^* E-mail: alexey.boyko@skolkovotech.ru
^** E-mail: i.oseledets@skoltech.ru

Поступила в редакцию 24.11.2020
После доработки 24.11.2020
Принята к публикации 14.01.2021

DOI: 10.31857/S0044466921050045

Полный текст (PDF)

Аннотация

Рассматривается задача стохастического оптимального управления общего вида с малым винеровским шумом. Данная задача аппроксимируется с помощью марковского процесса принятия решений. Решение уравнения Беллмана на функцию ценности вычисляется с помощью метода итерации ценности (VI) в формате малорангового тензорного поезда (ТТ-VI). Предложена модификация данного алгоритма (ТТ-QI): нелинейный оператор Беллмана итеративно применяется сначала с использованием малоранговых алгебраических операций, а затем с использованием алгоритма крестовой аппроксимации. Показана более низкая, чем в основном методе, сложность на одну итерацию в случае малых ТТ-рангов тензоров вероятностей перехода. На примере задач управления обратным маятником и машинами Дубинса показано ускорение времени расчета оптимального регулятора в 3–10 раз по сравнению с существующим методом. Библ. 13. Фиг. 6. Табл. 1.

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

1. ВВЕДЕНИЕ

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

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

Стохастические динамические системы могут быть представлены (с точностью до погрешности дискретизации) как марковский процесс принятия решений (см. [1], [2]). Формулировка задачи оптимального управления на языкe уравнения Беллмана также позволяет нахождение оптимального регулятора в случае произвольной нелинейной, разрывной или точечной награды. Это делает возможным переписать на беллмановский язык такие задачи, как терминальную задачу об оптимальном быстродействии, задачу о минимальном затраченном топливе (с произвольной, а не только квадратичной характеристикой), задачу максимизации вероятности попадания в целевую область.

Ключевой проблемой в таком подходе является тот факт, что решение задачи стохастического оптимального управления на языке Беллмана страдает от “проклятия размерности”, так как имеет крайне высокую асимптотику вычислительной сложности. Сложность по памяти растет как $O({{N}^{{{{d}_{s}} + {{d}_{a}}}}})$, где $N$ – количество узлов дискретизации сетки по одной координате, а ${{d}_{s}}$ и ${{d}_{a}}$ – размерности пространства состояний и пространства управляющих сигналов соответственно. Например, для простого беспилотного летательного аппарата ${{d}_{s}} + {{d}_{a}} = 16$.

Один из методов обойти данную проблему – использовать дифференцируемые функциональные аппроксиматоры, например нейронные сети. Такой подход изначально был развит Бертсекасом под названием нейродинамического программирования (см. [3]), а позже стал известен как (глубокое) машинное обучение с подкреплением. К настоящему моменту действие в этом направлении привело к значительным успехам в решении задач управления существенно нелинейными малоприводными системами, в том числе с неголономными связями (см. [4]). Пожалуй, самым выдающимся примером этого является синтез регулятора для движения бега в подробной многозвенной малоприводной математической модели человека, управляемой более чем 100 мышцами (см. [5]). Такой подход, однако, имеет свои недостатки. Главным образом, это очень высокие требования к количеству данных (симуляций Монте-Карло) и к вычислительным мощностям, а также необходимость вручную подбирать различные параметры оптимизатора и топологию нейросети и отсутствие гарантий сходимости оптимального управления к глобальному оптимуму даже после длительной оптимизации.

Другой способ обойти “проклятие размерности” уравнения Беллмана – это использовать малоранговые тензорные разложения. Эффективность такого подхода для решения задачи оптимального управления была продемонстрирована Н. Хоровицем (см. [6]) из Калифорнийского технологического института с использованием линеаризованного уравнения Гамильтона–Якоби–Беллмана для контрольно-аффинных систем совместно с каноническим тензорным разложением. Несколько иной метод был предложен А. Городетским из Массачуссетского технологического института с использованием нелинейного уравнения Беллмана для систем общего вида и ТТ-разложения (см. [7]) и его дифференцируемого обобщения (см. [8], [9]). В этих статьях были найдены решения различных нелинейных малоприводных задач управления единообразным подходом, в том числе задачи акробатического пилотажа квадрокоптером далеко за областью линейности, используя лишь однопроцессорную многоядерную рабочую станцию. В данной статье мы рассматриваем именно подход, описанный А. Городетским.

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

Рассмотрим стохастическую динамическую полностью наблюдаемую систему, задаваемую уравнением

(1)

$d{{s}_{i}}(t) = {{b}_{i}}(s,a)dt + {{\sigma }_{{ij}}}(s)dw.$

Система в каждый момент времени описывается вектором состояния $s \in \mathcal{S}$ и вектором действия $a \in \mathcal{A}$ (также иногда называемым управляющим сигналом). В качестве пространства состояний системы $\mathcal{S}$ и пространства управляющих сигналов $\mathcal{A}$ рассматриваются ${{d}_{s}}$ и ${{d}_{a}}$ – мерные области, порожденные произведением замкнутых отрезков

$\begin{gathered} s \in \mathcal{S} = [s_{{min}}^{{(1)}},s_{{max}}^{{(1)}}] \times \; \ldots \; \times [s_{{min}}^{{({{d}_{s}})}},s_{{max}}^{{({{d}_{s}})}}], \\ a \in \mathcal{A} = [a_{{min}}^{{(1)}},a_{{max}}^{{(1)}}] \times \; \ldots \; \times [a_{{min}}^{{({{d}_{a}})}},a_{{max}}^{{({{d}_{a}})}}], \\ \end{gathered} $

$dw$ – стандартный броуновский шум (винеровский процесс). Даны функции дрейфа $b(s,a):\mathcal{S} \times \mathcal{A} \to {{\mathbb{R}}^{{{{d}_{s}}}}}$ и диффузии ${{\sigma }_{{ij}}}(s):\mathcal{S} \to {{\mathbb{R}}^{{d_{s}^{2}}}}$. Для постановки задачи оптимального управления также зададимся моментальной функцией награды $r(s,a)$, зависящей от текущего состояния системы $s$ и текущего управляющего сигнала $a$, а также коэффициентом дисконтирования $\beta $. Функции $b(s,a)$, ${{\sigma }_{{ij}}}(s)$ полагаются ограниченными на своих областях определений.

На функцию $r(s,a)$ накладывается следующее ограничение:

(2)

$\left| {r(s,a)} \right| \leqslant C(1 + {{\left| s \right|}^{k}} + {{\left| a \right|}^{k}}),$

где $C \in {{\mathbb{R}}_{ + }}$ и $k \in \mathbb{N}$ – некоторые константы.

Потребуем выполнения еще одного условия: матрично-значная функция диффузии ${{\sigma }_{{ij}}}(s)$ является диагональной (${{\sigma }_{{ij}}} = 0$ если $i \ne j$). Такая ситуация имеет место, например, в случае присутствия нескоррелированного шума в датчиках, измеряющих текущее состояние системы $s$.

Ставится следующая задача: найти оптимальную детерминистическую политику (в зависимости от области науки также называемую оптимальным управлением, регулятором или стратегией) $\pi {\text{*}}(s):\mathcal{S} \to \mathcal{A}$. Оптимальность здесь подразумевается в следующем смысле: если данная политика используется для генерации сигнала управления в каждый момент времени $a(t) = \pi {\text{*}}\left( {s(t)} \right)$, то матожидание суммарной награды (в общем случае c дисконтированием) по реализациям всех возможных случайных траекторий за все доступное время будет максимальным:

(3)

$\pi {\kern 1pt} {\text{*}}(s) = \mathop {\arg \max }\limits_a \left[ {\mathbb{E}\int\limits_{t = 0}^T {\exp } ( - \beta t)r(s(t),a)dt} \right],$

где $\beta $ – коэффициент дисконтирования в непрерывном времени, $T$ – время конца эпизода (для дисконтированных нетерминальных задач управления $T = \infty $). Необходимым образом требуется, чтобы матожидание награды было ограничено: $\mathbb{E}\int_{t = 0}^T {{{e}^{{ - \beta t}}}} r(s(t),a)dt < \infty $.

3. АППРОКСИМАЦИЯ МАРКОВСКИМ ПРОЦЕССОМ ПРИНЯТИЯ РЕШЕНИЙ

Говорят, что задан марковский процесс принятия решений (МППР) в дискретном времени, если задан следующий кортеж:

(4)

$(\mathcal{S},\mathcal{A},T,R,\gamma ),$

где заданы конечное множество допустимых значений состояний $\mathcal{S}$, конечное множество допустимых значений действий (управляющих сигналов, управляющих решений) $\mathcal{A}$, известны все (условные) вероятности перехода из любого состояния $s$ в любое состояние $s{\text{'}}$ через любое действие $a$: $T(s,a,s{\text{'}}) = P(s{\text{'}}\,|\,s,a)$, а также функция награды для любого такого перехода $R(s,a,s{\text{'}})$ и коэффициент дисконтирования $\gamma \in (0,\;1]$. Если $\gamma = 1$, то МППР называется недисконтированным.

Для численного нахождения решения исходной задачи стохастического оптимального управления (1) дискретизуем ее. Общая теория дискретизации для таких задач может быть найдена в книгах В. Флеминга (см. [2]) и Г. Кушнера (см. [1]). Основная идея этой теории в том, чтобы спроектировать непрерывные пространства состояний и действий на сетку, а затем построить МППР, который был бы эквивалентен нашему стохастическому процессу управления (1) в смысле дрифта и диффузии за единицу времени:

$\mathop {lim}\limits_{h \to 0} \frac{{\mathbb{E}[{{s}_{{t + 1}}} - {{s}_{t}}\,|\,{{s}_{t}} = s,\;{{a}_{t}} = a]}}{{\Delta t}} = b(s,a),$

$\mathop {lim}\limits_{h \to 0} \frac{{{\text{Cov}}[{{s}_{{t + 1}}} - {{s}_{t}}\,|\,{{s}_{t}} = s,\;{{a}_{t}} = a]}}{{\Delta t}} = \sigma (s,a).$

Здесь $h = min({{h}_{i}})$, а ${{h}_{i}}$ – шаги дискретизации равномерной сетки вдоль $i$-й оси: ${{h}_{i}} = (s_{{{\text{max}}}}^{{(i)}} - s_{{{\text{min}}}}^{{(i)}}){\text{/(}}N_{s}^{{(i)}} - 1{\text{)}}$.

Мы будем использовать модифицированную версию схемы расщепления против потока, предложенной в [9], построенной по общей методике Г. Кушнера. Данная схема разрешает переходы с ненулевой вероятностью только между соседними узлами сетки, и при этом все диагональные переходы запрещены. С учетом того, что диффузионный член ${{\sigma }_{{ij}}}(s)$ диагонален, схема дискретизации принимает следующий вид:

$\begin{gathered} {{Q}^{h}} = \mathop {\max }\limits_{s,a} \left( {\sum\limits_i {\frac{{\left| {{{b}_{i}}(s,a)} \right|}}{{{{h}_{i}}}}} + \frac{{\sigma _{i}^{2}(s)}}{{h_{i}^{2}}}} \right), \\ \Delta {{t}^{h}} = 1{\text{/}}{{Q}^{h}}, \\ \end{gathered} $

${{b}^{ + }} = \left\{ \begin{gathered} b(s,a),\quad {\text{если}}\quad b(s,a) > 0, \hfill \\ 0\quad {\text{иначе,}} \hfill \\ \end{gathered} \right.$

(5)

${{b}^{ - }} = \left\{ \begin{gathered} - b(s,a),\quad {\text{если}}\quad b(s,a) < 0, \hfill \\ 0\quad {\text{иначе,}} \hfill \\ \end{gathered} \right.$

$\begin{gathered} T(s,a,s \pm {{e}_{i}}{{h}_{i}}) = \Delta {{t}^{h}}\left( {\frac{{b_{i}^{ \pm }(s,a)}}{{{{h}_{i}}}} + \frac{{\sigma _{i}^{2}(s)}}{{2h_{i}^{2}}}} \right), \\ T(s,a,s) = 1 - \sum\limits_{s'} {T(s,a,s{\kern 1pt} ' \ne s)} , \\ \end{gathered} $

$\begin{gathered} R(s,a,s{\kern 1pt} ') = r(s,a)\Delta {{t}^{h}}, \\ \gamma = {{e}^{{ - \beta \Delta {{t}^{h}}}}}, \\ \end{gathered} $

где $s \pm {{e}_{i}}{{h}_{i}}$ обозначает дискретное состояние (узел сетки), которое является соседним по отношению к узлу s вдоль i-й оси, в направлении увеличения или уменьшения индекса соответственно.

4. МАЛОРАНГОВОЕ РАЗЛОЖЕНИЕ В ТЕНЗОРНЫЙ ПОЕЗД

Рассмотрим $d$-мерный массив (тензор) $V({{i}_{1}},\; \ldots ,\;{{i}_{d}})$. Допустим, мы хотим найти его представление в следующем виде:

(6)

$V({{i}_{1}},{{i}_{2}},\; \ldots ,\;{{i}_{d}}) = \sum\limits_{{{m}_{1}}, \ldots ,{{m}_{{d - 1}}}}^{{{r}_{1}},{{r}_{2}}, \ldots ,{{r}_{{d - 1}}}} {{{G}^{{(1)}}}} ({{i}_{1}},{{m}_{1}}){{G}^{{(2)}}}({{m}_{2}},{{i}_{2}},{{m}_{2}}) \cdot \; \ldots \; \cdot {{G}^{{(d)}}}({{m}_{{d - 1}}},{{i}_{d}}).$

Такое разложение называется разложением в тензорный поезд или TT-разложением. Числа ${{r}_{1}},\; \ldots ,\;{{r}_{{d - 1}}}$ называются TT-рангами. Чтобы равенство (6) выполнялось точно, необходимо хранить в худшем случае $O(d{{N}^{3}})$ чисел из вместо $O({{N}^{d}})$. Зачастую имеет смысл найти такой тензорный поезд, что данное равенство соблюдается с некоторой погрешностью (в смысле среднеквадратичной ошибки), но TT-ранги малы ($\forall i:{{r}_{i}} < 100$). В таком случае можно ожидать, что функция на сетке будет сжата с требованием по памяти, равным $O(dN{{R}^{2}})$, где $R = max({{r}_{i}})$.

Важным фактом является то, что для достаточно широкого класса функций, спроецированных на равномерную сетку (более полное описание дано в [10]), если мы зафиксируем наибольшую допустимую ошибку аппроксимации, при измельчении сетки ТТ-ранги не будут расти, либо будут расти лишь логарифмически. Отметим, что в дополнение к сжатию в тензорных поездах сохраняется возможность выполнять алгебраические ($ + $, $ - $, $ * $, $ \otimes $, $ \circ $) и линейно-алгебраические операции (внешние и внутренние произведения, свертки, суммирования по индексам), операции взятия подмассива и доступа к индивидуальным элементам, а также применять произвольные функции к тензорным поездам с помощью алгоритмов из семейства многомерной крестовой аппроксимации (см. [11]).

5. УРАВНЕНИЕ БЕЛЛМАНА

Метод нахождения глобально-оптимальной политики решений для марковского процесса принятия решений был предложен Р. Беллманом в свом классическом труде [12]. Рассмотрим его основные положения.

Зададимся МППР ($\mathcal{S}$, $\mathcal{A}$, $T$, $R$, $\gamma $). Задачей является найти такую управляющую политику, которая, стартуя из любого состояния $s$, максимизирует матожидание (возможно, дисконтированной) награды за все последующее время:

$V(s) = \mathbb{E}{\kern 1pt} \sum\limits_{t = 0}^\infty {{{\gamma }^{t}}R({{s}_{t}},{{a}_{t}},{{s}_{{t + 1}}})} \to max.$

Величина $V(s)$ называется функцией ценности (value function) и является ключевым элементом исследований в стохастическом оптимальном управлении и науке о машинном обучении с подкреплением.

Один из главных результатов Беллмана заключался в следующем: глобально оптимальная политика (управление, регулятор) для МППР не зависит от предыстории, а зависит только от текущего состояния. Из этого следует, что для глобальной оптимальности политики достаточно, чтобы действие было локально оптимально, но не в смысле текущего приращения функции награды $R(s,a)$, а в смысле текущего приращения функции ценности $V(s)$.

Рекуррентные соотношения значений функции ценности в соседних состояниях, разделенных всего одним временным шагом, называются уравнениями Беллмана и имеют вид

(7)

$V{\text{*}}(s) = \mathop {max}\limits_a \left[ {R(s,a) + \gamma \sum\limits_{s{\text{'}}} {V{\text{*}}} (s{\text{'}})T(s,a,s{\text{'}})} \right].$

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

$R(s,a) = {{\mathbb{E}}_{{s{\text{'}}}}}R(s,a,s{\text{'}}) = \sum\limits_{s{\text{'}}} R (s,a,s{\text{'}})T(s,a,s{\text{'}}).$

Если рассмотреть правую часть уравнения (7) как оператор $\hat {\beta }$, действующий на функцию $V(s)$, то такой оператор в литературе называется оператором оптимальности Беллмана (Bellman Optimality Operator).

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

(8)

$\pi {\text{*}}(s) = {\text{arg}}\mathop {max}\limits_a \left[ {R(s,a) + \gamma \sum\limits_{s{\text{'}}} {V{\text{*}}} (s{\text{'}})T(s,a,s{\text{'}})} \right].$

Заметим, что ряд задач оптимального управления (называемых терминальными), например, задача оптимального быстродействия, требуют задания терминальной области. При попадании в терминальную область счетчик времени останавливается, и дальшейшее накопление награды невозможно. Чтобы математически ввести терминальную область ${{\mathcal{S}}_{{{\text{terminal}}}}} \subset \mathcal{S}$ на языке Беллмана, необходимо положить нулем все вероятности переходов для состояний из данной области, кроме переходов каждого состояния само в себя: $\forall (s,a) \in {{\mathcal{S}}_{{{\text{terminal}}}}} \times \mathcal{A}:T(s,a,s{\text{'}}) = {{\delta }_{{ss{\text{'}}}}}$, где ${{\delta }_{{ss{\text{'}}}}}$ – символ Кронекера. Также функция награды должна быть положена нулем для всех переходов, начинающихся из терминальной области: $\forall s \in {{\mathcal{S}}_{{{\text{terminal}}}}}:R(s,a) = 0$, $V(s) = 0$.

5.1. Итерация функции ценности в формате тензорного поезда

Уравнение Беллмана является уравнением на неподвижную точку для нелинейного оператора оптимальности Беллмана $\hat {\beta }$, и оно может быть решено методом простых итераций. Такой подход был предложен самим Р. Беллманом в его классической работе [12] и называется методом итерации функции ценности (Value Iteration, VI). Он использует то обстоятельство, что оператор $\hat {\beta }$ является сжимающим отображением почти во всех практически применимых случаях (в случае дисконтированных задач, либо в случах недисконтированных задач с терминальной областью, достижимой хотя бы одной политикой).

Для обоснования использования малорангового разложения рассмотрим количество занимаемой памяти для МППР, аппроксимирующего задачу стохастического оптимального управления (1). Поскольку используется равномерная сетка дискретизации вдоль каждой оси пространства состояний в ${{N}_{s}}$ узлов, результирующая сложность по памяти для хранения оптимальной функции ценности $V{\text{*}}(s)$ и оптимальной политики $\pi {\text{*}}(s)$ будет составлять $O(N_{s}^{{{{d}_{s}}}})$ для каждой из этих функций. Даже простые мобильные робототехнические аппараты, такие как БПЛА самолетного типа или квадрокоптеры, задаются ${{d}_{s}} = 12$-мерным вектором состояния и минимум ${{d}_{a}} = 3$-мерным вектором управляющих сигналов. В итоге при достаточно грубой дискретизации сетки в 100 узлов по каждой координате в задаче возникает минимум $N_{s}^{{{{d}_{s}}}} = {{10}^{{24}}}$ элементов, что делает прямое численное нахождение решения в беллмановском формализме невозможным. Если при поиске решения вероятности переходов $T(s,a,s{\text{'}})$ не вычисляются каждый раз, а хранятся, то дополнительно потребуется сохранить еще $(2{{d}_{s}} + 1)N_{s}^{{{{d}_{s}}}}N_{a}^{{{{d}_{a}}}}$ чисел (где ${{N}_{a}}$ – число узлов дискретизации вдоль осей пространства управляющих сигналов, а ${{d}_{a}}$ – размерность этого пространства).

В статье [7], предложенной в Массачуссетском технологическом институте А. Городетским, предлагалось хранить функцию ценности в виде малорангового тензорного поезда и применять оператор оптимальности Беллмана в рамках итерации функции ценности методом крестовой аппроксимации для тензорных поездов (см. [11]). В случае малых ТТ-рангов функции ценности ее сложность по памяти понижается с $O({{N}^{d}})$ до $O(dN{{r}^{2}})$, что делает задачу вычислительно решаемой даже на маломощном компьютере. Также в [7], [9] показано, что в ТТ-формате оператор $\hat {\beta }$ сохраняет свойство быть сжимающим отображением, и гарантии сходимости, выведенные для изначальной итерации функции ценности (см. [12]), сохраняются:

Data:$R(s,a),\gamma ,T(s,a,s{\text{'}})$
Result:$V{\text{*}}(s)$
while$\epsilon = \tfrac{{{{{\left\\| {{{V}^{{(k + 1)}}} - {{V}^{{(k)}}}} \right\\|}}_{2}}}}{{{{{\left\\| {{{V}^{{(k)}}}} \right\\|}}_{2}}}} < tol$do
	${{V}^{{k + 1}}} = TTCROSS((7),\epsilon )$
	$k = k + 1$
end

Алгоритм 1: TT-Value Iteration (TT-VI)

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

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

$\begin{gathered} {{R}_{V}} = max{\text{rank}}\,V(s,a), \\ {{R}_{P}} = \mathop {max}\limits_{i = 0,1, \ldots {{d}_{s}}, \pm } {\text{rank}}\,P_{i}^{ \pm }(s,a). \\ \end{gathered} $

В наивной имплементации TT-VI алгоритма крестовая аппроксимация требует $O(dNR_{V}^{3})$ на каждое применение оператора оптимальности Беллмана, а также требуется $O(dR_{V}^{2})$ операций для распаковывания значений из ТТ в каждой точке, что в итоге приводит к вычислительной сложности $O({{d}^{2}}NR_{V}^{5})$.

Однако вычисления в алгоритме TT-VI возможно оптимизировать. Для этого вспомним, что метод крестовой аппроксимации считывает значения функции в точках, которые сгруппированы в виде строк, столбцов и их многомерных обобщений (распорок) многомерной сетки. Оценим сложность такого алгоритма. Алгоритм крестовой аппроксимации требует $O(dR_{V}^{3})$ распорок, а извлечение одной распорки из ТТ-разложения требует $O(dR_{V}^{2})$ операций (что меньше, чем извлечение одного числа). Результирующая сложность алгоритма составляет $O({{d}^{2}}R_{V}^{5})$ операций. Далее в статье мы будем использовать именно эту оптимизированную версию как базовый алгоритм для сравнения.

5.2. Q-итерация в формате тензорного поезда

В данной статье мы предлагаем модификацию алгоритма TT-VI. Ее суть заключается в следующем: использовать стандарную итерацию функции ценности (алгоритм 1), но применять оператор оптимальности Беллмана в два этапа:

(9)

${{Q}^{k}}(s,a) = R(s,a) + \gamma \sum\limits_{s{\text{'}}} T (s,a,s{\text{'}}){{V}^{k}}(s{\text{'}}),$

(10)

${{V}^{{k + 1}}} = \mathop {max}\limits_a {{Q}^{k}}(s,a).$

Первый этап (9) включает в себя только стандартные тензорные операции, такие как сложения и суммирования по индексам, тогда как второй (10) требует взятия нелинейной функции (максимума), для чего мы используем крестовый метод, также как в TT-VI алгоритме. Так как у двух алгоритмов совпадают вторые этапы, мы опустим оценку сложности для них в сравнительном анализе.

Рассмотрим первый этап применения оператора $\hat {\beta }$, в частности суммирование по индексам с тензором $T$. Как следует из используемой нами дискретизации (5), из каждого состояния $s$ возможны только переходы в соседние состояния вдоль каждой оси $s{\text{'}} = s \pm {{e}_{i}}{{h}_{i}}$. Это означает, что сумма вероятностей по $s{\kern 1pt} '$ состоит из $2{{d}_{s}} + 1$ слагаемых: два на каждую координатную ось (вероятность перехода в направлении увеличения или уменьшения на 1) и еще одно на вероятность остаться на месте $s$:

$\sum\limits_{s{\text{'}}} T (s,a,s{\text{'}})V(s{\text{'}}) = {{P}_{0}}(s,a)V(s) + \sum\limits_{{\text{sign}} \in {\text{\{ }} + , - {\text{\} }}} {\sum\limits_{i = 1}^{{{d}_{s}}} {P_{i}^{{{\text{sign}}}}} } (s,a)V_{i}^{{ - {\text{sign}}}}(s),$

где $P_{i}^{{{\text{sign}}}}$ и ${{P}_{0}}$ – вероятности, посчитанные с помощью схемы дискретизации (5), а $V_{i}^{{ - {\text{sign}}}}(s) = {\text{circshift}}(V({{s}_{1}},\; \ldots ,\;{{s}_{{{{d}_{s}}}}}), - {\text{sign}},i)$ – функция ценности (на сетке), циклически сдвинутая вдоль $i$-й оси на 1 узел, в направлении $ - {\text{sign}}$.

Предложение 1. Циклический сдвиг по $k$-й координате функции на сетке, представленной в виде тензорного поезда, достигается циклическим сдвигом всего одного ($k$-го) тензорного ядра ${{G}^{{(k)}}}({{m}_{k}},{{i}_{k}},{{m}_{{k + 1}}})$. Он стоит $O(N{{R}^{2}})$ операций и не меняет тензорный ранг.

Доказательство. Из формулы ТТ-разложения видно, что при циклическом сдвиге любого из свободных индексов ${{i}_{k}} \to {{i}_{k}} \pm 1mod{{N}_{k}}$ равенство сохраняется:

$V({{i}_{1}},\; \ldots ,\;{{i}_{k}} \pm 1,\; \ldots ,\;{{i}_{d}}) = \sum\limits_{{{m}_{1}}, \ldots ,{{m}_{{d - 1}}}}^{{{r}_{1}}, \ldots ,{{r}_{{d - 1}}}} {{{G}^{{(1)}}}} ({{i}_{1}},{{m}_{1}}) \cdot \; \ldots \; \cdot {{G}^{{(k)}}}({{m}_{k}},{{i}_{k}} \pm 1mod{{N}_{k}},{{m}_{k}}) \cdot \; \ldots \; \cdot {{G}^{{(d)}}}({{m}_{{d - 1}}},{{i}_{d}}).$

Замена ${{G}^{{(k)}}}({{m}_{2}},{{i}_{k}},{{m}_{2}}) \to {{G}^{{(k)}}}({{m}_{2}},{{i}_{k}} \pm 1mod{{N}_{k}},{{m}_{2}})$ требует перемещения всех элементов массива ${{G}^{{(k)}}}$, которых насчитывается ${{r}_{k}}{{N}_{k}}{{r}_{{k + 1}}}$ для всех ядер кроме первого и последнего. Для первого и последнего ядра потребуется ${{N}_{1}}{{r}_{1}}$ и ${{N}_{d}}{{r}_{{d - 1}}}$ элементов соответственно, что в общем случае оценивается как $O(N{{R}^{2}})$. Очевидно, что такая операция не меняет размер ядра ${{G}^{{(k)}}}$, поэтому ранг ${{r}_{k}}$ сохраняется.

Чтобы вычислить поэлементные произведения $P(s,a)V(s)$ в формате тензорного поезда, необходимо добавить в функцию ценности “лишние” измерения, добавив новые (единичные) ядра: $P({{s}_{1}}\; \ldots \;{{s}_{{{{d}_{s}}}}},{{a}_{1}}\; \ldots \;{{a}_{{{{d}_{a}}}}}) \circ ((V(s)) \otimes {{1}_{{{{a}_{1}}}}}\; \ldots \; \otimes {{1}_{{{{a}_{{{{d}_{a}}}}}}}}).$

Предложение 2. Cложность первой половины (9) применения оператора $\hat {\beta }$ в TT-QI составляет

$O({{d}^{2}}NR_{V}^{3}R_{P}^{3}).$

Доказательство. Одно применение оператора оптимальности Беллмана в алгоритме TT-QI (до этапа применения алгоритма крестовой аппроксимации) использует следующие тензорные операции:

• $2{{d}_{s}}$ циклических сдвигов $V(s)$, каждый по $O(NR_{V}^{2})$, что в итоге составляет $O(dNR_{V}^{2})$,

• 1 акт добавления в $V(s)$ единичных ядер с лишними измерениями, что составляет $O({{d}_{a}})$,

• взятие $2{{d}_{s}} + 1$ поэлементных произведений от тензорных поездов, сложностью $O(dNR_{P}^{2}R_{V}^{2})$ каждый, что дает $O({{d}^{2}}NR_{P}^{2}R_{V}^{2})$ в сумме,

• взятие $2{{d}_{s}} + 1$ TT-SVD округлений, что дает $O(dNR_{P}^{3}R_{V}^{3})$ каждый, и $O({{d}^{2}}NR_{P}^{3}R_{V}^{3})$ в сумме.

Полная сложность первого этапа составляет в итоге

(11)

$O(dNR_{V}^{2} + {{d}_{a}} + {{d}^{2}}NR_{P}^{2}R_{V}^{2} + {{d}^{2}}NR_{P}^{3}R_{V}^{3}) = O({{d}^{2}}NR_{P}^{3}R_{V}^{3}).$

В случае существенно малых рангов ${{R}_{P}}$ тензоров вероятностей эта асимптотика более выгодна относительно $O({{d}^{2}}R_{V}^{5})$ в алгоритме TT-VI. Тензоры вероятностей ${{P}_{{0 {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} {{d}_{s}}}}}(s,a)$ вычисляются только один раз и их ранги ${{R}_{P}}$ постоянны во время итерационного процесса, что в итоге дает асимптотическую сложность $O(R_{V}^{3})$ в случае TT-QI вместо $O(R_{V}^{5})$ в случае TT-VI. В наших численных экспериментах это привело к существенному (в 5–6 раз) выигрышу в производительности.

6. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

Чтобы показать верность оптимальных политик (регуляторов), полученных на базе уравнения Беллмана, следует рассмотреть достаточно сложную систему, для которой решение не может быть получено более простыми методами из теории управления (такие как PID или линейно-квадратичные регуляторы).

Данный раздел содержит графики сходимости к решению уравнения Беллмана, а также результаты прямых симуляций траекторий, полученных под управлением оптимальной политики. Для интегрирования стохастических дифференциальных уравнений по Ито и расчета траекторий мы используем схему высокого порядка точности (1.0), предложенную Росслером (см. [13]).

6.1. Система 1: непериодический обратный маятник

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

(12)

$s = {{[\phi ,\dot {\phi }]}^{{\text{т}}}},$

(13)

$b(s,a) = {{\left[ {\dot {\phi },a - sin(\phi )} \right]}^{{\text{т}}}},$

(14)

$\sigma (s) = {\text{diag}}\left( {[{{{10}}^{{ - 2}}},{{{10}}^{{ - 2}}}]} \right).$

Первая модификация – существенно малая мощность управляющего момента на валу маятника (максимальный момент сил на вал, вызванный управляющим сигналом, равен лишь $30\% $ от максимального момента сил тяжести). Такая постановка задачи делает невозможным для регулятора достижение верхней точки любым способом, кроме как с помощью эксплуатации резонанса системы. Вторая модификация – это использование непериодического маятника, что также делает задачу управления сложнее.

Границы областей состояний и действий определены как

(15)

$s \in [ - 3\pi ,3\pi ] \times [ - \pi ,\pi ],$

(16)

$a \in [ - 0.3,0.3].$

Мы рассмотрим решения двух вариантов данной задачи: с квадратичной наградой и задачу оптимального быстродействия.

6.1.1. Задача управления с квадратичной наградой. Задачи управления с квадратичной наградой/ценой хорошо изучены в литературе, поэтому в данной статье мы будем использовать такую постановку как тестовую для сравнения производительности (см. табл. 1).

Таблица 1.

Сравнение производительности

Параметр	TT-VI (оптимизизирован-ный) : время, с	TT-QI : время, с
3-маятник, квадратичная награда, целевая ошибка $\delta = {{10}^{{ - 4}}}$	1959.9	536.7
3-маятник, задача оптимального быстродействия, целевая ошибка $\delta = 5 \times {{10}^{{ - 4}}}$	6024.2	2116.8
3-машина, квадратичная награда, целевая ошибка $\delta = 2 \times {{10}^{{ - 3}}}$	2952.2	293.4
4-машина, квадратичная награда, целевая ошибка $\delta = 5 \times {{10}^{{ - 3}}}$	3061.8	1109.3

Квадратичная награда задается следующим образом: ${{R}_{{QR}}}(s,a) = - {{(\phi - \pi )}^{2}} - 0.8{{(\phi )}^{2}} - 0.01{{a}^{2}}$ и коэффициент дисконтирования $\gamma = 0.999$.

6.1.2. Задача оптимального быстродействия. Для того чтобы сформулировать задачу оптимального быстродействия на беллмановском языке, необходимо положить награду для каждого перехода $(s\mathop \to \limits^a as{\text{'}})$ равной отрезку времени, которое требуется, чтобы этот переход совершить:

(17)

$R(s,a,s{\text{'}}) = \left\{ \begin{gathered} - \Delta t(s,a,s{\text{'}}),\quad {\text{если}}\quad s \notin {{\mathcal{S}}_{{{\text{terminal}}}}}, \hfill \\ 0,\quad {\text{если}}\quad s \in {{\mathcal{S}}_{{{\text{terminal}}}}}. \hfill \\ \end{gathered} \right.$

В таком случае функция ценности становится равна матожиданию времени прибывания в терминальную область:

(18)

$V(s) = \sum\limits_{t = 0}^\infty \mathbb{E} {\kern 1pt} R(s,a,s{\text{'}}) = \sum\limits_{t = 0}^\infty \mathbb{E} {\kern 1pt} \Delta t(s,a,s') = - \sum\limits_{t = 0}^\tau \mathbb{E} {\kern 1pt} \Delta t = - \tau .$

Терминальную область в данном случае мы полагаем малой окрестностью целевого состояния $(\pi ,0)$, что соответствует перевернутому маятнику с нулевой угловой скоростью:

(19)

${{\mathcal{S}}_{{{\text{terminal}}}}} = [\pi - 2{{h}_{0}},\pi + 2{{h}_{0}}] \times [ - 2{{h}_{1}},2{{h}_{1}}].$

Коэффициент дисконтирования положен равным единице ($\gamma = 1$), т.е. решается недисконтированная задача.

Решение данной задачи на сетке с дискретизацией $301 \times 151 \times 51$ дает оптимальную функцию ценности, оптимальный регулятор (см. фиг. 1) и соответствующим образом изменяет фазовый портрет динамической системы под действием регулятора, как показано на фиг. 2.

Фиг. 1.

Решение для задачи оптимального быстродействия непериодическим маятником: (а) – оптимальная функция ценности $V(s)$, (б) – оптимальный регулятор $\pi {\kern 1pt} {\text{*}}(s)$.

Фиг. 2.

Фазовые портреты динамической системы непериодического маятника: (а) – без управления ($a \equiv 0$), (б) – под управлением оптимального регулятора ($a = \pi {\kern 1pt} {\text{*}}(s)$).

6.2. Система 2: машина Дубинса

Другим известным примером малоприводной системы является машина (автомобиль) Дубинса.

6.2.1. Простая машина Дубинса. Рассмотрим простейшую модель колесного автомобиля, описываемого уравнениями

(20)

$s = {{[x,y,\phi ]}^{{\text{т}}}},$

(21)

$a = {{[u,\theta ]}^{{\text{т}}}},$

(22)

$b(s,a) = {{\left[ {ucos\phi ,\;usin\phi ,\;\frac{u}{L}{\text{tg}}\,\theta } \right]}^{{\text{т}}}},$

(23)

$\sigma (s) = \operatorname{diag} ({{10}^{{ - 3}}},{{10}^{{ - 3}}},{{10}^{{ - 3}}}).$

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

$\begin{gathered} \mathcal{S} = [ - 100,\;100] \times [ - 100,\;100] \times SO(2), \\ \mathcal{A} = [0,\;1] \times \left[ { - \frac{\pi }{3}, - \frac{\pi }{3}} \right], \\ \end{gathered} $

где $SO(2)$ – одномерная окружность. Функция награды задана следующим образом:

$r(s,a) = - {{10}^{{ - 4}}}{{x}^{2}} - {{10}^{{ - 4}}}{{y}^{2}} - {{10}^{{ - 5}}}{{\phi }^{2}} - {{10}^{{ - 3}}}{{u}^{2}}.$

Результаты тестов сходимости показаны ниже на фиг. 6а.

6.2.2. Машина Дубинса c инерцией. Данная модель машины является усложненной версией модели из п. 6.2.1. В ней сигнал управления не влияет на скорость напрямую, а лишь управляет ускорением. В этой модели автомобиль так же может ехать только вперед, а ускоряться может и вперед, и назад:

(24)

$s = {{[x,y,{v},\phi ]}^{{\text{т}}}},$

(25)

$a = {{[u,\theta ]}^{{\text{т}}}},$

(26)

$b(s,a) = {{\left[ {{v}cos\phi ,\;{v}sin\phi ,\;u,\;\frac{u}{L}{\text{tg}}\,\theta } \right]}^{{\text{т}}}},$

(27)

$\sigma (s) = \operatorname{diag} ({{10}^{{ - 3}}},{{10}^{{ - 3}}},{{10}^{{ - 3}}},{{10}^{{ - 3}}}),$

$\begin{gathered} \mathcal{S} = [ - 70,70] \times [ - 70,70] \times SO(2), \\ \mathcal{A} = [ - 2,1] \times \left[ { - \frac{\pi }{3}, - \frac{\pi }{3}} \right]. \\ \end{gathered} $

Мы рассмотрим два варианта награды для этой системы: $(A)$ – для задачи оптимального быстродействия ($r = {{r}_{A}}$, $\gamma = 1$), $(B)$ – для задачи с квадратичной наградой (${{r}_{B}}$, $\gamma = 0.999$):

(28)

${{r}_{A}}(s,a) = - \Delta t(s,a),$

(29)

${{r}_{B}}(s,a) = - {{10}^{{ - 4}}}{{x}^{2}} - {{10}^{{ - 4}}}{{y}^{2}} - {{10}^{{ - 5}}}{{\phi }^{2}} - {{10}^{{ - 3}}}{{u}^{2}}.$

На фиг. 3 видно несколько симулированных траекторий координат $(x,y)$ для задачи оптимального быстродействия машины Дубинса с инерцией (вариант $A$). Все траектории прибыли в терминальную область. Также видно, что траектории близко соответствуют теоретически предсказаннной форме кривых (путям Дубинса).

Фиг. 3.

Симулированные траектории $\left( {x,y} \right)$ машины Дубинса с инерцией.

7. СРАВНЕНИЕ ПРОИЗВОДИТЕЛЬНОСТИ

В данном разделе мы сравниваем производительность нашего алгоритма (TT-QI) с оптимизированной версией алгоритма TT-VI. На фиг. 4 изображена зависимость максимального ранга функций ценности ${{R}_{V}}$ от целевой ошибки сжатия $\varepsilon $ с помощью TT-SVD. Видно, что функции ценности для данных задач действительно с высокой точностью являются малоранговыми.

Фиг. 4.

Сжимаемость функции ценности в разных задачах.

Для экспериментов, приведенных ниже, функция ценности инициализировалась нулевыми малоранговыми тензорами. Вычисления проводились на рабочей станции с 3.5 ГГц 8-ядерным процессором производства AMD и 12 Гб оперативной памяти. На фиг. 5 и 6 видно, что оба алгоритма ведут себя практически идентично и имеют степенной закон сходимости, но отличаются по затраченному времени на константный множитель.

Фиг. 5.

График относительной ошибки решения уравнения Беллмана для алгоритмов TT-VI и TT-QI: (а) – 3-маятник (квадратичная награда), (б) – 3-маятник (задача оптимального быстродействия).

Фиг. 6.

График относительной ошибки решения уравнения Беллмана для алгоритмов TT-VI и TT-QI: (а) – простая машина Дубинса (квадратичная награда), (б) – машина Дубинса с инерцией (квадратичная награда).

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

На фиг. 5 и 6 построена относительная ошибка $\delta $, которая для $k$-й итерации определяется как

(30)

${{\delta }_{k}} = \frac{{{{{\left\| {{{V}_{k}} - {{V}_{{k - 1}}}} \right\|}}_{2}}}}{{{{{\left\| {{{V}_{{k - 1}}}} \right\|}}_{2}}}}.$

Сравнение времени решения с относительной ошибкой $\delta $ показано в табл. 1.

8. ЗАКЛЮЧЕНИЕ

Мы предложили модифицированный алгоритм для итерации функции ценности для задач стохастического оптимального управления в формате малорангового тензорного поезда. Если ТТ-ранги функции вероятностей перехода (после применения схемы расщепления) малы, а также ТТ-ранги функции награды малы, наш алгоритм имеет меньшую вычислительную сложность, чем существующий метод решения стохастического оптимального управления для общего случая неаффинных систем, существующий в литературе (см. [7]).

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

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

Авторы благодарны С. Долгову (университет Бата), С. Матвееву (ИВМ РАН, Сколтех) и Г. Овчинникову (Сколтех) за ценные обсуждения.

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

Kushner H., Dupuis P.G. Numerical methods for stochastic control problems in continuous time. Springer, 2013, V. 24.
Fleming W.H., Soner H.M. Controlled Markov Processes and Viscosity Solutions. Springer, 2006.
Bertsekas D.P., Tsitsiklis J.N. Neuro-Dynamic Programming, 1st ed. Athena Scientific, 1996.
Lillicrap T.P. et al. Continuous control with deep reinforcement learning // 4th Inter. Conf. Learn. Represent. ICLR, 2016.
Kidzinski et al. Learning to run challenge solutions: Adapting reinforcement learning methods for neuromusculoskeletal environments // Proceed. NIPS, 2017.
Horowitz M., Damle A., Burdick J.W. Linear Hamilton-Jacobi-Bellman equations in high dimensions // IEEE Conf. Decis. Control, 2014.
Gorodetsky A.A., Karaman S., Marzouk Y.M. Efficient high-dimensional stochastic optimal motion control using Tensor Train decomposition // Robotics: Sci. Syst. 2015.
Gorodetsky A.A., Karaman S., Marzouk Y.M. High-dimensional stochastic optimal control using continuous tensor decompositions // Inter. J. Robot. Res. 2018. № 37. Iss. 2–3.
Tal E., Gorodetsky A., Karaman A. Continuous Tensor Train-based dynamic programming for high-dimensional zero-sum differential games // Am. Control Conf. 2018.
Oseledets I.V., Tyrtyshnikov E.E. Breaking the curse of dimensionality, or how to use SVD in many dimensions // SIAM J. Sci. Comp. 2009. V. 31. № 5. P. 3744–3759.
Oseledets I.V., Tyrtyshnikov E.E. TT-cross approximation for multidimensional arrays // Lin. Alg. Appl. 2010. V. 432. № 1. P. 70–88.
Bellman R.E. Dynamic programming. Princeton Univ. Press, 1957.
Rossler A. Runge–Kutta methods for the strong approximation of solutions of stochastic differential equations // SIAM J. Numer. Anal. 2010. V. 3. № 48. P. 922–952.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал