Журнал вычислительной математики и математической физики, 2023, T. 63, № 1, стр. 165-174

Обучение порт-гамильтоновым системам: алгоритмы

Д. Лозиенко 1*, В. Сальников 1**, А. Фалез 1***

1 ЦНРС & Университет Ла-Рошели
17042 Ла Рошель, Пр. Мишеля Крепо, Франция

* E-mail: daria.loziienko1@univ-lr.fr
** E-mail: vladimir.salnikov@univ-lr.fr
*** E-mail: antoine.falaize@univ-lr.fr

Поступила в редакцию 04.08.2022
После доработки 21.08.2022
Принята к публикации 10.09.2022

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

Аннотация

Изучаются возможности определения структуры порт-гамильтоновых систем, соответствующих описывающим задачи механики системам обыкновенных дифференциальных уравнений по переменным без “определенной природы”. Предложенный алгоритм решает эту задачу в два этапа. Сначала с помощью методов машинного обучения восстанавливается структура связей системы, таким образом строится граф, характеризующийся выделенными подсистемами и взаимодействиями между ними. Затем этот граф дополняется построенными гамильтоновыми структурами каждой подсистемы, а также соответствующими портами. Описанный второй этап основывается на результатах из симплектической и пуассоновой геометрии, которые мы вкратце изложим. А конкретные решения строятся с использованием методов компьютерной алгебры и символьных вычислений. Представленный алгоритм позволяет расширить область применения порт-гамильтонова формализма до произвольных обыкновенных дифференциальных уравнений, потенциально вводя тем самым новое понятие нормальных форм обыкновенных дифференциальных уравнений. Библ. 18. Фиг. 1.

Ключевые слова: симплектические структуры, пуассонова геометрия, порт-гамильтоновы системы.

1. ВВЕДЕНИЕ. МОТИВАЦИЯ

Настоящая статья – часть большого проекта, связанного с направлением, которое мы называем “геометризация механики”. Сам проект состоит из нескольких этапов, связанных одной целью: изучить геометрический формализм, подходящий для описания внутренней структуры дифференциальных уравнений, моделирующих задачи механики, а также для вывода этих уравнений, исходя из физических свойств системы, используя геометрические структуры (см. обзор в [1], а также работы, на которые этот обзор ссылается).

Наиболее подробно изученный пример такого подхода касается естественных механических систем, описанных в канонических переменных (координаты и импульсы), в таком случае удобным формализмом будут гамильтоновы системы, определенные на симплектических многообразиях. Это описание, действительно, классическое и восходит к Лагранжу, Лежандру, Пуанкаре, хотя в их работах формализм не был написан современным языком Сурьё и Арнольда. Важно отметить, что эта логика работает “в обоих направлениях”. С одной стороны, для конкретной симплектической структуры, описывающей свойства фазового пространства системы, зная гамильтониан, можно явно выписать уравнения движения. Это особенно важно, так как сохранение симплектической структуры потоком гарантирует сохранение энергии; данный факт дает начало целой “индустрии” симплектических интеграторов (см. [2], [3]) – численных методов, сохраняющих геометрическую структуру. С другой стороны, зная уравнения движения, можно, естественно, при выполнении некоторых условий восстановить симплектическую структуру и гамильтониан. Мы обсудим подробнее это “обратное направление”.

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

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

В предыдущем абзаце мы использовали слово “сборка” неслучайно: оно подчеркивает, что в контексте ПГС есть “простое” направление: когда физическая система естественным образом явно разделена на подсистемы, между которыми известны законы взаимодействия, соответствующая ПГС выписывается явно (см. [9]). Сложности же возникают при решении обратной задачи: зная систему дифференциальных уравнений, восстановить структуру соответствующей ПГС – это и будет основной темой исследования в настоящей статье. В [10] сформулирован ряд открытых вопросов и задач, к решению которых разумно подходить с использованием методов компьютерной алгебры, в том числе и задача оптимального восстановления ПГС. Как мы обясним ниже, структура ПГС может быть задана так называемым размеченным графом – его топология восстанавливается методами машинного обучения, а “разметка” естественным образом определяется с помощью символьных вычислений и компьютерной алгебры.

