Customization and Optimization of the Upper Limb Rehabilitation Trajectory and Compliance Tracking Control
-
摘要: 卒中后早期的被动康复训练可以促进患者脑神经重塑甚至重获肢体运动能力。多数上肢康复机器人中，只将参数化规则曲线作为被动训练的轨迹，没能融合康复医师的经验及患者的个体特征，缺乏3维个性轨迹的光滑优化；此外，训练过程中的安全柔顺交互的问题未得到解决。针对以上问题，本文基于新型末端牵引式3自由度上肢康复机器人开展轨迹个性化定制优化以及跟踪控制研究。首先，在笛卡尔空间应用导纳算法及力补偿策略实现理疗师对轨迹的定制。然后，用道格拉斯–普克法压缩原始轨迹数据得到型值点，保留初始轨迹的拓扑形状。接下来，用非均匀有理B样条曲线（non-uniform rational b–spline，NURBS）进行插值，并融合动态切换概率和t变异改进蝴蝶优化算法（butterfly optimization algorithm，BOA）优化拟合的轨迹。最后，设计基于径向基（radial basis function，RBF）神经网络的滑模自适应以及速度前馈加比例–积分–微分（proportional–integral–derivative，PID）的控制策略实现末端轨迹跟踪，并增加导纳–阻抗控制形成基于力阈值的多模式柔顺跟踪控制，保证人机交互力较大时患者的安全。结果表明：定制机器人末端轨迹的拖动力在5 N以内；改进的BOA算法具有更高的收敛速度和精度，优化的轨迹曲率和更小；轨迹跟踪控制策略可以将跟踪误差控制在6 mm内；模式切换在0.3 s内响应，能顺应较大交互力，从而提高训练的安全性。Abstract: Early passive rehabilitation training after stroke can promote brain nerve remodeling and even regain motor ability. Most upper limb rehabilitation robots use parameterized regular curves as passive training trajectories, which fail to integrate the experience of rehabilitation doctors and the individual characteristics of patients, or lack smooth optimization of three-dimensional customized trajectories and safe interaction during training. Therefore, this research studied the trajectory customization, optimization and tracking control strategy based on the novel end-effector 3 degrees of freedom (DOFs) upper limb rehabilitation robot. Firstly, the admittance algorithm and force compensation strategy were applied in Cartesian space to realize the customization of the trajectory by the physiotherapist. Then, the original trajectory data were compressed by Douglas-Peucher method to obtain a few via points. The non-uniform rational B-spline curve (NURBS) was used for interpolation, and an improved butterfly optimization algorithm (BOA) based on dynamic switching probability and t mutation was used to minimize the trajectory curvature. Finally, the sliding mode adaptive control based on the radial basis function (RBF) network and the speed feedforward plus PID control strategy were designed to control the trajectory tracking. In addition, the admittance-impedance control was added, forming a multi-mode tracking controller, to make the training safe with the large human-robot interaction force. The results showed that the dragging force of trajectory customization is within 5N. The improved BOA algorithm has higher convergence speed and accuracy, and the curvature of the optimized trajectory is smaller. The control strategy can resist the interference within the threshold and keep the tracking error within 6mm. The mode can switch in 0.3s, to compliant to large interaction force and improve the safety of training .
-
随着全球老龄化及神经损伤造成的肢体运动功能障碍患者人数的逐步上升，如何满足庞大的康复需求、提高患者生活质量是亟需解决的关键问题。临床结果表明，在卒中后早期，被动康复训练能够促进神经重塑，帮助患肢肌肉重获自主收缩的能力[1]。相较于传统的手动康复方式，机器人辅助康复技术具有的重复精度高、定量化评估准确等优点，已经促使其成为世界范围内的研究热点[2]。传统康复理疗师拖拽病人的患肢进行康复训练的方式结合了理疗师的实践经验并可以实时交互，可以给患者更好的训练效果。而目前大多数康复机器人的训练轨迹是一些提前定义好的规则曲线，在训练范围上有一定的限制。另外，由于不能适应较大的突发人机交互力，在跟踪轨迹康复训练过程中存在危险或者训练直接中断过程的情况。日常上肢的活动通常是在参数多变的非结构环境中，所以有必要研究能自由定制被动康复训练轨迹且可以柔顺交互的上肢康复机器人。
目前，Morita等[3]使用基于位置的阻抗控制策略实现机器人末端轨迹定制，使用最小二乘法拟合原始的噪声数据来近似康复医师的意图，并用拟合后的轨迹作为被动训练的期望轨迹。游有鹏等[4]基于自行研制的直流电动机设计力矩控制及机器人连杆重量和关节摩擦力自测量补偿算法，实现无力传感器的轻松拖动。杨浩等[5]先在笛卡尔空间应用导纳算法，求得笛卡尔空间期望位置后再用逆解求关节位置来拖动下肢康复机器人，并未对个性化轨迹进行优化和跟踪控制。但是，康复轨迹的构建过程中，医生和病人的手同时与机器人的末端耦合，会因为病人、医生抖动或者环境噪声等原因，构建出不平滑的轨迹。Zeestraten等[6]基于统计模型进行轨迹学习和规划，但是计算效率较低，且生成的轨迹不连续。对轨迹进行光滑优化可以采用具有规范性和连续可微性等重要性质的B样条曲线[7]。Wan等[8]利用改进B样条曲线规划机器人的避障轨迹。董甲甲等[9]没有直接对笛卡尔空间的3维轨迹的优化，而是改进B样条的求解过程，用于6R机器人的关节轨迹拟合优化。天津大学的Mei等[10]将关节最小急动度和B样条曲线相结合，对6自由度高速并联机器人的末端轨迹进行了优化，重点是降低关节的急动度，拟合轨迹相对于原始轨迹有较大的变形。Ijspeert等[11]利用动态运动基元的方法来拟合期望轨迹，当目标点改变后能够快速的拟合出新的轨迹，泛化能力强。但是轨迹逼近的算法复杂，理解困难且需要确定的参数较多。另外，以上研究并没有对插值拟合的轨迹进一步光滑优化。
轨迹跟踪控制的方法研究也较多。吴青聪等[12]结合外骨骼上肢康复机器人设计模糊滑模控制器实现被动训练的位置跟踪，理疗师可以根据患者状况手动调整外环导纳参数，使得人机交互力参与到机器人末端位置的控制中。Mushage等[13]设计模糊神经网络和基于状态观测的误差自适应非线性控制器对5自由度上肢外骨骼的轨迹跟踪进行控制。罗定吉等[14]采用模糊比例–积分–微分（proportional–integral–derivative，PID）控制算法对下肢外骨骼进行轨迹跟踪控制。Shi等[15]将径向基（radial basis function，RBF）神经网络和PID结合起来进行下肢外骨骼机器人运动的轨迹跟踪控制。现有多数研究都以提高轨迹跟踪精度为目标，而康复训练中应以保证安全为前提，再提高轨迹跟踪精度。Zhang等[16]设计变刚度执行器和力矩限制机构用于上肢康复外骨骼来保证训练中的安全性。Wu等[17]设计基于RBF神经网络的干扰观测器来补偿外界干扰和误差来提高轨迹跟踪精度，同时用自适应变导纳控制策略使机器人具有柔顺性来适应不同状态的患者。Trigili等[18]利用串联弹性关节元件来实现康复机器人的柔顺性。Miao等[19]在双边末端牵引式上肢康复机器人上设计位置控制用于被动训练，并基于个性化工作空间设计自适应变参数控制实现柔顺性。Guo等[20]利用强化学习算法设计变导纳控制算法实现匹配患者下肢刚度特性的康复训练，也具有一定的柔顺性。李超等[21]基于患者病情设计模糊自适应的粒子群位置阻抗控制算法用于下肢的康复训练；粒子群算法优化导纳参数提高位置跟踪精度，模糊自适应算法实现一定的柔顺性。以上研究中，特殊柔性元件结构及建模较复杂而且动力学特性会存在误差。其他的柔顺控制算法虽然可以在线修改参数，但是算法略显复杂。卒中病人在早期运动能力是比较弱的甚至没有，经过一段时间的训练恢复后会进入主动训练阶段才需要参数的在线实时修改。
针对以上问题，作者基于自行研制的3自由度末端牵引式上肢康复机器人，在3维空间内个性化定制训练轨迹，并对其进行光滑优化。设计多模式控制器使机器人末端既能跟踪个性化训练轨迹，又能顺应较大干扰，且参数设置简单。当干扰降低后能够继续跟踪期望轨迹，保证训练的安全性和连续性，匹配应用需求。
1. 平台及方法
1.1 上肢康复机器人系统
本研究中使用的上肢康复机器人是自行研制的末端牵引式空间3自由度上肢康复机器人，主要是针对上肢肩关节和肘关节的康复训练，包括两个水平转动关节和1个垂直移动关节，结构简图如图1所示。其中Z1、Z2、Z3代表3个关节的轴线，q1、q2、q3为3个关节的位置变量。两个转动关节由两个配置绝对编码器的交流伺服电机驱动，垂直移动关节由双作用气缸驱动。在气缸末端安装拉线式位移传感器，在气缸和小臂的连接处安装3维力传感器。康复机器人的试验平台如图2所示。本系统的控制系统包括上位机（基于MATLAB软件进行系统程序的设计和编译）和下位机（执行编译后的控制程序）。将电机的驱动方式设置为力矩模式，通过控制输入电压设定电机输出力矩，通过脉冲闭环反馈控制其转动位置。而气缸通过比例调压阀控制气缸输出力，利用末端位移反馈控制其位置。因为机器人整体高度的限制，气缸的总行程为150 mm，可实现上肢关节在矢状面内一定范围的训练。
图 1 机器人结构简图Fig. 1 Structure diagram of the robot
下载:
全尺寸图片
图 2 康复机器人试验平台Fig. 2 Platform of the rehabilitation robot system下载:
全尺寸图片
1.2 轨迹定制
理疗师慢速拖动机器人时，忽略机器人的惯性力，需要克服患者加在末端的力
${{\boldsymbol{F}}_{\text{p}}}$ 以及机器人的关节静摩擦力${{\boldsymbol{F}}_{\text{f}}}$ ，如式（1）所示：$$ {{\boldsymbol{F}}_{\rm{t}}} \geqslant {{\boldsymbol{F}}_{\text{p}}} + {{\boldsymbol{F}}_{\text{f}}} $$ (1) 式中，
${{\boldsymbol{F}}_{\text{t}}}$ 代表理疗师施加在末端的力。对于几乎完全丧失上肢运动能力的软瘫期病人，假设${{\boldsymbol{F}}_{\rm{p}}} = 0$ 。检测并补偿转动关节的静摩擦力后，利用机器人静力学原理将交互力
${{\boldsymbol{F}}_{{{\rm{int}}} }}$ 转换为关节空间的力矩${\tau _{\rm{h}}}$ ，如式（2）所示：$$ {{\boldsymbol{\tau }}_{\rm{h}}} = {{\boldsymbol{J}}^{\rm{T}}}{{\boldsymbol{F}}_{{\text{int}}}} $$ (2) 式中，J为机器人的雅克比矩阵。然后，在关节空间中应用导纳算法产生期望关节位置输入，如式（3）所示：
$$ \Delta{\boldsymbol{ q}} = {\boldsymbol{k}}{{{\tau }}_{\rm{h}}} $$ (3) 式中，k表示导纳系数矩阵，
$\Delta {{\boldsymbol{q}}}$ 表示关节位置增量向量。由位置增量和实际的关节位置变量q可以得到所期望的关节位置变量qd，将其输入到电机驱动系统后由内部算法得到关节的驱动力矩${\mathbf{\tau }}$ ，从而实现转动关节拖动。对于气缸驱动的移动关节，末端初始的垂直负载和理疗师施加的垂直力（大小为
${F_{{\text{pz}}}}$ ）可由末端的力传感器检测。根据理疗师施加的力大小${F_{{\text{pz}}}}$ ，由式（4）计算得到气缸无杆腔的压力改变量$\Delta P$ 。气缸位移随无杆腔气体压力的改变而改变，实现移动关节的轻松拖动。由3个关节的位置传感器记录的位置点得到初始的定制轨迹。$$ {{F}_{{\text{pz}}}} = \Delta {P} \times {{A}_1} $$ (4) 式中，
${{\text{A}}_{\text{1}}}$ 为气缸无杆腔的有效面积。1.3 NURBS插值曲线优化
对初始轨迹点应用数据压缩算法[22]得到有限个型值点，其坐标用Pi表示， i=1,2,···,n。选择可以如实反映型值点按弦长分布情况且拟合精度高的积累弦长参数化法[23]对型值点进行参数化处理得到节点矢量
${\boldsymbol{U}}$ ，${\boldsymbol{U}} = [{u_0},{u_1}, \cdots ,{u_{n + k + 1}}]$ 。积累弦长参数化法表示如式（5）所示：
$$ \left\{ \begin{gathered} {u_0} = {u_1} = {u_2} = {u_3} = 0, \\ {u_{i + 3}} = {u_{i + 2}} + \left| {{{\boldsymbol{P}}_i} - {{\boldsymbol{P}}_{i - 1}}} \right|/\sum\limits_{i = 1}^n {\left| {{{\boldsymbol{P}}_i} - {{\boldsymbol{P}}_{i - 1}}} \right|}, \\ {u_{n + 3}} = {u_{n + 4}} = {u_{n + 5}} = {u_{n + 6}} = 1 \\ \end{gathered} \right. $$ (5) 再按 De Boor算法[24]计算
$\bar{k} $ 次规范B样条基函数${N_{i,k}}(u)$ ，如式（6）所示：$$ \left\{ \begin{array}{l}{N}_{i,0}(u)=\left\{ \begin{array}{l}1\text{，}\qquad {u}_{i}\le u \le {u}_{i+1};\\ 0\text{，}\qquad 其他\end{array}\right.\\ {N}_{i,\bar{k}}(u)=\dfrac{u-{u}_{i}}{{u}_{i+\bar{k}}-{u}_{i}}{N}_{i,{k}-1}(u)+\dfrac{{u}_{i+\bar{k}+1}-u}{{u}_{i+\bar{k}+1}-{u}_{i+1}}{N}_{i+1,\bar{k}-1}(u);\end{array} \right.$$ (6) 式中，定义
$\dfrac{0}{0}=0 $ 。将式（6）得到的基函数带入$\bar{k} $ 次NURBS曲线公式中得到插值后的曲线$p(u)$ ，如式（7）所示：$$ p(u) = \frac{{\displaystyle\sum\limits_{i = 0}^n {{\omega _i}} {d_i}{N_{i,\bar{k}}}(u)}}{{\displaystyle\sum\limits_{i = 0}^n {{\omega _i}} {N_{i,\bar{k}}}(u)}} $$ (7) 式中，
$ {w_i} $ 为权因子，$ {d_i} $ 为控制顶点的值。数据压缩过程设定的阈值不同，得到的关键型值点
${{\boldsymbol{P}}_i}$ 也不同，曲线的形状和光顺性也不同。采用智能优化算法，以曲线光顺性为目标，以压缩阈值为变量优化插值曲线，生成连续光顺的插值曲线作为训练轨迹$p(u)$ 。目标函数为：
$$ Obj = \mathop \sum \limits_{j = 0}^l {\kappa _{\rm{c}}}(j) $$ (8) 式中，l为插值曲线离散点的总个数，
${\kappa _{{\rm{c}}}}(j)$ 为曲线在点j处的曲率，其公式为：$$ {\kappa _{\text{c}}} = \frac{{||{p^\prime }(u) \times {p^{\prime \prime }}(u)||}}{{||{p^\prime }(u)|{|^3}}} $$ (9) 插值曲线点要限制在机器人末端的安全工作空间内。根据机器人的逆运动学模型，由笛卡尔空间的坐标点可以得到3个关节的位置变量q。则优化约束条件表示为：
$$ {{\boldsymbol{q}}_{{\rm{min}}}} \leqslant {\boldsymbol{q}} \leqslant {{\boldsymbol{q}}_{{\rm{max}}}} $$ (10) 式中，
${{\boldsymbol{q}}_{{\text{min}}}}$ 和${{\boldsymbol{q}}_{{\text{max}}}}$ 分别表示机器人3个关节的最小限位和最大限位。气缸的行程比较小，在拖动的时候很容易就会到达的极限位置且不会产生任何危险，所以将第3关节超出极限位置的数据点修改为极限位置点。所用的蝴蝶优化算法[25]为一种新型的元启发式群智能优化算法，该算法在工程实践的初步应用中取得了良好的效果。
该算法的搜索行为依据为香味大小f，如式（11）所示。
$$ f = c{I^\alpha } $$ (11) 式中：
$I$ 为刺激强度，与蝴蝶的适应度有关；$c$ 为感官因子，通常取值为0.01；$\alpha $ 为依赖于$f$ 的幂指数，通常取值为0.1。标准蝴蝶算法的位置更新公式如下：
$$ {x}_{s}^{t+1}=\left\{ \begin{array}{l}{x}_{s}^{t}+{f}_{s}({r}_{1}^{2}{g}^{*}-{x}_{s}^{t})\text{，}\text{ }{r}_{3}\le {{ P}}_{{\rm{static}}};\\ {x}_{s}^{t}+{f}_{s}({r}_{2}^{2}{x}_{j}^{t}-{x}_{k}^{t}),\text{ }{r}_{3} > {P}_{{\rm{static}}}\end{array}\right. $$ (12) 式中：
$x_s^t$ 为第s只蝴蝶在第t次迭代中的位置坐标；$x_j^t$ 、$x_k^t$ 为从解空间中随机选择的第k只和第j只蝴蝶的位置坐标；${r_1}$ ，${r_2}$ 和${r_3}$ 为[0,1]内的随机数；${f_s}$ 是由式（11）得到的第s只蝴蝶的香味大小；g*表示第t次迭代后的全局最优解；${P_{{\rm{static}}}}$ 表示静态切换概率。算法处于局部搜索，还是全局搜索阶段依据随机数${r_3}$ 是否小于Pstatic的条件（Pstatic=0.8）。有许多对BOA算法的改进方法[26-27]，其设计的主要思想是：迭代前期种群个体应该侧重全局寻优，在解空间内进行大范围搜索；后期侧重局部寻优，使算法快速收敛。本文中融合了矢量数据压缩算法和BOA算法，将高维的插值点优化问题转化为1维的压缩阈值优化问题，简化了优化算法的求解过程。所以，仅对标准BOA算法做部分改进而不融合其他优化算法，其中，动态切换概率
${P_{\rm{d}}}$ 表达式为：$$ {P_{\rm{d}}} = {\left(\frac{{{T_{{\rm{max}}}} - n_1}}{{{T_{{\rm{max}}}}}}\right)^3} $$ (13) 式中，
${T_{{\rm{max}}}}$ 表示最大迭代次数，n1表示当前的迭代次数。采用t变异是为了防止陷入局部最优，提高优化收敛速度。具体是指在每一代的迭代过程中用t分布的概率密度函数代替随机数
${r_i}^2，i = 1,2$ 。将标准t分布的概率密度函数与优化算法的迭代次数进行关联，概率密度函数f(t)表示为：
$$ f(t) = 2.5\dfrac{{\varGamma \left(\dfrac{{n_2 + 1}}{2}\right)}}{{\sqrt {n{\text{π}} } \varGamma \left(\dfrac{n}{2}\right)}}{\left(1 + \dfrac{{{{(2t/{T_{{\rm{max}}}})}^2}}}{n}\right)^{ - \frac{{n_2 + 1}}{2}}} $$ (14) 式中，n2表示t分布的自由度；
$\varGamma ()$ 是Gamma函数，其定义表示为：$$ \varGamma (x) = \int \nolimits_0^\infty {{\rm{e}}^{ - t}}{t^{x - 1}}{\rm{d}}t $$ (15) 改进后的BOA算法中个体的位置更新公式为：
$$ {x}_{s}^{t+1}=\left\{ \begin{array}{l}{x}_{s}^{t}+{f}_{s}(f(t){g}^{*}-{x}_{s}^{t})\text{，}\text{ }{r}_{3}\le {P}_{{\rm{d}}};\\ {x}_{s}^{t}+{f}_{s}(f(t){x}_{j}^{t}-{x}_{k}^{t}),\text{ }{r}_{3} > {P}_{{\rm{d}}}\end{array}\right. $$ (16) 1.4 轨迹跟踪控制算法
根据Brunstorm的卒中后康复6阶段理论[28]，在卒中后早期的软瘫期，病人需要接受被动康复训练。其作用是通过患者的活动引起患者大脑的某些区域的兴奋，从而对当前的患肢活动产生意识。为了实现康复医师规划的训练动作，被动康复训练在保证安全性的前提下以良好的伺服跟踪性能为主要目标。得到优化的插值曲线之后，将曲线点作为被动康复训练的期望位置，设计轨迹跟踪控制器来保证机器人良好的轨迹跟踪能力。
1.4.1 跟踪控制器设计
m个关节机器人的动力学方程为：
$$ {\boldsymbol{M}}({\boldsymbol{q}}){\boldsymbol{\ddot q}} + {\boldsymbol{C}}({\boldsymbol{q}},{\boldsymbol{\dot q}}){\boldsymbol{\dot q}} + {\boldsymbol{G}}({\boldsymbol{q}}) + {\boldsymbol{F}}({\boldsymbol{\dot q}}) + {{\boldsymbol{\tau }}_{{\rm{d}}}} = {\boldsymbol{\tau }} $$ (17) 式中，
${\boldsymbol{q}}$ 、${\boldsymbol{\dot q}}$ 、${\boldsymbol{\ddot q}}$ 分别表示机器人关节的实际位置、速度和加速度；${\boldsymbol{M}}({\boldsymbol{q}})$ 为m×m阶正定惯性矩阵；${\boldsymbol{C}}({\boldsymbol{q}},{\boldsymbol{\dot q}})$ 为m×m阶科氏矩阵；${\boldsymbol{G}}({\boldsymbol{q}})$ 为m×1阶重力项向量；${\boldsymbol{F}}({\boldsymbol{\dot q}})$ 为摩擦力矩阵；${{\boldsymbol{\tau }}_{{\rm{d}}}}$ 为未知外加干扰，且满足$\left\| {{{\boldsymbol{\tau }}_{{\rm{d}}}}} \right\| \leqslant {\boldsymbol{D}}$ ，其中D表示未知干扰的上界；${\boldsymbol{\tau }}$ 为控制输入。用${{\boldsymbol{q}}_{{\rm{d}}}}$ 、${{\boldsymbol{\dot q}}_{{\rm{d}}}}$ 、${{\boldsymbol{\ddot q}}_{{\rm{d}}}}$ 分别表示期望的关节位置、速度和加速度。定义关节位置跟踪误差e为：$$ {\boldsymbol{e}} = {{\boldsymbol{q}}_{{\rm{d}}}} - {\boldsymbol{q}} $$ (18) 定义滑模面为r，表示为：
$$ {\boldsymbol{r}} = {\boldsymbol{\dot e}} + {\boldsymbol{\varLambda e}} $$ (19) 式中，
$\;{\boldsymbol{\varLambda }} $ 为对角矩阵，$\;{\boldsymbol{\varLambda }} = {{\boldsymbol{\varLambda }}^{\text{T}}} > {\text{0}}$ ，则$$ {\boldsymbol{\dot q}} = - {\boldsymbol{r}} + {{\boldsymbol{\dot q}}_{{\rm{d}}}} + {\boldsymbol{\varLambda e}} $$ (20) 由式（17）到（20）计算得到：
$$ {\boldsymbol{M}}({\boldsymbol{q}}){\boldsymbol{\dot r}} = - {\boldsymbol{C}}({\boldsymbol{q}},{\boldsymbol{\dot q}}){\boldsymbol{r}} - {\boldsymbol{\tau }} + {\boldsymbol{N}}({\boldsymbol{x}}) + {{\boldsymbol{\tau }}_{{\rm{d}}}} $$ (21) 式中，
${\boldsymbol{N}}({\boldsymbol{x}})$ 为需要逼近的模型不确定项。采用RBF神经网络逼近未知的N(x)，网络算法为：
$$ {\boldsymbol{N}}({\boldsymbol{x}}) \triangleq {{\boldsymbol{W}}^ * }^{\text{T}}{\text{φ}}({\boldsymbol{x}}) + {\boldsymbol{\varepsilon }} $$ (22) 式中，x为网络的输入，
${\text{φ}}({\boldsymbol{x}})$ 为网络的高斯基函数输出，W*为网络的理想权值且有界，${\boldsymbol{\varepsilon }}$ 为网络的逼近误差，$\left\| {\boldsymbol{\varepsilon }} \right\| \leqslant {{\boldsymbol{\varepsilon }}_{\text{N}}}$ ，${\varepsilon _{\rm{N}}}$ 是逼近误差的上界。取神经网络的输入x为：
$$ {\boldsymbol{x}} = \left[ {\begin{array}{*{20}{l}} {{{\boldsymbol{e}}^{\rm{T}}}}&{{{{\boldsymbol{\dot e}}}^{{\rm{T}}}}}&{{\boldsymbol{q}}_{{\rm{d}}}^{{\rm{T}}}}&{{\boldsymbol{\dot q}}_{{\rm{d}}}^{{\rm{T}}}}&{{\boldsymbol{\ddot q}}_{{\rm{d}}}^{{\rm{T}}}} \end{array}} \right] $$ (23) 设计控制输入
${\mathbf{\boldsymbol{\tau }}}$ 为：$$ {\boldsymbol{\tau }} = {\boldsymbol{\hat N}}({\boldsymbol{x}}) + {{\boldsymbol{K}}_{\rm{v}}}{\boldsymbol{r}} - {\boldsymbol{v}} $$ (24) 式中：
${\boldsymbol{\hat N}}({\boldsymbol{x}})$ 为RBF网络对动力学模型中不确定项${\boldsymbol{N}}({\boldsymbol{x}})$ 的估计值，${\boldsymbol{\hat N}}({\boldsymbol{x}}) = {{\boldsymbol{\hat W}}^{\rm{T}}}{\text{φ}}({\boldsymbol{x}})$ ，其中${\boldsymbol{\hat W}}$ 为网络连接权值的估计值；${{\boldsymbol{K}}_{\rm{v}}} $ 为对角正定矩阵，${{\boldsymbol{K}}_{\rm{v}}} \in {\mathbb{R}^{n \times n}}$ ；v为鲁棒项。令${\boldsymbol{\widetilde W}} = {{\boldsymbol{W}}^*} - {\boldsymbol{\hat W}}$ ，并定义权值矩阵的自适应律为：$$ {\boldsymbol{\dot {\hat W}}} = \delta {\text{φ}}({\boldsymbol{x}}){{\boldsymbol{r}}^{\rm{T}}} $$ (25) 式中，
$\delta $ 是一个正常数。将式（24）带入到式（21）中，得：$$ {\boldsymbol{M}}({\boldsymbol{q}}){{\dot {\boldsymbol{r}}}} = - \left( {{{\boldsymbol{K}}_{{\rm{v}}}} + {\boldsymbol{C}}({\boldsymbol{q}},{\boldsymbol{\dot q}})} \right){\boldsymbol{r}} + {{\boldsymbol{\zeta }}_1} $$ (26) 式中，
${{\boldsymbol{\zeta }}_1} = {{\boldsymbol{\widetilde W}}^{\rm{T}}}{\text{φ}}({\boldsymbol{x}}) + \left( {{\boldsymbol{\varepsilon }} + {{\boldsymbol{\tau }}_{{\rm{d}}}}} \right) + {\boldsymbol{v}}$ 。1.4.2 稳定性及收敛性分析
为了证明所设计的控制策略的动态稳定性，定义Lyapunov函数V为：
$$ V = \frac{1}{2}{{\boldsymbol{r}}^{\rm{T}}}{\boldsymbol{M}}({\boldsymbol{q}}){\boldsymbol{r}} + \frac{1}{{2\delta }}{tr} \left( {{{{\boldsymbol{\widetilde W}}}^{\rm{T}}}{\boldsymbol{\widetilde W}}} \right) $$ (27) 式中，
${tr}({{\boldsymbol{\widetilde W}}^{\rm{T}}}{\boldsymbol{\widetilde W}})$ 表示求矩阵$ {{\boldsymbol{\widetilde W}}^{\rm{T}}}{\boldsymbol{\widetilde W}} $ 的迹。计算V对时间的微分$\dot V $ 可得：$$ \dot V = {{\boldsymbol{r}}^{\rm{T}}}{\boldsymbol{M}}({\boldsymbol{q}}){\boldsymbol{\dot r}} + \frac{1}{2}{{\boldsymbol{r}}^{\rm{T}}}{\boldsymbol{\dot M}}({\boldsymbol{q}}){\boldsymbol{r}} + \frac{1}{\delta }{tr} \left( {{{{\boldsymbol{\widetilde W}}}^{\rm{T}}}{\boldsymbol{\dot {\widetilde W}}}} \right) $$ (28) 已知
${\boldsymbol{\dot M}}({\boldsymbol{q}}) - 2{\boldsymbol{C}}({\boldsymbol{q}},{\boldsymbol{\dot q}}) = 0$ ，将式（25）、（26）带入式（28）后，取鲁棒项v为${\boldsymbol{v}} = - \left( {{{\boldsymbol{\varepsilon }}_{\rm{N}}} + {\boldsymbol{d}}} \right){{\rm{sgn}}} ({\boldsymbol{r}})$ ，得：$$ \begin{aligned}[b] \dot V =& - {{\boldsymbol{r}}^{\rm{T}}}{{\boldsymbol{K}}_{{\rm{v}}}}{\boldsymbol{r}} + tr(({{{\boldsymbol{\widetilde W}}}^{\rm{T}}})\left({\text{φ}}({\boldsymbol{x}}){{\boldsymbol{r}}^{\rm{T}}} - \frac{1}{\delta }{\boldsymbol{\dot{ \hat W}}})\right) + \\ & {{\boldsymbol{r}}^{\rm{T}}}({\boldsymbol{\varepsilon }} + {{\boldsymbol{\tau }}_{\rm{d}}}) - {{\boldsymbol{r}}^{\rm{T}}}({{\boldsymbol{\varepsilon }}_{\rm{N}}} + {\boldsymbol{d}}) = - {{\boldsymbol{r}}^{\rm{T}}}{{\boldsymbol{K}}_{{\rm{v}}}}{\boldsymbol{r}} +\\& {{\boldsymbol{r}}^{\rm{T}}}({\boldsymbol{\varepsilon }} + {{\boldsymbol{\tau }}_{\rm{d}}}) - {{\boldsymbol{r}}^{\rm{T}}}({{\boldsymbol{\varepsilon }}_{\rm{N}}} + {\boldsymbol{d}}) \leqslant 0 \end{aligned} $$ 根据LaSalle不变性原理，闭环系统渐进稳定。
1.4.3 多模式控制策略设计
为保证安全，在上面控制器的基础上增加导纳–阻抗控制形成多模式切换控制策略。控制模式切换逻辑如图3所示，其中，
${{\boldsymbol{F}}_{{\rm{int}}}}$ 代表人机交互力，E代表位置误差。在训练过程中，若${{\boldsymbol{F}}_{{\rm{int}}}}$ 超过设定阈值${{\boldsymbol{F}}_{{\text{th}}}}$ ，切换为导纳控制使机器人跟随交互力顺应移动。当${{\boldsymbol{F}}_{{\rm{int}}}}$ 小于${{\boldsymbol{F}}_{{\text{th}}}}$ 且E大于设定阈值${{\boldsymbol{e}}_{{\text{th}}}}$ 时，启用阻抗控制将末端拉回到期望轨迹上，当E也小于${{\boldsymbol{e}}_{{\text{th}}}}$ 时重新切换为基于RBF网络的滑模控制。
图 3 轨迹跟踪过程模式切换逻辑Fig. 3 Mode switching logic of trajectory tracking下载:
全尺寸图片
本文提出的多模式柔顺跟踪控制策略的原理图如图4所示。轨迹定制部分用绿色框表示，由此得到的初始轨迹数据（
${\boldsymbol{q}}$ ，${\boldsymbol{\dot q}}$ ，${\boldsymbol{\ddot q}}$ ）经过数据处理和插值优化（图4中橙色框表示部分）得到期望的轨迹数据（${{\boldsymbol{q}}_{{\rm{d}}}}$ ，${{\boldsymbol{\dot q}}_{{\rm{d}}}}$ ，${{\boldsymbol{\ddot q}}_{{\rm{d}}}}$ ），然后输入到所设计的控制策略（图4中蓝色框表示部分）中进行轨迹跟踪康复训练。图4中，kp为阻抗系数。
图 4 多模式控制原理图Fig. 4 Schematic diagram of multi-mode control strategy下载:
全尺寸图片
2. 仿真与试验
2.1 轨迹定制拖动力试验
根据第1.4节提出的轨迹定制算法进行机器人末端拖动力试验。X和Y方向的最大拖动力的结果如图5所示。图5中的力仅表示力的大小，可以看出改变水平面内位置的力在4 N以内。
图 5 X和Y方向上的拖动力Fig. 5 Dragging force of X and Y direction下载:
全尺寸图片
Z方向的初始负载补偿和拖动力试验结果如图6所示。由两条点画线划分成的3个区域（Ⅰ、Ⅱ、Ⅲ）表示不同初始负载（分别约为3、5、8 N）的情况。可以看出3种情况下气缸都能保持在行程中的任意某个位置，且改变气缸位置的拖动力在5 N以内，即Z方向也可以轻松拖动。
图 6 垂直负载补偿和拖动力Fig. 6 Load compensation and dragging force of Z direction下载:
全尺寸图片
2.2 改进BOA算法试验
标准BOA和本文改进的BOA迭代优化过程如图7所示。由图7可知，改进的BOA算法具有更快的优化速度和更高的收敛精度，经过100次的迭代得到最小曲率和为13.377 4，对应的最优矢量数据压缩阈值为14.817 3。
图 7 BOA算法的迭代优化过程Fig. 7 Iterative optimization process of the BOA algorithm下载:
全尺寸图片
选取压缩阈值为20生成的插值曲线作为和优化插值曲线的对比，如图8所示。其中，小方框标记为7个未优化压缩阈值得到的型值点，小圆圈为9个优化压缩阈值得到的型值点。两条曲线产生的曲率的对比数据汇总于表1中，可以看出，优化曲线的曲率和更小为13.377 4，曲率的变化范围更小，标准差为0.224 2。
图 8 任意空间轨迹的NURBS插值曲线对比Fig. 8 Comparison of NURBS interpolation curves for arbitrary spatial trajectory下载:
全尺寸图片
表 1 插值曲线的曲率数据对比Table 1 Comparison of The Curvature Data for Interpolation Curves曲线类型 最小曲率 最大曲率 标准差 曲率总和 对比曲线 1.4879×10–5 17.8556 0.5882 33.6585 优化曲线 2.3985×10–5 6.8028 0.2242 13.3774 2.3 自由状态跟踪实验
机器人的第3关节垂直于前两个转动关节，存在结构解耦关系，可以分开对其进行控制。对气缸的位置控制采用PID加速度前馈控制。
自由状态轨迹跟踪实验中，1位健康志愿者模拟卒中初期患者，训练过程中不主动施加力。空间轨迹和水平面内的轨迹的跟踪结果如图9、10所示。由图9、10可以看出，实际轨迹可以较好的跟踪末端期望轨迹。在X、Y、Z 3个方向上的具体的跟踪误差如图11所示，由图11可以看出3个方向的误差都在6 mm内，由此可知，机器人的末端具有较好的轨迹跟踪能力。
图 9 末端自由状态下轨迹跟踪Fig. 9 Trajectory tracking under the free state of the end下载:
全尺寸图片
图 10 末端水平面轨迹跟踪Fig. 10 Trajectory tracking in the horizontal plane下载:
全尺寸图片
图 11 末端自由状态轨迹跟踪误差Fig. 11 Trajectory tracking error under the free state of the end下载:
全尺寸图片
2.4 控制模式切换实验
因为移动关节由气缸驱动，总行程并不大且气压驱动具有柔顺性。为了简化控制程序而仅考虑X和Y方向上有较大交互力的时候的顺应性。
本实验中，由1位健康的志愿者来模拟被动训练过程中突发较大人机交互力的情景。设置模式切换的交互力阈值
${F_{{\text{th}}}} 为10{\text{ N}}$ ，水平面内往返1个周期的柔顺跟踪结果如图12所示。 X和Y方向的位置跟踪情况分别如图13、14所示。
图 12 切换模式下的水平面内轨迹跟踪对随机干扰的响应Fig. 12 Response of the trajectory tracking to the random interference in the horizontal plane with switching mode下载:
全尺寸图片
图 13 切换模式下的X方向位置对随机干扰的响应Fig. 13 Response of the X–direction position under switch mode to random interference下载:
全尺寸图片
图 14 切换模式下的Y方向位置对随机干扰的响应Fig. 14 Response of the Y–direction position under switch mode to random interference下载:
全尺寸图片
以图13中标注的位置点为例：在 X方向上，当t1=30.77 s时，FX超过了阈值，机器人快速顺应交互力而偏离期望轨迹；t2=31.47 s时，FX小于阈值，机器人也能在0.20 s（t2–t3=0.16 s）内响应条件的改变而开始回归到期望轨迹上。图14中，在Y方向上，t2=20.53 s时，切换为导纳控制顺应交互力，交互力小于阈值时能在0.30 s（t3–t2=0.28 s）内响应条件的改变。最终机器人末端的实际位置都能回到期望轨迹上，继续完成被动康复训练。
3. 结 论
本文基于自行研制的末端牵引式3自由度上肢康复机器人进行末端拖动轨迹插值优化和轨迹柔顺跟踪控制研究，主要工作如下：
1）利用机器人末端安装的3维力传感器和关节伺服电机对初始负载和关节的静摩擦力补偿，基于静力学原理将末端力转换到关节空间，再采用导纳算法实现机器人拖动的方法避免了建立机器人的逆解模型。实现康复训练轨迹的轻松拖动定制。
2）采用矢量数据压缩算法保留了原始轨迹的最远点并保持轨迹的拓扑形状，不会降低上肢关节的运动范围。利用NURBS曲线插值和改进BOA算法优化插值曲线的曲率和，提升训练轨迹的平滑和光顺性。
3）结合RBF网络的滑模控制算法和导纳算法形成模式切换的柔顺轨迹跟踪控制策略。保证小干扰下的良好轨迹跟踪训练和人机交互力较大时患者的安全。跟踪精度可满足康复训练的需求。模式切换的响应快速，切换的阈值可根据康复医师的经验或者其他评定手段提前确定，控制器的设计简单，便于实施应用。
未来将研究人机交互力的识别和分离技术，实时补偿患者产生的作用力，减少康复医生的拖拽力，改善应用体验。此外，还将研究训练轨迹的能量最优规划及对于恢复了一定肌力的患者的按需辅助主动康复策略。
-
图 1 机器人结构简图
Fig. 1 Structure diagram of the robot
下载:
全尺寸图片
图 2 康复机器人试验平台
Fig. 2 Platform of the rehabilitation robot system
下载:
全尺寸图片
图 3 轨迹跟踪过程模式切换逻辑
Fig. 3 Mode switching logic of trajectory tracking
下载:
全尺寸图片
图 4 多模式控制原理图
Fig. 4 Schematic diagram of multi-mode control strategy
下载:
全尺寸图片
图 5 X和Y方向上的拖动力
Fig. 5 Dragging force of X and Y direction
下载:
全尺寸图片
图 6 垂直负载补偿和拖动力
Fig. 6 Load compensation and dragging force of Z direction
下载:
全尺寸图片
图 7 BOA算法的迭代优化过程
Fig. 7 Iterative optimization process of the BOA algorithm
下载:
全尺寸图片
图 8 任意空间轨迹的NURBS插值曲线对比
Fig. 8 Comparison of NURBS interpolation curves for arbitrary spatial trajectory
下载:
全尺寸图片
图 9 末端自由状态下轨迹跟踪
Fig. 9 Trajectory tracking under the free state of the end
下载:
全尺寸图片
图 10 末端水平面轨迹跟踪
Fig. 10 Trajectory tracking in the horizontal plane
下载:
全尺寸图片
图 11 末端自由状态轨迹跟踪误差
Fig. 11 Trajectory tracking error under the free state of the end
下载:
全尺寸图片
图 12 切换模式下的水平面内轨迹跟踪对随机干扰的响应
Fig. 12 Response of the trajectory tracking to the random interference in the horizontal plane with switching mode
下载:
全尺寸图片
图 13 切换模式下的X方向位置对随机干扰的响应
Fig. 13 Response of the X–direction position under switch mode to random interference
下载:
全尺寸图片
图 14 切换模式下的Y方向位置对随机干扰的响应
Fig. 14 Response of the Y–direction position under switch mode to random interference
下载:
全尺寸图片
表 1 插值曲线的曲率数据对比
Table 1 Comparison of The Curvature Data for Interpolation Curves
曲线类型 最小曲率 最大曲率 标准差 曲率总和 对比曲线 1.4879×10–5 17.8556 0.5882 33.6585 优化曲线 2.3985×10–5 6.8028 0.2242 13.3774 -
[1] Gassert R,Dietz V.Rehabilitation robots for the treatment of sensorimotor deficits:A neurophysiological perspective[J].Journal of Neuroengineering and Rehabilitation,2018,15(1):1–15. doi: 10.1186/s12984-018-0383-x[PubMed [2] Lin D J,Finklestein S P,Cramer S C.New directions in treatments targeting stroke recovery[J].Stroke,2018,49(12):3107–3114. doi: 10.1161/STROKEAHA.118.021359[PubMed [3] Morita Y,Nagasaki M,Ukai H,et al.Development of rehabilitation training support system of upper limb motor function for personalized rehabilitation[C]//Proceedings of the 2008 IEEE International Conference on Robotics and Biomimetics.Bangkok:IEEE,2009:300−305. [4] 游有鹏,张宇,李成刚.面向直接示教的机器人零力控制[J].机械工程学报,2014,50(3):10–17. doi: 10.3901/JME.2014.03.010 You Youpeng,Zhang Yu,Li Chenggang.Force-free control for the direct teaching of robots[J].Journal of Mechanical Engineering,2014,50(3):10–17 doi: 10.3901/JME.2014.03.010 [5] 杨浩,韩建海,李向攀.卧式下肢康复训练机器人力觉拖动示教研究[J].机械设计与制造,2020,351(5):272–275. doi: 10.19356/j.cnki.1001-3997.2020.05.066 Yang Hao,Han Jianhai,Li Xiangpan.Research on drag and teach of horizontal lower limb rehabilitative robot[J].Machinery Design & Manufacture,2020,351(5):272–275 doi: 10.19356/j.cnki.1001-3997.2020.05.066 [6] Zeestraten M J A,Havoutis I,Silvério J,et al.An approach for imitation learning on Riemannian manifolds[J].IEEE Robotics and Automation Letters,2017,2(3):1240–1247. doi: 10.1109/LRA.2017.2657001[LinkOut [7] 冯峰,蒋维.一种带约束限制的三次B样条曲线矢量数据压缩算法[J].数学杂志,2021,41(3):247–256. doi: 10.13548/j.sxzz.2021.03.007 Feng Feng,Jiang Wei.A cubic b-spline-based vector data compression algorithm with boundary constraints[J].Journal of Mathematics,2021,41(3):247–256 doi: 10.13548/j.sxzz.2021.03.007 [8] Wan Neng,Xu Degang,Ye Huawen.Improved cubic B-spline curve method for path optimization of manipulator obstacle avoidance[C]//Proceedings of the 2018 Chinese Automation Congress (CAC).Xi’an:IEEE,2019:1471−1476. [9] 董甲甲,王太勇,董靖川,等.改进B样条曲线应用于6R机器人轨迹优化[J].中国机械工程,2018,29(2):193–200. doi: 10.3969/j.issn.1004-132X.2018.02.010 Dong Jiajia,Wang Taiyong,Dong Jingchuan,et al.Applications of improved B-spline curves to 6R robot trajectory optimization[J].China Mechanical Engineering,2018,29(2):193–200 doi: 10.3969/j.issn.1004-132X.2018.02.010 [10] Mei Jiangping,Zhang Fan,Zang Jiawei,et al.Trajectory optimization of the 6-degrees-of-freedom high-speed parallel robot based on B-spline curve[J].Science Progress,2020,103(1):36850419880115. doi: 10.1177/0036850419880115[PubMed [11] Ijspeert A J,Nakanishi J,Hoffmann H,et al.Dynamical movement primitives:Learning attractor models for motor behaviors[J].Neural Computation,2013,25(2):328–373. doi: 10.1162/NECO_a_00393[PubMed [12] 吴青聪,王兴松,吴洪涛,等.上肢康复外骨骼机器人的模糊滑模导纳控制[J].机器人,2018,40(4):457–465. doi: 10.13973/j.cnki.robot.18093 Wu Qingcong,Wang Xingsong,Wu Hongtao,et al.Fuzzy sliding mode admittance control of the upper limb rehabilitation exoskeleton robot[J].Robot,2018,40(4):457–465 doi: 10.13973/j.cnki.robot.18093 [13] Mushage B O,Chedjou J C,Kyamakya K.Fuzzy neural network and observer-based fault-tolerant adaptive nonlinear control of uncertain 5-DOF upper-limb exoskeleton robot for passive rehabilitation[J].Nonlinear Dynamics,2017,87(3):2021–2037. doi: 10.1007/s11071-016-3173-7[LinkOut [14] 罗定吉,高学山,李健,等.下肢外骨骼康复机器人轨迹跟踪控制[J].兵工自动化,2020,39(11):87–91. doi: 10.7690/bgzdh.2020.11.021 Luo Dingji,Gao Xueshan,Li Jian,et al.Trajectory tracking control of lower limb exoskeleton rehabilitation robot[J].Ordnance Industry Automation,2020,39(11):87–91 doi: 10.7690/bgzdh.2020.11.021 [15] Shi Jia,Xu Linsen,Cheng Gaoxin,et al.Trajectory tracking control based on RBF neural network of the lower limb rehabilitation robot[C]//Proceedings of the 2020 IEEE International Conference on Mechatronics and Automation (ICMA).Beijing:IEEE,2020:117−123. [16] Zhang Songyuan,Guo Shuxiang,Fu Yili,et al.Integrating compliant actuator and torque limiter mechanism for safe home-based upper-limb rehabilitation device design[J].Journal of Medical and Biological Engineering,2017,37(3):357–364. doi: 10.1007/s40846-017-0228-2[LinkOut [17] Wu Qingcong,Chen Bai,Wu Hongtao.Adaptive admittance control of an upper extremity rehabilitation robot with neural-network-based disturbance observer[J].IEEE Access,2019,7:123807–123819. doi: 10.1109/ACCESS.2019.2938566[LinkOut [18] Trigili E,Crea S,Moisè M,et al.Design and experimental characterization of a shoulder-elbow exoskeleton with compliant joints for post-stroke rehabilitation[J].IEEE/ASME Transactions on Mechatronics,2019,24(4):1485–1496. doi: 10.1109/TMECH.2019.2907465[LinkOut [19] Miao Qing,Peng Yuxin,Liu Li,et al.Subject-specific compliance control of an upper-limb bilateral robotic system[J].Robotics and Autonomous Systems,2020,126:103478. doi: 10.1016/j.robot.2020.103478[LinkOut [20] Guo Bingjing,Han Jianhai,Li Xiangpan,et al.Human–robot interactive control based on reinforcement learning for gait rehabilitation training robot[J].International Journal of Advanced Robotic Systems,2019,16(2):172988141983958. doi: 10.1177/1729881419839584 [21] 李超.外骨骼下肢康复机器人系统开发及其柔顺控制研究[D].杭州:浙江理工大学,2019. Li Chao.Research on the development of exoskeleton lower limb rehabilitation robot system and its compliance control[D].Hangzhou:Zhejiang Sci-Tech University,2019. [22] 毋东,王庆,王慧青.基于样条的矢量数据多尺度表达与压缩算法[J].计算机工程,2012,38(8):201–203. doi: 10.3969/j.issn.1000-3428.2012.08.066 Wu Dong,Wang Qing,Wang Huiqing.Multi-scale representation and compression algorithm for vector data based on spline[J].Computer Engineering,2012,38(8):201–203 doi: 10.3969/j.issn.1000-3428.2012.08.066 [23] 霍亚光,高扬,宋绪丁.不同参数化法对三次NURBS曲线拟合误差的影响[J].机电工程技术,2019,48(4):54–57. doi: 10.3969/j.issn.1009-9492.2019.04.018 Huo Yaguang,Gao Yang,Song Xuding.Effect of different parameterization methods on the cubic NURBS curve fitting errors[J].Mechanical & Electrical Engineering Technology,2019,48(4):54–57 doi: 10.3969/j.issn.1009-9492.2019.04.018 [24] Hashemian A,Hosseini S F.An integrated fitting and fairing approach for object reconstruction using smooth NURBS curves and surfaces[J].Computers & Mathematics With Applications,2018,76(7):1555–1575. doi: 10.1016/j.camwa.2018.07.007[LinkOut [25] 高文欣,刘升,肖子雅,等.柯西变异和自适应权重优化的蝴蝶算法[J].计算机工程与应用,2020,56(15):43–50. doi: 10.3778/j.issn.1002-8331.1907-0048 Gao Wenxin,Liu Sheng,Xiao Ziya,et al.Butterfly optimization algorithm based on cauchy variation and adaptive weight[J].Computer Engineering and Applications,2020,56(15):43–50 doi: 10.3778/j.issn.1002-8331.1907-0048 [26] 王依柔,张达敏.融合正弦余弦和无限折叠迭代混沌映射的蝴蝶优化算法[J].模式识别与人工智能,2020,33(7):660–669. doi: 10.16451/j.cnki.issn1003-6059.202007008 Wang Yirou,Zhang Damin.Butterfly optimization algorithm combining sine cosine and iterative chaotic map with infinite collapses[J].Pattern Recognition and Artificial Intelligence,2020,33(7):660–669 doi: 10.16451/j.cnki.issn1003-6059.202007008 [27] Tubishat M,Alswaitti M,Mirjalili S,et al.Dynamic butterfly optimization algorithm for feature selection[J].IEEE Access,2020,8:194303–194314. doi: 10.1109/ACCESS.2020.3033757[LinkOut [28] Naghdi S,Ansari N N,Mansouri K,et al.A neurophysiological and clinical study of Brunnstrom recovery stages in the upper limb following stroke[J].Brain Injury,2010,24(11):1372–1378. doi: 10.3109/02699052.2010.506860[PubMed
