CN101837164A - 功能性电刺激中pid参数的双源特征融合蚁群整定方法 - Google Patents

Abstract

本发明涉及以电脉冲刺激进行肢体康复的器械领域。为实现准确稳定实时地控制FES***地电流强度，有效地提高FES***准确性和稳定性，本发明采用的技术方案是：功能性电刺激中PID参数的双源特征融合蚁群整定方法，包括下列步骤：首先，利用助行过程的肌肉模型HRV预测膝关节角度；其次，利用蚁群算法整定PID参数，实时调控FES电流水平强度；采用蚁群算法对PID参数进行控制，如未达到预期目标继续寻优；在新的PID系数下计算***输出及其与肌肉模型HRV的偏差后再进入下一步蚁群算法的自学习与加权系数自调整；反复前一步骤，实现PID控制参数的自适应在线整定，并用于FES***。本发明主要应用于整定PID参数。

Description

功能性电刺激中PID参数的双源特征融合蚁群整定方法

技术领域

本发明涉及以电脉冲刺激进行肢体康复的器械领域，尤其是功能性电刺激中PID参数的双源特征融合蚁群整定方法。

背景技术

功能性电刺激(Functional Electrical Stimulation，FES)是通过电流脉冲序列来刺激肢体运动肌群及其外周神经，有效地恢复或重建截瘫患者的部分运动功能的技术。据统计，由于脊髓再生能力微弱，针对脊髓损伤瘫痪患者，目前尚未有可直接修复损伤的有效医治方法，实施功能康复训练是一有效的措施。脊髓损伤瘫痪患者人数逐年增多，功能康复训练是亟待需求的技术。20世纪60年代，Liberson首次成功地利用电刺激腓神经矫正了偏瘫患者足下垂的步态，开创了功能性电刺激用于运动和感觉功能康复治疗的新途径。目前，FES已经成为了恢复或重建截瘫患者的部分运动功能，是重要的康复治疗手段。然而，如何精密控制FES的触发时序和脉冲电流强度以保证电刺激作用效果能准确完成预定的功能动作仍是FES的技术关键。据统计，目前FES的触发控制的方式研究尚少，而且根据作用效果与预定动作偏差，用闭环控制来自动调整FES刺激强度和时序参数，从而大大提高了FES***的准确性和稳定性，但是现在有效的控制方法仍然在探索之中。

柄反作用矢量(handle reactions vector，HRV)是根据步行器帮助下的站立及行走的过程中，步行器提供给患者的效用实际上可以分为明确独立的3个部分：前后向的力推进，左右向的力平衡和上下向的力支持，这其实也可理解为患者为维持自身正常站立行走对外界所需的附加力学诉求提出的新概念，即是患者在站立行走过程中对步行器的作用合成简化为集中载荷，分别用手柄中点横截面形心处的两个力学矢量来表示，如图1所示，矢量在x，y，z轴上的方向分量合力大小可以分别表征患者借助步行器所获得的力推进，力平衡和力支持水平。其中，定义坐标系所设定的x轴正向为患者的右向，y轴正向为患者的前向，z轴正向为患者的上向。这样，HRV的定义公式也可以写为：

[HRV]＝[HRV₁，HRV_r]^T＝[F_lx，F_ly，F_lz，F_rx，F_ry，F_rz]^T                   (1)

目前，HRV被广泛地应用在监视在电刺激过程中病人行走时的状况，继而防止病人摔倒，造成二次伤害。本专利提出利用此参数预测膝关节角度，继而精密控制FES***的电流水平强度，保证电刺激作用效果能准确完成预定的功能动作，并且防止肌疲劳。

比例微积分(proportional-integral-differential，PID)是一种非常实用的反馈调节算法，它根据***检测或操作偏差，利用比例、积分、微分运算获得所需调节量以对***进行反馈控制，因其操作方便而广泛用于工程实践。尤其当被控***特性参数不明确或难以及时在线测定时，稳妥的闭环控制即可采用PID整定算法。面对肌肉的复杂性和时变性操作环境，由于PID的稳定性好、工作可靠，目前仍在功能性电刺激领域得到了广泛的应用。PID核心技术是精密确定其中比例、积分、微分系数，尤其在FES领域，对***稳定性要求极为严格，所以对PID参数选择尤为重要。PID控制要取得较好的控制效果，必须调整好比例、积分和微分三种控制作用，形成控制量中既相互配合又相互制约的关系。

