Прямое интегрирование дифференциально-алгебраических уравнений индекса 3 для динамики многозвенных механизмов в избыточных координатах с комплексной схемой Розенброка
Моделирование механизмов с замкнутыми кинематическими цепями в избыточных (максимальных) координатах даёт систему дифференциально-алгебраических уравнений (ДАУ) индекса 3 в формулировке Лагранжа первого рода. Обычное средство — продифференцировать уравнения связей на уровне положений, что понижает индекс, но разрушает алгебраическую связь, так что для контроля возникающего дрейфа становятся необходимыми стабилизация Баумгарте или проектирование на многообразие связей, а коэффициенты стабилизации приходится подбирать. В настоящей работе рассматривается противоположный выбор: прямое интегрирование формы индекса 3 одностадийной (одношаговой) комплексной схемой Розенброка семейства CROS, без какого-либо понижения индекса, стабилизации или проектирования. На плоских тестовых примерах (бенчмарках) динамики многозвенных механизмов, охватывающих пассивные и активные связи, мы находим, что связи на уровне положений удовлетворяются со вторым порядком по шагу интегрирования на каждом шаге по времени, без систематического дрейфа на протяжении шестидесятисекундных расчётов, и что не требуется ни одного параметра настройки. На том же механизме и с той же матрицей Якоби понижение индекса до первого, интегрируемое неявным методом, оставляет нарушения связей на два порядка величины большими и без обратной связи не сходится по шагу; стабилизация Баумгарте уменьшает нарушение лишь ценой подбора коэффициентов и привнесённой энергии, а проектирование достигает замыкания на уровне ошибок округления только за счёт дополнительной нелинейной коррекции на каждом шаге. Схема A- и L-устойчива, поэтому она допускает жёсткую динамику связей, которая вынуждает явные интеграторы с пониженным индексом использовать очень малые шаги. Поскольку вся правая часть представляет собой единую функцию невязки от состояния, её матрица Якоби доступна посредством автоматического дифференцирования, что устраняет наиболее подверженную ошибкам часть реализации для многозвенных механизмов. Множители Лагранжа получаются из того же решения: мы проверяем множитель при наличии привода против движущего момента в замкнутой форме. Мы также сообщаем о встреченном практическом ограничении: схема деградирует, когда матрица Якоби связей вырождена по рангу, что имеет место для некоторых конфигураций с закреплённым телом.
динамика многозвенных механизмов, дифференциально-алгебраические уравнения, ДАУ индекса 3, методы Розенброка, избыточные координаты, автоматическое дифференцирование, стабилизация связей
1 Введение
Моделирование многозвенных механизмов с замкнутыми кинематическими цепями — параллельных манипуляторов, рычажных механизмов, шагающих машин — наиболее естественно формулируется в избыточных (максимальных) координатах, где каждое твёрдое тело несёт свои собственные положение и ориентацию, а шарниры наложены как алгебраические связи. Альтернатива — вывод уравнений движения в наборе минимальных координат — требует аналитического исключения зависимых координат и даёт выражения, специфичные для конкретного механизма, которые приходится выводить заново всякий раз, когда меняется топология. Для параллельных механизмов с несколькими замкнутыми кинематическими цепями этот вывод составляет преобладающую часть трудозатрат на моделирование (Kuts и др. 2026).
Путь через избыточные координаты приводит, посредством формулировки Лагранжа первого рода, к системе дифференциально-алгебраических уравнений индекса 3:
\[ \begin{cases} \mathbf{M}\ddot{\mathbf{q}} + \mathbf{C}^{\mathsf T}(\mathbf{q})\,\boldsymbol{\lambda} = \mathbf{Q}(\mathbf{q},\dot{\mathbf{q}},t),\\[2pt] \boldsymbol{\Phi}(\mathbf{q},t) = \mathbf{0}, \end{cases} \]
где \(\mathbf{q}\in\mathbb{R}^{n}\) собирает декартовы координаты всех тел, \(\mathbf{M}\) — блочно-диагональная матрица масс, \(\mathbf{C}=\partial\boldsymbol{\Phi}/\partial\mathbf{q}\) — матрица Якоби связей, а \(\boldsymbol{\lambda}\) — множители Лагранжа, которые физически представляют собой реакции связей и усилия приводов.
ДАУ индекса 3 неудобны для стандартных интеграторов ОДУ (Gear и Petzold 1984; Brenan, Campbell, и Petzold 1996). Установившаяся практика состоит в том, чтобы понизить индекс: однократное или двукратное дифференцирование \(\boldsymbol{\Phi}=\mathbf{0}\) даёт связь на уровне скоростей или ускорений, и полученная система индекса 2 или индекса 1 дискретизируется напрямую. Платой является то, что связь по положению больше не обеспечивается, так что численное решение покидает многообразие связей. Стандартными являются два средства:
- стабилизация Баумгарте (Baumgarte 1972), которая заменяет \(\boldsymbol{\Phi}=\mathbf{0}\) демпфированной обратной связью \(\ddot{\boldsymbol{\Phi}}+2\alpha\dot{\boldsymbol{\Phi}}+\beta^{2}\boldsymbol{\Phi}=\mathbf{0}\) и требует коэффициентов, которые обеспечивают компромисс между нарушением связей и жёсткостью;
- проектирование или разделение координат (Eich 1993; Haug 1989), которое повторно ортогонализует состояние на многообразие связей после каждого шага, ценой дополнительного нелинейного решения на каждом шаге.
Оба средства добавляют механизмы, и оба взаимодействуют с контролем ошибки лежащего в основе интегратора.
В настоящей работе изучается прямая альтернатива, заключающаяся в прямом интегрировании системы индекса 3 в том виде, как она есть, одностадийной (одношаговой) комплексной схемой Розенброка семейства CROS (Al’shin и др. 2006), которая A- и L-устойчива и была разработана для жёстких и дифференциально-алгебраических систем.
2 Математическая формулировка
2.1 Уравнения движения в избыточных координатах
Рассмотрим \(n_b\) плоских твёрдых тел. Тело \(i\) имеет массу \(m_i\), момент инерции относительно центра масс \(J_i\) и координаты \(\mathbf{q}_i=(x_i,y_i,\theta_i)^{\mathsf T}\) в абсолютной системе отсчёта. Объединение всех тел даёт \(\mathbf{q}\in\mathbb{R}^{3n_b}\). Уравнения Ньютона–Эйлера без связей имеют вид
\[ \mathbf{M}\ddot{\mathbf{q}}=\mathbf{Q}_{\mathrm{ext}}(\mathbf{q},\dot{\mathbf{q}},t), \qquad \mathbf{M}=\mathrm{blockdiag}\bigl(\mathrm{diag}(m_i,m_i,J_i)\bigr), \]
где \(\mathbf{Q}_{\mathrm{ext}}\) собирает силы тяжести и приложенные силы. Шарниры представляют собой голономные связи \(\boldsymbol{\Phi}(\mathbf{q},t)=\mathbf{0}\in\mathbb{R}^{n_c}\); добавление их реакций через множители Лагранжа даёт приведённое выше ДАУ индекса 3. Поскольку силы связей не совершают работы на виртуальных перемещениях, \(\boldsymbol{\lambda}\) восстанавливается как часть решения.
2.2 Библиотека связей
Формулировка является общей относительно \(\boldsymbol{\Phi}\); необходимы лишь невязка и её матрица Якоби. В таблице 1 перечислены используемые здесь плоские шарниры. Локальные точки тела — это \(\mathbf{p}^{b}_c=(x_c,y_c)^{\mathsf T}\), а \(\mathbf{R}(\theta)\) — матрица поворота на плоскости.
| Шарнир | \(\boldsymbol{\Phi}\) | \(n_c\) |
|---|---|---|
| Закрепление к основанию | \(\mathbf{r}_1-\mathbf{r}^{\mathrm{ref}}\) и \(\theta_1-\theta^{\mathrm{ref}}\) | 3 |
| Вращательная пара (шарнир) | \(\mathbf{r}_i+\mathbf{R}_i\mathbf{p}^{i}_c-\mathbf{r}_j-\mathbf{R}_j\mathbf{p}^{j}_c\) | 2 |
| Поступательная пара (ползун) | \((\mathbf{r}_j-\mathbf{r}_i)\cdot\hat{\mathbf{n}}_i\) и \((\theta_j+\alpha_j)-(\theta_i+\alpha_i)\) | 2 |
| Вращательный позиционный двигатель | \((\theta_j-\theta_i)-u(t)\) | 1 |
| Линейный позиционный привод | \((\mathbf{r}_j-\mathbf{r}_i)\cdot\hat{\mathbf{t}}_i-u(t)\) | 1 |
Здесь \(\hat{\mathbf{n}}_i\) и \(\hat{\mathbf{t}}_i\) — единичные нормаль и касательная оси скольжения на теле \(i\), обе повёрнутые вместе с телом, а \(\alpha_i\) — постоянный локальный угол направления скольжения.
2.3 Активные связи
Позиционные приводы трактуются как дополнительные голономные связи, присоединённые к \(\boldsymbol{\Phi}\) с собственным множителем. Заданная траектория \(u(t)\) делает связь реономной, \(\Phi(\mathbf{q},t)=a(\mathbf{q})-u(t)=0\). Поэтому ограничения приводов и заданные профили входят точно так же, как шарниры, и тот же множитель несёт усилие привода. Соответствующий элемент матрицы Якоби по времени равен \(\partial\Phi/\partial t=-\dot u(t)\), который должен быть задан, чтобы дискретизация оставалась согласованной.
3 Численная схема
3.1 Шаг CROS
Запишем форму первого порядка
\[ \mathbf{M}_{\mathrm{aug}}\dot{\mathbf{z}}=\mathbf{F}(\mathbf{z}), \qquad \mathbf{z}=(\mathbf{q},\dot{\mathbf{q}},\boldsymbol{\lambda},t)^{\mathsf T}, \]
где \(\mathbf{M}_{\mathrm{aug}}\) расширяет матрицу масс единичным блоком для кинематических координат, нулями в строках множителей и единственной \(1\) для временной компоненты; невязка \(\mathbf{F}\) содержит уравнения количеств движения, уравнения связей \(\boldsymbol{\Phi}(\mathbf{q},t)\) и \(\dot t=1\). Шаг CROS имеет вид
\[ \boldsymbol{\zeta}=\Bigl(\mathbf{M}_{\mathrm{aug}} -\tfrac{1+\mathrm{i}}{2}\,h\,\mathbf{J}(\mathbf{z}_k)\Bigr)^{-1}\mathbf{F}(\mathbf{z}_k), \qquad \mathbf{z}_{k+1}=\mathbf{z}_k+h\,\mathrm{Re}\,\boldsymbol{\zeta}, \]
где \(\mathbf{J}=\partial\mathbf{F}/\partial\mathbf{z}\).
Здесь важны два свойства.
A- и L-устойчивость. Поскольку на каждом шаге берётся вещественная часть стадии, рекурсия остаётся вещественной, и функция устойчивости схемы на вещественной оси равна
\[ R(x)=1+\operatorname{Re}\frac{x}{1-\tfrac{1+\mathrm{i}}{2}x} =\frac{2}{x^{2}-2x+2}, \qquad R(x)\le 1\ \text{ при }\ x\le 0, \qquad R(-\infty)=0. \]
Схема, как и всё семейство CROS (Al’shin и др. 2006), A-устойчива и L-устойчива: бесконечно жёсткие затухающие моды не сохраняются, а подавляются полностью. Именно это снимает ограничение на шаг интегрирования, от которого страдает явный интегратор с пониженным индексом при жёсткой динамике связей ДАУ индекса 3.
Порядок на дифференциальной части. Разложение той же функции по степеням \(x\) даёт
\[ R(x)=1+x+\tfrac{1}{2}x^{2}+O(x^{4}), \]
то есть на вещественной оси рекурсия согласуется с \(\exp(x)\) вплоть до второго порядка (кубический член обращается в нуль тождественно). Это согласуется с тем, что мы наблюдаем на ДАУ: как показывают разделы 4 и 5, и ошибка траектории, и невязка связей сходятся со вторым порядком. Порядок дискретизации Розенброка для ДАУ определяется, кроме того, индексом и тем, как обрабатываются алгебраические строки (Lubich и Roche 1992; Hairer и Wanner 1996).
3.2 Автоматическое дифференцирование невязки
Поскольку вся правая часть собирается как единая функция невязки от \(\mathbf{z}\), матрица Якоби \(\mathbf{J}\) получается прямым (forward) автоматическим дифференцированием (Revels, Lubin, и Papamarkou 2016). Для невязки, построенной из тригонометрических выражений, подобных приведённым в таблице 1, это точно с рабочей точностью и устраняет необходимость вручную программировать \(\partial\mathbf{F}/\partial\mathbf{z}\), что для механизма с несколькими замкнутыми контурами является наиболее подверженной ошибкам частью реализации. Добавление нового шарнира означает написание только его невязки связи.
4 Результаты: выполнение связей
Все результаты получены с помощью описанной выше эталонной реализации. Эталонное решение, используемое для оценки ошибки, всегда строится другим путём: ДАУ индекса 3 проектировалось на ядро (нуль-пространство) матрицы Якоби связей, что давало ОДУ в минимальных координатах, которое интегрировалось классическим RK4 с шагом на два порядка величины меньшим. Этот эталон не имеет общего кода с ядром CROS.
Тестовый пример: массивное основание (тело-основание) и два однородных звена, \(m=1\) кг, \(L=1\) м, \(J=mL^{2}/12\), отпущенные из \(q_1=0.3\) и \(q_2=-1.0\) рад с нулевой скоростью. Механизм имеет \(n_b=3\) тела и \(n_c=7\) связей (три для закрепления к основанию, по две для каждого шарнира), так что \(n=6\cdot 3+7+1=26\).
4.1 Невязка связей при индексе 3
В таблице 2 приведена невязка уравнений связей \(\max|\boldsymbol{\Phi}(\mathbf{q})|\) вдоль траектории.
| Шарнир | \(h=10^{-2}\) | \(h=10^{-3}\) | \(h=10^{-4}\) |
|---|---|---|---|
| Закрепление к основанию (неподвижное) | \(6.41\cdot10^{-24}\) | \(1.11\cdot10^{-25}\) | \(1.16\cdot10^{-27}\) |
| Шарнир 1 | \(1.51\cdot10^{-3}\) | \(1.68\cdot10^{-5}\) | \(1.70\cdot10^{-7}\) |
| Шарнир 2 | \(6.01\cdot10^{-3}\) | \(7.70\cdot10^{-5}\) | \(7.67\cdot10^{-7}\) |
Два наблюдения. Во-первых, закрепление к основанию удовлетворяется с машинной точностью. Это не случайно: его строки связей содержат только единичные блоки, поэтому решение линейной системы воспроизводит их точно. Строки шарниров связывают положения и углы и удовлетворяются со вторым порядком, причём невязка убывает в примерно 90 и 100 раз на декаду уменьшения шага. Во-вторых, невязка шарнира не равна нулю, так что схема не обеспечивает связь точно; она обеспечивает её с точностью дискретизации.
В таблице 3 приведено геометрическое замыкание шарнира — расстояние между точкой шарнира, вычисленной по одному звену и по другому, — оценённое по «сырому» состоянию без использования множителей, на \(t\in[0,20]\) с.
| \(h\) | макс. замыкание, м | отношение |
|---|---|---|
| \(8\cdot10^{-3}\) | \(3.9631\cdot10^{-3}\) | – |
| \(4\cdot10^{-3}\) | \(1.1537\cdot10^{-3}\) | 3.44 |
| \(2\cdot10^{-3}\) | \(3.0145\cdot10^{-4}\) | 3.83 |
| \(10^{-3}\) | \(7.7147\cdot10^{-5}\) | 3.91 |
| \(5\cdot10^{-4}\) | \(1.9317\cdot10^{-5}\) | 3.99 |
Отношение сходится к 4, что является признаком второго порядка.
4.2 Отсутствие систематического дрейфа
Именно это свойство мотивирует весь подход. Та же система интегрировалась до 60 с при \(h=10^{-3}\), и нарушение исследовалось окно за окном (таблица 4).
| Временное окно, с | макс. \(\|\boldsymbol{\Phi}\|\), все шарниры | макс. геометрическое замыкание, м |
|---|---|---|
| 0 – 5 | \(6.1531\cdot10^{-5}\) | \(6.2629\cdot10^{-5}\) |
| 5 – 15 | \(7.6950\cdot10^{-5}\) | \(7.7147\cdot10^{-5}\) |
| 15 – 30 | \(7.4816\cdot10^{-5}\) | \(7.5192\cdot10^{-5}\) |
| 30 – 45 | \(6.5455\cdot10^{-5}\) | \(6.5458\cdot10^{-5}\) |
| 45 – 60 | \(6.1290\cdot10^{-5}\) | \(6.1328\cdot10^{-5}\) |
Нарушение возрастает до плато, определяемого шагом интегрирования, и затем остаётся на нём. Оно не растёт со временем, и состояние остаётся ограниченным (\(\max|\mathbf{z}|=187.5\) за 60 с). Следовательно, ни член Баумгарте, ни шаг проектирования не нужны, чтобы удерживать алгебраические связи под контролем.
4.3 Точность траектории
Положение рабочего органа (выходного звена) сравнивалось с независимым эталоном на \(t\in[0,4]\) с (таблица 5).
| \(h\) | макс. ошибка, м | отношение |
|---|---|---|
| \(8\cdot10^{-3}\) | \(4.7819\cdot10^{-2}\) | – |
| \(4\cdot10^{-3}\) | \(1.4376\cdot10^{-2}\) | 3.33 |
| \(2\cdot10^{-3}\) | \(6.3441\cdot10^{-3}\) | 2.27 |
| \(10^{-3}\) | \(2.0107\cdot10^{-3}\) | 3.16 |
| \(5\cdot10^{-4}\) | \(5.5519\cdot10^{-4}\) | 3.62 |
Отношения приближаются к 4, что снова указывает на второй порядок; разброс на промежуточных шагах объясняется тем, что максимум нормы достигается в разные моменты времени при разных шагах, что нормально для нелинейной двухзвенной цепи.
5 Результаты: сравнение с понижением индекса
Чтобы изолировать роль индекса, тот же самый механизм интегрировался с той же самой матрицей Якоби, полученной автоматическим дифференцированием, прямым путём и четырьмя способами понижения индекса:
- прямая форма индекса 3 со схемой CROS — метод настоящей работы;
- форма индекса 1, в которой каждая связь продифференцирована дважды и на уровне ускорений налагается \(\ddot{\boldsymbol{\Phi}}=\mathbf{0}\), интегрируемая неявной схемой
Trapezoid(SDIRK, A-устойчивая); - та же форма индекса 1, интегрируемая неявной схемой
Rodas4(метод Розенброка, L-устойчивый); - форма индекса 1 с обратной связью Баумгарте \(\ddot{\boldsymbol{\Phi}}+2\alpha\dot{\boldsymbol{\Phi}}+\beta^{2}\boldsymbol{\Phi}=\mathbf{0}\);
- форма индекса 1 с проектированием состояния на многообразие связей после каждого шага.
Понижение проводится именно до индекса 1: при однократном дифференцировании (индекс 2) связь по положению исчезает из формулировки целиком, и никакое уменьшение шага её не возвращает. Форма индекса 1 при этом остаётся ДАУ — множители Лагранжа входят в неё как алгебраические неизвестные, — поэтому явные методы для неё непригодны в принципе, и во всех вариантах используются неявные схемы. В реализации множители исключаются внутри правой части точечным решением системы Карруша–Куна–Таккера с матрицей приложения реакций \(\mathbf{B}=\partial\mathbf{F}/\partial\boldsymbol{\lambda}\), а \(\boldsymbol{\Phi}\), \(\dot{\boldsymbol{\Phi}}\) и конвективная часть \(\ddot{\boldsymbol{\Phi}}\) получаются автоматическим дифференцированием по направлению потока; таким образом, одна и та же невязка связей обслуживает и форму индекса 3, и форму индекса 1, а отдельная матрица Якоби связей не программируется нигде.
В таблице 6 приведено максимальное геометрическое замыкание второго шарнира на \(t\in[0,5]\) с.
| \(h\) | CROS, индекс 3 | индекс 1, Trapezoid | индекс 1, Rodas4 | индекс 1 + Баумгарте, Rodas4 | индекс 1 + проектирование |
|---|---|---|---|---|---|
| \(4\cdot10^{-3}\) | \(9.875\cdot10^{-4}\) | \(7.753\cdot10^{-4}\) | \(1.527\cdot10^{-2}\) | \(1.151\cdot10^{-3}\) | \(5.79\cdot10^{-14}\) |
| \(2\cdot10^{-3}\) | \(2.497\cdot10^{-4}\) | \(2.919\cdot10^{-2}\) | \(7.849\cdot10^{-3}\) | \(5.842\cdot10^{-4}\) | \(6.04\cdot10^{-14}\) |
| \(10^{-3}\) | \(6.263\cdot10^{-5}\) | \(4.810\cdot10^{-2}\) | \(4.003\cdot10^{-3}\) | \(2.963\cdot10^{-4}\) | \(5.62\cdot10^{-14}\) |
Прямая форма индекса 3 даёт второй порядок с отношениями, стремящимися к 4. Форма индекса 1 без обратной связи на два порядка хуже и, что важнее, не сходится по шагу: позиционная связь в ней не наложена, и нарушение определяется не ошибкой дискретизации, а возбуждением паразитных мод, которых в исходной форме индекса 3 просто нет.
Выбор неявного метода при этом не безразличен. A-устойчивый Trapezoid эти моды не гасит: при измельчении шага растёт число шагов, вместе с ним накапливается возбуждение, и нарушение растёт (\(7.8\cdot10^{-4}\to4.8\cdot10^{-2}\)). L-устойчивый Rodas4 их подавляет, и нарушение убывает с шагом монотонно (\(1.5\cdot10^{-2}\to4.0\cdot10^{-3}\)). Иначе говоря, пониженной форме нужен не просто неявный, а жёстко-устойчивый метод, тогда как прямая форма индекса 3 довольствуется A-устойчивой схемой CROS.
| \(2\alpha\) | \(\beta^{2}\) | Trapezoid: замыкание / дрейф энергии | Rodas4: замыкание / дрейф энергии |
|---|---|---|---|
| 0 | 0 | \(4.81\cdot10^{-2}\) / \(1.8\cdot10^{-2}\) | \(4.00\cdot10^{-3}\) / \(1.1\cdot10^{-3}\) |
| 2 | 1 | \(7.91\cdot10^{-3}\) / \(2.4\cdot10^{-2}\) | \(7.64\cdot10^{-4}\) / \(1.4\cdot10^{-3}\) |
| 10 | 25 | \(4.51\cdot10^{-3}\) / \(2.5\cdot10^{-2}\) | \(2.96\cdot10^{-4}\) / \(2.0\cdot10^{-3}\) |
| 40 | 400 | \(3.35\cdot10^{-3}\) / \(5.5\cdot10^{-2}\) | \(1.99\cdot10^{-4}\) / \(4.1\cdot10^{-3}\) |
| 100 | 2500 | \(2.02\cdot10^{-3}\) / \(5.6\cdot10^{-2}\) | \(1.09\cdot10^{-4}\) / \(5.8\cdot10^{-3}\) |
Обратная связь Баумгарте уменьшает нарушение связей монотонно с ростом коэффициентов, но вносит в систему работу, которой в исходной задаче нет: дрейф механической энергии растёт вместе с коэффициентами. Компромисс между нарушением связей и искажением динамики приходится выбирать подбором, и он зависит от шага, — именно этого выбора прямая форма индекса 3 не требует. С неявным методом расходимости, которая возникала бы у явной схемы при больших коэффициентах, не наблюдается: ценой больших коэффициентов становится не потеря устойчивости, а привнесённая энергия.
Проектирование даёт замыкание на уровне ошибок округления (\(\sim6\cdot10^{-14}\)) и не зависит от шага, но требует нелинейной коррекции состояния на каждом шаге, а сама коррекция тоже возмущает динамику: её дрейф энергии (\(5.3\cdot10^{-4}\) при \(h=10^{-3}\)) на порядок больше, чем у прямой схемы индекса 3 (\(5.1\cdot10^{-5}\)).
Прямая схема индекса 3 не требует ни коэффициентов, ни коррекций: ограниченное, контролируемое шагом нарушение порядка \(O(h^{2})\), получаемое за одно вычисление невязки и одно решение линейной системы на каждом шаге, и никакого параметра для настройки.
6 Результаты: замкнутая кинематическая цепь с приводом и пружиной
Второй пример — кривошипно-ползунный механизм с замкнутой кинематической цепью: массивное основание, кривошип, шатун и ползун. Основание закреплено, кривошип соединён с ним шарниром, шатун замыкает цепь между кривошипом и ползуном, ползун удерживается на прямолинейной направляющей, а вращательный позиционный привод задаёт абсолютный угол кривошипа \(\theta=\pi t\). В механизме четыре тела и двенадцать связей: три на закрепление основания, по две на каждый из трёх шарниров, две на поступательную пару и одна на привод. Подвижность механизма нулевая: движение полностью задано приводом, а содержательным результатом являются усилия. Вдоль направляющей действует линейная пружина жёсткости \(k=20\) Н/м — силовой элемент, не вносящий дополнительных связей.
Начальные скорости выбираются согласованными: состояние проектируется на многообразие \(\{\boldsymbol{\Phi}=\mathbf{0},\ \dot{\boldsymbol{\Phi}}=\mathbf{0}\}\), так что заданная угловая скорость кривошипа \(\omega=\pi\) рад/с выполняется с первого шага и разгонный переходный процесс в измерения не попадает.
В таблице 8 приведены максимальное геометрическое замыкание замкнутого контура, ошибка отслеживания заданного угла и невязка баланса мощности.
| \(h\) | макс. замыкание, м | отношение | \(\max|\theta-\pi t|\), рад | \(\Delta E-\int\lambda\omega\,dt\), Дж | относительная невязка |
|---|---|---|---|---|---|
| \(4\cdot10^{-3}\) | \(5.92\cdot10^{-5}\) | – | \(3.4\cdot10^{-13}\) | \(29.5\) | \(7.4\cdot10^{-1}\) |
| \(2\cdot10^{-3}\) | \(1.48\cdot10^{-5}\) | 4.00 | \(1.0\cdot10^{-12}\) | \(14.8\) | \(3.7\cdot10^{-1}\) |
| \(10^{-3}\) | \(3.71\cdot10^{-6}\) | 3.99 | \(1.0\cdot10^{-12}\) | \(7.4\) | \(1.9\cdot10^{-1}\) |
| \(5\cdot10^{-4}\) | \(9.59\cdot10^{-7}\) | 3.86 | \(1.6\cdot10^{-12}\) | \(3.7\) | \(9.2\cdot10^{-2}\) |
Замыкание контура убывает со вторым порядком (отношения 4.00, 3.99, 3.86), хотя цепь замкнута и содержит четыре тела: матрица Якоби строится автоматическим дифференцированием всей невязки целиком, поэтому добавление замыкающего контура, пружины и привода не меняет в реализации ни одного элемента. Заданный угол отслеживается на уровне ошибок округления (\(10^{-12}\) рад), поскольку он наложен алгебраической строкой.
Множители в замкнутой цепи проверялись двумя независимыми способами. Во-первых, в статике: при постоянном заданном угле кривошипа механизм приходит в состояние покоя, и множитель привода обязан равняться производной потенциальной энергии по углу кривошипа \(dU/d\theta\), взятой вдоль многообразия связей (гравитация плюс энергия пружины). Численно это выполняется с относительной ошибкой \(5.4\cdot10^{-3}\) на диапазоне \(\theta\in[0.4,\,2.2]\) рад. Во-вторых, в движении: работа привода \(\int\lambda\omega\,dt\) сходится к изменению механической энергии \(\Delta E\) (таблица 8), причём относительная невязка убывает примерно вдвое при каждом измельчении шага (\(0.74\to0.09\)). Обе проверки подтверждают, что множитель привода несёт движущий момент и в замкнутой цепи с силовым элементом.
Отдельно проверялась согласованность самих реакций: матрица приложения реакций \(\mathbf{B}=\partial\mathbf{F}/\partial\boldsymbol{\lambda}\) должна совпадать с \((\partial\boldsymbol{\Phi}/\partial\mathbf{q})^{\mathsf T}\), иначе силы связей совершают паразитную работу, а множители перестают быть физическими силами, хотя кинематика остаётся верной. Для используемой библиотеки связей это тождество выполняется с машинной точностью (проверено и для двухзвенной цепи, и для кривошипно-ползунного механизма); такая проверка дешева и должна быть частью отладки любой реализации. # Обсуждение и ограничения
Вырожденность по рангу. Подход деградирует, когда матрица Якоби невязки вырождена по рангу. Мы встретили это при построении тестового примера, в котором первое звено имело конечную массу и одновременно было закреплено неподвижным шарниром; в этой конфигурации матрица Якоби была вырождена на два, множители приобрели паразитную компоненту порядка \(1/h\), и состояние разошлось. Замена закреплённого первого звена массивным основанием (телом-основанием) полностью устранила проблему, и все приведённые выше результаты используют эту корректно поставленную конфигурацию. В примерах из библиотеки вырожденность также проявляется (например, ранг 23 из 26 для цепи основание–кривошип–коромысло), но она безобидна, поскольку вырожденные направления не возбуждаются. Практическое следствие состоит в том, что матрицу Якоби невязки следует контролировать при промышленном использовании и что закрепление лёгкого тела напрямую к основанию — это конфигурация, которой следует избегать.
Второй порядок, а не точное выполнение. Связи по положению удовлетворяются с точностью \(O(h^{2})\), а не с машинной точностью. Этого достаточно для устранения дрейфа, но это не то же самое, что точно спроектированное решение, и достижимый геометрический допуск определяется шагом интегрирования.
Плотная матрица Якоби и плотная линейная алгебра. Именно автоматическое дифференцирование делает реализацию короткой, но оно же делает матрицу Якоби плотной, так что стоимость шага растёт кубически с числом неизвестных. Аналитические и разреженные матрицы Якоби с разреженной факторизацией — очевидный следующий шаг; тем более что большинство невязок связей локальны для нескольких тел.
Фиксированный шаг. В настоящей реализации используется постоянный шаг, и управление шагом отсутствует, поэтому шаг приходится выбирать консервативно. Для механизмов с быстро меняющимися профилями приводов это практическое неудобство. Добавление управления шагом потребовало бы встроенной оценки ошибки для этой схемы на ДАУ, которая в изучаемой здесь реализации недоступна; для ДАУ такую оценку нельзя заимствовать без изменений из теории Розенброка для индекса 0, поскольку порядок дискретизации здесь определяется ещё и индексом, и трактовкой алгебраических строк (Lubich и Roche 1992; Hairer и Wanner 1996).
A- и L-устойчивость. Поскольку \(R(-\infty)=0\), бесконечно жёсткие затухающие моды подавляются полностью, и жёсткость связей не ограничивает шаг интегрирования со стороны устойчивости.
Область применения. Изложение двумерное, а библиотека связей охватывает голономные шарниры и позиционные приводы. Односторонние контакты, трение и неголономные связи находятся за пределами настоящей работы и потребовали бы аппарата ограничений-неравенств или комплементарности.
7 Заключение
Мы исследовали прямое интегрирование дифференциально-алгебраических уравнений индекса 3 динамики многозвенных механизмов в избыточных координатах одностадийной (одношаговой) комплексной схемой Розенброка семейства CROS, без понижения индекса, без стабилизации Баумгарте и без проектирования.
На плоских тестовых примерах (бенчмарках), охватывающих как пассивные, так и активные связи, алгебраические связи на уровне положений удовлетворяются со вторым порядком по шагу интегрирования на каждом шаге, а нарушение достигает зависящего от шага плато без систематического роста на протяжении 60 с. На том же механизме и с той же матрицей Якоби понижение индекса до первого, интегрируемое неявным методом, оставляло нарушения связей на два порядка величины большими и без обратной связи не сходилось по шагу, стабилизация Баумгарте уменьшала нарушение лишь ценой подбора коэффициентов и привнесённой энергии, а проектирование достигало замыкания на уровне ошибок округления лишь ценой дополнительной нелинейной коррекции на каждом шаге. Прямой путь через индекс 3 вообще не требовал параметра настройки.
Множители Лагранжа, полученные из того же решения, несут физический смысл: в замкнутой кинематической цепи с приводом и пружиной множитель привода совпадает с производной потенциальной энергии по углу кривошипа, а его работа — с изменением механической энергии механизма.
Практическая привлекательность подхода — его низкая стоимость реализации. Вся правая часть представляет собой единую функцию невязки от состояния, поэтому матрица Якоби следует из автоматического дифференцирования, добавление шарнира означает написание только его невязки, а ограничения приводов и заданные траектории входят так же, как шарниры. Для жёстких механизмов с замкнутыми контурами, где приемлема точность второго порядка при умеренном фиксированном шаге, это выгодный компромисс между трудозатратами на моделирование и достигаемой точностью. За пределами этой области — в частности там, где матрица Якоби невязки вырождена по рангу или где требуется разреженная и масштабируемая формулировка, — подход следует применять с осторожностью.
8 Благодарности
Автор благодарит разработчиков пакетов экосистемы Julia, использованных в эталонной реализации, в частности ForwardDiff.jl, StaticArrays.jl и Makie.jl, и признаёт вклад литературы по схеме CROS, без которой изучаемая здесь схема не была бы выявлена.