我国康复市场需求庞大。人口老龄化加速、慢病患者逐年增加、意外事故多发等导致运动功能障碍人数逐渐增加，急需大量康复器械[1]。作为康复医疗器械的一个细分领域，电刺激治疗仪通过采集、分析表面肌电信号和生物反馈训练，对人体施加电刺激以改善神经或肌肉功能，主要用于盆底康复[2]、神经康复[3-4](脑卒中、帕金森等)两个领域。20世纪60年代，Liberson开创了功能性电刺激(Functional Electrical Stimulation, FES)疗法用于矫正偏瘫患者足下垂的步态，能够有效恢复或重建截瘫患者部分运动功能[5]。电刺激治疗仪需要在专业人员监督下操作，电刺激信号大小应根据患者病情和治疗效果确定。根据患者病情及治疗效果实时精准控制FES的脉冲电流强度和触发时序来保证电刺激作用效果，进而准确完成预定功能是FES的技术关键。由于PID控制稳定性好、工作可靠，目前在功能性电刺激领域得到了广泛应用[6-7]。但传统PID参数整定过程相对复杂，需要丰富的经验和技巧，导致参数设定不够精确，影响电刺激效果；PID对于非线性、时变系统的适应性较差，而人体本身就是一个高度复杂的非线性系统，因此，直接使用PID进行功能性电刺激控制无法完全满足实际需要。此外，PID在处理快速响应系统时也存在一定局限性，而功能性电刺激往往需要快速、准确地响应患者的生理状态变化。

为解决传统PID控制算法在电刺激治疗仪应用中存在的局限性，提高控制效率，国内外学者提出了多种优化方法来确定PID控制参数。例如，利用遗传算法对PID控制参数进行整定，但需要根据不同控制对象，预先选定合理的控制参数[8];提出一种针对可穿戴设备电刺激系统的自适应PID算法，在大体方向上可以直接调节控制参数，不足之处在于得不到准确的控制参数[9]。长城建造算法(GWCA)是一种简单、鲁棒且具有竞争力的有效求解约束和无约束问题的方法。GWCA的灵感来自古长城建造过程中工人之间的竞争和淘汰机制，引入了一个劳工运动数学模型来模拟算法动态。与其他使用多个模型生成新解决方案的方法不同，GWCA在每次迭代中随机为每个工人分配一个预定义的运动模型。这种独特的方法展示了GWCA的动态性和高收敛性，从而在效率方面优于其他优化方法[10]。在GWCA中，距离因子H(t)表示第t次迭代中第i个工作者将建筑材料从当前位置通过不同工具运输到山顶，随着迭代次数增加，距离因子呈线性变化。由于功能性电刺激信号的自适应控制是非线性的，采用线性变化的距离因子对 PID控制参数寻优存在一定缺陷，导致GWCA各阶段的寻优不平衡，易陷入局部最优，使得PID控制参数不精确，影响PID控制算法对功能性电刺激信号的自适应控制。

为解决GWCA在求解时存在的问题，本文提出一种长城建造优化算法的改进算法(IGWCA)确定增量式PID算法的3个控制参数，旨在优化电刺激信号自适应控制响应速度、控制精度和稳定性，进而提升电刺激治疗仪的康复治疗效果。

1、基于IGWCA自适应增量式PID设计

1.1 增量式PID原理

增量式PID由比例单元、积分单元和微分单元3部分组成，通过比较实际输出值与期望输出值之间的差异来修正控制信号[11]。增量式PID算法数学模型：

式中：u(t)为PID的输出信号；e(t)为给定值与测量值之间的关系；Kp、Ki、Kd分别对应电刺激治疗仪功能性电刺激PID算法的比例参数、积分参数和微分参数。

将式(1)离散化后得到

式中：u(k)为k时刻电刺激信号输出；k为采样序列(k=0,1,2,…);e(k)为k时刻电刺激信号值与电刺激信号目标值的误差。

将式(2)代入到式(3)

Δu(k)=u(k)-u(k-1) , (3)

可进一步简化得到增量式PID离散化公式：

u(k)=Kp[e(k)-e(k-1)]+Kie(k)+Kd[e(k)-2e(k-1)+e(k-2)]。 (4)

由式(4)可知，当e(k)、e(k-1)、e(k-2)、Kp、Ki和Kd已知时，u(k)的计算和存储量相比传统PID算法大幅减少。

图1 电刺激自适应控制原理

图2 电刺激自适应控制流程

1.2 IGWCA自适应增量式PID控制原理