发明内容

为克服现有技术的不足，提供一种功能性电刺激中PID参数的双源特征融合蚁群整定方法，能够准确稳定实时地控制FES***地电流强度，有效地提高FES***准确性和稳定性，并获得可观的社会效益和经济效益。为达到上述目的，本发明采用的技术方案是：功能性电刺激中PID参数的双源特征融合蚁群整定方法，包括下列步骤：

首先，利用助行过程的肌肉模型HRV预测膝关节角度；

其次，利用蚁群算法整定PID参数，实时调控FES电流水平强度，其整定流程为：首

先根据PID的三个决策变量K_p、K_i和K_d取值范围的上下界，确定包括蚁群群体规模、搜索空间维数的参数，并对其进行编码，然后利用通过实际关节角度与肌肉模型HRV输出关节角度的相应关系作为适度评价函数计算的适应度值；采用蚁群算法对PID参数进行控制，即确定蚁群算法的参数设置，利用蚂蚁随机搜索使其变量优化PID的K_p、K_i和K_d三个系数，利用适应度函数调节蚂蚁每次寻索路径以及判断是否达到预设目标，如达到预设目标，计算最终最佳的位置即得PID的K_p、K_i和K_d三个系数，如未达到预期目标继续寻优，直到达到预设目标；在新的PID系数下计算***输出yout及其与肌肉模型HRV的偏差后再进入下一步蚁群算法的自学习与加权系数自调整；

反复前一步骤，最终实现PID控制参数的自适应在线整定，并用于FES***。

所述实际关节角度与肌肉模型HRV输出关节角度的相应关系为：

L＝M×HRV^-1                           (3)

M表示膝关节角度，HRV表示使用者施加在步行器上力的柄反作用矢量，L表示HRV与M之间的关系，采用偏最小二乘回归的方法求解L：

设有m个HRV变量HRV1，…，HRVm，p个M变量，M1，…，Mp，共i(i＝1，…，n)个观测值的数据集，T、U分别为从HRV变量与M变量中提取的成分，从原始变量集中提取第一对成分T1、U1的线性组合为：

T₁＝ω₁₁HRV₁+…+ω_1mHRV_m＝ω₁′HRV                                (4)

        U₁＝v₁₁M₁+…+v_1pM_p＝v₁′M                                 (5)

其中ω₁＝(ω₁₁，…，ω_1m )′为模型效应权重，v₁＝v₁₁，…，v_1p)′为M变量权重，将上述提取第一成分的要求转化为求条件极值问题：

其中t1、u1为由样本求得的第一对成分的得分向量，HRV0、M0为初始变量，利用拉格朗日乘子法，上述问题转化为求单位向量ω₁和v₁，使θ₁＝ω₁′HRV₀′M₀v₁最大，即求矩阵HRV₀′M₀M₀′HRV₀的特征值和特征向量，其最大特征值为θ₁ ²，相应的单位特征向量就是所求的解ω₁，而v₁由公式

得到；

其次建立初始变量对T1的方程

$\left\{\begin{array}{c}{\mathrm{HRV}}_0=t_1{\mathrm{\α}}_1^{\′}+E_1\\ M_0=t_1{\mathrm{\β}}_1^{\′}+F_1\end{array}\right.---\left(7\right)$

其中t1意义同前，α₁′＝(α₁₁，…，α_1m)，β₁′＝β₁₁，…，β_1p)为仅一个M量t1时的参数向量，E1、F1分别为n×m和n×p残差阵，按照普通最小二乘法可求得系数向量α₁和β₁，其中α₁成为模型效应载荷量；

如提取的第一成分不能达到回归模型的精度，运用残差阵E1、F1代替X0、Y0，重复提取成分，依次类推，假设最终提取了r个成分，HRV0、M0对r个成分的回归方程为：