В следующих разделах мы разобьем задачу восстановления структуры ПГС на более явные подзадачи и объясним идеи решения каждой из них. Вначале мы напомним классический гамильтонов формализм (разд. 2), затем введем конкретную терминологию для ПГС (разд. 3) и опишем желаемый результат решения задачи вместе с некоторыми методами его достижения (разд. 4). При написании данной статьи мы осознали, что при программной реализации каждого из шагов метода возникает некоторое количество нюансов. Но чтобы не перегружать статью техническими подробностями, мы решили вынести вопросы программирования и бенчмаркинга, а также сложные примеры в отдельный текст (см. [11]), сосредоточившись здесь на теоретических основах и математическом аппарате.

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

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

Самая простая гамильтонова система имеет вид

(1)
${\mathbf{\dot {q}}} = \frac{{\partial H}}{{\partial {\mathbf{p}}}},\quad {\mathbf{\dot {p}}} = - \frac{{\partial H}}{{\partial {\mathbf{q}}}},$
где ${\mathbf{q}}$ – векторная переменная размера $n$, обозначающая (обобщенные) координаты системы, а ${\mathbf{p}}$ – соответствующие импульсы. Легко проверить, что поток такой системы сохраняет поверхности уровня $H({\mathbf{q}},{\mathbf{p}})$, что физически обозначает сохранение энергии. Даже для этого простого случая “обратная” задача имеет смысл: для произвольной системы
${\mathbf{\dot {q}}} = {\mathbf{f}}({\mathbf{q}},{\mathbf{p}}),\quad {\mathbf{\dot {p}}} = {\mathbf{g}}({\mathbf{q}},{\mathbf{p}})$
можно ли найти функцию $H$, определяющую эквивалентную гамильтонову систему? Ответ основан на условии совместности
$\frac{{\partial {\mathbf{f}}}}{{\partial {\mathbf{q}}}} = - \frac{{\partial {\mathbf{g}}}}{{\partial {\mathbf{p}}}},$
которое, очевидно, необходимо и (по крайней мере, локально) достаточно. Но чтобы сделать ответ конструктивным, требуется некоторая явная операция (в данном случае – интегрирование).

Хотя мы не упомянули симплектическую структуру в уравнении (1), она там неявно присутствует: систему можно переписать в более общем виде

(2)
${\mathbf{\dot {x}}} = J({\mathbf{x}})\frac{{\partial H}}{{\partial {\mathbf{x}}}},$
где $J({\mathbf{x}})$ – кососимметричная невырожденная матрица, удовлетворяющая некоторому дифференциальному условию. (Условие замкнутости соответствующей дифференциальной формы, см. упрощенное изложение терминологии в [10] или более подробное построение общей теории в [12].) Пример (1) соответствует ${\mathbf{x}} = ({\mathbf{q}},{\mathbf{p}})$, а
(3)
$J = {{J}_{D}} = \left( {\begin{array}{*{20}{c}} 0&{{{I}_{n}}} \\ { - {{I}_{n}}}&0 \end{array}} \right),$
т.е. гамильтоновой системе на каноническом симплектическом пространстве, обозначаемом $T{\kern 1pt} *{\kern 1pt} Q$. Теорема Дарбу (см., например, [13]) утверждает, что любая симплектическая структура может быть приведена к такому каноническому виду. В заданной точке это простой результат из линейной алгебры, но в ее окрестности конструкция более сложная: чтобы выбрать подходящие координаты, нужно интегрировать некоторый поток.

