Human-in-the-Loop Optimization of Exoskeleton Assistance Via Online Simulation of Metabolic Cost

Human-in-the-Loop Optimization of Exoskeleton Assistance Via Online Simulation of Metabolic Cost

文章来源https://ieeexplore.ieee.org/abstract/document/9698243

本文介绍了一种人在环的外骨骼控制优化算法，主要内容包括：

1. 通过OpenSim来估计人体新陈代谢速率$V$。
2. 通过贝叶斯优化算法来优化控制算法中的参数$p$，使之最小化人体形成代谢速率$V$。

一. 实验硬件部分

1. 带有位置与力传感的跑步机&动捕系统

如上图所示，在实验过程中记录受试者走路时脚掌与跑步机界面的力与力矩${\mathrm{f}}_g$：
$${\mathrm{f}}_g=(f_x,f_y,f_z,0,T_y,0)\in {\mathbb{R}}^6,$$
其中$y$为垂直于接触界面的方向，$T_y$代表力矩。实验同时记录了接触时压力中心的位置：
$${\mathrm{c}}_p=(c_x,c_y,c_z)\in {\mathbb{R}}^3.$$
同时通过动捕系统获取Markers的位置信息$x_o$
$$x_o=(p_x,p_y.p_z)\in {\mathbb{R}}^3.$$

2. Active Pelvis Orthosis (APO)外骨骼

外骨骼的参数如下所示：

二. 外骨骼控制算法(四种模式)

1. Transparent Mode

通过底层的PD控制器，使得外骨骼的助力${\tau }_{des}=0$，即模仿穿戴者没有穿戴外骨骼时的场景。

2. Adaptive Oscillator Mode

$${\tau }_{des}=K_v(\theta (\varphi -\delta \varphi )-\theta (\varphi )),$$
其中$K_v,\delta \theta$为可调整参数，分别代表虚拟刚度矩阵(Virtual Stiffness)以及相位差(文章用相位$\theta \in [0,100]$来表征当前行走处在的不同阶段，例如0代表左脚刚刚着地，右脚刚准备抬起；50代表右脚抬起到最高处)。$\theta (\cdot )$是将相位$\theta$映射到关节角度的函数。

这种控制算法相当于估计下一时刻关节角度位置，然后给予合适大小的力矩帮助穿戴者到达该关节角度位置。

3. Generic Mode

