Метод оптимального управления экзоскелетом нижних конечностей с упругими элементами

Введение. Экзоскелеты нижних конечностей (ЭНК) вызывают всё больший интерес, обусловленный необходимостью решения глобальных проблем здравоохранения — старения населения и роста числа нервно-мышечных травм [1]. Эти устройства призваны обеспечить эффективные решения для поддержки и улучшения двигательных функций человека, например, для помощи при ходьбе [2], реабилитации и компенсации потери равновесия [3], тем самым повышая независимость и качество жизни.

ЭНК часто проектируются без упругих элементов из‑за возросшей сложности процесса стабилизации во время движения и влияния дополнительных факторов, возникающих при использовании упругих элементов [4]. С другой стороны, для этих экзоскелетов могут требоваться упругие элементы, которые улучшают способность конструкции адаптироваться к неровным поверхностям. Упругие элементы могут быть установлены в области голеностопного сустава [5] или использоваться для его полной замены [6].

Управление ЭНК является сложной задачей. Ранее были предложены различные подходы к её решению. Разработка эффективного метода управления зависит от множества факторов, поддерживая постоянную актуальность исследований в этой области. Для управления ЭНК существуют разные методы: адаптивный метод управления [1][7], робастные (устойчивые) методы [8][9] и метод оптимального управления [10][11]. Несмотря на эффективность первых двух методов, третий является наиболее удачным. Метод оптимального управления учитывает не только повышение устойчивости и эффективности управления ЭНК в динамических и непредсказуемых условиях, но и позволяет снизить энергозатраты и расход ресурсов системы управления [10][11]. Однако интеграция упругих элементов в стопы ЭНК приводит к осложнению управления ими. Кроме того, учет внешних возмущений при управлении ЭНК сопровождается новыми проблемами в обеспечении динамической устойчивости системы управления ЭНК.

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

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

Для достижения поставленной цели были поставлены следующие задачи:

• разработка математической модели перевёрнутого маятника с упругими элементами в стопах;
• разработка метода оптимального управления ЭНК с упругими элементами;
• учет влияния внешних возмущений (гауссовского белого шума);
• проведение численного моделирования в среде Wolfram Mathematica для двух типов управлений: терминальное управление и управление с обратной связью;
• анализ полученных результатов исследования переходных процессов основных параметров динамики ЭНК с упругими элементами.

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

Характеристики упругих элементов и параметры модели следующие: k_y, k_z — коэффициенты жёсткости горизонтальной и вертикальной пружин, соответственно. Их значения — k_y ∈ [ 500, 1 500] Н/м [13], k_z ∈ [ 7 000, 20 000] Н/м [14]. c_y, c_z — коэффициенты горизонтального и вертикального демпфирования, соответственно, c_y ∈ [ 30, 200 Нс/м], c_z ∈ [ 500–2 000 Нс/м] [15]. — момент инерции экзоскелета относительно центра масс, кгм²; m = 70 — масса экзоскелета с пациентом, кг; h = 1 — длина экзоскелета до центра масс, м; μ_y = 0,7 и μ_z = 0 — коэффициенты трения; N_y = 0 и N_z = m g — нормальные силы, Н; g = 9,8 — ускорение свободного падения, м/с².

Методы. Динамика ЭНК описывается уравнениями Лагранжа второго рода в общем виде²:

(1)

где L = T – V — функция Лагранжа; T — кинетическая энергия системы; V — потенциальная энергия системы; Q — обобщенные силы; q = (θ, y_b, z_b)^T — вектор обобщенных координат; θ — угол поворота звена экзоскелета, отсчитываемый от горизонтальной поверхности (параллельной опорной плоскости) в направлении против часовой стрелки; y_b, z_b — горизонтальное и вертикальное перемещение основания.

Координаты центра масс маятника определяются следующими уравнениями:

(2)

Кинетическая энергия маятника задается по формуле:

(3)

где v — скорость центр масс, равен .

Потенциальная энергия систем описываются следующим уравнением:

(4)

Обобщённые силы демпфирования, применяемые с помощью функции диссипации Рэлея для моделирования линейного демпфирования, записываются в виде [15]:

(5)

Модель трения Кулона задаётся следующими уравнениями [15]:

(6)

где α >> 1 — параметр регуляризации (для аппроксимации разрывной функции sign(v)), в данной работе выбран α = 100.

После промежуточных вычислений, уравнения (2–6) принимают вид:

(7)

где — внешние обобщенные силы.

В окрестности вертикального равновесия (θ ≈ 0) справедливы допущения: sin θ ≈ θ, cos θ = 1, z_b ≈ –mg/k_z. Уравнения (7) линеаризуются с учетом малости и при малых скоростях:

(8)

При предположениях о малых углах вертикальное перемещение z становится динамически независимым от горизонтального y и углового θ состояний. Его поведение сводится к гармоническому осциллятору, который может анализироваться отдельно [12].

