1 引言

我国和世界上许多国家一样，正在步入老龄化.在老龄人群中有大量的脑血管疾病或神经系统疾病患者，这类患者多数伴有偏瘫症状^[1].近年由于患心脑血管疾病，特别是脑卒中使中老年患者出现偏瘫的人数不断增多，并且在年龄上呈现年轻化趋势^[2].世卫组织数据显示，2010年中国因脑卒中死亡人数列全球首位，发病率以每年8.7%的速度增长.每年因中风而死的中国人占所有死亡人数的22.45%^[3].一些医疗机构对患者治疗和恢复情况进行长期的观察和分析，证明了强制运动对上肢康复的效果确实有显著作用^[4].

将机器人技术与康复医学结合，用机器人代替理疗师来辅助偏瘫患者进行康复训练，已经成为国内外相关领域的研究热点^[5-6].近十年来，有很多康复机器人系统涌现出来，用来帮助医生和病人进行康复训练^[7-8].最著名的如美国麻省理工学院开发的MIT-MANUS上肢康复机器人^[9]，瑞士Hocoma公司的Locomat步态康复机器人等^[10].目前，机器人辅助康复训练的有效性已经得到临床验证^[11-12].但是现有的康复训练系统在进行主动训练时实现困难，所以患者大部分的训练都是被动的，而且训练方案相对单调；另外，在康复运动训练过程中，患者手臂的训练轨迹以及施加在患者手臂上的力很难保证一致性.因此在训练方式和控制策略上，上肢康复机器人研究需要有所创新和突破.

根据患者病情与身体状况不同，康复训练主要分成两个阶段，第一阶段称为被动运动阶段，该阶段患者上肢无运动能力，只能依靠机器人手臂带动患者上肢运动：第二阶段为自主运动阶段，此阶段患者上肢已有一定的运动能力，但是仍存在一定的运动障碍，需要适当的助力^[8].其中，被动训练通过常规的机器人轨迹跟踪控制即可实现，但是由于缺乏患者的主动参与，无法充分地刺激神经回路，训练效果相比主动康复训练要差^[13].脑卒中上肢偏瘫患者在主动康复训练阶段的患者主动力是康复机器人研究中十分重要的一个方面，由于主动训练过程中存在康复机器人与运动功能受损的患者手臂之间的相互作用，而患者是有自主运动意图的训练对象，因此对机器人和患者手臂之间的主动训练控制策略的研究是十分重要的.

本文的目标是设计一类主动训练控制策略，能为患者提供一个既舒适又安全的康复训练环境，防止患肢产生震颤、痉挛等不安全且非正常的肌肉反应，避免患者受到二次伤害.本文首先建立了带有患者交互力的上肢康复机器人主动训练系统模型.然后根据交互力设计一类非线性滚动时域控制器，以人工免疫优化算法为滚动优化策略，计算各个采样时刻的最优解，保证控制器的可行性和稳定性，最后对本文提出的控制策略进行了仿真验证.

本文研究对象为五自由度上肢外骨骼康复机器人，该上肢康复机器人的5个自由度能够模拟人体上肢5个主要关节的运动，分别是肩部外展/内收(关节1) 和俯/仰(关节2)、肘部俯/仰(关节3)、腕部旋内/旋外(关节4) 和腕部俯/仰(关节5). 5个关节处均设有力矩传感器，可实时采集患者与机械臂之间的交互力.康复机械臂的上臂、前臂长度可调，可实现左、右手穿戴，满足患者双侧康复的需求.另外，机器人还可实现单关节运动以及复杂的多关节联合运动，例如进行洗脸、提裤等符合人体运动规律的动作，满足患者日常生活的要求.五自由度上肢外骨骼机器人的原理示意图如图 1所示.

图 1 五自由度上肢外骨骼康复机器人机械结构示意图 Figure 1 Schematic diagram of the 5 DOF upper limb rehabilitation robot

图选项

脑卒中病程分为急性期(也称软瘫期或者弛缓期)、恢复期和后遗症期^[14].

脑卒中急性期持续时间一般为2~4周，在发病后的最初几天，根据患者的实际病情，一般以临床抢救为主，但同时提倡尽早介入康复治疗，早期康复治疗对患者的最终康复有很大的促进作用.病发一周后，患者进入亚急性期(也称痉挛期或硬瘫期)，血压、脉搏、呼吸等生命体征已基本稳定，意识清醒，多数能够理解医护人员的语言，并能配合康复治疗.因此应鼓励患者发挥自身的主观能动作用，积极进行康复功能训练.在此期间，肌肉痉挛加重，出现一些低级中枢控制下的原始运动形式.此时应通过正确的康复技术来降低肌张力、尽快使肌肉痉挛减轻.痉挛期内应鼓励患者进行自我被动运动，即患者利用健侧的力量活动患肢的运动.如果患者的能力允许，则还应鼓励患者进行肢体的主动运动，即依靠患肢自身力量进行的运动.

