Границы погрешности линеаризации с параллельным ускорением на GPU для систем реального времени
Робастное оптимальное управление нелинейными и нейросетевыми динамическими системами

Аннотация

В данной статье исследуется робастное оптимальное управление в реальном времени для неопределенных нелинейных систем, где линейно-изменяющиеся во времени (LTV) аппроксимации делают планирование выполнимым, но требуют надежных границ ошибки линеаризации (LEB) для гарантии робастного выполнения ограничений. Мы разрабатываем точные, дифференцируемые, параллельные на GPU LEB для LTV-аппроксимаций нелинейной динамики и динамики нейронных сетей (NN). Для аналитической динамики мы вводим основанные на путях границы Гессиана, которые являются более точными, чем стандартные интервальные методы. Для динамики NN мы выводим сертифицированные LEB, используя аффинные релаксации, сгенерированные верификатором NN, и локальные поправки Якобиана. Мы адаптируем решатель робастного управления на основе LTV с использованием системного синтеза (system-level synthesis), параллельный на GPU, для совместимости с этими LEB, расширяя его для работы с правообратимыми матрицами возмущений и множествами возмущений с ненулевым центром для точного распространения зонотопической неопределенности. Наш метод, GPUSLS-LEO, обеспечивает онлайн-оптимизацию робастных стратегий обратной связи, которые учитывают ошибку линеаризации, создавая точные, формально верифицированные достижимые трубки. На сложной нелинейной динамике и динамике NN размерностью до 168 состояний наш метод позволяет вычислять стратегии робастного управления на GPU с частотой до 67 Гц, сокращая время решения и консервативность по сравнению с базовыми методами при сохранении формальных гарантий и производительности в реальном времени..

I Введение

Безопасное нелинейное управление в реальном времени в условиях неопределенности имеет решающее значение для отказоустойчивой работы роботов, дронов и космических аппаратов. Это стимулирует разработку методов решения задачи робастного нелинейного оптимального управления (RNOCP), целью которой является поиск стратегии управления, гарантирующей выполнение ограничений при наихудших возмущениях. На практике решатели RNOCP обычно оптимизируют (1) номинальную траекторию и (2) стабилизирующий контроллер обратной связи для ее отслеживания при наличии возмущений. Робастное выполнение ограничений обеспечивается путем ужесточения номинальных ограничений с использованием аппроксимации сверху достижимой трубки ошибки слежения замкнутой системы, что гарантирует нахождение динамики замкнутой системы в безопасном подмножестве пространства состояний и входов. Однако точный анализ достижимости посредством распространения ошибки через нелинейную динамику [1, 2, 3], как правило, невыполним в реальном времени. Чтобы решить эту проблему, [4, 5, 6, 7] аппроксимируют нелинейную динамику вдоль номинальной траектории с помощью линейной нестационарной (LTV) системы и вычисляют достижимые множества, используя LTV-динамику. Синтез на системном уровне (SLS) [8, 4] является одной из таких структур робастного управления LTV, которую можно решать в реальном времени на больших горизонтах для многомерных систем [9, 10].

Несмотря на высокую эффективность, методы на основе LTV гарантируют соблюдение ограничений устойчивости только при наличии границы ошибки линеаризации, которая отражает отклонение между истинной нелинейной динамикой и ее LTV-аппроксимацией. Существующие методы часто используют глобальные границы [4, 11], что приводит к консервативным (широким) трубкам достижимости. Более точные локальные границы вокруг номинальной траектории часто возможны, но должны вычисляться эффективно для планирования в реальном времени и быть достаточно информативными, чтобы направлять оптимизатор к безопасным траекториям с минимальной ошибкой линеаризации. Таким образом, задача состоит в разработке границ, которые одновременно являются (1) точными, (2) вычислимыми в реальном времени и (3) дифференцируемыми, чтобы решатель RNOCP мог использовать их для изменения формы номинальной траектории и сужения трубок достижимости. Поскольку существующие методы не удовлетворяют всем трем требованиям, многие решатели RNOCP на основе LTV вынуждены игнорировать ошибку линеаризации, жертвуя формальными гарантиями ради производительности [12, 13].

⟦Fig. 1: (a): Устойчивые трубки из GPUSLS-LEO для нейронной системы толкателя T и возмущенные развертки для траектории вращения и толчка против случайных возмущений. (b): Развертка MPC с использованием нашего метода на реальной динамике, успешно перемещающая T к цели. (c): Устойчивые трубки положения $x$ и $y$, демонстрирующие узкие трубки и все смоделированные развертки, остающиеся внутри трубок.

Чтобы устранить этот разрыв, мы предлагаем семейство точных дифференцируемых границ ошибки линеаризации (LEBs) для решателей RNOCP с нелинейной динамикой, описываемой аналитическими функциями и нейронными сетями (NNs), с реализацией параллельных вычислений на GPU в JAX. Чтобы сделать эти LEBs совместимыми с решателем SLS реального времени из [10, 9], мы расширяем решатель: 1) для поддержки правообратимых матриц возмущений, что позволяет выполнять прямое зонотопическое распространение LEBs, избегая консервативных сверхаппроксимаций из [9]; и 2) для обработки множеств возмущений, не центрированных в нуле. Наконец, мы эмпирически показываем, что эти границы обеспечивают вычисление формально верифицированных решений RNOCP в реальном времени на GPU с ограниченным консерватизмом для систем с числом состояний до 168 и временем решения от 15 мс. Наш вклад заключается в следующем::

• ⟦ Набор точных, доказательно обоснованных методов оценки ошибки линеаризации, включая интервальные границы на основе путей для аналитических систем и границы на основе верификаторов NN для динамики NN.
• ⟦ Новая формулировка нелинейного SLS, включающая ошибку линеаризации посредством точного зонотопного распространения, включая обработку правообратимых матриц возмущений и множеств возмущений с ненулевым центром, что гарантирует робастное выполнение ограничений для исходной нелинейной системы.
• ⟦ Ускоренный графическим процессором решатель нелинейных SLS, GPUSLS-LEO, который использует формулировку этих границ линеаризации в JAX для обеспечения дифференцируемой, параллельной на графическом процессоре оценки и робастного нелинейного оптимального управления на основе градиента.
• ⟦ Оценка на аналитических и NN моделях динамики с поддержанием частоты управления в реальном времени (до 67 Гц) для систем размерностью до 168D и задач с $\approx 2\times 10^{5}$ переменными решения.

II Связанные работы

Хотя существуют методы анализа достижимости без линеаризации [1, 2, 3], прямое распространение ошибки через нелинейную динамику, как правило, является вычислительно неразрешимым в режиме реального времени, особенно при размерности более $\approx$10 состояний [3, 1]. Таким образом, решатели RNOCP часто полагаются на линеаризованное распространение ошибки для повышения эффективности [14, 13], однако это вносит ошибку линеаризации, которую трудно строго ограничить, и которую часто игнорируют [14, 15, 13, 10], недооценивают посредством выборки [11, 4] или обрабатывают, предполагая наличие оракула границ линеаризации [16, 17, 18, 19, 20]. Решатели RNOCP, которые явно вычисляют LEB, обычно полагаются на глобальные оценки Липшица [11] или интервальную арифметику [7, 5, 6, 21]. Интервальная арифметика использовалась как для вычисления достижимых множеств [22, 23], так и для поддержки анализа достижимости путем построения границ линейного дифференциального включения (LDI) для нелинейной динамики [24]. В частности, в [25] используется ограничение якобиана для получения более точных LDI для анализа достижимости, но этот подход не применяется для проектирования управления или прямого вычисления LEB. С вычислительной точки зрения существующие методы ограничения линеаризации являются автономными [4], медленными или недифференцируемыми [7], что ограничивает их использование при планировании. Чтобы устранить эти пробелы, наш метод использует интервальную арифметику и JAX для вычисления точных дифференцируемых LEB для непосредственного использования в оптимальном управлении. Более того, мы также адаптируем инструменты верификации neural-network (NN) [26] для вывода аффинных границ ошибки линеаризации для динамики NN..

Поскольку вычисление LEB может занимать много времени [4], их часто генерируют один раз в автономном режиме [4, 11], однако они оказываются неточными, если мы хотим ограничить ошибку вокруг одной траектории, как в RMPC. Таким образом, специалисты часто пренебрегают ошибкой линеаризации, получая эвристические достижимые трубки, которые жертвуют формальными гарантиями ради производительности [12]. Чтобы решить эту проблему, мы показываем, что точные локальные LEB с параллельным ускорением на GPU могут вычисляться в режиме реального времени, обеспечивая при этом как быстродействие для многомерных систем, так и строгие гарантии..

III Предварительные сведения и постановка задачи

Обозначения: Мы обозначаем $I_{n}\in\mathbb{R}^{n\times n}$ как единичную матрицу, норму Фробениуса матрицы $A\in\mathbb{R}^{m\times n}$ как $\|A\|_{\mathcal{F}}:=\sqrt{\text{Trace}(A^{\top}A)}$. Пусть $[N]:=\{0,...,N-1\}$, и $[M,N]:=\{M,...,N\}$ для $M,N\in\mathbb{N}$. Для матрицы $A\in\mathbb{R}^{n\times p}$ мы определяем построчную $\ell_{1}$-норму $\|A\|_{1,\text{r}}\in\mathbb{R}^{n}$ как $\|A\|_{1,\text{r}}\coloneqq\left[\|A_{1,:}\|_{1},\ldots,\|A_{n,:}\|_{1}\right]^{\top}$. Для $a\in\mathbb{R}^{n}$ мы определяем $\text{diag}(a)\in\mathbb{R}^{n\times n}$ как её диагональную матрицу. Для заданной $A\in\mathbb{R}^{m\times n}$ с линейно независимыми строками мы определяем $A^{\ddagger}\in\mathbb{R}^{n\times m}$ так, что $AA^{\ddagger}=I_{m}$, где $A^{\ddagger}:=A^{\top}(AA^{\top})^{-1}$. Мы обозначаем сумму Минковского двух множеств $A$ и $B$ как $A\oplus B:=\{a+b\mid a\in A,b\in B\}$. Для функции $p(x)$ с $x\in\mathbb{R}^{n}$ (или явно $p\left((x_{1},\dots,x_{n})\right)$) мы определяем $\partial_{j}p(\alpha):=\frac{\partial p}{\partial x_{j}}|_{x=\alpha}$ и $\partial_{jk}p(\alpha):=\frac{\partial^{2}p}{\partial x_{j}\partial x_{k}}|_{x=\alpha}$ как первую и вторую частные производные по $j$-му и $k$-му компонентам векторного входа. Мы определяем $\mathbb{I}_{n}$ как множество всех $n$-мерных интервалов. Пусть $\Delta,\hat{\Delta}\in\mathbb{I}_{n}$, тогда $\Delta\cdot\hat{\Delta}:=\{\delta\cdot\hat{\delta}\mid\delta\in\Delta,\hat{\delta}\in\hat{\Delta}\}$⟧.

Определения: Мы рассматриваем неопределенную нелинейную динамику

$$x_{k+1}=f(x_{k},u_{k})+E(x_{k})w_{k},$$ (1)

где $x_{k}\in\mathcal{X}\subseteq\mathbb{R}^{n_{x}}$ — состояние системы на временном шаге $k$, $u_{k}\in\mathcal{U}\subseteq\mathbb{R}^{n_{u}}$ обозначает управляющее воздействие, $f:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{x}}$ — функцию динамики, $E:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{x}\times n_{x}}$ — функцию масштабирования возмущений, а $w_{k}\in\mathcal{E}_{n_{x}}\coloneqq\{w\in\mathbb{R}^{n_{x}},\|w\|_{\infty}\leq 1\}$ — возмущение, нормализованное таким образом, чтобы оно находилось внутри единичного $\ell_{\infty}$-шара.

Мы стремимся разработать контроллер $\pi(\cdot)$, который решает следующую задачу робастного нелинейного оптимального управления (RNOCP):

$$\min_{\pi(\cdot)}\quad J\left(\bar{x},\pi(\cdot)\right)$$ (2a)

s.t. $$x_{k+1}=f(x_{k},u_{k})+E(x_{k})w_{k},\quad\forall k\in[N],$$ (2b)

$$x_{0}=\bar{x}_{0},$$ (2c)

$$u_{k}=\pi_{k}(x_{0},\dots,x_{k}),\quad\forall k\in[N],$$ (2d)

$$g_{k}(x_{k},u_{k})\leq 0,\quad\forall w_{k}\in\mathcal{E}_{n_{x}},\quad\forall k\in[N],$$ (2e)

$$g^{f}(x_{N})\leq 0,\quad\forall w_{k}\in\mathcal{E}_{n_{x}},\quad\forall k\in[N].$$ (2f)

где $\pi:=\{\pi_{i}\}_{i=0}^{N-1}$ представляет собой последовательность причинно-следственных стратегий управления, $\bar{x}_{0}\in\mathbb{R}^{n_{x}}$ — начальное состояние, функции $g_{k}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{c}}$ обозначают поэтапные ограничения на состояния и входы, а $g^{f}:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{f}}$ — терминальное ограничение. Мы определяем замкнутое достижимое множество $\Omega_{k}:=(\Omega_{k}^{x},\Omega_{k}^{u})$ состояний и управлений как

$$\small\hskip-7.0pt\Omega_{k}:=\left\{(x_{k},u_{k})\;\middle|\;\begin{aligned} &x_{j+1}=f(x_{j},u_{j})+E(x_{j})w_{j},\ \ \forall j\in[k],\\ &\hskip 18.49988pt\forall w_{j}\in\mathcal{E}_{n_{x}},\hskip 9.24994pt\forall j\in[k],\\ &u_{j}=\pi_{j}(x_{0},\ldots,x_{j}),\hskip 9.24994pt\forall j\in[k+1].\end{aligned}\right\}.$$ (3)

Поскольку (2) является неразрешимой бесконечномерной задачей, решатели RNOCP часто упрощают её путем вычисления (A) номинальной траектории состояния-входа и (B) следящего контроллера вокруг этой траектории. В данной статье мы используем один из таких методов, GPUSLS [10], который адаптирует [9] для решения в реальном времени на GPU. Затем мы расширяем [10], чтобы учесть ошибку линеаризации в разд. V. Наш подход использует SLS, который мы теперь рассмотрим.

III-A Синтез на системном уровне (SLS)

Рассмотрим следующую неопределенную LTV динамику

$$x_{k+1}=A_{k}x_{k}+B_{k}u_{k}+E_{k}w_{k}.$$ (4)