Систему (2) можно также рассматривать для любой кососимметричной матрицы $J$, не обязательно невырожденной. Если она снова удовлетворяет некоторому дифференциальному условию, то задает так называемую структуру Пуассона. Тождеству Якоби для соответствующей скобки Пуассона (см. [10]). Согласно теореме Вайнштена о расщеплении (см. [13]), пространство со структурой Пуассона представляет расслоение, слои которого симплектические. Иными словами, максимальное пространство, ограничение $J$ на которое невырождено, допускает координаты Дарбу $({\mathbf{q}},{\mathbf{p}})$, а в трансверсальном направлении координаты $({\mathbf{y}})$ можно выбрать специальным образом, независимо от остальных:

$J = {{J}_{W}} = \left( {\begin{array}{*{20}{c}} {{{J}_{D}}}&0 \\ 0&{\Phi ({\mathbf{y}})} \end{array}} \right),$
где ${{J}_{D}}$ – блочная матрица вида (3), а $\Phi ({\mathbf{y}})$ – кососимметричная матрица, коэффициенты которой не зависят от координат $({\mathbf{q}},{\mathbf{p}})$, а также все ${{\Phi }_{{ij}}}({\mathbf{0}}) = 0$.

Мы заметили, что в примере (1) поток системы сохраняет поверхности уровня энергии. Это свойство выполняется и для системы (2), оно следует, по сути, из кососимметричности матрицы $J$. Верно даже более сильное утверждение: сохраняются не только уровни энергии, но и вся геометрическая структура фазового пространства – симплектическая форма или соответственно пуассоново расслоение. Это более общее утверждение чуть сложнее “обратить”: на самом деле, мы хотим по заданному потоку системы подобрать гамильтониан, который порождает его для некоторой симплектической или пуассоновой структуры. Важно также заметить, что сохранение геометрической структуры некоторым потоком само по себе не гарантирует, что он является гамильтоновым. На самом деле, разница между множествами сохраняющих структуру и гамильтоновых векторных полей характеризуется некоторой когомологией фазового пространства. Не вдаваясь в геометрические подробности, дадим лишь корректные названия этим объектам. В симплектическом случае речь идет о классе первой когомологии де Рама, который характеризует топологию пространства и вполне вычислим. В случае структуры Пуассона класс тоже называется первой когомологией Пуассона, и он вычислим для большинства важных структур. Мотивированному читателю мы рекомендуем [13], где даны подробные определения и приведены примеры. Таким образом, разумно анализировать поток системы в два этапа: сначала убедиться, что он сохраняет структуру, а затем проверить, есть ли препятствия к его гамильтоновости. Подобный подход использовался в [4] для построения пуассоновых интеграторов.

3. ПОРТ-ГАМИЛЬТОНОВ ФОРМАЛИЗМ ДЛЯ ВЗАИМОДЕЙСТВУЮЩИХ И ОТКРЫТЫХ СИСТЕМ

В этом разделе мы приводим идеи порт-гамильтонова формализма и напоминаем или определяем соответствующую терминологию. Как мы уже упоминали, основная мысль состоит в изменении гамильтоновых уравнений добавлением негамильтоновых слагаемых. Например, систему (1) можно обобщить до

(4)
${\mathbf{\dot {q}}} = \frac{{\partial H}}{{\partial {\mathbf{p}}}},\quad {\mathbf{\dot {p}}} = - \frac{{\partial H}}{{\partial {\mathbf{q}}}} + {\mathbf{F}},$
где ${\mathbf{F}}$ – некоторая сила, природу которой мы пока не уточняем. Подобным образом в неканоническом симплектическом или пуассоновом случае (2) система теперь может принимать вид
(5)
${\mathbf{\dot {x}}} = (J({\mathbf{x}}) - R({\mathbf{x}}))\frac{{\partial H}}{{\partial {\mathbf{x}}}} + w({\mathbf{x}}){\mathbf{u}},$
где $J({\mathbf{x}})$ – по-прежнему кососимметричная матрица, а $R({\mathbf{x}}),w({\mathbf{x}})$ и ${\mathbf{u}}$ – новые слагаемые. С качественной точки зрения $R$ соответствует внутренним силам в системе, в то время как $w$ и ${\mathbf{u}}$ отвечают за взаимодействие с “внешним миром”. Мы различаем $w$ и ${\mathbf{u}}$, чтобы подчеркнуть, что это взаимодействие может разделяться на части, зависящие от состояния системы и на полностью контролируемые средой. Мы записываем соответствующее слагаемое в виде произведения для простоты изложения с очевидным возможным обобщением. В [6], [7] эти слагаемые различаются более подробно, исходя из некоторой физической терминологии: они могут описывать накопление энергии, диссипацию, управление и т.д. Но для данной статьи уровень детализации (5) достаточен для объяснения предложенного подхода.