在主动康复训练过程中，患肢所产生的主动交互力是时变的，即使患者处于非痉挛期，发生肌肉痉挛也是不可避免的，且其不确定性可能会导致机器人系统的不稳定.本文通过各关节处的力矩传感器采集交互力信息，并将患者肌肉痉挛力矩作为一项扰动考虑到机器人动力学模型中^[15]，则上肢康复机器人动力学模型可描述为

(1)

其中，分别为关节角度，速度和加速度矢量，M(q)∈R^n×n为广义惯性矩阵，为哥氏力和离心力矩阵，G(q)∈Rⁿ为重力矩阵，B∈Rⁿ为康复过程中存在的患者未知有界肌肉痉挛力矩扰动，T_int为人机交互作用力矩，T_u∈Rⁿ为控制力矩矩阵.影响到康复训练系统的肌肉痉挛力矩模型定义为

(2)

其中，b_i=β_i(t-T₀)f(t)，i=1，2，3，4，5，表示由肌肉痉挛造成的机器人系统动态的偏移量. β_i(t-T₀)用来刻画此偏移力矩的随时间变化的特征，其中假设痉挛发生在未知时刻T₀，f(t)表示痉挛发生的概率和痉挛大小与训练时间成非线性关系，具体地说，我们把痉挛造成的力矩偏移随时间变化的特征表示为

(3)

其中，σ_i为大于0的正数，表示未知的变化率.较小值的σ_i刻画了慢变的肌肉收缩痉挛过程，也称初期痉挛：对于较大值的σ_i，时间特征β_i接近为阶跃函数，它表示了突发痉挛.本文中，假设偏移力矩向量和它们的时间导数都是有界的.

首先定义新的状态变量x₁=q_d-q为位置跟踪误差，其导数项，则式(1) 转化为

(4)

由此可得系统状态空间表达式

(5)

式中，δ(t，x)代表不确定性和患者痉挛力矩的非线性干扰.令控制量u为

通过以上定义，可将上肢康复机器人系统线性化，从而得到带有人机交互抵抗力矩和痉挛力矩不确定性的五自由度上肢康复机器人系统的跟踪误差动态方程如下：

(6)

在每个采样时刻t_k，将式(6) 固定线性化为

(7)

假设 给定的期望轨迹q_d，期望速度与期望加速度均有界.

滚动时域跟踪控制器设计如下：

(8)

其中，A_cl，k=A_k+B_2kK_k，C_cl，k=C_k+D_2kK_k，当且仅当存在一个对称矩阵P使得下面的不等式成立^[16]：

(9)

上式等价于P-A_cl^TPA_cl-γ^-1C_cl^TC_cl > 0，即A_cl^TPA_cl-P < 0.因此，A_cl是内部稳定的.

令Q^-1=P，进而得到式(7) 满足如下不等式

(10)

取状态反馈增益K_k=Y_kQ_k^-1.式(6) 的控制输入需要满足如下约束：

(11)

根据文[17]中L₂增益的定义，当式(6) 的初值x(t₀)=0时，对于大于0的常数γ，如果下面的不等式成立：

(12)

则称式(6) 所表示的系统有标准的L₂增益γ.

引理1 在每个采样时刻t_k，如果式(6) 满足如下不等式

(13)

可保证式(6) 满足L₂增益小于，其中.

引理2 定义椭圆域ε是一个半径为r_c的椭圆形区域，记为ε(Q_k^-1，r_c)={x^TQ_k^-1x≤r_c，x∈Rⁿ}.在每个采样时刻t_k，如果式(4) 满足如下不等式

(14)

且在采样间隔[t_k，t_k+1)内，系统的误差x(t)在椭圆域ε内，即x(t)∈ε(Q_k^-1，r_c)，则系统满足控制约束(11).其中.

引理3 假设(a)不确定性扰动有界，即 < ∞；(b)采样时刻t_k的误差x(t_k)在ε(Q_k^-1，r_c)={x∈Rⁿ|x^TQ_k^-1x≤r_c}内；那么如果采样时间非常短，系统在任意采样时刻[t_k，t_k+1)的误差也在相同的ε内，即x(t)∈ε(Q_k^-1，r_c).

证明 由引理3的条件(b)可得

(15)

由于扰动有界及γ_k < ∞得

(16)

进而

(17)

由式(10)、式(15) 和式(17) 可得，即在任意时刻t∈[t_k，t_k+1)的状态x(t)也在同一个椭圆域ε(Q_k^-1，r_c)内.

对于采样周期较长的情况，如果外部干扰满足：

(18)