$${\tau }_{des}=\{\begin{array}{ll}-{\tau }_e,& 0\le x<17,\\ {\tau }_f,& 40\le x<64,\\ -{\tau }_e,& 92\le x\le 100,\\ 0,& \text{else}.\end{array}$$
将人体行走分为固定的四个阶段，每个阶段给予固定的助力(缺乏个性化)。

4. Human-in-the-Loop Mode（文章的卖点所在）

与3.Generic Mode中相同，文章将人体行走分为了五个阶段，并根据事先设定好的函数形式提供助力。不同的是，这种模式下五个阶段的分界线$p_i,i=1,2,3,4$变为了可调节参数。(注意只有$p_i$是可调节参数，${\tau }_e,{\tau }_f$仍为固定大小的值)。

文章在接下来的部分讲解了如何利用OpenSim的肌肉骨骼模型来优化$p_i$，以达到最小化人体新陈代谢速率的效果。

三. 通过OpenSim进行仿真，估计人体新陈代谢速率

文章在OpenSim中建立了Gait2354人体肌肉骨骼模型与IUVO建立的外骨骼CAD模型。外骨骼与人体肌肉骨骼模型之间通过固定点进行连接。模型中包括了许多Markers，与** 一.实验硬件部分**中的Markers点位置一一对应。

1.通过模型计算关节力矩

由模型中的Markers点与实验中的Markers点位置一一对应，可以通过最小化以下函数求解当前各关节角度$q$:
$$\underset{q}{\min }\sum _{i=1}^m{\omega }_i||x_{oi}-x_i(q)|{|}^2,$$
其中$m$为marker点的数量，$x_{oi}$为实验中通过动捕系统测得的第$i$个Marker的三维坐标，$x_i(\cdot )$为OpenSim的函数，其功能是计算前向运动学，即根据当前的关节角度$q$计算当前Marker点的位置。通过最小化以上函数，就可以求得当前各关节角度$q$。

在求得关节角度后，再根据逆动力学模型计算当前关节力矩$\tau$:
$$\tau =M(q)\ddot{q}+C(q,\dot{q})+g(q)+J(q{)}^T{f}_{\text{ext}},$$
其中$M,J$分别为模型的质量矩阵(mass matrix)以及Jacobian，$C,g$为科里奥利力项和重力项，$f_{\text{ext}}$为外界力。这些函数均为OpenSim内部函数。

2.通过关节力矩计算肌肉激活度

每一个关节力矩都由不同的肌肉-肌腱模块提供的力矩组成：
$${\tau }_j=\sum _{m=1}^c{\varphi }_m^T{C}_m,$$
其中$C_m$为力臂，${\varphi }_m$为肌肉力。而肌肉力通过以下的函数计算得到:
$${\varphi }_m=H({\alpha }_m,{\lambda }_m,v_m,{\pi }_m),$$
其中$H$为OpenSim内置黑箱函数，${\alpha }_m$为肌肉激活度，${\lambda }_m,v_m$为当前肌肉骨骼模型参数，${\pi }_m$为肌肉的参数。

通过联立以上两式，可以得到肌肉激活度${\alpha }_m$。注意到以上的问题一般都为超定方程，无法求得单一解，一般情况下都需要增加一目标函数以求得单一解。文章选择的目标为最小化肌肉激活度:
$$\min \sum _{m=1}^n({\alpha }_m{)}^2.$$

3.通过肌肉激活度计算肌肉能量(与新陈代谢速率相关)

通过黑箱模型，可以计算得到第m块肌肉运用的能量比例${\dot{ϵ}}_m$:
$${\dot{ϵ}}_m=U({\alpha }_m,{\lambda }_m,v_m,{\pi }_m),$$
肌肉总能量比例为
$$\dot{ϵ}=\sum _{m=1}^n{\dot{ϵ}}_m$$

四. 利用贝叶斯优化更新模型参数

在实验进行当中，实时记录markers的位置信息(用于计算关节角度)以及外部力信息(用于计算关节力矩):
P = { p 1 , p 2 , . . . , p n } , F = { f 1 , f 2 , . . . , f n } . P = \{p_1, p_2, ..., p_n\},\\ F = \{ f_1, f_2, ..., f_n\}. P={p1​,p2​,...,pn​},F={f1​,f2​,...,fn​}.
根据以上信息，可以依次得到关节角度信息$Q$：
Q = { q 1 , q 2 , . . . , q n } . Q = \{ q_1, q_2, ..., q_n\}. Q={q1​,q2​,...,qn​}.
肌肉激活度信息$A$:
A = { α 1 , α 2 , . . . , α n } . A = \{ \alpha_1, \alpha_2, ..., \alpha_n\}. A={α1​,α2​,...,αn​}.
以及肌肉能量比例信息$E$。
E = { ϵ 1 , ϵ 2 , . . . , ϵ n } . E = \{ \epsilon_1, \epsilon_2, ..., \epsilon_n\}. E={ϵ1​,ϵ2​,...,ϵn​}.
根据肌肉能量比例信息$E$，可以计算得到肌肉的能量消耗$V_i,V_s$:
$$\begin{gathered} V_i=\frac{1}{M(t_i-t_{i-1})}{\int }_{t_{i-1}}^{t_i}\dot{{ϵ}_i}(t)dt, \\ {V}_s=\frac{1}{n}\sum _{i=1}^n{V}_i. \end{gathered}$$
假设外骨骼控制算法中的参数$p_i,i=1,2,3,4$与$V_s$满足如下关系：
$$V_s=f_s^C(p_1,p_2,p_3,p_4).$$
则根据贝叶斯优化，可以求得$p_i^{\ast },i=1,2,3,4$使得V_s最小。