Поиск идеальной траектории в пространстве состояний

Задачи оптимального управления окружают нас повсюду: как посадить ракету SpaceX на морскую платформу с минимальным расходом топлива? Как управлять химическим реактором для получения максимальной прибыли за заданное время? Как алгоритмическому трейдеру продать огромный пакет акций, минимально обвалив рыночную цену? В математике эти задачи сводятся к поиску управляющей функции (закона управления), которая минимизирует определенный функционал затрат (интеграл потерь).

Существуют два главных математических пути решения таких задач. Первый подход — это принцип максимума Понтрягина, который приводит к сложнейшим краевым задачам для обыкновенных дифференциальных уравнений (ОДУ). Второй, не менее мощный подход — метод динамического программирования, разработанный американским математиком Ричардом Беллманом в 1950-х годах. Принцип оптимальности Беллмана философски гласит: «Оптимальное поведение обладает тем свойством, что каковы бы ни были начальное состояние и решение в начальный момент, последующие решения должны составлять оптимальное поведение относительно состояния, получающегося в результате первого решения». В непрерывном времени этот принцип приводит к знаменитому уравнению Гамильтона-Якоби-Беллмана (HJB).

Сложность УЧП и вязкостные решения

Уравнение HJB — это нелинейное дифференциальное уравнение в частных производных (УЧП) первого или второго порядка (если в системе есть стохастический шум). Его неизвестной является функция ценности (Value function) — минимальная ожидаемая стоимость достижения цели из любого текущего состояния. Решив уравнение HJB, мы получаем закон управления в виде обратной связи (closed-loop control), что критически важно для робастности реальных систем управления, подверженных помехам.

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

Численные схемы: метод Годунова и Upwind

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

Для стабильного решения необходимо применять специализированные односторонние разностные схемы (Upwind schemes) первого порядка. Идея, заимствованная из вычислительной газовой динамики (схемы Годунова и Лакса-Фридрихса), заключается в том, что разностная производная должна браться строго против направления потока информации. Кроме того, уравнение HJB требует решения задачи минимизации (вычисления Гамильтониана) на каждом узле сетки, что делает алгоритм крайне ресурсоемким. Но главная преграда — это «проклятие размерности» Беллмана. Размерность сетки для УЧП равна количеству переменных состояния системы. Для 6-мерного дрона-квадрокоптера построение полной сетки невозможно ни на одном суперкомпьютере мира. В последние годы для пробивания этого барьера начали активно применяться методы машинного обучения, где уравнения HJB решаются с использованием глубоких нейронных сетей (Physics-Informed Neural Networks), открывая новую эру в робототехнике и теории управления.

Когда прошлое непрерывно влияет на будущее

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

Такие физические процессы с «исторической памятью» адекватно описываются интегральными уравнениями Вольтерры. Главной отличительной чертой уравнения Вольтерры (в отличие от уравнений Фредгольма) является то, что верхний предел интеграла не является фиксированным числом, а представляет собой саму независимую переменную (чаще всего время t). То есть мы интегрируем искомую функцию от момента зарождения системы (t=0) до текущего момента времени. Искомая функция может стоять только под знаком интеграла (уравнения первого рода) или присутствовать и вне его (уравнения второго рода).

Численные квадратуры и сведение к рекуррентной системе

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

Идея метода квадратур интуитивно понятна. Мы вводим на оси времени сетку дискретных узлов. Затем мы заменяем непрерывный интеграл из уравнения Вольтерры его конечно-разностным приближением (например, по формуле трапеций или методу Симпсона). Благодаря тому, что верхний предел интеграла переменный, мы получаем не единую гигантскую матрицу (как у Фредгольма), а удобную треугольную систему линейных алгебраических уравнений (СЛАУ).

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

Методы коллокации и аналоги Рунге-Кутты

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

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

Ограничения полиномиального разложения Тейлора

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

Однако у рядов Тейлора есть фундаментальные недостатки. Во-первых, их радиус сходимости строго ограничен особыми точками функции (полюсами) в комплексной плоскости. Если функция имеет разрыв или уходит в бесконечность, ряд Тейлора «взорвется» и не сможет ее описать. Во-вторых, даже внутри области сходимости ряд Тейлора работает хорошо только в очень узкой окрестности точки разложения. Попытка экстраполировать полином Тейлора далеко от центра приводит к феномену катастрофического расхождения: полином устремляется в бесконечность, игнорируя реальное поведение функции (например, асимптотическое стремление к константе).

Могущество дробно-рациональных функций Паде

В 1890 году французский математик Анри Паде (ученик знаменитого Эрмита) формализовал метод, который кардинально превосходит полиномы Тейлора по качеству аппроксимации. Аппроксимация Паде предлагает искать приближение функции не в виде одного многочлена, а в виде рациональной функции — отношения двух полиномов (числителя степени M и знаменателя степени N).

Идея построения аппроксимации Паде поразительно элегантна. Мы приравниваем исходную функцию (или ее известный ряд Тейлора) к искомой дроби. Затем умножаем обе части на знаменатель дроби. Приравнивая коэффициенты при одинаковых степенях независимой переменной слева и справа, мы получаем систему линейных алгебраических уравнений. Решив эту простую СЛАУ, мы находим искомые коэффициенты полиномов в числителе и знаменателе. Обозначается такая аппроксимация обычно как [M/N].

Преодоление полюсов и аналитическое продолжение

Наличие полинома в знаменателе наделяет аппроксимацию Паде математической «суперспособностью». Если исходная физическая функция имеет полюс (то есть обращается в бесконечность, как, например, функция тангенса при 90 градусах), ряд Тейлора ломается. Аппроксимация Паде, напротив, блестяще моделирует этот физический разрыв: знаменатель дроби просто сам обращается в ноль в этой точке, идеально повторяя топологию исходной функции!

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