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

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

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

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

一. 实验硬件部分

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

如上图所示，在实验过程中记录受试者走路时脚掌与跑步机界面的力与力矩fg\boldsymbol{\mathrm{f}}_gfg：
fg=(fx,fy,fz,0,Ty,0)∈R6,\boldsymbol{\mathrm{f}}_g = (f_x,f_y,f_z,0,T_y,0) \in \mathbb{R}^6, fg=(fx,fy,fz,0,Ty,0)∈R6,
其中yyy为垂直于接触界面的方向，TyT_yTy代表力矩。实验同时记录了接触时压力中心的位置：
cp=(cx,cy,cz)∈R3.\boldsymbol{\mathrm{c}}_p = (c_x,c_y,c_z) \in \mathbb{R}^3. cp=(cx,cy,cz)∈R3.
同时通过动捕系统获取Markers的位置信息xox_{o}xo
xo=(px,py.pz)∈R3.x_o = (p_x,p_y.p_z) \in \mathbb{R}^3. xo=(px,py.pz)∈R3.

2. Active Pelvis Orthosis (APO)外骨骼

外骨骼的参数如下所示：

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

1. Transparent Mode

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

2. Adaptive Oscillator Mode

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

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

3. Generic Mode

τdes={−τe,0≤x<17,τf,40≤x<64,−τe,92≤x≤100,0,else.\tau_{des} = \begin{cases} -\tau_e, & 0\le x<17,\\ \tau_f, & 40 \le x <64, \\ -\tau_e, & 92 \le x \le 100, \\ 0, & \text{else}. \end{cases} τdes=⎩⎨⎧−τe,τf,−τe,0,0≤x<17,40≤x<64,92≤x≤100,else.
将人体行走分为固定的四个阶段，每个阶段给予固定的助力(缺乏个性化)。

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

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

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

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

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

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

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

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

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

每一个关节力矩都由不同的肌肉-肌腱模块提供的力矩组成：
τj=∑m=1cϕmTCm,\tau_j = \sum_{m=1} ^c{\phi_m^TC_m}, τj=m=1∑cϕmTCm,
其中CmC_mCm为力臂，ϕm\phi_mϕm为肌肉力。而肌肉力通过以下的函数计算得到:
ϕm=H(αm,λm,vm,πm),\phi_m = H(\alpha_m,\lambda_m, v_m,\pi_m), ϕm=H(αm,λm,vm,πm),
其中HHH为OpenSim内置黑箱函数，αm\alpha_mαm为肌肉激活度，λm,vm\lambda_m, v_mλm,vm为当前肌肉骨骼模型参数，πm\pi_mπm为肌肉的参数。

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

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

通过黑箱模型，可以计算得到第m块肌肉运用的能量比例ϵ˙m\dot\epsilon_mϵ˙m:
ϵ˙m=U(αm,λm,vm,πm),\dot\epsilon_m = U(\alpha_m,\lambda_m, v_m,\pi_m), ϵ˙m=U(αm,λm,vm,πm),
肌肉总能量比例为
ϵ˙=∑m=1nϵ˙m\dot\epsilon = \sum_{m=1}^n\dot\epsilon_m ϵ˙=m=1∑nϵ˙m

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

在实验进行当中，实时记录markers的位置信息(用于计算关节角度)以及外部力信息(用于计算关节力矩):
P={p1,p2,...,pn},F={f1,f2,...,fn}.P = \{p_1, p_2, ..., p_n\},\\ F = \{ f_1, f_2, ..., f_n\}. P={p1,p2,...,pn},F={f1,f2,...,fn}.
根据以上信息，可以依次得到关节角度信息QQQ：
Q={q1,q2,...,qn}.Q = \{ q_1, q_2, ..., q_n\}. Q={q1,q2,...,qn}.
肌肉激活度信息AAA:
A={α1,α2,...,αn}.A = \{ \alpha_1, \alpha_2, ..., \alpha_n\}. A={α1,α2,...,αn}.
以及肌肉能量比例信息EEE。
E={ϵ1,ϵ2,...,ϵn}.E = \{ \epsilon_1, \epsilon_2, ..., \epsilon_n\}. E={ϵ1,ϵ2,...,ϵn}.
根据肌肉能量比例信息EEE，可以计算得到肌肉的能量消耗Vi,VsV_i,V_sVi,Vs:
Vi=1M(ti−ti−1)∫ti−1tiϵi˙(t)dt,Vs=1n∑i=1nVi.V_i = \frac{1}{M(t_i-t_{i-1})}\int_{t_{i-1}}^{t_i}{\dot{\epsilon_i}(t)dt},\\ V_s = \frac{1}{n}\sum_{i=1}^{n}{V_i}. Vi=M(ti−ti−1)1∫ti−1tiϵi˙(t)dt,Vs=n1i=1∑nVi.
假设外骨骼控制算法中的参数pi,i=1,2,3,4p_i, i=1,2,3,4pi,i=1,2,3,4与VsV_sVs满足如下关系：
Vs=fsC(p1,p2,p3,p4).V_s = f_s^C(p_1,p_2,p_3,p_4). Vs=fsC(p1,p2,p3,p4).
则根据贝叶斯优化，可以求得pi∗,i=1,2,3,4p_i^*, i=1,2,3,4pi∗,i=1,2,3,4使得V_s最小。