SLS — это концепция проектирования управления, которая выполняет оптимизацию по причинно-следственным регуляторам обратной связи по возмущению

$$u_{k}=v_{k}+\sum_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{u}}_{k,j}E_{j}w_{j},$$ (5)

с номинальным управлением $v_{k}\in\mathbb{R}^{n_{u}}$, то есть мы назначаем матрицу обратной связи по возмущению $\mathbf{\Phi}^{\mathrm{u}}_{k,j}\in\mathbb{R}^{n_{u}\times n_{x}}$ для каждого возмущения $E_{j}w_{j}$ и управление $u_{k}$ для $k>j$. Для динамики LTV (4) можно показать с помощью алгебраических преобразований [8], что результирующая последовательность состояний замкнутой системы может быть выражена как

$$x_{k}=z_{k}+\sum_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{x}}_{k,j}E_{j}w_{j},~z_{0}=\bar{x}_{0},$$ (6)

где $z_{k}\in\mathbb{R}^{n_{x}}$ — номинальное состояние, а замкнутый отклик $\mathbf{\Phi}^{\mathrm{x}}_{k,j}\in\mathbb{R}^{n_{x}\times n_{x}}$ отражает влияние возмущения $w_{j}$ на состояние $x_{k}$. Начиная с $\mathbf{\Phi}^{\mathrm{x}}_{j+1,j}=I_{n_{x}}$, SLS распространяет возмущение через $\mathbf{\Phi}^{\mathrm{x}}_{k+1,j}=A_{k}\mathbf{\Phi}^{\mathrm{x}}_{k,j}+B_{k}\mathbf{\Phi}^{\mathrm{u}}_{k,j}$, для всех $j\in[N]$ и $k\in[j+1,N-1]$. Поскольку (5)–(6) являются истинной траекторией замкнутой системы при заданной последовательности возмущений, точные достижимые множества замкнутой системы представляют собой $\Omega_{x}^{k}=z_{k}+\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{x}}_{k,j}E_{j}\mathcal{E}_{n_{x}}$ и $\Omega_{k}^{u}=v_{k}+\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{u}}_{k,j}E_{j}\mathcal{E}_{n_{x}}$. Чтобы применить методы разд. III-A к нелинейной динамике (1), мы аппроксимируем (1) с помощью номинальной динамики $z_{k+1}=f(z_{k},v_{k})$ и LTV-динамики ошибки вдоль номинальной траектории состояния-управления $\mathbf{z}:=\{z_{k}\}_{k=0}^{N}\in\mathbb{R}^{n_{x}\times(N+1)}$, $\mathbf{v}:=\{v_{k}\}_{k=0}^{N-1}\in\mathbb{R}^{n_{u}\times N}$. Для номинального $(z_{k},v_{k})\in\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}$ мы строим LTV-аппроксимацию вида (4), определяя

$$A_{k}:=\left.\frac{\partial f}{\partial x}\right|_{(z_{k},v_{k})},B_{k}:=\left.\frac{\partial f}{\partial u}\right|_{(z_{k},v_{k})},E_{k}:=E(z_{k}).$$ (7)

Обозначим $\delta z_{k}:=x_{k}-z_{k}$ и $\delta v_{k}:=u_{k}-v_{k}$. Тогда остаточную ошибку $d_{k}(x_{k},u_{k},w_{k})$ между (1) и её аппроксимацией Тейлора первого порядка относительно $(z_{k},v_{k})$, то есть $f(z_{k},v_{k})+A_{k}\delta z_{k}+B_{k}\delta v_{k}$, можно записать как,

$$d_{k}(x_{k}, u_{k},w_{k}):=E_{k}w_{k}+\underbrace{r_{k}^{\mathrm{lin}}+r_{k}^{\mathrm{dist}}}_{\hskip 16.38895pt:=r(x_{k},u_{k},w_{k})},$$ (8a)

$$r_{k}^{\mathrm{dist}} :=\bigl(E(x_{k})-E_{k}\bigr)w_{k},$$ (8b)

$$r_{k}^{\mathrm{lin}} :=f(x_{k},u_{k})-f(z_{k},v_{k})-A_{k}\delta z_{k}-B_{k}\delta v_{k}.$$ (8c)

В частности, (8c) учитывает остаточный член Тейлора нелинейной динамики $f$, в то время как (8b) учитывает ошибку, вызванную зависимостью матрицы возмущений $E$ от состояния. Учитывая данную постановку, мы изучаем следующие задачи.

III-B Постановка задачи

Задача 1 (Граница ошибки линеаризации). Для всех $k\in[N]$ найти множественнозначную аппроксимацию сверху $\bar{\mathcal{E}}_{k}\subseteq\mathbb{R}^{n_{x}}$, а также соответствующую масштабирующую матрицу $\bar{R}_{k}\in\mathbb{R}^{n_{x}\times n_{x}}$ и смещение $\bar{c}_{k}\in\mathbb{R}^{n_{x}}$, где

$$r_{k}(x_{k},u_{k},w_{k}) \in\bar{\mathcal{E}}_{k}\subseteq\big(\bar{R}_{k}\mathcal{E}_{n_{x}}\oplus\bar{c}_{k}\big),$$ (9)

$$\forall(x_{k},u_{k})\overset{\eqref{eq:reach}}{\in}\Omega_{k},\ \forall w_{k}\in\mathcal{E}_{n_{x}}.$$

Задача 2 (RNOCP в реальном времени). Используя $\bar{\mathcal{E}}_{k}$ для учета расхождения между (1) и его LTV-аппроксимацией, эффективно решите RNOCP (2), используя суррогатную модель LTV (4) для надежного удовлетворения ограничений (2e) и (2f).

IV Границы погрешности линеаризации с параллельным ускорением на GPU

Для упрощения изложения мы рассматриваем векторнозначную функцию $p(\bm{\alpha}):\mathbb{R}^{n}\to\mathbb{R}^{m}$ от расширенного состояния $\bm{\alpha}:=[x^{\top},u^{\top}]^{\top}$. Пусть $\mathcal{A}\subset\mathbb{R}^{n}$ — компактное множество, а $\tilde{\bm{\alpha}}\in\mathcal{A}$ — номинальная точка. Ошибка линеаризации, то есть остаточный член Лагранжа разложения в ряд Тейлора первого порядка для $p$ в точке $\tilde{\bm{\alpha}}$, имеет вид,

$$\bm{r}(\bm{\alpha};\tilde{\bm{\alpha}}):=p(\bm{\alpha})-p(\tilde{\bm{\alpha}})-J(\tilde{\bm{\alpha}})(\bm{\alpha}-\tilde{\bm{\alpha}}),\quad\bm{\alpha}\in\mathcal{A}.$$ (10)

Для частного случая, когда $p(\bm{\alpha})=f(x,u)$, (10) совпадает с (8c). Наша цель — вычислить корректную покомпонентную интервальную оболочку $\mathcal{R}\in\mathbb{I}_{m}$ для $\bm{r}(\mathcal{A};\tilde{\bm{\alpha}})$, то есть $\bm{r}(\mathcal{A};\tilde{\bm{\alpha}})\subseteq\mathcal{R}$.

IV-A Границы погрешности линеаризации для аналитической динамики

Согласно теореме Тейлора в форме Лагранжа, для каждого компонента $\ell\in\{1,\dots,m\}$ существует $\xi\in\mathcal{A}$ такое, что

$$\bm{r}_{\ell}(\bm{\alpha};\tilde{\bm{\alpha}})=\tfrac{1}{2}\delta^{\top}H_{\ell}(\xi)\delta,\quad\delta:=\bm{\alpha}-\tilde{\bm{\alpha}}.$$ (11)

где $H_{\ell}(\xi)$ является $\ell$-м элементом Гессиана $H$ в точке $\xi$. Для оценки интервального включения $\bm{r}_{\ell}(\mathcal{A};\tilde{\bm{\alpha}})$ в предыдущих работах [7, 5, 6, 21] отдельно вычисляются $\delta^{\top},H_{\ell}(\xi)$ и $\delta$ в интервальной арифметике, а затем выполняются интервальные произведения, как в (12), что дает включение $\mathcal{R}_{\ell,\mathrm{classic}}$⟧:

$$\bm{r}_{\ell}(\mathcal{A};\tilde{\bm{\alpha}})\subseteq\mathcal{R}_{\ell,\mathrm{classic}}:=\tfrac{1}{2}\Delta^{\top}[H_{\ell}](\mathcal{A})\Delta,$$ (12)

где $\Delta:=\mathcal{A}\ominus\tilde{\bm{\alpha}}\subset\mathbb{R}^{n}$ обозначает множество смещений, а $[H_{\ell}](\mathcal{A})$ представляет собой интервальное приближение Гессиана на $\mathcal{A}$ с элементами $[H_{\ell,jk}](\mathcal{A})$ (где $H_{\ell,jk}(\bm{\alpha}):=\partial_{jk}p_{\ell}(\bm{\alpha})$). На практике $\tilde{\bm{\alpha}}$ обычно выбирается в центре $\mathcal{A}$ [7], что приводит к симметричному $\Delta$ для более узкого приближения $\bm{r}_{\ell}(\mathcal{A};\tilde{\bm{\alpha}})$⟧.

Сохраняя строгую квадратичную форму,  (11) можно сначала разложить покомпонентно, а затем выполнить интервальную арифметику над этим разложением. Это дает (13), где интервал $\Delta_{j}$ обозначает $j$-й срез $\Delta$ после применения симметрии интервального гессиана. Когда $\Delta$ является симметричным, $\mathcal{R}_{\ell,\mathrm{expand}}\subseteq\mathcal{R}_{\ell,\mathrm{classic}}$ гарантируется, поскольку (13) сохраняет неотрицательность $\delta_{k}^{2}$ для всех $\delta_{k}\in\Delta_{k}$ и для всех $k\in\{1,\ldots,n\}$:

$$\bm{r}_{\ell}(\mathcal{A};\tilde{\bm{\alpha}})\subseteq\;\mathcal{R}_{\ell,\mathrm{expand}}= \textstyle\sum_{\begin{subarray}{c}1\leq j<k\leq n\end{subarray}}[H_{\ell,kj}](\mathcal{A})\,\Delta_{j}\cdot\Delta_{k}$$

$$+ \tfrac{1}{2}\textstyle\sum_{k=1}^{n}[H_{\ell,kk}](\mathcal{A})\,\Delta_{k}^{2}.$$ (13)

Однако классические границы в 12 и 13 оценивают все интервальные оболочки элементов Гессиана по всему множеству $\mathcal{A}$, что приводит к консервативности из-за эффектов зависимости и наматывания. Чтобы смягчить это, мы используем покоординатный путь от $\tilde{\bm{\alpha}}$ до $\bm{\alpha}$, который позволяет получить более точные границы для членов второго порядка. Для любого $\bm{\alpha}\in\mathcal{A}$ определим последовательность промежуточных точек

$${\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(k)}:=(\bm{\alpha}_{1},\dots,\bm{\alpha}_{k},\tilde{\bm{\alpha}}_{k+1},\dots,\tilde{\bm{\alpha}}_{n}),\quad k\in[n+1],$$ (14)

с ${\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(0)}=\tilde{\bm{\alpha}}$ и ${\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(n)}=\bm{\alpha}$. Это определяет покоординатный путь [24, 25], который инкрементально перемещается от $\tilde{\bm{\alpha}}$ к $\bm{\alpha}$ вдоль каждой координаты. Соответственно, определим префиксные множества $\mathcal{A}^{(k)}(\tilde{\bm{\alpha}}):=\{\bm{\alpha}\in\mathcal{A}\mid\bm{\alpha}_{\ell}=\tilde{\bm{\alpha}}_{\ell},\forall\ell\in[k+1,n]\}$.

Theorem 1 (Граница Гессиана на основе путей)

Пусть $p:\mathbb{R}^{n}\to\mathbb{R}^{m}$ будет дважды непрерывно дифференцируемой на интервале $\mathcal{A}\subset\mathbb{R}^{n}$, и пусть $\tilde{\bm{\alpha}}\in\mathcal{A}$. Определим множество смещений $\Delta:=\mathcal{A}\ominus\tilde{\bm{\alpha}}$. Для каждого $\bm{\alpha}\in\mathcal{A}$ определим покоординатный путь $\bm{\alpha}^{(k)}$ и соответствующие множества префиксов $\mathcal{A}^{(k)}(\tilde{\bm{\alpha}})$. Тогда для каждого компонента $\ell\in\{1,\dots,m\}$ ошибка линеаризации удовлетворяет

$$\hskip-6.0pt\bm{r}_{\ell}(\mathcal{A};\tilde{\bm{\alpha}})\subseteq\mathcal{R}_{\ell,\mathrm{path}} \hskip-3.0pt:=\hskip-3.0pt\textstyle\sum_{1\leq j<k\leq n}[H_{\ell,kj}](\mathcal{A}^{(j)}(\tilde{\bm{\alpha}}))\,\Delta_{j}\cdot\Delta_{k}$$

$$+\tfrac{1}{2}\textstyle\sum_{k=1}^{n}[H_{\ell,kk}](\mathcal{A}^{(k)}(\tilde{\bm{\alpha}}))\,\Delta_{k}^{2}.$$ (15)

Доказательство:

Зафиксируем $\bm{\alpha}\in\mathcal{A}$ и определим $\delta:=\bm{\alpha}-\tilde{\bm{\alpha}}$. Используя метод телескопической суммы вдоль покомпонентного пути, получаем $p_{\ell}(\bm{\alpha})-p_{\ell}(\tilde{\bm{\alpha}})=\textstyle\sum_{k=1}^{n}\big(p_{\ell}({\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(k)})-p_{\ell}({\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(k-1)})\big)$. Для каждого $k$ применение теоремы Тейлора в одномерном случае дает

$$p_{\ell}({\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(k)})-p_{\ell}({\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(k-1)})=\partial_{k}p_{\ell}({\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(k-1)})\delta_{k}+\tfrac{1}{2}\partial_{kk}p_{\ell}(\xi_{k})\,\delta_{k}^{2},$$ (16)