Определяющее свойство ПГС – возможность их “соединять”, что объясняет терминологию “портов”, как в электрических схемах. Иными словами, если рассмотреть несколько изначально независимых систем вида (5) и обеспечить их взаимодействие через соответствующие порты (слагаемые $w{\mathbf{u}}$), то результат будет системой, большей, чем каждая из составляющих, но снова имеющей вид (5). Легко убедиться, что в таком случае глобальная структура этой системы будет довольно специфическая. Матрицы $J$ и $R$ будут блочного вида, характеризуя внутреннюю динамику подсистем, в то время как слагаемые $w{\mathbf{u}}$ будут в точности обозначать взаимодействие между блоками.

ПГС удобно описывать графами. Обычно их строят из элементарных узлов, как на фиг. 1. Вершины графа соответствуют внутренним (гамильтоновым) составляющим системы, а ребра – портам (или наоборот в двойственном представлении). Некоторые из ребер идут к другим подсистемам, а некоторые заканчиваются виртуальными вершинами (т.е. идут как бы в пустоту), что соответствует диссипации или управлению. В терминах [7] портам соответствуют переменные потока (flux) и напряжения (effort), что позволяет снабдить ребра графа метками, отвечающими за негамильтонову часть системы (5). Эти переменные естественным образом двойственны друг другу, что делает связанные системы переопределенными. Таким образом, переменные будут связаны некоторыми алгебраическими соотношениями, которые определяют упомянутые ранее структуры Дирака, а физически они соответствуют балансу полной энергии системы.

Фиг. 1.

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

Такое представление с помощью графов и используется явно, чтобы “собрать” ПГС из блоков, а затем выбрать оптимальную стратегию для численного моделирования. Этот подход был успешно реализован в программном коде с использованием распределенных вычислений (см. [9]) и применен для некоторых сложных физических систем (см. [14], [15]). Ниже мы снова рассмотрим обратную задачу по системе дифференциальных уравнений в наиболее общем виде:

(6)
${\mathbf{\dot {x}}} = {\mathbf{f}}({\mathbf{x}}),$
и построим ее порт-гамильтоново представление (5).

4. АЛГОРИТМЫ И ПРИМЕРЫ

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

Конкретный вопрос, на самом деле, несколько более деликатен, чем просто алгоритм перехода от (6) к (5). Заметим, что решение данной задачи совершенно не обязано быть единственным: мы уже упоминали, что, если соединить несколько ПГС, получится новая ПГС – понятно, что в решении обратной задачи нет никаких причин полагать, что такого слияния не происходит. Более того, для некоторых слагаемых возможен выбор: отнести их к гамильтонианам или к портам – это очень важная “степень свободы”, позволяющая значительно упростить некоторые процедуры. С учетом всех этих замечаний, нужно потребовать, чтобы процесс был в некотором смысле оптимальным, иначе можно построить совершенно неразумные решения. Рассмотрим, например, все правые части уравнений системы (6) как независимые порты – формально это решение задачи, но с большим количеством вырождений и полностью игнорирующее гамильтонову составляющую. Или, как в другом крайнем случае, если система большая, но изолированная от среды, скорее всего, ее можно полностью переписать в гамильтоновом виде, естественно, полностью теряя информацию о внутренних связях.

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

4.1. Определение структуры графа

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

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