为提高GWCA的寻优精度和跳出全局次优解的能力，通过改进距离因子H(t)建立IGWCA的寻优数学模型，在电刺激治疗仪PID算法控制参数寻优过程中，满足功能性电刺激信号的非线性变化，更好地自适应寻优过程，提高控制参数寻优精度。为避免GWCA陷入局部最优或者使其在短时间内跳出局部最优，引入“次优解聚集”策略改进GWCA的寻优策略，从而进一步提高电刺激治疗仪PID算法控制参数的精度。

如图1所示，基于 IGWCA优化的电刺激自适应控制系统主要包括电刺激信号目标值输入模块、电刺激信号实时输入模块、电刺激信号误差计算模块、改进长城建造算法模块、增量式PID模块和电刺激信号输出模块。工作原理：1)电刺激治疗仪将目标电刺激信号值输入电刺激信号目标值输入模块。2)利用电刺激电传感器实时采集皮肤电刺激信号，利用电刺激治疗仪控制器的模数转换电路将模拟信号转换为数字信号，得到皮肤电刺激信号值。3)在电刺激信号误差计算模块中计算当前皮肤电刺激信号值与电刺激信号目标值的误差e(k),并将误差e(k)输入增量式PID算法模块计算得到当前电刺激信号输出值u(k)。4)经电刺激信号输出模块反复执行上述步骤，快速达到电刺激信号目标值。其中，增量式PID算法数学模型为通过目标函数引导、利用 IGWCA优化后的增量式PID。利用 IGWCA优化电刺激治疗仪控制系统的PID算法控制功能性电刺激信号，可以提高功能性电刺激的自适应性响应和稳定性，从而改善电刺激治疗仪的康复治疗效果，具体实施步骤见图2。

2、IGWCA算法数学模型

GWCA算法搜索过程首先从创建一个随机的工作者群体(即一组Kp、Ki和Kd参数候选解)开始。工作者群体包括工程师、士兵和劳动者3种角色，不同工人通过不同方式到达山顶建造长城。在迭代过程中，模拟工程师、士兵和劳动者建造长城的不同方式及所在阶层位置建立数学模型。工人根据最佳工程师位置、最佳士兵位置和最佳劳动者位置，通过数学模型计算领导者的位置，即为PID最佳控制参数。在GWCA中，工程师、士兵、劳动者3种工作模式之间的随机切换控制着GWCA的全局搜索和局部开发。

2.1 初始化阶段