для некоторого $\xi_{k}\in[{\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(k-1)},{\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(k)}]\subseteq\mathcal{A}^{(k)}(\tilde{\bm{\alpha}})$. Суммирование (16) по $k$ и вычитание $\nabla p_{\ell}(\tilde{\bm{\alpha}})^{\top}\delta$ дает $p_{\ell}(\bm{\alpha})-p_{\ell}(\tilde{\bm{\alpha}})-\nabla p_{\ell}(\tilde{\bm{\alpha}})^{\top}\delta$ в левой части, что равно $\bm{r}_{\ell}(\bm{\alpha};\tilde{\bm{\alpha}})$, а правая часть находится в $\textstyle\sum_{k=1}^{n}\big(\partial_{k}p_{\ell}({\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(k-1)})-\partial_{k}p_{\ell}(\tilde{\bm{\alpha}})\big)\delta_{k}+\tfrac{1}{2}\textstyle\sum_{k=1}^{n}\partial_{kk}p_{\ell}(\xi_{k})\,\delta_{k}^{2}$. Разлагая разность градиентов вдоль того же пути для $k$, получаем $\partial_{k}p_{\ell}({\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(k-1)})-\partial_{k}p_{\ell}(\tilde{\bm{\alpha}})=\sum_{j=1}^{k-1}\big(\partial_{k}p_{\ell}({\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(j)})-\partial_{k}p_{\ell}({\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(j-1)})\big)$. Согласно теореме о среднем значении, для каждого $j<k$ существует $\eta_{j,k}\in[{\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(j-1)},{\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(j)}]\subseteq\mathcal{A}^{(j)}$ такое, что $\partial_{k}p_{\ell}({\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(j)})-\partial_{k}p_{\ell}({\bm{\alpha}}_{\tilde{\bm{\alpha}}}^{(j-1)})=\partial_{kj}p_{\ell}(\eta_{j,k})\,\delta_{j}$. Подстановка дает $\bm{r}_{\ell}(\bm{\alpha};\tilde{\bm{\alpha}})=\sum_{1\leq j<k\leq n}\partial_{kj}p_{\ell}(\eta_{j,k})\,\delta_{j}\delta_{k}+\tfrac{1}{2}\sum_{k=1}^{n}\partial_{kk}p_{\ell}(\xi_{k})\,\delta_{k}^{2}$. Поскольку $\eta_{j,k}\in\mathcal{A}^{(j)}(\tilde{\bm{\alpha}})$, $\xi_{k}\in\mathcal{A}^{(k)}(\tilde{\bm{\alpha}})$ и $\delta\in\Delta$, взятие интервальных оболочек дает (15). ∎

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

Corollary 1 (Путевой Гессиан никогда не бывает слабее)

Когда $\Delta$ является симметричной, гарантируется, что охват в (15) не менее точен, чем (13) и (12). То есть,

$$\mathcal{R}_{\ell,\mathrm{path}}\subseteq\mathcal{R}_{\ell,\mathrm{expand}}\subseteq\mathcal{R}_{\ell,\mathrm{classic}}.$$ (17)

Доказательство:

Для каждого $j\in\{1,\dots,n\}$, множество префиксов удовлетворяет $\mathcal{A}^{(j)}\subseteq\mathcal{A}$ по построению. Затем, для каждого элемента Гессиана,

$$[H_{\ell,kj}](\mathcal{A}^{(j)})\subseteq[H_{\ell,kj}](\mathcal{A}),[H_{\ell,kk}](\mathcal{A}^{(k)})\subseteq[H_{\ell,kk}](\mathcal{A}).$$

Умножение на $\Delta_{j}\Delta_{k}$ и $\Delta_{k}^{2}$ сохраняет включение, а суммирование полученных членов дает $\mathcal{R}_{\ell,\mathrm{path}}\subseteq\mathcal{R}_{\ell,\mathrm{expand}}$. ∎

В отличие от $\mathcal{R}_{\ell,\mathrm{classic}}$, $\mathcal{R}_{\ell,\mathrm{expand}}$ и $\mathcal{R}_{\ell,\mathrm{path}}$ сохраняют покомпонентное расширение для получения точных интервалов со смещенным центром (посредством квадратичной формы, а не симметричных границ). Мы используем повышенную точность в контексте RNOCP (разд. V).

IV-B Границы погрешности линеаризации для нейронной динамики

Для динамики NN $f$ вычисление точного интервала гессиана может быть затруднительным. Вместо этого мы используем CROWN [26, 27] для получения аффинных границ NN на множестве входных данных в виде параллелепипеда и преобразуем их в интервальное включение $\bm{r}(\bm{\alpha};\tilde{\bm{\alpha}})$⟧.

Theorem 2 (Границы CROWN (Thm. 3.2 [26]))

Для заданной NN, $p:\mathbb{R}^{n}\to\mathbb{R}^{m}$, и ограниченного входного сигнала $\bm{\alpha}\in\mathcal{A}$, где $\mathcal{A}$ является интервалом, существуют линейные нижние/верхние границы для $p(\bm{\alpha})$:

$$\mathbf{\underline{W}}\bm{\alpha}+\mathbf{\underline{b}}\leq p(\bm{\alpha})\leq\mathbf{\overline{W}}\bm{\alpha}+\mathbf{\overline{b}},$$ (18)

где $\mathbf{\underline{W}},\mathbf{\overline{W}}\in\mathbb{R}^{m\times n}$ и $\mathbf{\underline{b}},\mathbf{\overline{b}}\in\mathbb{R}^{m}$⟧.

Граница в 18 не ограничивает напрямую $\bm{r}(\bm{\alpha};\tilde{\bm{\alpha}})$, поскольку наклоны $\mathbf{\underline{W}},\mathbf{\overline{W}}$, как правило, не совпадают с якобианом $J(\tilde{\bm{\alpha}})$. Мы решаем эту проблему посредством дальнейшей релаксации в теореме 3.⟧.

Theorem 3 (Граница ошибки линеаризации от CROWN)

Пусть $p:\mathbb{R}^{n}\to\mathbb{R}^{m}$ — это NN, $\mathcal{A}\subset\mathbb{R}^{n}$ — это параллелепипед, $\tilde{\bm{\alpha}}\in\mathcal{A}$, $J(\tilde{\bm{\alpha}})$ — это якобиан $p$ в точке $\tilde{\bm{\alpha}}$, и пусть $(\mathbf{\underline{W}},\mathbf{\underline{b}},\mathbf{\overline{W}},\mathbf{\overline{b}})$ удовлетворяет границам CROWN (18) на $\mathcal{A}$. Тогда $\forall\bm{\alpha}\in\mathcal{A}$⟧:

$$\underline{\bm{r}}(\bm{\alpha})\leq\bm{r}(\bm{\alpha};\tilde{\bm{\alpha}})\leq\overline{\bm{r}}(\bm{\alpha}),$$ (19)

$$\underline{\bm{r}}(\bm{\alpha}) :=(\mathbf{\underline{W}}-J(\tilde{\bm{\alpha}}))\bm{\alpha}+J(\tilde{\bm{\alpha}})\tilde{\bm{\alpha}}-p(\tilde{\bm{\alpha}})+\mathbf{\underline{b}},$$ (20)

$$\overline{\bm{r}}(\bm{\alpha}) :=(\mathbf{\overline{W}}-J(\tilde{\bm{\alpha}}))\bm{\alpha}+J(\tilde{\bm{\alpha}})\tilde{\bm{\alpha}}-p(\tilde{\bm{\alpha}})+\mathbf{\overline{b}}.$$

Затем интервальное замыкание $\bm{r}(\mathcal{A};\tilde{\bm{\alpha}})\subseteq[\underline{\bm{r}},\overline{\bm{r}}]$ задается как

$$\underline{\bm{r}}_{i}=\min_{\bm{\alpha}\in\mathcal{A}}\;\underline{\bm{r}}_{i}(\bm{\alpha}),\quad\overline{\bm{r}}_{i}=\max_{\bm{\alpha}\in\mathcal{A}}\;\overline{\bm{r}}_{i}(\bm{\alpha}),\quad\forall i\in[1,m],$$ (21)

Доказательство:

Подставляя нижнюю границу CROWN 18 в выражение для $\bm{r}(\bm{\alpha};\tilde{\bm{\alpha}})$ 10, мы получаем $\bm{r}(\bm{\alpha};\tilde{\bm{\alpha}})\geq\mathbf{\underline{W}}\bm{\alpha}+\mathbf{\underline{b}}-p(\tilde{\bm{\alpha}})-J(\tilde{\bm{\alpha}})(\bm{\alpha}-\tilde{\bm{\alpha}})$, что упрощается до $\bm{r}(\bm{\alpha};\tilde{\bm{\alpha}})\geq(\mathbf{\underline{W}}-J(\tilde{\bm{\alpha}}))\bm{\alpha}+J(\tilde{\bm{\alpha}})\tilde{\bm{\alpha}}-p(\tilde{\bm{\alpha}})+\mathbf{\underline{b}}=\underline{\bm{r}}(\bm{\alpha}).$ Аналогичным образом мы применяем верхнюю границу CROWN. Таким образом, $\forall\bm{\alpha}\in\mathcal{A}$: $\underline{\bm{r}}(\bm{\alpha})\leq\bm{r}(\bm{\alpha};\tilde{\bm{\alpha}})\leq\overline{\bm{r}}(\bm{\alpha})$. Для границ, справедливых на $\mathcal{A}$, мы берем поэлементный минимум $\underline{\bm{r}}(\bm{\alpha})$ и максимум $\overline{\bm{r}}(\bm{\alpha})$ по $\mathcal{A}$, как в 21). ∎

Теорема 3 согласует наклоны CROWN и наклон локальной линеаризации $J(\tilde{\bm{\alpha}})$, тем самым создавая корректную границу остатка вокруг номинальной точки. Поскольку $\underline{\bm{r}}(\bm{\alpha})$ и $\overline{\bm{r}}(\bm{\alpha})$ являются аффинными относительно $\bm{\alpha}$, задача (21) не требует больших вычислительных затрат, когда $\mathcal{A}$ представляет собой параллелепипед..

IV-C Зонотопные границы ⟦K0⟧ $r_{k}(x_{k},u_{k},w_{k})$

Для получения границы зонотопа для $r_{k}(x_{k},u_{k},w_{k})$ (9) мы вычисляем оболочку $[\underline{r}_{k}^{\mathrm{lin}},\overline{r}_{k}^{\mathrm{lin}}]$ ошибки линеаризации $r_{k}^{\mathrm{lin}}$, как в разд. IV-A или IV-B, и оболочку $[\underline{r}_{k}^{\mathrm{dist}},\overline{r}_{k}^{\mathrm{dist}}]$ ошибки возмущения $r_{k}^{\mathrm{dist}}$ с помощью интервальной арифметики в (22), где $[E](\Omega_{k}^{x})$ — это интервальная оболочка $E$ на $\Omega_{k}^{x}$⟧:

$$r_{k}^{\mathrm{dist}}=\bigl(E(x_{k})-E_{k}\bigr)w_{k}\in\bigl([E](\Omega_{k}^{x})\ominus E_{k}\bigr)\mathcal{E}_{n_{x}}.$$ (22)

Наконец, мы получаем $r_{k}(x_{k},u_{k},w_{k})\in[\underline{r}_{k},\overline{r}_{k}]$, где $\underline{r}_{k}=\underline{r}_{k}^{\mathrm{lin}}+\underline{r}_{k}^{\mathrm{dist}}$ и $\overline{r}_{k}=\overline{r}_{k}^{\mathrm{lin}}+\overline{r}_{k}^{\mathrm{dist}}$. Его зонотопная граница задается выражением $\big(\bar{R}_{k}\mathcal{E}_{n_{x}}\oplus\bar{c}_{k}\big)$, где $\bar{c}_{k}=\tfrac{1}{2}(\underline{r}_{k}+\overline{r}_{k})$ и $\bar{R}_{k}=\tfrac{1}{2}\mathrm{diag}(\overline{r}_{k}-\underline{r}_{k})$. Мы обозначаем это отображение компактно как

$$\hskip-1.0pt\bar{R}_{k}\hskip-1.0pt=\hskip-1.0pt\Gamma_{r}(\Omega_{k},z_{k},v_{k};f,E),\hskip 6.0pt\bar{c}_{k}\hskip-1.0pt=\hskip-1.0pt\Gamma_{c}(\Omega_{k},z_{k},v_{k};f,E).$$ (23)

IV-D Дифференцируемая параллельная реализация границ

Мы реализуем наши LEBs в JAX, при этом PBH использует immrax для интервальной арифметики, а jax_verify предоставляет границы на основе CROWN для динамики NN. JAX обеспечивает дифференцируемость и параллелизацию на GPU, что позволяет выполнять уточнение номинальной точки $\tilde{\bm{\alpha}}$ на основе градиента и масштабируемое разбиение входной области для снижения консервативности для SLS.

V Использование границ ошибки линеаризации в GPUSLS

Мы обсуждаем использование наших LEBs в SLS (разд. V-A) и описываем эффективную реализацию на базе GPU (разд. V-B⟧).

V-A Формулировка нелинейного SLS

Следуя [9, App. A], SLS может быть расширен на нелинейные системы путем планирования номинальной траектории для динамики без возмущений $z_{k+1}=f(z_{k},v_{k})$ и моделирования динамики ошибки слежения через LTV-систему из (4). Для приближенного решения (2) мы модифицируем формулировку нелинейного SLS из [9, App. A] так, чтобы она была совместима с 1) множествами возмущений с ненулевым центром и 2) LEB, приведенными в разд. IV-A-IV-B⟧:

$$\hskip-13.0pt\min_{\begin{subarray}{c}\mathbf{z},\,\mathbf{v},\mathbf{\Phi}\end{subarray}}\ \ J(\mathbf{z},\mathbf{v})+\tilde{H}_{0}(\mathbf{\Phi},\bm{\mathcal{D}})$$ (24a)

s.t. $$z_{k+1}=f(z_{k},v_{k}),\quad\forall k\in[N],\quad z_{0}=\bar{x}_{0},$$ (24b)

$$\mathbf{\Phi}^{\mathrm{x}}_{k+1,j}=A_{k}\mathbf{\Phi}^{\mathrm{x}}_{k,j}+B_{k}\mathbf{\Phi}^{\mathrm{u}}_{k,j},\ \ \mathbf{\Phi}^{\mathrm{x}}_{j+1,j}=I_{n_{x}},$$ (24c)

$$\quad\quad\forall j\in[N],\hskip 5.0pt\forall k\in[j+1,N-1],$$

$$\bar{E}_{k}:=\mathcal{Z}_{r}(z_{k},v_{k},\bm{\mu}_{k},\bm{\sigma}_{k}),\quad\forall k\in[N],$$ (24d)

$$\bar{c}_{k}:=\mathcal{Z}_{c}(z_{k},v_{k},\bm{\mu}_{k},\bm{\sigma}_{k}),\quad\forall k\in[N],$$ (24e)