$\left\{\begin{array}{c}{\mathrm{HRV}}_0=t_1{\mathrm{\α}}_1^{\′}+...+t_r{\mathrm{\α}}_r^{\′}+E_r\\ M_0=t_1{\mathrm{\β}}_1^{\′}+...+t_r{{\mathrm{\β}}_r}^{\′}+F_r\end{array}\right.---\left(8\right)$

把第一步分析所得HRV量中提取成分Tk(k＝1，…，r)线性组合带入M量对r个成分建立的回归方程，即t_r＝ω_k1HRV₁+…+ω_kmHRV_m代入M_j＝t₁β_1j+…+t_rβ_rj(j＝1，…，p)，即得标准化变量的回归方程M_j＝α_j1HRV₁+…+α_jmHRV_m；

最后根据式(3)，即可求出L。

所述进行编码是，根据膝关节角度和电流模式，以及实际误差情况设定PID的三个参数为5位有效数字，其中K_p小数点前2位，小数点后3位；K_i和K_d小数点前1位，小数点后4位，具体编码如以下公式：

K_p＝y_1，j×10¹+y_2，j×10⁰+y_3，j×10^-1+y_4，j×10^-2+y_5，j×10^-3

K_i＝y_6，j×10⁰+y_7，j×10^-1+y_8，j×10^-2+y_9，j×10^-3+y_10，j×10^-4

K_d＝y_6，j×10⁰+y_7，j×10^-1+y_8，j×10^-2+y_9，j×10^-3+y_10，j×10^-4

并且根据输出的超调等误差设计适应函数具体函数为：

fit＝α1×σ+β1×t+c×error                              (9)

其中，σ为超调量，t为上升时间，error为输出关节角度与预设关节角度的相对误差，α1＝0.1，β1＝0.8，c＝2；

蚁群算法整定PID控制器参数的具体工作流程为：

Step1：参数初始化，令时间t＝0和循环次数N_max＝0，设置最大循环次数N_cmax，将m个蚂蚁置于起始点；

Step2：设置蚂蚁个数和循环次数；

Step3：蚂蚁随机搜索，在一次爬行结束后，决定哪些特征变量被选中作为实际输入变量，修改禁忌表指针，即选择好之后将蚂蚁移动到新的元素，并把该元素移动该蚂蚁个体的禁忌表中；

Step4：计算蚂蚁个体据状态转移概率公式计算的概率，根据此概率选择元素；

Step5：若蚂蚁元素未遍历完，转Step3，否则为Step6；

Step6：更新信息素浓度划分正确率高的特征变量的信息素浓度得到增强，下次搜索时会以更大的概率被选中；

Step7：满足结束调节，整定结束。

本发明的特点在于：利用助行器的HRV变化预测膝关节角度变化，然后通过蚁群群算法优化PID的比例系数、微分系数以及积分系数，继而控制FES***的电流脉冲强度，有效地提高了FES***准确性和稳定性。

附图说明

图1柄反作用矢量(HRV)定义示意图。

图2基于HRV的FES***结构框图。

图3蚁群群算法整定PID参数控制方法的结构框图。

图4助行功能性电刺激中的人体模型。

图5蚁群算法整定PID参数编码示意图。

图6实验场景。

图7蚁群群算法整定的PID控制追踪结果。

具体实施方式

本发明的主旨是提出一种新的FES的精密控制方法，利用步行器的HRV参数预测的膝关节角度与实际的膝关节角度预测的关节角度的误差，通过蚁群算法优化PID的比例系数、积分系数以及微分系数，继而准确稳定实时地控制FES***地电流强度。该项发明可有效地提高FES***准确性和稳定性，并获得可观的社会效益和经济效益。

基于HRV的功能性电刺激助行中的精密控制新技术的应用的结构如图2所示，其工作流程为：首先，利用助行过程的HRV预测膝关节角度，其次，利用蚁群算法整定PID参数，实时调控FES电流水平强度。其整定结构示意图如图3所示，为：首先根据PID的三个决策变量K_p、K_i和K_d取值范围的上下界，确定蚁群群体规模、搜索空间维数等参数，并对其进行编码，然后利用通过实际关节角度与肌肉模型输出关节角度的相应关系作为适度评价函数计算的适应度值，以及确定蚁群算法的参数设置，利用蚂蚁随机收索使其变量优化PID的K_p、K_i和K_d三个系数，利用适应度函数调节蚂蚁每次寻索路径以及判断是否达到预设目标。如达到预设目标，计算最终最佳的位置即得PID的K_p、K_i和K_d三个系数，如未达到预期目标继续寻优，直到达到预设目标。在新的PID系数下计算***输出yout及其与肌肉模型的偏差后再进入下一步蚁群算法的自学习与加权系数自调整。反复此过程，最终实现PID控制参数的自适应在线整定，并用于FES***。

1HRV预测膝关节角度模型

助行过程中，当使用者在功能性电刺激作用下，抬腿迈步时，为了支持身体稳定，使用者在步行器上所施加的力则有所不同，因为关节的大小不同会使人体重心处于不同位置，则克服重力所施加的力也不同，同时人体所处的平面位置也有所改变，为了位置倾翻则在平面上所施加的力也有所变化，因此，关节角度和使用者对步行器所施加的力有一定的关系，如图4所示。

M＝L·HRV+wPW                                                    (1)

其中，M表示膝关节角度，HRV表示使用者施加在步行器上力的柄反作用矢量，L表示HRV与M之间的关系，w表示系数，W表示上臂、躯干和下肢的重心，P表示三重心与M之间的关系。

实际中，由于步行器的作用，人体重心移动较小，膝关节角度则可表示成

M＝L·HRV                                                            (2)

其中，M表示膝关节角度，HRV表示使用者施加在步行器上力的柄反作用矢量，L表示HRV与M之间的关系。根据式2所示，确定L就可以利用HRV取出相应时刻的膝关节角度。

L＝M×HRV^-1                                                      (3)

本发明求解L时，采用了偏最小二乘回归的方法。

设有m个HRV变量HRV1，…，HRVm，p个M变量，M1，…，Mp，共i(i＝1，…，n)个观测值的数据集。T、U分别为从HRV变量与M变量中提取的成分，这里提取的成分通常称为偏最小二乘因子。

从原始变量集中提取第一对成分T1、U1的线性组合为：

T₁＝ω₁₁HRV₁+…+ω_1mHRV_m＝ω₁′HRV                               (4)

U₁＝v₁₁M₁+…+v_1pM_p＝v₁′M                                        (5)

其中ω₁＝(ω₁₁，…，ω_1m)′为模型效应权重，v₁(v₁₁，…，v_1p)′为M变量权重。为保证T1、U1各自尽可能多地提取所在变量组的变异信息，同时保证两者之间的相关程度达到最大，据成分的协方差可由相应成分的得分向量的内积来计算的性质，上述提取第一成分的要求转化为求条件极值问题。

其中t1、u1为由样本求得的第一对成分的得分向量，HRV0、M0为初始变量。利用拉格朗日乘子法，上述问题转化为求单位向量ω₁和v₁，使θ₁＝ω₁′HRV₀′M₀v₁最大，即求矩阵HRV₀′M₀M₀′HRV₀的特征值和特征向量，其最大特征值为θ₁ ²，相应的单位特征向量就是所求的解ω₁，而v₁由公式

得到。

其次建立初始变量对T1的方程

$\left\{\begin{array}{c}{\mathrm{HRV}}_0=t_1{\mathrm{\α}}_1^{\′}+E_1\\ M_0=t_1{\mathrm{\β}}_1^{\′}+F_1\end{array}\right.---\left(7\right)$

其中t1意义同前，α₁′＝(α₁₁，…，α_1m)，β₁′＝β₁₁，…，β_1p)为仅一个M量t1时的参数向量，E1、F1分别为n×m和n×p残差阵。按照普通最小二乘法可求得系数向量α₁和β₁，其中α₁成为模型效应载荷量。

如提取的第一成分不能达到回归模型的精度，运用残差阵E1、F1代替X0、Y0，重复提取成分，依次类推。假设最终提取了r个成分，HRV0、M0对r个成分的回归方程为：

$\left\{\begin{array}{c}{\mathrm{HRV}}_0=t_1{\mathrm{\α}}_1^{\′}+...+t_r{\mathrm{\α}}_r^{\′}+E_r\\ M_0=t_1{\mathrm{\β}}_1^{\′}+...+t_r{{\mathrm{\β}}_r}^{\′}+F_r\end{array}\right.---\left(8\right)$

把第一步分析所得HRV量中提取成分Tk(k＝1，…，r)线性组合带入M量对r个成分建立的回归方程，即t_r＝ω_k1HRV₁+…+ω_kmHRV_m代入M_j＝t₁β_1j+…+t_rβ_rj(j＝1，…，p)，即得标准化变量的回归方程M_j＝α_j1HRV₁+…+α_jmHRV_m。

最后根据式3，即可求出L。

蚁群算法对PID参数的控制，首先对PID的三个参数进行编码，根据膝关节角度和电流模式，以及实际误差等情况设定PID的三个参数为5位有效数字，其中K_p小数点前2位，小数点后3位；K_i和K_d小数点前1位，小数点后4位。其具体编码示意图如图5所示。

如以下公式

K_p＝y_1，j×10¹+y_2，j×10⁰+y_3，j×10^-1+y_4，j×10^-2+y_5，j×10^-3

K_i＝y_6，j×10⁰+y_7，j×10^-1+y_8，j×10^-2+y_9，j×10^-3+y_10，j×10^-4

K_d＝y_6，j×10⁰+y_7，j×10^-1+y_8，j×10^-2+y_9，j×10^-3+y_10，j×10^-4

并且根据输出的超调等误差设计适应函数具体函数为：

fit＝α1×σ+β1×t+c×error                                  (9)

其中，σ为超调量，t为上升时间，error为输出关节角度与预设关节角度的相对误差，α1＝0.1，β1＝0.8，c＝2。

蚁群算法整定PID控制器参数的具体工作流程为：

Step1：参数初始化。令时间t＝0和循环次数N_max＝0，设置最大循环次数N_cmax，将m个蚂蚁置于起始点。

Step2：设置蚂蚁个数和循环次数

Step3：蚂蚁随机搜索，在一次爬行结束后，决定哪些特征变量被选中作为实际输入变量，修改禁忌表指针，即选择好之后将蚂蚁移动到新的元素，并把该元素移动该蚂蚁个体的禁忌表中

Step4：计算蚂蚁个体据状态转移概率公式计算的概率，根据此概率选择元素

Step5：若蚂蚁元素未遍历完，转Step3，否则为Step6

Step6：更新信息素浓度划分正确率高的特征变量的信息素浓度得到增强，下次搜索时会以更大的概率被选中

Step7：满足结束调节，整定结束。

2蚁群群算法整定PID参数的控制

PID由比例单元P、积分单元I和微分单元D三部分组成，根据***的误差，通过设定的K_p、K_i和K_d三个参数对***进行控制。

$\mathrm{yout}\left(t\right)=K_p\mathrm{error}\left(t\right)+K_i\underset{j=0}{\overset{t}{\mathrm{\Σ}}}\mathrm{error}\left(j\right)+K_d[\mathrm{error}\left(t\right)-\mathrm{error}(t-1)]---\left(9\right)$

其中K_p是比例系数，K_i是积分系数，K_d是微分系数，error为预设输出与实际输出的偏差，u(t)为PID的输出，同时又是受控***的输入。

由PID输出公式(1)可以得到

$u(t-1)=K_p\mathrm{error}(t-1)+K_i\underset{j=0}{\overset{t-1}{\mathrm{\Σ}}}\mathrm{error}\left(j\right)+K_d[\mathrm{error}(t-1)-\mathrm{error}(t-2)]---\left(10\right)$

根据：

Δu(t)＝u(t)-u(t-1)

＝K_p(error(t)-error(t-1))+K_ierror(t)+K_d(error(t)-2error(t-1)+error(t-2))

……………………………………………………………(11)

有：

u(t)＝Δu(t)+u(t-1)＝

u(t-1)+K_p(error(t)-error(t-1))+K_ierror(t)+K_d(error(t)-2error(t-1)+error(t-2))

                                                             ………………(12)

本发明采用蚁群算法进行PID控制参数的自适应优化，把PID的三个参数作为一个组合，利用蚁群算法寻优来解决这个组合问题。蚁群算法是一种源于大自然生物世界的新型仿生算法，用蚁群算法求解最优化问题时，首先将最优化问题转化为了求解最短路径问题。每只蚂蚁从初始接点N₀₀，N₀₁...N_0n出发，顺序走过N₁，N₂...，的其中一子结点，直到终结点N_k0、N_k1...N_kn组成路径(N_0t N_1t...N_kt)，t∈[0，1，2…9]。其路径可代表一个二进制的可行解。每次蚂蚁访问城市时有以下的特征：

状态转化规则：蚁群算法使用的状态转化规则为基于TSP问题提出的随机比例规则，它给出位于城市i的蚂蚁k选择移动到城市j的概率，

其中τ_ij(i，j)为(i，j)的适应度，η_ij(i，j)＝(10-|y(i)-y(i)^*|)/10，y(i)蚁群搜索时在i处的值，y(i)^*为上次搜索时在i处的值。α为残留信息的相对重要程度、β为期望值的相对重要程度。

在蚁群算法中，选择方式为

其中，q为均匀分布在[0，1]上的一个随机数，q₀为[0，1]上的参变量。

全局更新规则：蚂蚁算法有不同的更新算法，蚁群***采用的全局更新原则，只允许全局最优解的蚂蚁释放信息素，这样是为了使蚂蚁的搜索主要集中在当前循环为止所找出的最好路径的邻域。

τ_ij(i，j)←(1-ρ)□τ_ij(i，j)+ρ·Δτ_ij(i，j)                       (14)

其中ρ为信息数挥发系数，L_gb为目前为止找到的全局最优路径

局部更新信息：每只蚂蚁建立一个解的过程中也有进行信息数素迹的更新

τ_ij(i，j)←(1-γ)□τ_ij(i，j)+γ·Δτ_ij(i，j)                       (16)

其中γ∈[0，1]。

蚁群算法对PID参数的控制，首先对PID的三个参数进行编码，根据膝关节角度和电流模式，以及实际误差等情况设定PID的三个参数为5位有效数字，其中K_p小数点前2位，小数点后3位；K_i和K_d小数点前1位，小数点后4位。其具体编码示意图如图5所示。

如以下公式

K_p＝y_1，j×10¹+y_2，j×10⁰+y_3，j×10^-1+y_4，j×10^-2+y_5，j×10^-3

K_i＝y_6，j×10⁰+y_7，j×10^-1+y_8，j×10^-2+y_9，j×10^-3+y_10，j×10^-4

K_d＝y_6，j ×10⁰+y_7，j×10^-1+y_8，j ×10^-2+y_9，j ×10^-3+y_10，j×10^-4

蚁群算法整定PID控制器参数的具体工作流程为：

Step1：参数初始化。令时间t＝0和循环次数N_max＝0，设置最大循环次数N_cmax，将m个蚂蚁置于起始点。

Step2：设置蚂蚁个数和循环次数

Step3：蚂蚁随机搜索，在一次爬行结束后，决定哪些特征变量被选中作为实际输入变量，修改禁忌表指针，即选择好之后将蚂蚁移动到新的元素，并把该元素移动该蚂蚁个体的禁忌表中

Step4：计算蚂蚁个体据状态转移概率公式计算的概率，根据此概率选择元素

Step5：若蚂蚁元素未遍历完，转Step3，否则为Step6

Step6：更新信息素浓度划分正确率高的特征变量的信息素浓度得到增强，下次搜索时会以更大的概率被选中

Step7：满足结束调节，整定结束。

3实验方案

实验装置采用无线传输的步行器***和美国SIGMEDICS公司生产的Parastep功能性电刺激***，该***包含微处理器和刺激脉冲发生电路，含六条刺激通道，电池供电。实验内容为：利用FES***对下肢相关肌群进行刺激，使受试者按照预定的动作，同时记录施加在步行器上的HRV首先通过安装在步行器上的12导联应变片(BX350-6AA)电桥网络转化成的电压信号以及膝关节角度运动轨迹。要求受试者身体健康，无下肢肌肉、骨骼疾患，无神经疾患及严重心肺疾患。实验时受试者安坐于步行器前，

将刺激电极固定于相应的部位，未施加电刺激时，受试者保持轻松。FES实验场景如图5所示。电刺激脉冲序列采用经典的Lilly波形，脉冲频率为25Hz、脉宽150μs，脉冲电流在0～120m范围内可调。实验中，实时地记录HRV以

及可通过改变脉冲电流大小来调整刺激强度以改变由刺激产生的膝关节角度。实验前，设定期望的膝关节角度运动轨迹，实验中利用角度测量计实时检测膝关节张角变化。实验数据采样率为128Hz，数据记录时长为60s。

有益效果

蚁群算法整定PID参数的新算法对FES脉冲电流幅值进行测算和调整，使FES作用所产生的膝关节角度运动贴近预期的运动轨迹。图7为蚁群算法整定的PID控制追踪结果。图中红线表示预期运动轨迹、蓝线为实际输出关节角度。X轴为时间，Y轴为膝关节运动角度。为更清楚地观察蚁群算法整定PID的控制误差，如蚁群算法整定PID控制下预设输入膝关节角度与实际膝关节角度的相对误差所示，则可以看出误差均在5％之内，可以达到精确的控制。

本发明的主旨是提出一种新的FES的精密控制方法，利用步行器的HRV参数预测的膝关节角度与实际的膝关节角度预测的关节角度的误差，通过蚁群算法优化PID的比例系数、积分系数以及微分系数，继而准确稳定实时地控制FES***地电流强度。该项发明可有效地提高FES***准确性和稳定性，并获得可观的社会效益和经济效益。最佳实施方案拟采用专利转让、技术合作或产品开发。

Claims (4)

1.一种功能性电刺激中PID参数的双源特征融合蚁群整定方法，其特征是，包括下列步骤：

首先，利用助行过程的肌肉模型HRV预测膝关节角度；

其次，利用蚁群算法整定PID参数，实时调控FES电流水平强度，其整定流程为：首先根据PID的三个决策变量K_p、K_i和K_d取值范围的上下界，确定包括蚁群群体规模、搜索空间维数的参数，并对其进行编码，然后利用通过实际关节角度与肌肉模型HRV输出关节角度的相应关系作为适度评价函数计算的适应度值；采用蚁群算法对PID参数进行控制，即确定蚁群算法的参数设置，利用蚂蚁随机搜索使其变量优化PID的K_p、K_i和K_d三个系数，利用适应度函数调节蚂蚁每次寻索路径以及判断是否达到预设目标，如达到预设目标，计算最终最佳的位置即得PID的K_p、K_i和K_d三个系数，如未达到预期目标继续寻优，直到达到预设目标；在新的PID系数下计算***输出yout及其与肌肉模型HRV的偏差后再进入下一步蚁群算法的自学习与加权系数自调整；

反复前一步骤，最终实现PID控制参数的自适应在线整定，并用于FES***。

2.根据权利要求1所述的一种功能性电刺激中PID参数的双源特征融合蚁群整定方法，其特征是，实际关节角度与肌肉模型HRV输出关节角度的相应关系为：

L＝M×HRV^-1                           (3)

M表示膝关节角度，HRV表示使用者施加在步行器上力的柄反作用矢量，L表示HRV与M之间的关系，采用偏最小二乘回归的方法求解L：

设有m个HRV变量HRV1，…，HRVm，p个M变量，M1，…，Mp，共i(i＝1，…，n)个观测值的数据集，T、U分别为从HRV变量与M变量中提取的成分，从原始变量集中提取第一对成分T1、U1的线性组合为：

T₁＝ω₁₁HRV₁+…+ω_1mHRV_m＝ω′₁HRV    (4)

U₁＝v₁₁M₁+…+v_1pM_p＝v′₁M             (5)

其中ω₁＝(ω₁₁，…，ω_1m)′为模型效应权重，v₁＝(v₁₁，…，v_1p)′为M变量权重，将上述提取第一成分的要求转化为求条件极值问题：

其中t1、u1为由样本求得的第一对成分的得分向量，HRVO、MO为初始变量，利用拉格朗日乘子法，上述问题转化为求单位向量ω₁和v₁，使θ₁＝ω′₁HRV′₀M₀v₁最大，即求矩阵HRV′₀M₀M′₀HRV₀的特征值和特征向量，其最大特征值为θ₁ ²，相应的单位特征向量就是所求的解ω₁，而v₁由公式得到；

其次建立初始变量对T1的方程

$\left\{\begin{array}{c}{\mathrm{HRV}}_0=t_1{\mathrm{\α}}_1^{\′}+E_1\\ M_0=t_1{\mathrm{\β}}_1^{\′}+F_1\end{array}\right.---\left(7\right)$

其中t1意义同前，α′₁＝(α₁₁，…，α_1m)，β′₁＝(β₁₁，…，β_1p)为仅一个M量t1时的参数向量，E1、F1分别为n×m和n×p残差阵，按照普通最小二乘法可求得系数向量α₁和β₁，其中α₁成为模型效应载荷量；

如提取的第一成分不能达到回归模型的精度，运用残差阵E1、F1代替XO、YO，重复提取成分，依次类推，假设最终提取了r个成分，HRVO、MO对r个成分的回归方程为：

$\left\{\begin{array}{c}{\mathrm{HRV}}_0=t_1{\mathrm{\α}}_1^{\′}+...+t_r{\mathrm{\α}}_r^{\′}+E_r\\ M_0=t_1{\mathrm{\β}}_1^{\′}+...+t_r{{\mathrm{\β}}_r}^{\′}+F_r\end{array}\right.---\left(8\right)$

把第一步分析所得HRV量中提取成分Tk(k＝1，…，r)线性组合带入M量对r个成分建立的回归方程，即t_r＝ω_k1HRV₁+…+ω_kmHRV_m代入M_j＝t₁β_1j+…+t_rβ_rj(j＝1，…，p)，即得标准化变量的回归方程M_j＝α_j1HRV₁+…+α_jmHRV_m；

最后根据式(3)，即可求出L。

3.根据权利要求1所述的一种功能性电刺激中PID参数的双源特征融合蚁群整定方法，其特征是，所述进行编码是，根据膝关节角度和电流模式，以及实际误差情况设定PID的三个参数为5位有效数字，其中K_p小数点前2位，小数点后3位；K_i和K_d小数点前1位，小数点后4位，具体编码如以下公式：

K_i＝y_6，j×10⁰+y_7，j×10^-1+y_8，j×10^-2+y_9，j×10^-3+y_10，j×10^-4

K_p＝y_1，j×10¹+y_2，j×10⁰+y_3，j×10^-1+y_4，j×10^-2+y_5，j×10^-3

K_d＝y_6，j×10⁰+y_7，j×10^-1+y_8，j×10^-2+y_9，j×10^-3+y_10，j×10^-4

并且根据输出的超调等误差设计适应函数具体函数为：

fit＝α1×σ+β1×t+c×error    (9)

其中，σ为超调量，t为上升时间，error为输出关节角度与预设关节角度的相对误差，α1＝0.1，β1＝0.8，c＝2。

4.根据权利要求1所述的一种功能性电刺激中PID参数的双源特征融合蚁群整定方法，其特征是，蚁群算法整定PID控制器参数的具体工作流程为：

Step1：参数初始化，令时间t＝0和循环次数N_max＝0，设置最大循环次数N_cmax，将m个蚂蚁置于起始点；

Step2：设置蚂蚁个数和循环次数；

Step3：蚂蚁随机搜索，在一次爬行结束后，决定哪些特征变量被选中作为实际输入变量，修改禁忌表指针，即选择好之后将蚂蚁移动到新的元素，并把该元素移动该蚂蚁个体的禁忌表中；

Step4：计算蚂蚁个体据状态转移概率公式计算的概率，根据此概率选择元素；

Step5：若蚂蚁元素未遍历完，转Step3，否则为Step6；

Step6：更新信息素浓度划分正确率高的特征变量的信息素浓度得到增强，下次搜索时会以更大的概率被选中；

Step7：满足结束调节，整定结束。

Images