则仍能由式(10) 和式(15) 得到x^TQ_k^-1x≤r_c.

引理4 给定γ > 0，假设式(6) 在平衡点处的L₂增益小于γ；A是渐近稳定的.那么存在平衡点的邻域使得非线性系统(1) 跟踪误差局部满足L₂增益小于γ.

定理 如果式(6) 的解唯一并且充分光滑，不确定性和患者痉挛力矩的扰动有界，式(6) 是零状态可观测的，采样间隔充分小或干扰满足式(18) 和满足引理1的条件；在每个采样时刻均有优化解，即(γ_k，Q_k，Y_k)满足如下不等式

(19)

(20)

那么对于式(1) 所表示的康复机器人系统：

(a)对于非持续性干扰是渐近稳定的；

(b)跟踪误差局部L₂增益小于，这里定义为

(21)

证明 假设在每个采样时刻均有优化解，则由式(11) 可得

(22)

由于V(x)=x^TQ^-1x，式(20) 可写成

(23)

由式(21) 和引理1的式(14) 可得式(6) 满足下式

(24)

假设干扰是有界的，那么对于有界的γ(t)，则有.当t→∞，可得

(25)

由式(25) 可知由零状态可观得到闭环系统是渐近稳定的.性质(a)得证.

由定义(23)，式(24) 变为

(26)

即式(6) 在平衡点(0，0) 处的L₂增益小于，应用引理4可知存在平衡点的领域使得康复机器人系统的跟踪误差局部L₂增益小于，性质(b)得证.

对于非线性康复机器人系统滚动优化控制而言，优化问题的约束在每个采样时刻是变化的.而且，患者的干扰信息和康复情况是不可预知的.因此，出现突发干扰时，就不能保证优化问题的可行性.本文采用如下最大—最小(minimax)优化问题替代最小化问题^[18]，以保证控制策略的可行性，达到实时对抗患者扰动的目的.

其中，U是所有允许控制函数的集合，W是所有允许干扰函数的集合.本节应用拉格朗日松弛法^[18]将约束minimax问题转化为最小化问题.

假设允许干扰函数的能量为α，即允许干扰函数集合为，则定义如下拉格朗日函数：

其中，γ²为拉格朗日算子.那么，对任意的允许干扰，可以得出：

(27)

对于给定的γ，定义储能函数

(28)

进而，对于任意的有界干扰，式(6) 满足如下不等式：

由式(27) 和式(28) 可得

(29)

由式(29) 可得

因此minimax问题转化为如下最小化问题：

上式等价于

(30)

其中，ω₁和ω₂为加权系数.由上述推导可知，r是作为的上界引入的，所以r和γ相同，都是反映系统性能的参数，所以二者之间不存在冲突.因此，ω₁和ω₂的选择没有固定的要求.

文[19]讨论了不考虑外部干扰的情况，即取ω₁=1，ω₂=0的情况.文[20]讨论了有外部干扰的情况，即取ω₁=0，ω₂=1的情况.本节所设计的滚动优化策略的基本思想是保留ω₁和ω₂为可调节参数，将r作为一个独立变量，引入一个正数η≥0，令r < r_c(1+η)，在每个采样时刻t_k，通过求解优化问题(30)，得到(γ_k，Q_k，Y_k)，如果无解，则需要通过人工免疫优化算法来计算η来优化椭圆域ε，直到其包含x(t_k)，从而使优化问题有解.计算K_k=Y_kQ_k^-1，最终得到控制量u=K_kx.改进的人工免疫算法流程如图 2所示.

图 2 改进免疫算法流程图 Figure 2 Flow chart of improved immune algorithm

图选项

本节对参数η进行二进制编码.首先设定参数变量的编码长度为l，得到二进制字符串，构成免疫遗传算法的个体.

定义η_min和η_max分别为参数η的下界和上界，编码精度为.参数η所对应的二进制编码为

(31)

其中，Round(·)表示四舍五入对η取整，对应字符串长度若小于l，则在首位补0.

聚集密度表示相似抗体的稀松程度，既能表示出每个抗体之间的联系，也能表征出种群的分散性和多样性，在种群进化过程中起到调节种群的作用.抗体相似性可由抗体编码的Euclidean距离表示，在一个规模为N的种群中，第i个抗体X_i和第j个抗体X_j之间的Euclidean距离为

(32)

其中，X_i^(k)和X_j^(k)分别表示抗体抗体X_i和抗体X_j的第k位上的基因.进而，得到抗体X_i的聚集密度为

(33)

在传统的进化算法(例如GA)中，每一次进化之后，父代个体会立刻被子代个体取代.尽管这样可以保证种群的质量，但是种群的多样性和分散性可能会恶化.本节提出一种基于小生境选择的种群更新策略，通过引入一个临时种群Q′来存储交叉和变异操作产生的新个体，然后在集合Q∪Q′中选择质量高的个体，作为下一次进化的父代种群，记Q为种群，η为种群中的个体.具体步骤如下：