Анализ разных алгоритмов показал, что задача отлично подходит для методов машинного обучения, что объясняет название статьи. С одной стороны, для обучения и тестов мы можем генерировать любое необходимое количество примеров ПГС с заданной топологией связей (см. [9]). С другой, мы четко определили набор параметров, которые нужно оптимизировать. Подобные методы известны под названием алгоритмы свертки (pooling operators), это довольно простые нейронные сети с небольшим количеством слоев.

В данной статье мы упоминаем машинное обучение в режиме “черного ящика” для полноты конечного алгоритма. Технические подробности, а также сравнение производительности методов будет подробно описано в [11]. Тем не менее уже здесь важно сделать пару замечаний для четкого понимания алгоритма. Исходный граф строится явно по правым частям системы дифференциальных уравнений (6): в первом приближении все переменные ${{x}_{i}}$ объявляются независимыми вершинами, их пары связаны (направленным) ребром, если одна их них присутствует в правой части уравнения для другой. А на конечном этапе группы вершин с большим количеством попарных связей объединяются в узлы, они соответствуют гамильтоновым подсистемам, а связи между узлами объявляются портами. Некоторая “тонкая настройка” возможна на разных этапах этой оптимизации, а именно, слагаемые определенного вида в правых частях дифференциальных уравнений можно насильно учитывать внутри или снаружи узлов, мы вернемся к этой идее ниже.

4.2. Разметка графа

Результат предыдущего этапа – граф, который по идее меньше, чем исходный, он описывает скопления переменных (узлы). А слагаемые в дифференциальных уравнениях разделяются на две группы: внутренние и внешние (отвечающие за взаимодействие узлов).

Разметка внутренней части – наиболее сложный этап с математической точки зрения. Изолированные (небольшие) системы дифференциальных уравнений нужно разделить на структурную часть (симплектическая/пуассонова геометрия для гамильтоновых систем, которую мы определили в разд. 2) и “порт-” часть (добавленная в разд. 3).

Даже если это, возможно, было понятно уже в разд. 2, подчеркнем еще раз, что обратная задача в общем случае не решается явно методами компьютерной алгебры и символьных вычислений. Основная трудность заключается в приведении структуры (симплектической или Пуассона) к нормальной форме, оно основано на теоремах существования, которые доказываются интегрированием определенных гамильтоновых потоков. Понятно, что возможность сделать это символьными методами – очень сильное условие, сравнимое с существованием явного аналитического решения системы. Для некоторых приложений можно решить задачу приближенно для коротких траекторий, например, с помощью преобразования Коссра (аналог уравнения Гамильтона–Якоби, см. [4]). Но такой подход позволяет построить эквивалентную систему только локально, что недостаточно для ее дальнейшего использования в моделировании.

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

Более конкретно, процедура выполняется в несколько этапов. На первом составляется каталог простых симплектических и пуассоновых структур, т.е. матриц $J$ подходящих размеров с заданным вырождением. Это могут быть симплектические структуры, записанные в координатах Дарбу, все структуры Пуассона постоянного ранга, линейные и квадратичные структуры Пуассона с некоторыми особенностями. Список можно расширить более нелинейными системами, но даже упомянутые структуры покрывают достаточный геометрический спектр. На втором этапе вычисляются соответствующие первые когомологии (см. разд. 2 и [13]). Для приведенных примеров они будут не очень большими, а иногда и вообще тривиальными, для некоторых структур они известны явно. На третьем этапе из правых частей уравнений выбираются слагаемые, сохраняющие построенные структуры. Те из них, которые не лежат в нетривиальном классе когомологии, интегрируются до гамильтониана. Этот последний шаг, в отличие от нормальных форм структуры, явный и может быть выполнен с помощью символьных вычислений. И даже если он оказывается технически сложным, класс гамильтонианов всегда может быть ограничен до полиномов или порожден другим конкретным набором элементарных функций. Более “академический” вариант последнего этапа – построить (один раз) список порождающих векторных полей, сохраняющих упомянутые выше структуры. Этот результат может иметь независимые геометрические приложения. А затем идентифицировать соответствующие слагаемые в правых частях уравнений. По выполнении этого этапа будет восстановлена гамильтонова часть уравнений и перечислены оставшиеся неидентифицированные слагаемые.

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