Для синтеза управления представим систему в канонической форме пространства состояний с учетом внешних возмущений и шумов измерений. Задав вектор состояния , запишем систему в форме Коши:

(9)

где A — матрица состояния; B — матрица управления; u(t) — вектор управляющих сигналов; D — матрица возмущений; w(t) — вектор внешних возмущений; C — матрица измерений; v(t) — вектор внешних возмущений.

Основная задача регулятора заключается в переводе динамической системы из начального состояния x(t₀) в заданное конечное состояние x(t₁) за определенное время t₁.

Грамиан управляемости W характеризует способность системы достигать произвольных состояний за конечное время T и определяется по формуле [12]:

(10)

где Ф(t₀, t) = e^{A(t 0 – t)} — матричная экспонента.

Система является полностью управляемой на [t₀, t₁] тогда и только тогда, когда W(t₁, t₀) обратима. Если Грамиан обратим, любое состояние x_T ∈ ℝⁿ может быть достигнуто с помощью соответствующего управления u(t) [12].

Управление, которое минимизирует квадратичную функцию энергозатрат и переводит систему из состояния x(t₀) в состояние x(t₁), имеет вид [12]:

(11)

После промежуточных вычислений квадратичную функцию энергозатрат можно рассчитать с помощью выражения (12):

(12)

Управление, заданное уравнением (11), является программным (разомкнутым) и зависит от времени. Для обеспечения асимптотической стабилизации системы в начало координат при произвольных начальных условиях необходим синтез закона управления замкнутого контура, основанного на обратной связи по состоянию в реальном времени. Задача управления формулируется в виде минимизации квадратичной функции энергозатрат [12]:

(13)

где R — положительно определенная матрица; Q — положительно полуопределенная матрица.

Оптимальное управление, обеспечивающее минимизацию энергозатрат, может быть реализовано в виде отрицательной обратной связи с переменными параметрами:

(14)

где k(t) = R^–1B^TP(t); k(t) — матрица коэффициента обратной связи; P(t) — решение уравнения Риккати.

Для нахождения P(t) при t → ∞ необходимо решить алгебраическое уравнение Риккати, которое задается следующим образом:

(15)

Система с управляющим сигналом u(t) описывается уравнением:

(16)

На практике измерение всех состояний невозможно из-за ограниченного числа датчиков, шумов и погрешностей измерений [12]. Для оценки вектора состояния по входам и выходам используется оптимальный фильтр Калмана.

Система с управляющим сигналом u(t) и фильтр Калмана для определения вектора состояния имеет следующий вид [12]:

(17)

где — оценочное состояние; C — выход наблюдателя; L — матрица усиления Калмана.

При наличии внешних возмущений w(t) и шума измерений v(t) с нулевыми математическими ожиданиями M[v(t)] = 0, M[w(t)] = 0 и ковариационными матрицами s_w и s_v, соответственно, матрица коэффициента усиления L определяется из решения алгебраического уравнения Риккати (18) и соответственно найдется стационарный фильтр Калмана:

(18)

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

(19)

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

(20)

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

Результаты исследования. В работе были исследованы переходные процессы управляющего момента, силы, смещения и угла при переходе ЭНК с упругими элементами из неустойчивого положения в положение вертикального равновесия. Результаты численного моделирования приведены на рис. 2–6 при трех начальных условиях: Y₁ = [ 1/10, 0, 0 ,0], Y₂ = [ 1/10, 0, 1/10 ,0], Y₃ = [ 1/10 0 1/10 0 0 1/10 0 1/20 0].

Результаты терминального управления рассчитаны для времени стабилизации 0,6 с. На рис. 2 представлены кривые управляющих моментов разомкнутой системы при начальных условиях Y₁ и целевом положении в состоянии вертикального равновесия. Синяя кривая соответствует управляющему моменту M, оранжевая — управляющей силе F_y.

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

На рис. 3 представлены кривые переходных процессов основных параметров динамической системы ЭНК (θ, y_b) при начальных условиях Y₁. Синяя кривая представляет собой переходный процесс угла отклонения θ, а оранжевая — изменение координаты по оси ординат y_b в течение заданного времени стабилизации 0,6 с при терминальном управлении.

Результаты моделирования на рис. 4–6 получены с применением типа управления с обратной связью, которое позволяет провести поиск оптимальной промежуточной траектории движения. Для исследования влияния жесткости упругих элементов (изменения коэффициента жесткости пружины k_y) на устойчивость динамической системы ЭНК при минимизации энергозатрат были рассмотрены переходные процессы основных параметров динамической системы ЭНК (рис. 4) при начальных условиях Y₂ и следующих параметрах системы с обратной связью: c_y = 100, k_z = 20 000, Q = eye(4) ⋅ 10³, R = eye(2).