$$\bm{\sigma}_{k}=h^{\sigma}_{k}(\mathbf{\Phi},\bm{\mathcal{D}}),~\bm{\mu}_{k}=h^{\mu}_{k}(\mathbf{\Phi},\bm{\mathcal{D}}),\,\,\forall k\in[N],$$ (24f)

$$g_{k}(z_{k},v_{k})+h_{k}(\mathbf{\Phi},\bm{\mathcal{D}})\leq 0,\qquad\forall k\in[N],\hskip-5.0pt$$ (24g)

$$g^{f}(z_{N})+h^{f}(\mathbf{\Phi},\bm{\mathcal{D}})\leq 0,$$ (24h)

где $A_{k},B_{k}$ — это линеаризованная динамика (7) в момент времени $k$, а функции $\mathcal{Z}_{r}(\cdot),\mathcal{Z}_{c}(\cdot)$ (24d)–(24e) представляют собой суммы Минковского внешнего возмущения и ошибки линеаризации как функции номинальных значений ($z_{k}$,$v_{k}$), смещения центра достижимой трубки $\bm{\mu}_{k}\in\mathbb{R}^{n_{x}+n_{u}}$ и ширин $\bm{\sigma}_{k}\in\mathbb{R}^{n_{x}+n_{u}}$ соответственно [12, 4] и $\bm{\mathcal{D}}:=\{(\bar{c}_{j},\bar{E}_{j})\}_{j=0}^{N-1}$. Формально, для всех $k\in[N]$ в (24d)–(24e), $\bar{c}_{k}\in\mathbb{R}^{n_{x}}$ и $\bar{E}_{k}\in\mathbb{R}^{n_{x}\times 2n_{x}}$ задаются как

$$\hskip-4.0pt\mathcal{A}_{k} =\left([z_{k}^{\top},v_{k}^{\top}]^{\top}+\bm{\mu}_{k}\right)\oplus\textrm{diag}(\bm{\sigma}_{k})\mathcal{E}_{n_{x}+n_{u}},$$ (25a)

$$\hskip-4.0pt\bar{E}_{k} :=[E(z_{k})\ \ \bar{R}_{k}]=\underbrace{\begin{bmatrix}E(z_{k})\ \ \Gamma_{r}\!\left(\mathcal{A}_{k},z_{k},v_{k};f,E\right)\end{bmatrix}}_{\mathcal{Z}_{r}(z_{k},v_{k},\bm{\mu}_{k},\bm{\sigma}_{k})},$$ (25b)

$$\hskip-4.0pt\bar{c}_{k} :=\underbrace{\Gamma_{c}\!\left(\mathcal{A}_{k},z_{k},v_{k};f,E\right)}_{\mathcal{Z}_{c}(z_{k},v_{k},\bm{\mu}_{k},\bm{\sigma}_{k})},$$ (25c)

где $\bar{c}_{k}$ — это смещение, а $\bar{E}_{k}$ — это порождающая матрица зонотопа для суммы Минковского внешнего возмущения и ошибки линеаризации: $E(z_{k})\mathcal{E}_{n_{x}}\oplus(\bar{c}_{k}\oplus\bar{R}_{k}\mathcal{E}_{n_{x}})\equiv\bar{c}_{k}\oplus\bar{E}_{k}\mathcal{E}_{2n_{x}}$. Используя (25b) и (25c), получаем аппроксимацию сверху достижимой трубки:,

$$\bar{\Omega}_{k}^{x}:=z_{k}+\textstyle\sum_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{c}_{j}+\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{E}_{j}\mathcal{E}_{2n_{x}},$$ (26)

$$\bar{\Omega}_{k}^{u}:=v_{k}+\textstyle\sum_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{u}}_{k,j}\bar{c}_{j}+\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{u}}_{k,j}\bar{E}_{j}\mathcal{E}_{2n_{x}}.$$

Устойчивое удовлетворение ограничений обеспечивается за счет ужесточения ограничений $h_{k}(\mathbf{\Phi},\bm{\mathcal{D}})$ и $h^{f}(\mathbf{\Phi},\bm{\mathcal{D}})$. Вкратце, эти ужесточения отражают распространение единично-нормализованных возмущений в (26), которые затем масштабируются границами возмущений $\bar{E}_{j}$ и смещаются на величину $\bar{c}_{j}$. Для простоты мы выводим ужесточения, предполагая линейные ограничения, т.е. $g_{k}(z_{k},v_{k})=G_{k}\begin{bmatrix}z_{k},\ v_{k}\end{bmatrix}^{\top}$ и $g^{f}(z_{N})=G^{f}z_{N}$, для $G_{k}\in\mathbb{R}^{n_{c}\times(n_{x}+n_{u})},G_{N}\in\mathbb{R}^{n_{c}\times n_{x}}$. Допустимые ужесточения могут быть получены для нелинейных $g_{k}$ и $g^{f}$ путем включения дополнительных членов ошибки [12, Eq. 22]. Формально мы определяем центр трубки (27) и ширину (28), объединенные (терминальные) ужесточения ограничений (29)-(30) и суррогатную стоимость отклика на системном уровне (31) как,

$$h_{k}^{\mu}(\mathbf{\Phi},\bm{\mathcal{D}}) =\textstyle\sum_{j=0}^{k-1}\mathbf{\Phi}_{k,j}\bar{c}_{j}$$ (27)

$$h_{k}^{\sigma}(\mathbf{\Phi},\bm{\mathcal{D}}) =\textstyle\sum_{j=0}^{k-1}\big\|\mathbf{\Phi}_{k,j}\bar{E}_{j}\big\|_{1,\mathrm{r}}$$ (28)

$$h_{k}(\mathbf{\Phi},\bm{\mathcal{D}}) =\textstyle\sum_{j=0}^{k-1}(G_{k}\mathbf{\Phi}_{k,j}\bar{c}_{j}+\big\|G_{k}\mathbf{\Phi}_{k,j}\bar{E}_{j}\big\|_{1,\mathrm{r}})$$ (29)

$$h^{f}(\mathbf{\Phi},\bm{\mathcal{D}}) =\textstyle\sum_{j=0}^{N-1}(G_{N}\mathbf{\Phi}^{\mathrm{x}}_{N,j}\bar{c}_{j}+\big\|G_{N}\mathbf{\Phi}^{\mathrm{x}}_{N,j}\bar{E}_{j}\big\|_{1,\mathrm{r}})\hskip-2.0pt$$ (30)

$$\tilde{H}_{0}(\mathbf{\Phi},\bm{\mathcal{D}}) =\textstyle\sum_{j=0}^{N-1}\!(\textstyle\sum_{k=j}^{N-1}\!(\|\tilde{Q}^{1/2}\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{E}_{j}\|_{\mathcal{F}}^{2}$$ (31)

$$\hskip 11.38109pt+\|\tilde{R}^{1/2}\mathbf{\Phi}^{\mathrm{u}}_{k,j}\bar{E}_{j}\|_{\mathcal{F}}^{2})+\|\tilde{Q}_{N}^{1/2}\mathbf{\Phi}^{\mathrm{x}}_{N,j}\bar{E}_{j}\|_{\mathcal{F}}^{2})\,,$$

где $\small{\mathbf{\Phi}_{k,j}:=\begin{bmatrix}{\mathbf{\Phi}^{\mathrm{x}}_{k,j}}^{\top}&\hskip-5.69054pt{\mathbf{\Phi}^{\mathrm{u}}_{k,j}}^{\top}\end{bmatrix}^{\top}}$, $\mathbf{\Phi}$ собирает $\mathbf{\Phi}_{k,j}$ для всех $k,j\in[N]$, и $\tilde{Q},\tilde{Q}_{N}\in\mathbb{S}_{++}^{n_{x}},\tilde{R}\in\mathbb{S}_{++}^{n_{u}}$. Сужения ограничений и ширины трубок в (28)-(30) выведены для наихудшего случая возмущений с использованием свойства двойственной нормы [28]. Функционал стоимости $\tilde{H}_{0}(\mathbf{\Phi},\bm{\mathcal{D}})$ (31) штрафует неопределенность (ширины трубок (28)). Следуя [9, 10], мы решаем (24) посредством итерационного процесса, который чередуется между поиском (A) номинальной траектории и (B) робастного контроллера. Чтобы решить (A), мы находим $(\mathbf{z},\mathbf{v})$ путем решения NOCP с суженными ограничениями (32⟧),

$$\hskip-5.0pt\min_{\begin{subarray}{c}\mathbf{z},\,\mathbf{v}\end{subarray}}\quad J(\mathbf{z},\mathbf{v})$$ (32a)

s.t. $$z_{k+1}=f(z_{k},v_{k}),\quad\forall k\in[N],\qquad z_{0}=\bar{x}_{0},$$ (32b)

$$g_{k}(z_{k},v_{k})+h_{k}(\mathbf{\Phi},\bm{\mathcal{D}})\leq 0,\quad\forall k\in[N],$$ (32c)

$$g^{f}(z_{N})+h^{f}(\mathbf{\Phi},\bm{\mathcal{D}})\leq 0.$$ (32d)

Затем (32) повторяет структуру (2) и может быть решена с использованием метода параллельных вычислений на GPU из [10, Sec. IV]. После решения (A) (32) мы проводим линеаризацию вокруг полученного $(\mathbf{z},\mathbf{v})$ и находим робастный контроллер (B) путем поиска $\mathbf{\Phi}$, который оптимизирует (31):

$$\hskip-3.0pt\min_{\mathbf{\Phi}^{\mathrm{x}},\mathbf{\Phi}^{\mathrm{u}}}\ \sum_{j=0}^{N-1}\sum_{k=j}^{N-1}\|\mathcal{Q}_{k,j}\mathbf{\Phi}_{k,j}\bar{E}_{j}\|^{2}_{\mathcal{F}}+\|\mathcal{Q}_{N,j}\mathbf{\Phi}^{\mathrm{x}}_{N,j}\bar{E}_{j}\|^{2}_{\mathcal{F}}$$ (33a)

s.t. $$\mathbf{\Phi}^{\mathrm{x}}_{k+1,j}=A_{k}\mathbf{\Phi}^{\mathrm{x}}_{k,j}+B_{k}\mathbf{\Phi}^{\mathrm{u}}_{k,j},$$ (33b)

$$\mathbf{\Phi}^{\mathrm{x}}_{j+1,j}=I_{n_{x}},~\forall j\in[N],~\forall k\in[j+1,N-1],\hskip-4.0pt$$ (33c)

где $\mathcal{Q}_{k,j},\mathcal{Q}_{N,j}$ определены в соответствии с [9, Eq. 20, 23] (см. прил. E) для обеспечения согласованности между (32) и (33). Примечательно, что умножение справа на $\bar{E}_{j}$ препятствует представлению функции стоимости в (33a) в виде стандартной целевой функции LQR в $\mathbf{\Phi}^{x}$ и $\mathbf{\Phi}^{u}$. Таким образом, она не может быть решена напрямую с использованием методов на основе уравнения Риккати из [9, 10]. Вместо этого мы решаем (34), которая имеет стоимость LQR и может быть решена с помощью эффективных решателей [9, 10]⟧:

$$\hskip-8.0pt\min_{\mathbf{\Phi}^{\mathrm{x}},\mathbf{\Phi}^{\mathrm{u}}}\ \sum_{j=0}^{N-1}\hskip-3.0pt\sum_{k=j}^{N-1}\|\mathcal{Q}_{k,j}\mathbf{\Phi}_{k,j}\|^{2}_{\mathcal{F}}+\|\mathcal{Q}_{N,j}\mathbf{\Phi}^{\mathrm{x}}_{N,j}\|^{2}_{\mathcal{F}}$$ (34a)

s.t. $$\mathbf{\Phi}^{\mathrm{x}}_{k+1,j}=A_{k}\mathbf{\Phi}^{\mathrm{x}}_{k,j}+B_{k}\mathbf{\Phi}^{\mathrm{u}}_{k,j},$$ (34b)

$$\mathbf{\Phi}^{\mathrm{x}}_{j+1,j}=I_{n_{x}},~\forall j\in[N],~\forall k\in[j+1,N-1].$$ (34c)

Заметим, что (34) совпадает с (33), если $\bar{E}_{j}=I$ для всех $j\in[N]$. В предложении 1 мы доказываем, что множества минимизаторов (33) и (34) эквивалентны, если $\bar{E}_{j}$ является правообратимым для всех $j\in[N]$ (см. прил. C для доказательства).

Proposition 1

Пусть $(\mathbf{\tilde{\Phi}}^{\mathrm{x}},\mathbf{\tilde{\Phi}}^{\mathrm{u}})$ является оптимизатором (33), а $(\mathbf{\hat{\Phi}}^{\mathrm{x}},\mathbf{\hat{\Phi}}^{\mathrm{u}})$ является оптимизатором (34). Если $\bar{E}_{j}$ является правообратимым для всех $j\in[N]$, то есть существует $\bar{E}_{j}^{\ddagger}$ такое, что $\bar{E}_{j}\bar{E}_{j}^{\ddagger}=I_{n_{x}}$, тогда два оптимизатора совпадают: $\mathbf{\tilde{\Phi}}^{\mathrm{x}}=\mathbf{\hat{\Phi}}^{\mathrm{x}},\mathbf{\tilde{\Phi}}^{\mathrm{u}}=\mathbf{\hat{\Phi}}^{\mathrm{u}}$.

При реализации метода из разд. IV мы добавляем небольшое положительное диагональное дополнение к $\bar{R}_{k}$, чтобы она всегда была обратимой, и, таким образом, $\bar{E}_{j}$ (25) всегда была правообратимой для всех $j\in[N]$. После решения (34) мы вычисляем члены ужесточения ограничений в (27)–(30). Это обновляет центры достижимых трубок $\bm{\mu}_{k}$ и ширины $\bm{\sigma}_{k}$: $\bm{\sigma}_{k}\leftarrow h_{k}^{\sigma}(\mathbf{\Phi},\bm{\mathcal{D}})~\text{and}~\bm{\mu}_{k}\leftarrow h_{k}^{\mu}(\mathbf{\Phi},\bm{\mathcal{D}}),$ которые используются в качестве ужесточений ограничений на следующей итерации оптимизации номинальной траектории (32). Этот процесс повторяется до сходимости (как в [9, Alg. 1]); мы завершаем работу после фиксированного числа итераций или если изменение итератов оказывается ниже порогового значения (см. прил. B для алгоритмического блока).