4.3. Алгоритм

Сформулируем теперь все замечания из предыдущих подразделов в виде единого алгоритма.

Входные данные: система дифференциальных уравнений в виде ${\mathbf{\dot {x}}} = {\mathbf{f}}({\mathbf{x}})$:

Шаг 1. По правым частям системы восстановить структуру графа ПГС.

В результате переменные ${{x}_{i}}$ разбиты на группы (вершины графа), а слагаемые в правых частях разделены на внутренние (зависящие только от переменных “своей” группы) и внешние (остальные).

Шаг 2. Для каждой отдельной группы

a) построить (или выбрать из заранее построенного каталога) все симплектические и пуассоновы структуры подходящих размеров;

б) из внутренних переменных выбрать максимальную комбинацию слагаемых, сохраняющих одну из структур из предыдущего пункта;

в) если когомология выбранной структуры не тривиальна, проверить выбранную комбинацию на принадлежность тривиальному классу;

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

В результате для каждой группы из внутренних переменных выделены компоненты Гамильтонова потока.

Шаг 3. Идентифицировать порты:

a) слагаемые, не выбранные в шаге 2б) объявить внутренними портами, сопоставить им виртуальные вершины;

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

В результате построен размеченный граф, как на фиг. 1.

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

4.4. Примеры

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

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

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

$\ddot {x} + a\dot {x} + x = 0.$
Она переписывается в форме системы (6):
${{\dot {x}}_{1}} = {{x}_{2}},\quad {{\dot {x}}_{2}} = - {{x}_{1}} - a{{x}_{2}},$
которая, в свою очередь, имеет вид (4), т.е. (5) для
$J = {{J}_{D}} = \left( {\begin{array}{*{20}{c}} 0&1 \\ { - 1}&0 \end{array}} \right),\quad H = \frac{1}{2}(x_{1}^{2} + x_{2}^{2}),$
а оставшемуся слагаемому $ - a{{x}_{2}}$ соответствует диссипативный порт.

Твердое тело с неподвижной точкой. Данная система описывает (см. [12]) динамику твердого тела в терминах компонент кинетического момента ${\mathbf{M}} \in {{\mathbb{R}}^{3}}$ под действием внешней силы ${\mathbf{N}} \in {{\mathbb{R}}^{3}}$:

${{\dot {M}}_{1}} = {{a}_{1}}{{M}_{2}}{{M}_{3}} + {{N}_{1}},\quad {{\dot {M}}_{2}} = {{a}_{2}}{{M}_{3}}{{M}_{1}} + {{N}_{2}},\quad {{\dot {M}}_{3}} = {{a}_{3}}{{M}_{1}}{{M}_{2}} + {{N}_{3}}.$
Здесь сила ${\mathbf{N}}$ задана явно, а коэффициенты ${{a}_{i}}$ зависят от геометрии тела. Если ${{I}_{1}},\;{{I}_{2}},\;{{I}_{3}}$ – моменты инерции тела, то
${{a}_{1}} = \frac{{{{I}_{2}} - {{I}_{3}}}}{{{{I}_{2}}{{I}_{3}}}},\quad {{a}_{2}} = \frac{{{{I}_{3}} - {{I}_{1}}}}{{{{I}_{3}}{{I}_{1}}}},\quad {{a}_{3}} = \frac{{{{I}_{1}} - {{I}_{2}}}}{{{{I}_{1}}{{I}_{2}}}}.$
Из соображений размерности ясно, что матрица $J$ не может быть невырожденной, поэтому разумно искать решение среди структур Пуассона. Действительно, линейная структура
$J = \left( {\begin{array}{*{20}{c}} 0&{ - {{M}_{3}}}&{{{M}_{2}}} \\ {{{M}_{3}}}&0&{ - {{M}_{1}}} \\ { - {{M}_{2}}}&{{{M}_{1}}}&0 \end{array}} \right)$
с квадратичным гамильтонианом
$H = \frac{{M_{1}^{2}}}{{2{{I}_{1}}}} + \frac{{M_{2}^{2}}}{{2{{I}_{2}}}} + \frac{{M_{3}^{2}}}{{2{{I}_{3}}}}$
восстанавливают консервативную часть системы, а сила ${\mathbf{N}}$ определяет порт.