Оранжевая линия представляет собой зависимость y_b, а зеленая линия — θ. Значения коэффициента жесткости были выбраны следующие: k_y = [ 700, 1 000, 1 500] Н/м. Переходные процессы представлены соответственно на рис. 4 a, б, в.

Для исследования влияния демпфирования упругих элементов (изменения коэффициента демпфирования c_y) на устойчивость динамической системы ЭНК при минимизации энергозатрат были рассмотрены переходные процессы основных параметров динамической системы ЭНК (θ, y_b) (рис. 5) при начальных условиях Y₂ и k_y =1 000 с обратной связью. Оранжевая линия представляет собой зависимость y_b, а зеленая линия — θ. Значения коэффициента демпфирования составляют c_y = [ 0, 50, 100] Нс/м, а переходные процессы изображены соответственно на рис. 5 a, б и в.

Это исследование также учитывает влияние шума путем добавления фильтра Калмана. На рис. 6 показаны кривые переходного процесса при начальных условиях Y₃ = [ 1/10 0 1/10 0 0 1/10 0 1/20 0] и ковариационные матрицы s_w = eye(1) и s_v = 0,1 eye(2). Оранжевая линия иллюстрирует переходный процесс y_b, а зеленая линия — переходный процесс θ.

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

Из таблицы 1 видно, что по мере увеличения времени, необходимого для достижения целевого устойчивого состояния, энергозатраты пропорционально уменьшаются. Данная зависимость согласуется с динамикой управляющих сигналов: с увеличением допустимого времени стабилизации амплитуды управляющих моментов M и F_y уменьшаются, что, согласно формулам (12)–(14), снижает значения энергозатрат. При этом можно сделать вывод о том, что увеличение времени стабилизации (в пределах 0,1–0,4 секунды) снижает энергозатраты на 98 %.

Кривые на рис. 3 демонстрируют переходный процесс длительностью 0,6 с до достижения положения вертикального равновесия, что согласуется с принципами терминального управления (достижение целевого устойчивого состояния за определенное время стабилизации).

Кривые на рис. 4 иллюстрируют влияние изменения коэффициента жесткости на переходные процессы основных параметров динамической системы ЭНК (θ, y_b). Изменение координаты по оси ординат y_b при низком значении коэффициента жесткости характеризуется значительным перерегулированием и медленной стабилизацией, а рост коэффициента жесткости исключает перерегулирование и ускоряет стабилизацию системы. Напротив, изменения угла θ демонстрирует обратную зависимость: при низком значении коэффициента жесткости наблюдается плавный переходный процесс без перерегулирования и короткое время стабилизации системы, в то время как увеличение значения коэффициента жесткости до 1 500 вызывает перерегулирование по углу, несмотря на ускоренную стабилизацию. Это противоречие подчеркивает конкурирующую динамику между y_b и θ. При коэффициенте жесткости k_y =1 000 (рис. 4 б) достигается оптимальный компромисс: минимизация перерегулирования по координате y_b, сохранение стабильности угла θ и быстрая сходимость. Данный режим обеспечивает сбалансированную работу системы.

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

Кривые на рис. 5 показывают влияние коэффициента демпфирования. При низком значении коэффициента демпфирования наблюдается значительная вибрация системы перед выходом на установившийся режим, сопровождаемая перерегулированием. Увеличение c_y снижает амплитуду колебаний и устраняет перерегулирование, однако приводит к росту энергозатрат. С целью обеспечения баланса между устойчивостью и значением квадратичной функции энергозатрат было выбрано значение коэффициента c_y = 100, что гарантирует устойчивость системы управления ЭНК с упругими элементами.

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

Полученные результаты показывают следующую картину. Во-первых, применение открытого управления как базового этапа позволяет минимизировать энергозатраты на стадии вывода системы из неустойчивого состояния в целевое вертикальное равновесие в заданное время стабилизации. Во-вторых, переход к замкнутому контурному управлению на основе обратной связи по состоянию обеспечивает асимптотическую устойчивость и устойчивость переходных режимов за счет решения уравнения Риккати и реализации эффективного фильтра Калмана. В-третьих, численное моделирование выявило важное влияние параметров упругих элементов на динамику системы: увеличение жесткости пружин уменьшает перерегулирование по углу θ и ускоряет стабилизацию, но может приводить к сложностям на пути к устойчивой траектории; увеличение демпфирования снижает колебания и уменьшает перерегулирование, однако повышает энергозатраты. Найден баланс между стабильностью и энергозатратами, который достигается при конкретных значениях k_y и c_y (в примерах k_y = 1 000 Н/м, c_y = 100 Н·с/м).

Исходя из предыдущего обсуждения и анализа полученных в статье результатов можно сказать, что разработанному методу оптимального управления ЭНК с упругими элементами удалось обеспечить устойчивость системы управления данного экзоскелета за определенное время стабилизации и при минимальных энергозатратах.

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

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