Step 1 由式(32) 计算种群Q的每个个体距离；

Step 2 随机选取两个父代个体，进行交叉和变异操作，生成一个子代个体，将该子代个体存入新种群Q′中；

Step 3 将Q∪Q′中的解按照其目标函数从小到大的顺序排序，即按照解的质量从高到低的顺序排序，得到一个解的序列Q_2N=(η₁，η₂，η₃，…，η_2N)；设置m=0，被删除解集合C为空集；

Step 4 定义小生境半径α，根据种群Q∪Q′中解的平均距离D₀来自适应计算α；

Step 5 如果解的序列Q_2N为空集，则转Step 9.否则，将Q_2N中的第一个解η₁加入到新种群Q_new中.这样可以保证当前的历史最优解一定会进入到下一代的新种群中.再将η₁从Q_2N中删除，更新Q∪Q′中解序列的编号，更新m=m+1；

Step 6 将Q_2N中与η₁距离小于或等于α的解全部删除，再将这些被删除的解插入到被删除解的集合C中，然后更新Q_2N中解序列的编号；

Step 7 如果m < N，则转至Step 5.否则，已经得到新种群；

Step 8 如果当前的新种群中解的数量小于N，从被删除解的集合C中随机选择解加入到新种群中，直到新种群中解的数量达到N；

Step 9：种群更新结束.

以上过程可以看出，通过计算抗体距离，并采用小生境方法来选择新种群中的解，不仅可以保证当前的历史最好解能够进入新种群，而且在每次插入新的解后，会将其周围抗体浓度较大的解删除，从而保证了种群分散性，提高了种群的质量.

假设康复机器人连杆质量均匀分布，系统参数的具体数值如表 1所示.

表 1 上肢康复机器人杆件参数 Table 1 Parameters of upper limb rehabilitation robot

表选项

系统的初值取为x(0)=[0.7，0.2，-0.5，0.4，0.3]^T，.将偏移力矩函数设计为：b₁=β₁(t-T₀)f₁(t)=(1-e^-(t-T₀)/10)，b₂=b₃=b₄=b₅=0.人工免疫算法参数设计为：N=100，交叉概率p_c=0.5，变异概率p_m=0.3.参数α=0.25D₀，ω₁=0.4，ω₂=0.6.设计最大迭代次数为50.

图 3和图 4为5个关节的角度和速度响应曲线，其中x₁~x₅为各关节的角度跟踪误差，单位为rad. 为各关节的速度跟踪误差，单位为rad/s.可以看出系统存在有界干扰时，系统的震荡幅度较小，角度误差和速度误差均可以保证在可接受范围内，这表明未知扰动和不确定性痉挛并不会影响系统的稳定性.系统的控制输入响应曲线如图 5所示，尽管控制量u有小范围跳跃，也可以保证系统的连续性.椭圆域优化结果如图 6所示，可以清晰地看到椭圆域半径r_c的优化过程，保证了状态x(t_k)始终在椭圆域内，其对应的目标函数值如图 7所示，说明采用人工免疫优化算法作为滚动优化策略是可行且有效的.性能指标γ的优化结果如图 8所示.

图 3 康复机器人主动训练系统角度误差响应曲线 Figure 3 Response curve of the angle error in rehabilitation robot active training system

图选项

图 4 康复机器人主动训练系统角速度误差响应曲线 Figure 4 Response curve of the velocity error in rehabilitation robot active training system

图选项

图 5 康复机器人主动训练系统控制输入 Figure 5 Control input of rehabilitation robot active training system

图选项

图 6 康复机器人主动训练系统椭圆域优化过程 Figure 6 Elliptic domain optimization process of rehabilitation robot active training system

图选项

图 7 滚动时域控制器期望最优值 Figure 7 The expected cost of receding horizon controller

图选项

图 8 滚动优化性能指标γ Figure 8 Performance index γ of receding horizon optimization

图选项

本文针对上肢偏瘫患者在主动康复训练阶段存在未知扰动和痉挛力矩不确定性问题，设计了一类基于人工免疫优化算法的滚动时域控制器，以提高人机交互系统的可靠性.通过构建椭圆域，采用人工免疫算法优化其半径，计算各个采样时刻的可行解，可以保证该滚动时域控制器的可行性和稳定性，即可使不确定性、未知扰动得到适当的处理而不会影响主动训练的目的.实验结果表明，在五自由度上肢康复机器人系统中，所提出的滚动时域控制策略是可行的.康复训练的安全性和更精确的轨迹跟踪控制将作为未来深入研究的重点.