Обобщение системы Лотки–Вольтерра. Данная система обобщает классическую модель динамики популяций на случай произвольной размерности:

${{\dot {x}}_{i}} = {{\varepsilon }_{i}}{{x}_{i}} + \sum\limits_{j = 1}^n \,{{A}_{{ij}}}{{x}_{i}}{{x}_{j}},\quad i = 1, \ldots ,n.$
Для нечетного $n = 2k + 1$ в отсутствие линейных правых частей (${{\varepsilon }_{i}} = 0$) она называется системой Богоявленского–Ито, если ${{A}_{{ij}}} = 1$ для $i + 1 \leqslant j \leqslant \min (i + k,n)$, и ${{A}_{{ij}}} = - 1$ для ${\text{min}}(i + k,n) \leqslant j \leqslant n$. Такие системы подробно изучены в [16]. В рамках описанного здесь подхода для каждого фиксированного значения $n$ легко проверить, что система гамильтонова для квадратичной структуры Пуассона ${{J}_{{ij}}}(x) = {{A}_{{ij}}}{{x}_{i}}{{x}_{j}}$ и линейного гамильтониана $H = {{x}_{1}} + \ldots + {{x}_{n}}$. Линейные слагаемые ${{\varepsilon }_{i}}{{x}_{i}}$ могут быть добавлены как диссипативные порты.

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

$a\ddot {y} + b\dot {y} + cy = dq,$
$\ddot {q} + k({{q}^{2}} - 1)\dot {q} + lq = m\ddot {y},$
$y$ и $q$ суть динамические переменные, отвечающие соответственно за положение тела и за интенсивность завихрений среды; остальные латинские буквы обозначают некоторые постоянные параметры системы. Перепишем систему в виде ${\mathbf{\dot {x}}} = {\mathbf{f}}$:
${{\dot {x}}_{1}} = {{x}_{2}},$
${{\dot {x}}_{2}} = - \frac{b}{a}{{x}_{2}} - \frac{c}{a}{{x}_{1}} + \frac{d}{a}{{x}_{3}},$
${{\dot {x}}_{3}} = {{x}_{4}},$
${{\dot {x}}_{4}} = - k(x_{3}^{2} - 1){{x}_{4}} - l{{x}_{3}} + m{{\dot {x}}_{2}}.$
Заметим, что система выше, строго говоря, еще не приведена к виду ${\mathbf{\dot {x}}} = {\mathbf{f}}({\mathbf{x}})$ из-за зависимости правой части четвертого уравнения от ${{\dot {x}}_{2}}$. Когда мы изучали ее в [18], ${{\dot {x}}_{2}}$ заменялась на соответствующие правые части в явном виде. Для текущего подхода это делать не обязательно: достаточно обозначить $m{{\dot {x}}_{2}}$ новой переменной $u({{x}_{1}},{{x}_{2}},{{x}_{3}})$, не забывая, от чего она зависит для построения графа.

Первое приближение графа на шаге 1 алгоритма задается матрицей инцидентности

$\left( {\begin{array}{*{20}{c}} 0&1&0&0 \\ 1&1&1&0 \\ 0&0&0&1 \\ 1&1&1&1 \end{array}} \right).$
Неудачная его свертка вернет один узел из всех четырех переменных, построенный вокруг ${{x}_{4}}$, а разумная – два узла: $({{x}_{1}},{{x}_{2}})$ и $({{x}_{3}},{{x}_{4}})$, соединенные между собой. Тогда на шаге 2 внутри каждого узла гамильтонова структура соответствует гармоническому осциллятору с диссипативным портом, как в самом первом примере этого раздела. Взаимодействие между узлами, естественно, задается портом $u$ и двойственным ему.