Новизна формулировки: Здесь мы суммируем изменения, внесенные в предыдущие формулировки нелинейного SLS [9, App. A] [4]⟧:

1. ⟦ Зонотопы позволяют выполнять точное вычисление сужения (27)–(30) путем конкатенации генераторов внешнего возмущения и ошибки линеаризации, в отличие от предыдущих эллипсоидальных методов [9, App. A], однако это приводит к правообратимым $\bar{E}_{j}$, несовместимым с существующими решателями [9, 10]. Мы решаем эту проблему путем установки $\mathbf{\Phi}^{\mathrm{x}}_{j+1,j}=I_{n_{x}}$ и включения масштабирования возмущения непосредственно в сужения и функцию стоимости (31) с помощью совместимой переформулировки (34), обоснованной предл. 1⟧.
2. ⟦ Наконец, мы расширяем [4, 9] для обработки возмущений с ненулевым средним значением, что позволяет напрямую использовать интервальные границы из (23) без дополнительного консерватизма, связанного с центрированием относительно нуля..

V-B Эффективная реализация SLS на GPU

Мы предлагаем GPUSLS-Linearization Error Optimization (GPUSLS-LEO), который опирается на GPUSLS [10] для явного учета распространения ошибки линеаризации. Мы сначала заменяем эллипсоидальное представление возмущений из [10] на зонотоп, модифицируя сужения ограничений и стоимостные слагаемые посредством (29)-(31). При каждом обновлении контроллера мы модифицируем GPUSLS, включая зонотопическую границу ошибки линеаризации посредством (25). Затем мы модифицируем инициализацию $\mathbf{\Phi}^{\text{x}}$ посредством (34c).

Все остальные шаги в точности следуют процедуре GPUSLS. В частности, GPUSLS-LEO решает (24) итеративно посредством последовательного квадратичного программирования (SQP). На каждой итерации GPUSLS-LEO аппроксимирует стоимость квадратичной функцией и линеаризует $x_{k+1}=f(x_{k},u_{k})$ вокруг текущего $(\textbf{z},\textbf{v})$ для получения $A_{k},B_{k}$ (4). Затем он решает (A) локальную QP-аппроксимацию (32) для обновления $(\textbf{z},\textbf{v})$, (B) (34) для получения $\mathbf{\Phi}$, и обновляет (C) сужения посредством (24f). Поскольку наши LEB дифференцируемы, мы можем штрафовать накопление ошибки во время шага (A) путем модификации номинальной целевой функции (32a) до $\bar{J}(\textbf{z},\textbf{v}):=J(\textbf{z},\textbf{v})+\lambda J_{\textrm{err}}(\textbf{z},\textbf{v})$, где

$$\hskip-2.0ptJ_{\textrm{err}}(\textbf{z},\textbf{v})=\sum_{j=0}^{N-1}\sum_{k=j}^{N-1}\|\mathbf{\Phi}_{k,j}\mathcal{Z}_{r}(z_{j},v_{j},\bm{\mu}_{j},\bm{\sigma}_{j})\|_{\mathcal{F}}^{2},$$ (35)

штрафует распространяемую остаточную ошибку (8a). Здесь $\mathbf{\Phi},\bm{\mu}_{j},\bm{\sigma}_{j}$ взяты из предыдущей итерации, а $\lambda$ является весовым коэффициентом. Примечательно, что (35) смещает $(\textbf{z},\textbf{v})$ в сторону областей с низкими возмущениями от распространяемых экзогенных ошибок и ошибок линеаризации. Наконец, мы показываем, что путем решения (24), GPUSLS-LEO обеспечивает удержание в достижимых трубках (26), т.е. $x_{k}\in\bar{\Omega}_{k}^{x}$ и $u_{k}\in\bar{\Omega}_{k}^{u}$, для всех $k\in[N]$. Формально мы имеем (доказательство в приложении D).⟧):

Theorem 4

Пусть $\mathbf{z}$, $\mathbf{v}$ и $\mathbf{\Phi}$ являются допустимым решением для (24). Тогда (26) гарантированно дает верхнюю аппроксимацию истинного достижимого множества замкнутой системы (3), то есть $\Omega_{k}\subseteq\bar{\Omega}_{k}$ для всех $k\in[N]$.

VI Результаты

Мы оцениваем точность LEB (разд. VI-A) и робастные контроллеры, полученные с помощью GPUSLS-LEO (разд. VI-B). Тесты времени выполнения используют NVIDIA RTX 4090, за исключением экспериментов с несколькими квадрокоптерами и длинным горизонтом, для которых требуется NVIDIA H200. Определения систем приведены в прил. A⟧.

VI-A Ошибка линеаризации

Мы сравниваем LEB, полученную с помощью нашего метода на основе траекторий Hessian (PBH), для систем спутника ($n_{x}=7$, $n_{u}=3$) и квадрокоптера ($n_{x}=12$, $n_{u}=4$) (прил. A) с несколькими базовыми методами. Глобальная выборка [4] оценивает кривизну наихудшего случая в автономном режиме (занимая более 7 минут [4]) по глобальному множеству $\mathcal{X}\times\mathcal{U}$, что дает эмпирическую аппроксимацию сверху. Интервальный Hessian (IH) (классический) [6, 7] вычисляет (12) на локальном множестве $\mathcal{T}\subset(\mathcal{X}\times\mathcal{U})$. IH (классический)-CORA обозначает стандартную реализацию IH в CORA [6, 29]. Случайная выборка использует 10 тысяч образцов из локального множества и дает аппроксимацию снизу. Мы выбираем номинальные точки $(z,v)$ из $\mathcal{X}\times\mathcal{U}$ и определяем локальные множества $\mathcal{T}:=(z,v)\oplus\epsilon\mathcal{E}_{n_{x}+n_{u}}$ с помощью $\epsilon=0.01,0.05,0.1,0.15,0.2$. Рис. 2 показывает среднюю ширину интервала (23). Для спутника, где нелинейность умеренная, все методы аппроксимации сверху остаются близкими к нижней границе случайной выборки по мере роста $\epsilon$, хотя PBH все равно является наиболее точным. Для квадрокоптера PBH неизменно дает наиболее точные границы. Напротив, глобальная выборка дает неточные результаты при малых $\epsilon$, что говорит о том, что глобальная кривизна не учитывает локальное поведение, в то время как интервальный Hessian становится чрезмерно консервативным при $\epsilon=0.20$. Таблица I сравнивает время выполнения в режиме онлайн. PBH выполняется менее чем за 1 мс, что позволяет использовать его в реальном времени. Несмотря на то, что он медленнее глобальной границы, которая эффективно вычисляется онлайн после медленных автономных вычислений, PBH на 60% и 27% (для спутника и квадрокоптера соответственно) быстрее плотной случайной выборки и имеет пренебрежимо малое время выполнения по сравнению с CORA IH.

⟦Fig. 2: (а): ⟦BROKEN_K1⟧ для спутниковой системы. (б): ⟦BROKEN_K1⟧ для системы квадрокоптера. В обеих системах наш метод достигает наиболее точной аппроксимации сверху..

⟦TABLE I: Сравнение времени выполнения (в мс) методов границ линеаризации.

System	Random	Global	IH-classic	IH-classic-CORA	Path-based
Satellite	0.894	0.143	0.218	11570	0.361
Quadrotor	1.278	0.156	0.453	227200	0.936

⟦TABLE II: Ошибка линеаризации при номинальной оптимизации с помощью PGD.

System	Initial	Sampling	PGD	PGD better (%)
Satellite	0.2680	0.2365	0.2241	100
Quadrotor	0.02253	0.01639	0.01590	100

Чтобы количественно оценить преимущество дифференцируемости в нашей реализации Гессиана на основе путей, мы ищем номинальную точку в пределах фиксированного радиуса возмущения $\epsilon$, которая минимизирует размер LEB. Начиная со случайных номинальных точек, мы запускаем проекционный градиентный спуск (PGD) в пределах локальной окрестности радиуса $\epsilon_{\text{opt}}$, чтобы уменьшить ширину интервала ошибки. Мы сравниваем PGD со случайной выборкой номинальных точек при том же временном бюджете, где оба метода используют PBH для оценки ошибки линеаризации. Таблица II показывает, что в 100% тестовых случаев PGD дает меньшую ширину интервала ошибки, чем случайная выборка для обеих систем, что демонстрирует ценность градиентной оптимизации для поиска номинальных точек с более узкими LEB, что может быть полезно для последующего SLS..

VI-B Использование границ линеаризации для RNOCP

В этом разделе мы оцениваем GPUSLS-LEO на наборе аналитических и нейронных динамических систем в сравнении с различными базовыми методами RNOCP. Мы дополнительно оцениваем влияние ключевых абляций, включая эффект оптимизации по градиентам ошибки линеаризации и использование зонотопов с ненулевым центром. Все сравниваемые методы оцениваются при одинаковом наборе возмущений для заданной динамической системы и, где это применимо, запускаются максимум на 100 итераций SQP..

⟦Fig. 3: Спутник (7D). (а): Размеры трубок для нашего метода, GPUSLS-LEO, в сравнении с базовыми методами и вариантом без ошибки линеаризации. Наш метод является наименее консервативным при минимальном увеличении размера трубки. (b): Трубки для первых четырех размерностей состояния; наши трубки охватывают 100% выборочных прогонов..

VI-B1 Спутник

Мы оцениваем наш метод на задаче управления ориентацией 7D спутника, где целью является стабилизация системы в желаемой ориентации при наличии ограниченных возмущений. Мы используем экспериментальную установку из [4] и проводим сравнение с их подходом. На рис. 3 (a) мы сравниваем ширины трубок $\bm{\sigma}_{k}$ для: GPUSLS-LEO, GPUSLS-E (эллипсоидальное распространение, как в [9, App. A]), GPUSLS (без распространения ошибки линеаризации) [10] и NL-SLS ( [4]). Наш метод позволяет получить наиболее узкие трубки среди всех базовых подходов. В среднем наши трубки на 20% меньше, чем у NL-SLS, который опирается на глобальную максимальную границу Гессиана и поэтому является локально чрезмерно консервативным вдоль траектории. По сравнению с GPUSLS-E наш метод позволяет получить трубки на 45% меньше, так как он включает аппроксимацию суммы Минковского для эллипсоидов сверху, в то время как наша комбинация на основе зонотопов является точной. Мы наблюдаем, что внедрение нашего конвейера LEB не вносит существенного консерватизма, поскольку ширина трубок увеличивается лишь на 4% по сравнению с GPUSLS. На рис. 3 (b) мы показываем результаты 20 случайных и 64 состязательных возмущений для состояний кватернионов, демонстрируя 100% удержание в рамках нашего метода..

⟦Fig. 4: Планарный квадрокоптер (6D). (a): Сравнение площади трубки с CCM и immrax, показывающее, что наш метод менее консервативен. (b): Трубки и траектории для GPUSLS-LEO и CCM. Оба метода обеспечивают 100% безопасность, но консервативность CCM вынуждает систему выбирать субоптимальный путь.

⟦Fig. 5: Планарный квадрокоптер (6D). (a): Система, перемещающаяся через поле препятствий по траектории длиной 20 м с горизонтом планирования 7500. (b): Ширина трубок вдоль горизонта. Несмотря на длинный горизонт, трубки остаются ограниченными благодаря оптимизации размеров трубок контроллером GPUSLS-LEO..

VI-B2 Планарный квадрокоптер

Мы оцениваем наш метод на 6D планарном квадрокоптере в сравнении с метриками сжатия управления (CCM) [30]. Оба метода должны перемещаться через плотное поле препятствий на горизонте в 40 шагов. Рис. 4(a) показывает, что наш метод достигает на 91% меньшей площади трубки, чем CCM для позиционных состояний $(x,y)$. Это происходит потому, что CCM подтверждает безопасность во всем пространстве состояний, что вносит присущую консервативность. GPUSLS-LEO, с другой стороны, подтверждает безопасность вдоль конкретной траектории, обеспечивая более узкие, зависящие от траектории границы. Мы также оцениваем замкнутую систему с помощью инструмента интервального анализа достижимости immrax [24], показывая, что рассчитанная им площадь трубки расходится через несколько шагов..

На рис. 4(b) мы сравниваем рассчитанные траектории и трубки для нашего метода и CCM. Наш подход менее консервативен, что позволяет выбрать более прямой, оптимальный для задачи путь, в то время как CCM вынужден выбирать более длинный маршрут. На рис. 5(a) мы показываем траекторию и трубки на длинном горизонте в 7500 шагов, что включает оптимизацию траектории состояния-управления с $1.95\times 10^{5}$ переменными решения. Благодаря точности наших LEB и оптимизации контроллера с обратной связью в SLS, наши трубки остаются хорошо ограниченными и конечными на протяжении всего горизонта (рис. 5(b)), что демонстрирует масштабируемость метода для длинных горизонтов..

⟦Fig. 6: Система из 14 связанных трехмерных квадрокоптеров (168D). (a): Вычисленные траектории и трубки в поле препятствий. Даже при высокой размерности наш метод способен вычислять точные трубки. (b): Трубки и 100 запусков для состояний $p_{y}$ и $v_{y}$ для 4 квадрокоптеров. 100% запусков безопасно остаются внутри трубок.

VI-B3 Связанный трехмерный квадрокоптер

Чтобы продемонстрировать масштабируемость для высокоразмерных систем, мы оцениваем наш метод на связанной системе из 14 трехмерных квадрокоптеров (168D) на горизонте 25. Перед системой стоит задача навигации через поле препятствий, при этом квадрокоптеры связаны силой пружины, направленной к центроиду группы. Мы показываем запланированные траектории и трубки на рис. 6(a) и наблюдаем, что все реализации остаются внутри узких трубок, несмотря на высокоразмерное пространство состояний, что демонстрирует масштабируемость метода по отношению к размерности состояния. На рис. 6(b) мы строим графики трубок и соответствующие им реализации для $p_{y}$ и $v_{y}$ для 100 состязательных возмущений.

⟦Fig. 7: Нейронный T-толкатель (5D). (a): Сравнение размеров трубки для различных абляций нашего метода, показывающее, что трубки с ненулевым центром и градиенты ошибки линеаризации уменьшают ширину трубок. (b): Прогоны нашего метода; все траектории остаются внутри вычисленных робастных трубок. (c): Прогоны GPUSLS. Без формального учета ошибки линеаризации прогоны выходят за пределы трубки.

VI-B4 Нейронный T-толкатель