在GWCA中，将工人的位置视为候选方案(Lacn),其中包含一些表示为决策变量的元素(Po

式中：d为每个解决方案中决策变量的数量；n为被认为是解决方案的候选工人数量。

在优化初始阶段，Po

由逻辑混沌映射生成混沌序列，可根据优化问题定义的内容随机确定决策变量的值。Po

在搜索空间中随机确定的初始位置：

为工人所在的初始位置；cxlj为逻辑混沌映射生成的混沌序列；LBj和UBj分别为第j个决策变量工作者个体位置的下界和上界。

功能性电刺激信号的自适应控制问题是通过PID算法实现调整实时电刺激信号的灵敏度和稳定性达到最佳鲁棒性，将这一过程转换成待优化数学模型，并将该模型作为IGWCA寻优目标函数，表达式为

式中：J(k)为k时刻的目标函数值。

2.2 改进距离因子模型

为适应功能性电刺激信号的非线性变化，引入布朗运动，增加距离因子H(t)的随机性，改进后的距离因子H

式中：Hmax为改进的距离因子最大值，取1;n为当前阶段对应的工作者数量；f(X

)为第t次迭代第i个工作者位置的适应度；f(X

)为第t-1次迭代最佳工作者位置的适应度，即第t-1次迭代适应度最小值；布朗函数β=randn(1,D),randn是从标准正态分布中随机生成的维数为1×D数组中的一个随机值；X

是第t次迭代第i个工作者位置；X

是第t-1次迭代最佳工作者位置；Z

之间的欧式距离，表达式为

式中：D为GWCA的问题维度，取3;y=1,2,3表示第y个问题维度。

2.3 次优解聚集策略模型

为避免GWCA陷入局部最优或在陷入局部最优后能够在短时间内跳出局部最优，提高功能性电刺激PID算法的控制精度，引入“次优解聚集”策略改进GWCA的寻优策略。次优解聚集策略通过计算工作者个体在局部最优值附近的密度来实现，并提出全局次优解聚集指标C

式中：N为工作者最大规模；pbest为全局最优解；UB和LB分别为工作者个体位置上界和下界；MaxIter为最大迭代次数；

为控制聚集范围。

利用式(10)计算每个工作者个体全局次优解聚集指标，并判断是否超过N/2个工作者个体的全局次优解聚集指标大于0.5。若否，则结束次优解聚集策略；若是，则说明GWCA在对功能性电刺激PID算法控制参数寻优时陷入了局部最优，利用式(6)跳出局部最优。

为第t+1次迭代第i个工作者位置；X

为第t次迭代最佳工作者位置；w为在-1到1之间生成的随机数；c为由标准柯西分布产生的随机数。

2.4 工程师运动数学模型

在IGWCA局部开发阶段，引入改进距离因子H

后的工程师速度因子v

式中：T为模拟工程师工具产生的推力；TL为工程师工具磨损水平，取值从1线性下降到0;m为岩石质量；g为重力加速度；θ为[0,80]的一个随机数，表示第i个工作者位置的斜率与地平线之间的夹角；gampdf为伽马分布的概率密度函数，P和Q为伽玛概率分布函数中的两个参数；t为当前迭代次数；C(t)为第t次迭代中山路地形的复杂程度，数学模型为

式中：Cmax为C(t)的最大值；Cmin为C(t)的最小值。

伽马分布概率密度函数表达式为[12]

工程师熟练使用工具来提高搬运石头的效率，工程师与效率最高的工作者(最佳位置工作者)位置竞争，利用改进的速度因子建立工程师位置更新数学模型：

为第t+1次迭代第i个工程师的位置；E

为第t次迭代最佳工程师的位置；E

为第t次迭代第i个工程师的位置；X

为第t次迭代最佳工作者的位置；R1、R2取值为0到1的随机数，表示工程师的随机权重；γ取0或1。

2.5 士兵运动数学模型

在IGWCA局部开发阶段，引入改进距离因子H

后的士兵速度因子v

模拟士兵前往最近工作者那里监督工作的行为，使用符号函数确定士兵的下一个位置，同时引入改进的士兵速度因子建立士兵的位置更新数学模型：

为第t次迭代最佳士兵的位置；Poide为离第t个工作者最近士兵的位置；f(Side)为Side的适应度；f(S

)为第t次迭代第i个士兵的适应度；R3、R4为取值0到1的随机数；sign为符号函数，当sign(f(Side)-f(S

))=1时，士兵接近领导者位置，当sign(f(Side)-f(S

2.6 劳动者运动数学模型

劳动者受到士兵和工程师约束，根据士兵和工程师的位置进行位置更新。劳动者位置更新数学模型：

为第t+1次迭代第i个劳动者的位置；L

为第t次迭代第i个劳动者的位置；L

为第t次迭代最佳劳动者的位置；R5取值为0到1的随机数。

2.7 替换和重组的数学模型

在长城建造期间，当工人失去了工作能力或效率低下时，代理人会招募新的工人来替换他们，以确保工程进度。将利用式(10)、式(12)和式(13)获得的工作者位置代入评价函数，淘汰具有较高适应度的工作者。代理替换工作者的数学模型：

式中：ceil是一个整数函数，返回大于或等于指定表达式的最小整数；R6为取值0到1的随机数；EL为被淘汰工作者的数量；PN为新招募工作者的位置；PR是在淘汰后仍保留工作者的位置；PE是淘汰工作者的位置；p为定期淘汰工作者的比例，取10%。

2.8 工作者位置更新流程

利用改进的距离因子和“次优解聚集”策略优化增量式PID算法，得到IGWCA工作者位置更新模型。实施步骤见图3。IGWCA工作者更新数学模型优化增量式PID算法Kp、Ki和Kd控制参数的步骤见图4。

图3 工作者位置更新流程

图4 IGWCA优化增量式PID控制参数

3、功能性电刺激仿真分析

目前，临床上使用最多的电刺激治疗仪功能性电刺激形式为恒流刺激，常见波形为低频双相脉冲波形。在实际电刺激治疗中，频率、脉宽、电流大小等都会影响功能性电刺激的刺激效果。通过控制变量法[13],在脉宽、频率等不变的情况下，以电流为电刺激信号输出目标，通过3种算法仿真分析验证IGWCA在功能性电刺激信号自适应控制方面的优越性。

3.1 算法寻优仿真分析

在Matlab中设计GWCA和IGWCA,设初始工作者最大规模N=90,最大迭代次数MaxIter=50;伽马分布两个函数参数P=2,Q=1.3,工作者个体位置上界UB=[300,200,300],个体位置下界LB=[0,0,0]。将电流信号问题转化为待优化数学模型，寻优过程中GWCA和IGWCA对电刺激治疗仪增量式PID算法控制参数的寻优过程见图5。

图5 控制参数寻优过程

图6 两种算法优化目标函数适应度对比

由图5可知：GWCA对增量式PID算法控制参数寻优最佳值为Kp=10.21、Ki=0.37、Kd=1.20;IGWCA对增量式PID算法控制参数寻优最佳值为Kp=52.36、Ki=0.51、Kd=5.15。

GWCA和IGWCA两种算法对目标函数的适应度见图6。由图6可知，随着迭代次数增加，两种算法计算目标函数适应度皆减少。不同之处在于GWCA的适应度在整个迭代过程中下降较慢，大约迭代15次后趋于稳定，维持在0.5左右；IGWCA在迭代10次时达到较低的适应度。由此可知，IGWCA的适应度下降速度更快，且适应度更低，说明IGWCA寻优速度和收敛精度更高，优化功能性电刺激信号的自适应控制性能更好。

3.2 Simulink控制器模型搭建及仿真分析

在Matlab/Simulink[14]中搭建电刺激治疗仪功能性电刺激自适应控制仿真模型(图7),仿真分析标准PID、GWCA-PID、IGWCA-PID 3种控制算法对电刺激治疗仪电流信号自适应系统响应的影响，结果见图8。

图7 电刺激自适应控制仿真模型

图8 不同控制算法作用下电刺激控制效果比较

由图8可知：1)在初始阶段，3种控制算法均迅速提升了系统输出，以尝试达到设定的电流目标值。2)在初始阶段，IGWCA-PID表现出最快的响应速度，并在大约1.2秒时达到峰值；随后，系统输出开始稳定，IGWCA-PID显示出较高稳定性和较小超调。3)标准PID和GWCA-PID在达到稳态前经历了较明显波动，在4 s后，所有控制算法均基本达到稳态，但IGWCA-PID的输出响应更接近目标值1 mA,显示出更优的稳定性和准确性。4)与其他两种控制算法相比，IGWCA-PID在整个仿真过程中显示出更小的波动性和更少的调节时间。