5. ОБСУЖДЕНИЕ. ПЕРСПЕКТИВЫ

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

Алгоритм исходно задумывался для “упорядочивания” взаимодействующих мультифизических систем, но он, конечно, может использоваться и более широко. Например, для искусственных моделей с нефизическими переменными (т.е. не обязательно интерпретируемыми в терминах состояния механической системы), подход все равно позволит построить некоторое разложение и выявить внутреннюю структуру. Пример такой ситуации для упрощенной модели взаимодействия жидкости с препятствием рассмотрен выше. Некоторые подобные примеры также будут изучены в [11].

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

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

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

  1. Salnikov V., Hamdouni A., Loziienko D. Generalized and graded geometry for mechanics: a comprehensive introduction // Math. and Mech. of Complex Syst. 2021. V. 9. № 1. P. 59–75.

  2. Verlet L. Computer “Experiments” on Classical Fluids // Phys. Rev. 1967. V. 159. P. 98–103.

  3. Yoshida H. Construction of higher order symplectic integrators // Phys. Lett. A. 1990. V. 150. № 5–7. P. 262–268.

  4. Cosserat O. Symplectic groupoids for Poisson integrators, Preprint: arXiv:2205.04838.

  5. Paynter H.M. Analysis and Design of Engineering Systems, MIT Press, Cambridge, Massachusetts, 1961.

  6. Maschke B.M., van der Schaft A.J., Breedveld P.C. An intrinsic Hamiltonian formulation of network dynamics: nonstandard Poisson structures and gyrators // J. Franklin Inst. 1992. V. 329. № 5. P. 923–966.

  7. van der Schaft A. Port-Hamiltonian systems: an introductory survey // Proceed. of the Inter. Congress of Math. 2006. V. III. P. 1339–1365, Madrid.

  8. Cosserat O., Laurent-Gengoux C., Kotov A., Ryvkin L., Salnikov V. On Dirac structures admitting a variational approach, Preprint: arXiv:2109.00313.

  9. Falaize A. Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios: Approche par réseau de composants et formulation hamiltonienne à ports, PhD thesis, Télécomm. et Électronique de Paris, Université Pierre et Marie Curie, 2016.

  10. Сальников В.Н., Хамдуни А. Дифференциальная геометрия и механика – источник задач для компьютерной алгебры // Программирование. 2020. № 2. С. 60–66.

  11. Salnikov V., Falaize A., Loziienko D. Learning port-Hamiltonian systems – applications, готовится к печати.

  12. Арнольд В.И. Математические методы классической механики. М.: Наука, 1974.

  13. Cannas Da Silva A., Weinstein A. Geometric Models for Noncommutative Algebras, Am. Math. Soc. 2000.

  14. Falaize A., Hélie T. Passive guaranteed simulation of analog audio circuits: A port-hamiltonian approach // Applied Science, Applied Acoustics, special issue Audio Signal Process. 2016. V. 6. P. 273.

  15. Falaize A., Hélie T. Passive simulation of the nonlinear port-hamiltonian modeling of a rhodes piano // J. of Sound and Vibrat. 2016. V. 390. P. 289–309.

  16. Evripidou C.A., Kassotakis P., Vanhaecke P. Integrable deformations of the Bogoyavlenskij-Itoh Lotka-Volterra systems // J. of Regular and Chaotic Dynam. 2017. V. 22 P. 721–739.

  17. Leclercq T., de Langre E. Vortex-induced vibrations of cylinders bent by the flow // J. of Fluids and Structur. 2018. V. 80. P. 77–93.

  18. Salnikov V., Hamdouni A. Geometric integrators in mechanics – the need for computer algebra tools // Proceed. of the Third Inter. Conf. “Computer algebra”, 40–46, 2019, Moscow, Russia.

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