Мы демонстрируем, что наша формулировка позволяет учитывать распространение ошибки линеаризации через обученную нейронную динамику. Мы обучаем модель на основе MLP для системы T-pusher и ставим перед ней задачу вращения и перемещения объекта в целевую конфигурацию (Рис. 1). Мы оцениваем наш метод на нескольких абляциях (Рис. 7(a)), включая GPUSLS, GSL (сокращение для GPUSLS-LEO), GSL-C (центрированные в нуле зонотопы) и GSL ($\lambda=1$), который включает штраф за ошибку линеаризации, как определено в V-B. По сравнению с GPUSLS, GPUSLS-LEO создает трубки, которые в среднем на 99% больше. Однако Рис. 7(b) показывает, что наш метод достигает 100% охвата всех развертываний, в то время как Рис. 7(c) иллюстрирует, как траектории из GPUSLS выходят за пределы трубки, демонстрируя необходимость формального учета ошибки линеаризации..

На рис. 7(a) показано, что добавление штрафов за ошибку линеаризации позволяет снизить консервативность, при этом $\lambda=1$ уменьшает ширину трубок на 8%. Мы также демонстрируем важность зонотопов с ненулевым центром, поскольку использование зонотопов с центром в нуле приводит к увеличению трубок на 25%. Мы представляем результаты прогонов и отклонения от возмущений на рис. 1(a). На рис. 1(b) мы демонстрируем возможности метода для замкнутого MPC на симулированной динамике контактов, где толкатель успешно выполняет задачу. Отметим, что мы формально не учитываем ошибку обучения; однако она может быть включена, как в [31, 32].

⟦TABLE III: Распределение времени выполнения (мс) на итерацию по системам.

System	GPUSLS	Lin. Error Grad.	Remainder	Total	DT
Satellite	11	3	$<1$	15	1000
Planar Quadrotor	18	1	$<1$	20	150
Coupled Quadrotor (96D)	40	–	48	88	100
Neural T-Dynamics	630	100	32	762	1000

⟦Fig. 8: Связанный квадрокоптер. Оценка времени выполнения каждой итерации при увеличении количества связанных квадрокоптеров. GPUSLS демонстрирует благоприятное масштабирование в зависимости от размерности состояния, увеличиваясь приблизительно логарифмически относительно пространства состояний и управления. Метод Path-Based Hessians (PBH) обеспечивает низкое время выполнения при малых размерностях, но масштабируется значительно хуже по мере увеличения размерности состояния..

VI-B5 Время выполнения

Мы оцениваем время выполнения нашего метода для всех систем. В табл. III мы приводим распределение времени выполнения различных компонентов нашей методологии на итерацию в миллисекундах. Во всех системах, за исключением связанного квадрокоптера, основная вычислительная нагрузка приходится на GPUSLS, что показывает, что предложенные нами модификации не приводят к резкому увеличению времени выполнения. Мы также показываем, что все общие времена выполнения меньше шага дискретизации DT, что демонстрирует жизнеспособность нашего метода в рамках схемы MPC с итерациями в реальном времени [10].

Наконец, мы оцениваем масштабируемость метода при увеличении размерностей состояния и управления с использованием возрастающего числа связанных квадрокоптеров. Рис. 8 показывает, что GPUSLS требует большего начального времени выполнения, но благоприятно масштабируется с размерностью состояния и управления благодаря логарифмической зависимости [10]. PBH, с другой стороны, демонстрирует низкое начальное время выполнения, но может замедляться в более высоких размерностях из-за стоимости вычислений Гессиана..

VII Заключение

В данной статье представлена структура для вычисления точных, дифференцируемых и корректных границ ошибок линеаризации (LEBs) для нелинейных аналитических и NN динамических систем, а также их использования в GPUSLS-LEO — решателе RNOCP с параллелизацией на GPU в реальном времени, основанном на SLS. Данный подход снижает консервативность, превосходя современные инструменты анализа достижимости, при этом сохраняя формальные гарантии и обеспечивая производительность на системах с количеством состояний до 168, частотой управления в реальном времени (до 67 Гц) и задачах с $\approx 2\times 10^{5}$ переменными решения. В будущих работах будет исследовано, как наши LEBs могут быть расширены для ограничения ошибок линеаризации при планировании с учетом контактов [28] и нелинейном управлении с обратной связью по выходу, в том числе посредством расширений методов на основе SLS, описанных в [33, 34, 35].

В приложении мы сначала обсуждаем детали динамики систем, использованной в экспериментальных результатах (App. A). Затем мы приводим блок алгоритма, описывающий, как последовательное квадратичное программирование используется для решения GPUSLS-LEO (App. B). Далее мы обсуждаем и доказываем Proposition 1 и вспомогательную лемму в App. C. Затем мы приводим доказательства для остальных теоретических результатов, которые были опущены в основной рукописи в App. D. Наконец, мы более подробно обсуждаем переменную множителя Лагранжа $\tau$, используемую в нашей структуре в App. E⟧.

Appendix A Определения динамических систем

В этом разделе мы приводим обзор определений динамических систем с дискретным временем для каждой системы. Динамика с дискретным временем получена посредством дискретизации методом Эйлера вперед с шагом по времени $\Delta t$, $x_{k+1}=x_{k}+\Delta t(\dot{x}_{k})$.

A-A Спутник

Мы используем ту же динамику спутника и экспериментальную установку, что и [4]. В частности, мы определяем динамику как,

$$\dot{x} =\begin{bmatrix}\zeta(\omega)q\\ I_{S}^{-1}(v-w\times(I_{S}\omega))\end{bmatrix}$$ (36a)

$$\zeta(\omega) :=\frac{1}{2}\begin{bmatrix}0&-\omega_{1}&-\omega_{2}&-\omega_{3}\\ \omega_{1}&0&\omega_{3}&-\omega_{2}\\ \omega_{2}&-\omega_{3}&0&\omega_{1}\\ \omega_{3}&\omega_{2}&-\omega_{1}&0\end{bmatrix}$$ (36b)

с состояниями $x:=(q,\omega)\in\mathbb{R}^{7}$, где $q\in\mathbb{R}^{4}$ — кватернион ориентации, $\omega\in\mathbb{R}^{3}$ — угловая скорость вращения, $v\in\mathbb{R}^{3}$ — входной управляющий момент, а $I_{S}=\mathrm{diag}(5,2,1)$ — симметричная матрица инерции. Мы определяем постоянное масштабирование возмущения $E=5\cdot 10^{-3}\cdot\mathrm{diag}(0,0,0,0,1,1,1)$.

A-B Плоский квадрокоптер

Мы рассматриваем планарный квадрокоптер со следующей динамикой

$$\dot{x} =\begin{bmatrix}v_{x}\\ v_{y}\\ \dot{\phi}\\ -\frac{1}{m}(u_{1}+u_{2})\sin(\phi)\\ \frac{1}{m}(u_{1}+u_{2})\cos(\phi)-g\\ \frac{L}{J}(u_{2}-u_{1})\end{bmatrix}$$ (37a)

Состояние определяется как $x:=(p_{x},p_{y},\phi,v_{x},v_{y},\dot{\phi})\in\mathbb{R}^{6}$, а вход как $u:=(u_{1},u_{2})\in\mathbb{R}^{2}$. Положение обозначается как $(p_{x},p_{y})$, $\phi$ — угол тангажа, $(v_{x},v_{y})$ — поступательные скорости, а $\dot{\phi}$ — угловая скорость. Входы $(u_{1},u_{2})$ представляют собой тягу отдельных роторов. Мы определяем массу $m=2.0576$, ускорение свободного падения $g=9.81$, длину плеча $L=0.25$ и момент инерции $J=0.01$. Мы определяем масштабирование возмущения как $E=5\cdot 10^{-2}\cdot\text{diag}(0,0,0,1,1,0)$.

A-C Связанный трехмерный квадрокоптер

Мы рассматриваем систему из $N$ квадрокоптеров. Для заданного квадрокоптера $i$ мы определяем его динамику как

$$\dot{\xi}_{i} =\begin{bmatrix}\dot{p}_{x,i}\\ \dot{p}_{y,i}\\ \dot{p}_{z,i}\\ \dot{\phi}_{i}\\ \dot{\theta}_{i}\\ \dot{\psi}_{i}\\ \dot{v}_{x,i}\\ \dot{v}_{y,i}\\ \dot{v}_{z,i}\\ \dot{p}_{i}\\ \dot{q}_{i}\\ \dot{r}_{i}\end{bmatrix}=\begin{bmatrix}v_{x,i}\\ v_{y,i}\\ v_{z,i}\\ p_{i}+(q_{i}\sin\phi_{i}+r_{i}\cos\phi_{i})\tan\theta_{i}\\ q_{i}\cos\phi_{i}-r_{i}\sin\phi_{i}\\ \frac{q_{i}\sin\phi_{i}+r_{i}\cos\phi_{i}}{\cos\theta_{i}}\\ \frac{T_{i}}{m}\left(\cos\psi_{i}\sin\theta_{i}\cos\phi_{i}+\sin\psi_{i}\sin\phi_{i}\right)\\ \frac{T_{i}}{m}\left(\sin\psi_{i}\sin\theta_{i}\cos\phi_{i}-\cos\psi_{i}\sin\phi_{i}\right)\\ \frac{T_{i}}{m}\left(\cos\theta_{i}\cos\phi_{i}\right)-g\\ \frac{J_{y}-J_{z}}{J_{x}}q_{i}r_{i}+\frac{\tau_{\phi,i}}{J_{x}}\\ \frac{J_{z}-J_{x}}{J_{y}}p_{i}r_{i}+\frac{\tau_{\theta,i}}{J_{y}}\\ \frac{J_{x}-J_{y}}{J_{z}}p_{i}q_{i}+\frac{\tau_{\psi,i}}{J_{z}}\end{bmatrix}.$$ (38a)

Состояние квадрокоптера $i$ определяется как $\xi_{i}:=(p_{x,i},p_{y,i},p_{z,i},\phi_{i},\theta_{i},\psi_{i},v_{x,i},v_{y,i},v_{z,i},p_{i},q_{i},r_{i})\in\mathbb{R}^{12}$, где $(p_{x,i},p_{y,i},p_{z,i})$ обозначает положение, $(\phi_{i},\theta_{i},\psi_{i})$ — углы крена, тангажа и рыскания, $(v_{x,i},v_{y,i},v_{z,i})$ — поступательные скорости, а $(p_{i},q_{i},r_{i})$ обозначают угловые скорости корпуса. Входные данные определяются как $\upsilon_{i}:=(T_{i},\tau_{\phi,i},\tau_{\theta,i},\tau_{\psi,i})\in\mathbb{R}^{4}$, где $T_{i}$ — суммарная тяга, а $(\tau_{\phi,i},\tau_{\theta,i},\tau_{\psi,i})$ — моменты сил, действующие на корпус..

Для системы из $N$ связанных квадрокоптеров мы объединяем состояния и входные воздействия, получая общие векторы состояния и управления,

$$x:=(\xi_{1},\ldots,\xi_{N})\in\mathbb{R}^{12N},\quad u:=(\upsilon_{1},\ldots,\upsilon_{N})\in\mathbb{R}^{4N}.$$ (39)

Мы дополнительно включаем центроидную пружинно-демпферную связь в поступательную динамику. Мы определяем положение и скорость центроида как

$$\bar{p}_{k}=\frac{1}{N}\sum_{i=1}^{N}p_{i,k},\quad\bar{v}_{k}=\frac{1}{N}\sum_{i=1}^{N}v_{i,k}.$$ (40)

Ускорение связи, примененное к квадрокоптеру $i$, составляет

$$a_{i,k}^{\text{coup}}=-\alpha N(p_{i,k}-\bar{p}_{k})-\gamma(v_{i,k}-\bar{v}_{k})$$ (41)

где $\alpha$ и $\gamma$ представляют собой коэффициенты связи, при этом $\alpha$ определяет жесткость пружинного члена, а $\gamma$ управляет коэффициентом демпфирования. В наших экспериментах мы задаем $\alpha=0.75$ и $\gamma=0.25$. Таким образом, динамика в дискретном времени принимает вид

$$\xi_{i,k+1}=\xi_{i,k}+\Delta t\begin{bmatrix}\dot{p}_{x,i,k}\\ \dot{p}_{y,i,k}\\ \dot{p}_{z,i,k}\\ \dot{\phi}_{i,k}\\ \dot{\theta}_{i,k}\\ \dot{\psi}_{i,k}\\ \dot{v}_{i,k}+a^{\mathrm{coup}}_{i,k}\\ \dot{p}_{i,k}\\ \dot{q}_{i,k}\\ \dot{r}_{i,k}\end{bmatrix}.$$ (42)

Мы используем матрицу масштабирования возмущений $E=2\cdot 10^{-2}\cdot I_{12N}$

Appendix B Блок алгоритма

⟦Algorithm 1 GPUSLS-LEO

1: $\bar{x}_{0}$, $J$, $f(x,u)$, $E(x)$, $g_{k}(x,u)$, $g^{f}(x)$

2: в то время как SQP не сошелся выполнять

3:   $\mathbf{A}$, $\mathbf{B}$, $\mathbf{G}$ $\leftarrow\textsc{SQP: Linearize System}$

4:   $\mathbf{Q},\mathbf{M},\mathbf{R},\mathbf{q},\mathbf{r}\leftarrow\textsc{SQP: Quadraticize Cost}$

5:   в то время как SLS не сошелся выполнять

6:    $\mathbf{z}$, $\mathbf{v}$ $\leftarrow\textsc{Optimize Nominal Trajectory}$

7:    $\bm{\sigma},\bm{\mu}\leftarrow\textsc{Compute interval}$

8:    $\mathbf{\bar{E}}$, $\mathbf{\bar{c}}$ $\leftarrow\textsc{Compute residual error zonotope}$

9:    $\bm{\tau},\bm{\mathcal{Q}}\leftarrow\textsc{Update duals and cost}$

10:    $\mathbf{\Phi}^{\textrm{x}},\mathbf{\Phi}^{\textrm{u}}\leftarrow\textsc{Optimize Controller}$

11:    $h(\mathbf{\Phi},\bm{\mathcal{D}}),h^{f}(\mathbf{\Phi},\bm{\mathcal{D}})\leftarrow\textsc{Update Tightenings}$

12:   конец в то время как

13: конец в то время как

14: вернуть $\mathbf{z},\mathbf{v},\bm{\Phi}^{x},\bm{\Phi}^{u}$

Здесь $\mathbf{Q},\mathbf{M},\mathbf{R},\mathbf{q},\mathbf{r}$ соответствует квадратичной аппроксимации номинальной стоимости.