综上可知，IGWCA-PID与传统PID及GWCA-PID相比响应速度更快、稳定性高、超调量小，输出响应与目标值之间的误差最小，显示出较强鲁棒性。

4、结论

本文通过改进距离因子和引入“次优解聚集”策略改进GWCA的寻优数学模型和寻优策略，提高了GWCA全局搜索和局部开发的寻优平衡性。“次优解聚集”策略对算法寻优进行监控，能够在电刺激治疗仪功能性电刺激PID算法控制参数寻优时及时跳出局部次优解局限。在Simulink中搭建自适应控制模型并进行仿真分析，结果表明：与标准PID和GWCA两种控制算法相比，IGWCA响应速度更快，在1.2秒时达到设定电流峰值且具有较高的稳定性和较小超调。同等条件下，IGWCA还具有更高的寻优速度和收敛精度，可以更快地达到目标电流信号值，从而保证了更为稳定的康复训练效果。

本文通过改进算法仿真研究，验证了IGWCA可有效提高电刺激治疗仪功能性电刺激PID控制参数的控制精度及电刺激控制的灵敏度，下一步将开展临床数据分析和有效性分析。

参考文献:

[1]王金翔.多通道功能性电刺激康复仪的系统设计[D].济南:济南大学,2022.

[2]肖丽丽.生物反馈电刺激治疗仪在产妇产后康复中的应用效果[J].中国医疗器械信息,2023,29(10):145-147.

[5]蓝宁,肖志雄,聂开宝,等.功能性电刺激的原理､设计与应用(一)[J].中国康复理论与实践,1997,3(4):151-154.

[8]郭新奇,李晴,郝林琳.遗传算法整定PID控制器在下肢电刺激康复训练中的应用[J].中国医学装备,2011,8(10):70-74.

[9]刘银华,付宝乐.针对可穿戴设备电刺激系统的自适应PID算法:202310822455.0 [P].2023-10-31.

[11]宋安伟,王良,刘艳梅,等.基于增量式PID算法的移动机器人定位系统设计[J].湖南文理学院学报(自然科学版),2024,36(2):19-26.

[13]滕飞,边远,平冰宇,等.基于控制变量法的车床主轴转速分析[J].通化师范学院学报,2022,43(6):49-54.

[14]杨雪春.基于Simulink的柴油机可变涡流控制系统仿真建模研究[J].通化师范学院学报,2021,42(6):10-17.

基金资助:吉林省中医药管理局中医药科技项目(2024154);吉林省教育厅科学技术研究项目(JJKH20240080KJ);北华大学博士科研启动专项基金;

文章来源:祝燕,苗蓓蓓,刘洪波.基于改进长城建造算法的电刺激治疗仪自适应增量式PID控制[J].北华大学学报(自然科学版),2025,26(01):133-140.