Appendix C О решении для $\mathbf{\Phi}$

Решение (33) дает оптимальные значения для $\mathbf{\Phi}^{\mathrm{x}}$ и $\mathbf{\Phi}^{\mathrm{u}}$, которые штрафуют размер сужений (31). Однако это приводит к появлению в (33a) членов, которые масштабируют матрицы отклика $\mathbf{\Phi}^{\mathrm{x}}$ и $\mathbf{\Phi}^{\mathrm{u}}$ на $\bar{E}_{j}$ в функции стоимости (33a). Это по-прежнему квадратичная стоимость относительно $\mathbf{\Phi}$, но она больше не сохраняет структуру функции стоимости LQR. Чтобы решить эту проблему, мы вместо этого решаем (34), которая имеет структуру LQR и может быть решена с помощью рекурсий Риккати. Для полноты изложения мы приводим исходную задачу (33), переформулированную задачу (34) и эквивалентное утверждение оптимизатора (Prop. 1) ниже.:

Исходная задача:

$$\hskip-3.0pt\min_{\mathbf{\Phi}^{\mathrm{x}},\mathbf{\Phi}^{\mathrm{u}}}\ \sum_{j=0}^{N-1}\sum_{k=j}^{N-1}\|\mathcal{Q}_{k,j}\mathbf{\Phi}_{k,j}\bar{E}_{j}\|^{2}_{\mathcal{F}}+\|\mathcal{Q}_{N,j}\mathbf{\Phi}^{\mathrm{x}}_{N,j}\bar{E}_{j}\|^{2}_{\mathcal{F}}$$ (33a)

s.t. $$\mathbf{\Phi}^{\mathrm{x}}_{k+1,j}=A_{k}\mathbf{\Phi}^{\mathrm{x}}_{k,j}+B_{k}\mathbf{\Phi}^{\mathrm{u}}_{k,j},$$ (33b)

$$\mathbf{\Phi}^{\mathrm{x}}_{j+1,j}=I_{n_{x}},~\forall j\in[N],~\forall k\in[j+1,N-1],\hskip-4.0pt$$ (33c)

Переформулированная задача:

$$\hskip-5.0pt\min_{\mathbf{\Phi}^{\mathrm{x}},\mathbf{\Phi}^{\mathrm{u}}}\ \textstyle\sum_{j=0}^{N-1}\hskip-3.0pt\sum_{k=j}^{N-1}\|\mathcal{Q}_{k,j}\mathbf{\Phi}_{k,j}\|^{2}_{\mathcal{F}}+\|\mathcal{Q}_{N,j}\mathbf{\Phi}^{\mathrm{x}}_{N,j}\|^{2}_{\mathcal{F}}$$ (34a)

s.t. $$\mathbf{\Phi}^{\mathrm{x}}_{k+1,j}=A_{k}\mathbf{\Phi}^{\mathrm{x}}_{k,j}+B_{k}\mathbf{\Phi}^{\mathrm{u}}_{k,j},$$ (34b)

$$\mathbf{\Phi}^{\mathrm{x}}_{j+1,j}=I_{n_{x}},~\forall j\in[N],~\forall k\in[j+1,N-1]$$ (34c)

В предложении 1 мы доказываем, что множества оптимизаторов (33) и (34) эквивалентны..

Proposition 1 (Эквивалентные оптимизаторы)

Пусть $(\mathbf{\tilde{\Phi}}^{\mathrm{x}},\mathbf{\tilde{\Phi}}^{\mathrm{u}})$ будет оптимизатором (33), а $(\mathbf{\hat{\Phi}}^{\mathrm{x}},\mathbf{\hat{\Phi}}^{\mathrm{u}})$ — оптимизатором переформулированной задачи (34). Предположим, что $\bar{E}_{j}$ является правообратимым для всех $j\in[N]$, то есть существует $\bar{E}_{j}^{\ddagger}$ такое, что $\bar{E}_{j}\bar{E}_{j}^{\ddagger}=I_{n_{x}}$. Тогда эти два оптимизатора совпадают: $\mathbf{\tilde{\Phi}}^{\mathrm{x}}=\mathbf{\hat{\Phi}}^{\mathrm{x}},\mathbf{\tilde{\Phi}}^{\mathrm{u}}=\mathbf{\hat{\Phi}}^{\mathrm{u}}$.

Доказательство:

Заметьте, что (33) может быть переписано в (43) с помощью замены переменной $\mathbf{\Psi}_{k,j}=\mathbf{\Phi}_{k,j}\bar{E}_{j}$, поскольку при использовании (33b) и (33c) мы имеем, что,

$$\mathbf{\Phi}^{\mathrm{x}}_{k+1,j}\bar{E}_{j}=A_{k}\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{E}_{j}+B_{k}\mathbf{\Phi}^{\mathrm{u}}_{k,j}\bar{E}_{j} \iff\text{Eq. \eqref{eq:fastsls_qp_x_prop_var}},$$

and                                           

$$\mathbf{\Phi}^{\mathrm{x}}_{j+1,j}\bar{E}_{j}=I_{n_{x}}\bar{E}_{j} \iff\text{Eq. \eqref{eq:fastsls_qp_x_init_var}.}$$

$$\hskip-3.0pt\min_{\mathbf{\Psi}^{\mathrm{x}},\mathbf{\Psi}^{\mathrm{u}}}\ \sum_{j=0}^{N-1}\hskip-3.0pt\sum_{k=j}^{N-1}\|\mathcal{Q}_{k,j}\mathbf{\Psi}_{k,j}\|^{2}_{\mathcal{F}}+\|\mathcal{Q}_{N,j}\mathbf{\Psi}^{\mathrm{x}}_{N,j}\|^{2}_{\mathcal{F}}$$ (43a)

s.t. $$\mathbf{\Psi}^{\mathrm{x}}_{k+1,j}=A_{k}\mathbf{\Psi}^{\mathrm{x}}_{k,j}+B_{k}\mathbf{\Psi}^{\mathrm{u}}_{k,j},$$ (43b)

$$\mathbf{\Psi}^{\mathrm{x}}_{j+1,j}=\bar{E}_{j},~\forall j\in[N],~\forall k\in[j+1,N-1]$$ (43c)

Из [10], (43) может быть решено посредством параллельных рекурсий Риккати,

$$\mathcal{G}_{N,j} =\mathcal{Q}^{\text{x}}_{N,j},\quad\quad\mathcal{K}_{k,j}=-\mathcal{G}_{k,j}\mathcal{B}_{k,j},$$ (44)

$$\mathcal{P}_{k,j} =\mathcal{Q}^{\text{x}}_{k,j}+\big(A_{k}\big)^{\top}\mathcal{P}_{k+1,j}A_{k}+\big(\mathcal{K}_{k,j}\big)^{\top}\mathcal{B}_{k,j},$$

$$\mathcal{G}_{k,j} =\big(\mathcal{Q}^{\text{u}}_{k,j}+\big(B_{k}\big)^{\top}\mathcal{P}_{k+1,j}B_{k}\big)^{-1},$$

$$\mathcal{B}_{k,j} =\mathcal{Q}^{\text{ux}}_{k,j}+\big(B_{k}\big)^{\top}\mathcal{P}_{k+1,j}A_{k},$$

и параллельные прямые распространения

$$\mathbf{\Psi}^{\mathrm{x}}_{j+1,j}=\bar{E}_{j},\qquad\mathbf{\Psi}^{\mathrm{u}}_{k,j}=\mathcal{K}_{k,j}\mathbf{\Psi}^{\mathrm{x}}_{k,j},$$ (45)

$$\mathbf{\Psi}^{\mathrm{x}}_{k+1,j}=(A_{k}+B_{k}\mathcal{K}_{k,j})\mathbf{\Psi}^{\mathrm{x}}_{k,j}.$$

Затем, поскольку $\bar{E}_{j}$ является правообратимой, а (43) представляет ограничение системного уровня (SLC), принятое в лемме 1, коэффициенты усиления контроллера $\mathbf{\Psi}$ подчиняются,

$$\mathbf{\Psi}^{\mathrm{x}}_{j+1,j}=\bar{E}_{j},~~~~~~~~~\mathbf{\Psi}^{\mathrm{u}}_{k,j}\bar{E}_{j}^{\ddagger}=\mathcal{K}_{k,j}\mathbf{\Psi}^{\mathrm{x}}_{k,j}\bar{E}_{j}^{\ddagger}$$ (46)

$$\mathbf{\Psi}^{\mathrm{x}}_{k+1,j}\bar{E}_{j}^{\ddagger}=\left(A_{k}+B_{k}\mathcal{K}_{k,j}\right)\mathbf{\Psi}^{\mathrm{x}}_{k,j}\bar{E}_{j}^{\ddagger}.$$

Затем, путем раскрытия рекурсии (46), мы получаем для $k>j$⟧,

$$\mathbf{\Psi}_{k,j}\bar{E}_{j}^{\ddagger}=\begin{bmatrix}I_{n_{x}}\\ \mathcal{K}_{k,j}\end{bmatrix}\left(\prod_{m=j+1}^{k-1}\left(A_{m}+B_{m}\mathcal{K}_{m,j}\right)\right)\mathbf{\Psi}^{\mathrm{x}}_{j+1,j}\bar{E}_{j}^{\ddagger},$$ (47)

в то время как остальные коэффициенты усиления ($\Psi_{k,j}$, где $k\leq j$) равны нулю, при этом мы определяем упорядоченное матричное произведение $\prod_{i=1}^{n}R_{i}:=R_{n}R_{n-1}\cdots R_{1}$.

Для восстановления оптимальных коэффициентов усиления для (33) мы используем свойство правой обратимости $\bar{E}_{j}$ и подставляем $\mathbf{\Psi}_{k,j}\bar{E}^{\ddagger}_{j}=\mathbf{\tilde{\Phi}}_{k,j}$ в (47⟧):

$$\mathbf{\tilde{\Phi}}_{k,j} =\begin{bmatrix}I_{n_{x}}\\ \mathcal{K}_{k,j}\end{bmatrix}\left(\prod_{m=j+1}^{k-1}\left(A_{m}+B_{m}\mathcal{K}_{m,j}\right)\right)\mathbf{\tilde{\Phi}}^{\mathrm{x}}_{j+1,j}$$

$$=\begin{bmatrix}I_{n_{x}}\\ \mathcal{K}_{k,j}\end{bmatrix}\left(\prod_{m=j+1}^{k-1}\left(A_{m}+B_{m}\mathcal{K}_{m,j}\right)\right)I_{n_{x}}.$$ (48)

Второе равенство следует из того, что $\mathbf{\tilde{\Phi}}^{\mathrm{x}}_{j+1,j}$ удовлетворяет (33c). Теперь заметим, что (34) имеет те же рекурсии Риккати, что и (44), что приводит к оптимальному коэффициенту усиления, структурированному аналогично (47), где $E_{j}^{\ddagger}$ заменено на $I_{n_{x}}$, $\mathbf{\hat{\Phi}}$:

$$\mathbf{\hat{\Phi}}_{k,j}=\begin{bmatrix}I_{n_{x}}\\ \mathcal{K}_{k,j}\end{bmatrix}\left(\prod_{m=j+1}^{k-1}\left(A_{m}+B_{m}\mathcal{K}_{m,j}\right)\right)\mathbf{\hat{\Phi}}^{\mathrm{x}}_{j+1,j},$$ (49)

где $\mathbf{\hat{\Phi}}^{\mathrm{x}}_{j+1,j}=I_{n_{x}}$, как ограничено (34c). Следовательно $\forall\,k,j\in[N]~\mathbf{\tilde{\Phi}}_{k,j}=\mathbf{\hat{\Phi}}_{k,j}$, таким образом $\mathbf{\tilde{\Phi}}=\mathbf{\hat{\Phi}}$. ∎

Lemma 1

Пусть $\bar{\bm{E}}:=\textrm{blk\_diag}(I_{n_{x}},\bar{E}_{0},\dots,\bar{E}_{n-1})$, где для всех $i\in[N],\bar{E}_{i}\in\mathbb{R}^{n_{x}\times n_{w}^{(i)}}$ является правообратимым, и $n_{w}^{(i)}\geq n_{x}$. Обеспечение неопределенной динамики LTV

$$x_{k+1}=A_{k}x_{k}+B_{k}u_{k}+\bar{E}_{k}w_{k},~w_{k}\in\mathcal{E}_{n_{w}^{(k)}},$$ (50)

эквивалентно соблюдению ограничения на системном уровне (51), то есть (50)$\iff$(51),

$$\begin{bmatrix}(I_{Nn_{x}}-\bm{Z}\bm{A})&-\bm{Z}\bm{B}\end{bmatrix}\begin{bmatrix}\mathbf{\Phi}^{\mathrm{x}}\\ \mathbf{\Phi}^{\mathrm{u}}\end{bmatrix}=\bar{\bm{E}}.$$ (51)

Более того, отображение коэффициентов усиления контроллера для SLC в (51) представляет собой,

$$\mathcal{K}=(\mathbf{\Phi}^{\mathrm{u}}\bm{E}^{\ddagger})(\mathbf{\Phi}^{\mathrm{x}}\bm{\bar{E}}^{\ddagger})^{-1}=(\mathbf{\Phi}^{\mathrm{u}}\bm{\bar{E}}^{\top})(\mathbf{\Phi}^{\mathrm{x}}\bm{\bar{E}}^{\top})^{-1},$$ (52)

где $\bm{\bar{E}}^{\ddagger}=\textrm{blk\_diag}(I,\bar{E}_{0}^{\ddagger},\dots,\bar{E}_{n-1}^{\ddagger})$. Наконец, при заданном $\mathcal{K}$, прямое распространение коэффициентов усиления возмущений подчиняется:

$$\mathbf{\Phi}^{\mathrm{x}}_{j+1,j}=\bar{E}_{j},~~~~~~~~~\mathbf{\Phi}^{\mathrm{u}}_{k,j}\bar{E}_{j}^{\ddagger}=\mathcal{K}_{k,j}\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{E}_{j}^{\ddagger}$$ (53)

$$\mathbf{\Phi}^{\mathrm{x}}_{k+1,j}\bar{E}_{j}^{\ddagger}=\left(A_{k}+B_{k}\mathcal{K}_{k,j}\right)\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{E}_{j}^{\ddagger}.$$

Доказательство:

Во-первых, пусть $\delta_{k}=E_{k}w_{k}$ будет таким, что

$$x_{k+1}=A_{k}x_{k}+B_{k}u_{k}+\delta_{k}.$$ (54)

Теперь пусть $\bm{x}=[x_{0}^{\top},\dots,x_{n-1}^{\top}]^{\top},~\bm{u}=[u_{0}^{\top},\dots,u_{n-1}^{\top}]^{\top},~\text{and}~\bm{\delta}=[\delta_{0}^{\top},\dots,\delta_{n-1}^{\top}]^{\top}$. Тогда, определяя отображения коэффициентов усиления $\mathbf{\hat{\Phi}}^{\mathrm{x}}$ и $\mathbf{\hat{\Phi}}^{\mathrm{u}}$ таким образом, что,

$$\begin{bmatrix}\bm{x}\\ \bm{u}\end{bmatrix}=\begin{bmatrix}\mathbf{\hat{\Phi}}^{\mathrm{x}}\\ \mathbf{\hat{\Phi}}^{\mathrm{u}}\end{bmatrix}\bm{\delta},$$ (55)

мы имеем из [8]:

$$\begin{bmatrix}(I_{Nn_{x}}-\bm{Z}\bm{A})&-\bm{Z}\bm{B}\end{bmatrix}\begin{bmatrix}\mathbf{\hat{\Phi}}^{\mathrm{x}}\\ \mathbf{\hat{\Phi}}^{\mathrm{u}}\end{bmatrix}=I_{Nn_{x}}$$ (56)

Определяя $\bm{w}=[w_{0}^{\top},\dots,w_{n-1}^{\top}]^{\top}$, мы наблюдаем:

$$\bm{\delta}=\bm{\bar{E}}\bm{w}$$ (57)

Затем отображение коэффициентов усиления от $\bm{w}$ к $\bm{x}$ и $\bm{u}$ через $\mathbf{\Phi}^{\mathrm{x}}$ и $\mathbf{\Phi}^{\mathrm{u}}$ следует из (55) и (57):

$$\begin{bmatrix}\bm{x}\\ \bm{u}\end{bmatrix}=\begin{bmatrix}\mathbf{\Phi}^{\mathrm{x}}\\ \mathbf{\Phi}^{\mathrm{u}}\end{bmatrix}\bm{w}=\begin{bmatrix}\mathbf{\hat{\Phi}}^{\mathrm{x}}\bm{\bar{E}}\\ \mathbf{\hat{\Phi}}^{\mathrm{u}}\bm{\bar{E}}\end{bmatrix}\bm{w}=\begin{bmatrix}\mathbf{\hat{\Phi}}^{\mathrm{x}}\\ \mathbf{\hat{\Phi}}^{\mathrm{u}}\end{bmatrix}\bm{\bar{E}}\bm{w}.$$ (58)

Это устанавливает, что $\mathbf{\Phi}^{\mathrm{x}}=\hat{\mathbf{\Phi}}^{\mathrm{x}}\bm{\bar{E}}$ и $\mathbf{\Phi}^{\mathrm{u}}=\hat{\mathbf{\Phi}}^{\mathrm{u}}\bm{\bar{E}}$. Следовательно, умножая (56) на $\bm{\bar{E}}$, мы приходим к (51). Таким образом, мы показали $\text{\eqref{eq:slc}}\implies\text{\eqref{eq:ltv_proof}}$. Чтобы показать эквивалентность между (50) и (51), мы можем показать $\text{\eqref{eq:ltv_proof}}\implies\text{\eqref{eq:slc}}$, путем прямого восстановления матриц $A_{k},B_{k},\bar{E}_{k}$ $\forall k\in[N]$ из (51) для получения (50⟧).

Чтобы доказать (52), мы замечаем, что $K=\mathbf{\hat{\Phi}}^{\mathrm{u}}(\mathbf{\hat{\Phi}}^{\mathrm{x}})^{-1}$, [8]. Затем, используя подстановки $\mathbf{\Phi}^{\mathrm{x}}\bm{\bar{E}}^{\ddagger}=\mathbf{\hat{\Phi}}^{\mathrm{x}}$ и $\mathbf{\Phi}^{\mathrm{u}}\bm{\bar{E}}^{\ddagger}=\mathbf{\hat{\Phi}}^{\mathrm{u}}$ из (58), мы приходим к первому равенству в (52). Второе равенство возникает из алгебраических преобразований с использованием определения $\bm{\bar{E}}^{\ddagger}=\bm{\bar{E}}^{\top}(\bm{\bar{E}}\bm{\bar{E}}^{\top})^{-1}$. Наконец, преобразуя (52), чтобы увидеть $\mathcal{K}(\mathbf{\Phi}^{\mathrm{x}}\bm{\bar{E}}^{\ddagger})=\mathbf{\Phi}^{\mathrm{u}}\bm{\bar{E}}^{\ddagger}$, и раскрывая ограничение системного уровня в (51), мы получаем, что

$$\mathbf{\Phi}^{\mathrm{x}}_{j+1,j}=\bar{E}_{j},~~~~~~~~~\mathbf{\Phi}^{\mathrm{u}}_{k,j}\bar{E}_{j}^{\ddagger}=\mathcal{K}_{k,j}\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{E}_{j}^{\ddagger}$$ (59)

$$\mathbf{\Phi}^{\mathrm{x}}_{k+1,j}\bar{E}_{j}^{\ddagger}=\left(A_{k}+B_{k}\mathcal{K}_{k,j}\right)\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{E}_{j}^{\ddagger}.$$

∎

Appendix D Доказательства

Theorem 4

Пусть $\mathbf{z}$, $\mathbf{v}$ и $\mathbf{\Phi}$ являются допустимым решением для (24). Тогда (26) гарантированно дает верхнюю аппроксимацию истинного достижимого множества замкнутой системы (3), то есть $\Omega_{k}\subseteq\bar{\Omega}_{k}$ для всех $k\in[N]$.

Доказательство:

Сначала мы докажем, что 26 является корректным при условии допустимых возмущений. Пусть $\mathcal{W}_{k}$ определяет возмущение $\forall k\in[N]$. Тогда, исходя из SLS [4], мы имеем, что аппроксимация сверху достижимого множества представляет собой,

$$\bar{\Omega}_{k}^{x}:=z_{k}+\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{x}}_{k,j}\mathcal{W}_{j},$$ (60)

$$\bar{\Omega}_{k}^{u}:=v_{k}+\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{u}}_{k,j}\mathcal{W}_{j}.$$

Теперь предположим, что $\mathcal{W}_{k}$ может быть выражено как зонотоп, $\bar{c}_{k}\oplus\bar{E}_{k}\mathcal{E}_{2n_{x}}$. Выполнив эту подстановку в (60), мы приходим к,

$$\bar{\Omega}_{k}^{x}:=z_{k}+\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{x}}_{k,j}(\bar{c}_{k}\oplus\bar{E}_{k}\mathcal{E}_{2n_{x}}),$$ (61)

$$\bar{\Omega}_{k}^{u}:=v_{k}+\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{u}}_{k,j}(\bar{c}_{k}\oplus\bar{E}_{k}\mathcal{E}_{2n_{x}}).$$

Используя ассоциативность и коммутативность суммы Минковского, 61 можно переписать в,

$$\bar{\Omega}_{k}^{x}:=z_{k}+\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{c}_{j}\oplus\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{E}_{j}\mathcal{E}_{2n_{x}},$$ (62)

$$\bar{\Omega}_{k}^{u}:=v_{k}+\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{u}}_{k,j}\bar{c}_{j}\oplus\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{u}}_{k,j}\bar{E}_{j}\mathcal{E}_{2n_{x}}.$$

Затем, вспоминая, что $\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{c}_{j}$ является суммой Минковского одноэлементных множеств, мы можем заменить его на $\textstyle\sum_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{x}}_{k,j}\bar{c}_{j}$ (и аналогично поступить для $\textstyle\bigoplus_{j=0}^{k-1}\mathbf{\Phi}^{\mathrm{u}}_{k,j}\bar{c}_{j}$) и прийти к множеству достижимости возмущений (26). Теперь мы покажем, что границы возмущений $\mathcal{Z}_{r}(z_{j},v_{j},\bm{\mu}_{j},\bm{\sigma}_{j})$, $\mathcal{Z}_{c}(z_{j},v_{j},\bm{\mu}_{j},\bm{\sigma}_{j})$ в (24d)-(24e) являются корректными аппроксимациями сверху $\mathcal{\bar{E}}_{k}$ (9) и подтвердим, что замкнутые множества достижимости $\bar{\Omega}^{x}_{k}$ и $\bar{\Omega}^{u}_{k}$ аппроксимируют сверху истинное множество достижимости $\Omega^{x}_{k}$ и $\Omega^{u}_{k}$ (3⟧).

Базовый случай: возмущение на временном шаге $k=0$ в точности равно $E(z_{0})=E(\bar{x}_{0})$, что учитывается в (25), и, таким образом, (24d)-(24e) обеспечивает верхнюю аппроксимацию границ возмущения.

Индуктивный шаг: Предположим, что для всех временных шагов $j\in[k]$, $d_{j}(x_{j},u_{j},w_{j})\in E_{j}\mathcal{E}_{n_{x}}\oplus\bar{R}_{j}\mathcal{E}_{n_{x}}\oplus\bar{c}_{j}$. Тогда мы имеем переаппроксимированные достижимые множества, $\Omega_{k}^{x}\subseteq\bar{\Omega}^{x}_{k}$ и $\Omega_{k}^{u}\subseteq\bar{\Omega}^{u}_{k}$, на временном шаге $k$. Определим $\bar{\Omega}_{k}:=\{[x^{\top},u^{\top}]^{\top}\mid x\in\bar{\Omega}^{x}_{k},u\in\bar{\Omega}^{u}_{k}\}$. Затем заметим, что $\bar{\Omega}_{k}\subseteq\mathcal{A}_{k}:=\left([z_{k}^{\top},v_{k}^{\top}]^{\top}+\bm{\mu}_{k}\right)\oplus\textrm{diag}(\bm{\sigma}_{k})\mathcal{E}_{n_{x}+n_{u}}$, как определено в (25a). Тогда, по построению, для всех $[x_{k}^{\top},u_{k}^{\top}]^{\top}\in\bar{\Omega}_{k}$ и для всех $w_{k}\in\mathcal{E}_{n_{x}}$, мы имеем, что $r(x_{k},u_{k},w_{k})\in\bar{c}_{k}\oplus\bar{R}_{k}\mathcal{E}_{n_{x}}$, где $(\bar{c}_{k},\bar{R}_{k})$ вычисляется с использованием $\mathcal{A}_{k}$ в (25). Аналогично, номинальное возмущение $E_{k}w_{k}\in E(z_{k})\mathcal{E}_{n_{x}}$. Таким образом, полное возмущение (8a), $r_{k}^{\textrm{lin}}+r_{k}^{\mathrm{dist}}+E_{k}w_{k}=d(x_{k},u_{k},w_{k})\in\bar{c}_{k}\oplus\bar{R}_{k}\mathcal{E}_{n_{x}}\oplus E(z_{k})\mathcal{E}_{n_{x}}$. Конкатенация двух матриц-генераторов дает (25b). Таким образом, граница возмущения на временном шаге $k$ в (24d)-(24e) является переаппроксимацией и содержит истинную остаточную ошибку, т.е. $d_{k}(x_{k},u_{k},w_{k})$ (8a).

После установления допустимых границ возмущений для всех шагов $k\in[N]$, SLS гарантирует, что аппроксимация сверху достижимого множества замкнутой системы (26), $\bar{\Omega}^{x}_{k}$ и $\bar{\Omega}^{u}_{k}$ являются корректными, то есть $\Omega^{x}_{k}\subset\bar{\Omega}^{x}_{k}$ и $\Omega^{u}_{k}\subset\bar{\Omega}^{u}_{k}$ для всех $k\in[N]$⟧. ∎

Appendix E Двойственные задачи GPUSLS $\tau$

Как обсуждалось в разд. V-A, задача RNOCP решается путем итеративного чередования между решением задачи оптимизации номинальной траектории и обновлением контроллера. Поскольку эти две задачи оптимизации решаются раздельно, мы используем двойственную переменную $\tau$ для обеспечения согласованности между ними, как обсуждалось в [10, App. B]. В связи с нашей переформулировкой мы изменяем процедуру следующим образом. Как и в [10, App. B], мы определяем вспомогательные члены $\beta_{k,j}\in\mathbb{R}^{n_{c}}$ и $\beta_{N,j}\in\mathbb{R}^{n_{f}}$, где

$$\beta_{k,j} =\left[G_{k}\mathbf{\Phi}_{k,j}\bar{c}_{j}+\|G_{k}\mathbf{\Phi}_{k,j}\bar{E}_{j}\|_{1,\text{r}}\right]^{\circ 2}$$ (63a)

$$\forall j\in[N],\qquad\forall k\in[j,N]$$ (63b)

$$\beta_{N,j} =\left[G_{N}\mathbf{\Phi}^{\text{x}}_{N,j}\bar{c}_{j}+\|G_{N}\mathbf{\Phi}^{\text{x}}_{N,j}\bar{E}_{j}\|_{1,\text{r}}\right]^{\circ 2}$$ (63c)

$$\forall j\in[N]$$ (63d)

где $[\cdot]^{\circ 2}$ обозначает поэлементное возведение в квадрат. Затем мы можем вычислить наши двойственные переменные $\tau$ с помощью

$$\tau_{k,j}=\frac{\bar{\lambda}_{k}}{\sqrt{\beta_{k,j}+\epsilon}}\quad\forall j\in[N],\quad\forall k\in[j,N]$$ (64)

где мы обозначаем $\bar{\lambda}_{k}$ как двойственную переменную оптимизации номинальной траектории (соответствующую двойственной переменной $\lambda_{k}$ в [10]), введенную здесь во избежание конфликта обозначений. Наконец, мы можем определить наши стоимостные слагаемые как

$$\mathcal{Q}_{k,j}=\left(\text{diag}\left(\sqrt{\tau_{k,j}}\right)G_{k},\begin{bmatrix}\tilde{Q}^{1/2}&0\\ 0&\tilde{R}^{1/2}\end{bmatrix}\right),$$ (65)

$$\mathcal{Q}_{N,j}=\left(\text{diag}\left(\sqrt{\tau_{N,j}}\right)G_{N},\tilde{Q}_{N}^{1/2}\right).$$

Следуя [9] и [10], данная формулировка сохраняет двойственную согласованность и приводит к эквивалентной переформулировке RNOCP в рамках предложенной декомпозиции.