CN113786556A - 足下垂功能性电刺激康复系统变长度迭代学习控制方法 - Google Patents

Abstract

本发明公开了足下垂功能性电刺激康复系统变长度迭代学习控制方法，涉及优化迭代学习控制领域。该方法首先将批次长度变化下重复运行的足下垂功能性电刺激康复系统转换为批次长度相等的时间序列输入输出提升矩阵模型，然后针对批次长度变化问题提出了一种基于多集合逐次投影的优化迭代学习控制设计框架，并基于该框架设计了一种适用于批次长度变化情况的因果前馈加反馈优化迭代学习控制算法。此外，基于多集合逐次投影框架，证明了所设计的优化迭代学习控制算法的收敛性。该方法可以解决批次长度变化情况下足下垂功能性电刺激康复系统的跟踪控制问题，从而实现对期望轨迹的高精度跟踪。

Description

足下垂功能性电刺激康复系统变长度迭代学习控制方法

技术领域

本发明涉及优化迭代学习控制领域，尤其是足下垂功能性电刺激康复系统变长度迭代学习控制方法。

背景技术

足下垂病症往往产生于中风，是一种由于踝关节背屈受限或无能，导致步行时腿部摆动阶段足尖上抬角度不足或无法产生角度，从而影响正常行走的一种现象。足下垂患者行走时会产生非正常步态，例如足部前端拖地，增加了行走过程中跌倒的风险，严重影响了日常生活，甚至危及了患者的生命安全。因此，足下垂康复成为了众多辅助康复研究人员研究的热点问题。

功能性电刺激(Functional electrical stimulation，FES)是一种利用安全的低频脉冲电信号来刺激由于中风等原因而导致功能障碍的目标肌肉，使得目标肌肉再次兴奋，从而再次获得功能恢复的理疗方法。自二十世纪六十年代FES第一次用于辅助卒中患者步态以来，该方法已被证明在运动功能恢复方面十分有效。系统控制策略是FES康复系统中的关键技术，当足下垂患者使用基于FES的康复系统时，为了保证行走时的安全性，通常要求每次电刺激信号至少持续到患者足部与地面的第一次接触。许多FES康复系统使用一个预先定义好的刺激曲线，当检测到患者足部离开地面时，就施加该刺激信号，从而实现足下垂患者腿部摆动阶段踝关节合适的背屈。然而，足下垂患者进行康复行走时，其步幅会由于肌肉疲劳、患者体能消耗或者外界环境干扰等因素而产生相应的变化，预先设定好FES输入曲线的设计思路对于患者康复的效率较低，甚至会影响患者的健康与人身安全。事实上，利用FES辅助足下垂患者恢复是一个不断重复的过程，因而可以利用迭代学习控制(Iterative learning control，ILC)方法来精确控制足下垂患者踝关节背屈角度，从而提高步态控制的精度。

然而，传统ILC方法要求每个批次的时间长度严格一致，而足下垂患者每一步的持续时间并不能保持严格一致。特别是当患者足部与地面第一次接触时，其踝关节背屈角度就不再拥有学习价值，可以看作失去了ILC过程中的输出信号，因而出现了批次长度变化问题。

批次长度变化问题的常用解决思路是设定一个期望的批次长度，在实际应用中，该期望批次长度一般设为可能出现的最大长度。若实际批次长度小于最大长度，则将该批次缺失的误差信息用设计的信息来补齐，例如可以采用0来补充，那么批次长度变化问题就可以用传统的ILC方法来解决。一些研究为了解决输出信息缺失而导致的批次长度变化问题，利用高阶ILC思想，使用输出缺失处以往批次存在的误差来参与设计ILC更新律，从而增强对于批次长度变化的鲁棒性。然而该类ILC方法在处理批次长度变化问题时，由于批次长度的随机性，其收敛速度并不能得到良好保证。

优化ILC方法是一种将优化思想与ILC相结合的高性能控制方法，针对上述批次长度可变问题收敛速率较差问题，优化ILC设计应运而生。

发明内容

本发明人针对上述问题及技术需求，提出了一种足下垂功能性电刺激康复系统变长度迭代学习控制方法，利用优化的思想来解决ILC实现过程中批次长度变化导致的收敛速度较差问题。通过引入多集合逐次投影的思想提出了一种针对批次长度变化问题的优化ILC设计框架，并在此框架下设计了一种适用于批次长度变化情况的因果前馈加反馈ILC算法，提升了系统对于批次长度可变情况的收敛速度。并且，由于引入了反馈机制，该方法能够有效降低非重复扰动对于系统跟踪效果的影响。

本发明的技术方案如下：

一种足下垂康复的功能性电刺激变长度迭代学习优化控制方法，该方法包括：

第一步、确定足下垂FES康复系统的输入输出模型：

由于足下垂患者病症的多样性，一种生物阻抗测量系统被用于实时测量踝关节背屈角度，以确定康复系统的输入输出模型；

足下垂FES康复系统的输出信号为踝关节背屈角度，记为y，由生物阻抗测量系统测量得到；输入信号为FES信号，记为u，具有固定幅值和频率的调制脉冲宽度，由外部可控刺激装置产生，作用于足下垂患者的胫骨前肌；由于踝关节背屈在腿部摆动阶段的运动幅度在20度以内，因而腿部摆动时输入信号与输出信号之间的关系由如下二阶渐近稳定离散时间传递函数来描述：

其中，t表示时间，q^-1定义为向后位移一步算子，即q^-1u(t)＝u(t-1)，q^-a表示有a个采样间隔的死区时间；y(t)、u(t)和d(t)分别表示t时刻的输出信号、输入信号和扰动信号；Φ(·)和Γ(·)是关于q^-1的实值多项式，根据式(1)，令Φ(q^-1)＝1+φ₁q^-1+φ₂q^-2，

其中

φ₁、φ₂、

是数据驱动模型，用二阶离散传递函数来近似；此外，u(t)应限制在幅值约束

内，其中u表示FES信号作用于目标肌肉使得目标肌肉单位开始运动的最小脉冲宽度，

表示患者所能承受的最大刺激强度；

第二步、构建足下垂FES康复系统的离散状态空间方程：

针对第一步确定的足下垂FES康复系统，将其从离散传递函数形式改写为离散状态空间方程，以适用于之后设计的基于模型的ILC方法；忽略死区时间，即令a＝0，则离散时间传递函数(1)改写为：

定义

则经反变换得到θ(t)与输入信号u(t)的关系为：

θ(t)+φ₁θ(t-1)+φ₂θ(t-2)＝u(t) (3)

定义状态向量为x(t)＝[θ(t-2) θ(t-1)]^T，输入输出变量分别为FES信号u(t)和踝关节背屈角度y(t)；用下标k表示批次，令每次摆动期间从足部与地面脱离到再次接触为一个批次，每个批次的实际运行时间长度为N_k，则建立如下足下垂FES康复重复过程的离散状态空间方程：

其中：

对于该重复过程，存在如下两点假设：一是该重复过程各批次的初始状态保持严格一致，即x_k(0)＝x₀，其中x₀是一个常数向量；二是CB≠0；

第三步、建立足下垂功能性电刺激康复系统变批次长度的输入输出提升矩阵模型：

考虑系统沿批次轴的性能，针对式(4)形式的足下垂FES康复重复过程，将其状态空间方程表达式转换为时间序列的输入输出提升矩阵模型，输入输出提升矩阵模型要求输出的时间长度一致；针对足下垂FES康复系统各批次的实际运行时间长度不一致的问题，依据实际情况事先设定好最大时间长度N及最小时间长度N_-，并令实际批次运行时间长度N_k在该范围内变化，即N_k∈{N_-,N_-+1,...,N}；由于是离散情形，一定存在有限个批次长度，定义出现的批次长度的个数为J，满足J＝N-N_-+1；于是，具有相同输出长度的输入输出提升矩阵模型为：

y_k＝Gu_k+d_k (5)

其中：

d_k＝[d_k(1),d_k(2),K,d_k(N)]^T

u_k＝[u_k(0),u_k(1),K,u_k(N-1)]^T

y_k＝[y_k(1),y_k(2),K,y_k(N)]^T

上述变量中，G是时间序列上的输入输出传递矩阵，d_k是外部干扰对输出的影响；输入Hilbert空间l₂[0,N-1]及输出Hilbert空间l₂[1,N]由如下内积及相关的诱导范数定义：

其中，u,v∈l₂[0,N-1]，y,e∈l₂[1,N]，权值R和Q为正常数；定义期望输出y_d∈l₂[1,N]为：

y_d＝[y_d(1) y_d(2) … y_d(N)]^T (8)

尽管该提升矩阵模型的输出长度均为N，但各批次实际输出信号仍然有所缺失，并且若仍采用足部落地后的踝关节背屈角度测量信号值作为输出，则会影响ILC控制效果，甚至使系统不可控；因此，将每次足部与地面接触后的踝关节角度直接设置为期望信号，即

y_k(t)＝y_d(t),t∈[N_k+1,N] (9)

在这种情况下，当N_k＜N时，e_k≠y_d-y_k，其中e_k为跟踪误差序列，于是引入如下随机矩阵F_k来消除该不等关系：

即使在[N_k+1,N]时刻测量到的踝关节背屈角度信号值对于ILC学习过程是无益的，将跟踪误差序列利用如下期望输出序列与实际输出序列来表示：

e_k＝F_k(y_d-y_k)∈l₂[0,N] (11)

第四步、针对批次长度变化问题提出优化ILC设计框架：

批次长度变化下优化ILC轨迹跟踪问题的设计目标是迭代地找到一个最优控制输入u_∞，使得跟踪误差序列收敛到零，等价于迭代地在属于Hilbert空间的多个闭子空间M_j和M_J+1的交集中寻找点(0,u_∞)；闭子空间M_j和M_J+1的定义如下：

M_j＝{(e,u)∈H:e＝F_j(y_d-Gu-d)}∈{M₁,M₂,...,M_J} (12)

M_J+1＝{(e,u)∈H:e＝0} (13)

其中，闭子空间M_j表示批次长度变化的系统动态，M_j包含的闭子空间代表了离散情况下每个批次能够到达的实际运行时间长度，闭子空间M_J+1表示跟踪需求；d∈l₂[1,N]表示一个扰动量，随机矩阵F_j的定义如下：

根据康复系统的能控性要求，假设多个闭子空间M_j和M_J+1在Hilbert空间中存在交集，即

Hilbert空间H的定义如下：

H＝l₂[0,N-1]×l₂[1,N] (15)

其中，×表示两个空间的笛卡尔积；定义的Hilbert空间由跟踪误差序列及输入信号的Hilbert空间组成，其内积和相关的诱导范数由(6)和(7)导出，如下所示：

定义投影算子P_M(z)表示z在集合M上的投影，定义{j_k}_k≥0是一个由批次k确定、在集合{1,2,...,J}中取值的整数序列；通过选择一个初始点z₀∈H，定义Hilbert空间中点序列{z_k}_k≥0满足如下关系：

当使用上述定义的Hilbert空间中的点z来表示式(12)及(13)所代表的闭子空间中的点时，即z＝(e,u)时，设计不同的取值序列{j_k}_k≥0，以不同的顺序对式(12)及(13)所代表的多个闭子空间进行投影，利用逐次投影思想，形成针对批次长度变化问题的优化ILC设计框架；

第五步、设计针对批次长度可变的优化ILC轨迹跟踪算法：

根据第四步提出的优化ILC设计框架，针对批次长度变化问题，只需设计对于多个闭子空间合理的投影顺序，即导出优化ILC轨迹跟踪算法；一种投影顺序的设计如下：

即选定初始点z₀∈H后，先对{M₁,M₂,...,M_J}中的一个闭子空间进行投影，具体投影空间的选定取决于当前实际运行时间长度，然后再对M_J+1进行投影，依次进行；在此过程中，序列{z_k}_k≥0中任意相邻两点的距离均由式(16)的内积进行定义；根据式(18)下Hilbert空间H中距离的变化，设计如下性能指标：

性能指标满足实际工业问题不同的实际要求，此处考虑了跟踪误差序列及控制信号批次间变化；对于具有性能指标(19)的ILC问题，采用如下范数优化ILC律解决：

u_k+1＝u_k+G^*(I+GG^*)^-1e_k (20)

其中，I表示适当维数的单位矩阵，G^*表示输入输出传递矩阵G在Hilbert空间中的伴随矩阵；

由于伴随矩阵的性质，式(20)所示的范数优化ILC律实际上是非因果的，即现实中无法直接实现，因此采用一种因果前馈加反馈框架进行实现；式(20)中存在如下等价形式：

e_k+1＝(I+GG^*)^-1e_k (21)

根据伴随算子的定义，存在如下等式：

<e,Gu>_Q＝e^TQGR^-1Ru＝<R^-1G^TQe,u>_R＝<G^*e,u>_R (22)

其中：

结合式(21)和(22)，控制律(20)改写为：

u_k+1＝u_k+R^-1G^TQe_k+1 (23)

即：

定义共态向量p_k+1(t)为：

则式(24)表示为：

u_k+1(t)＝u_k(t)+R^-1B^Tp_k+1(t) (26)

若确定边界条件p_k+1(N)＝0，则共态向量p_k+1(t)由如下递推关系计算：

p_k+1(t)＝A^Tp_k+1(t+1)+C^TQe_k+1(t+1),t∈[0,N-1] (27)

假设康复系统是可观测的，则式(27)的递推关系存在如下一种实现方式，即共态向量p_k+1(t)还表示为：

p_k+1(t)＝-K(t)[I_n+BR^-1B^TK(t)]^-1A[x_k+1(t)-x_k(t)]+ξ_k+1(t) (28)

其中，I_n表示n维单位矩阵，K(t)为状态反馈矩阵，ξ_k+1(t)为一个前馈项；为了得到状态反馈矩阵K(t)和前馈项ξ_k+1(t)的计算方式，首先结合式(4)、(27)和(28)得到如下等式：

然后将式(28)和(29)代入式(27)，得到：

其中，为了简化表述，使用f₁(·)和f₂(·)代表括号内存在项的一种多项式；若令f₁(·)和f₂(·)均等于0，则式(30)成立且与康复系统的状态无关，采取此步骤，得到状态反馈矩阵K(t)和前馈项ξ_k+1(t)的计算形式如下：

K(t)＝A^TK(t+1)[I_n+BR^-1B^TK(t+1)]^-1A+C^TQC (31)

ξ_k+1(t)＝[I_n+K(t)BR^-1B^T]^-1[A^Tξ_k+1(t+1)+C^TQe_k(t+1)] (32)

并且，将边界条件p_k+1(N)＝0代入式(28)，当t＝0时，若式(28)成立且与康复系统状态无关，则得到状态反馈矩阵和前馈项的边界条件为K(N)＝0和ξ_k+1(N)＝0；

因此给定批次长度变化系统初始输入u₀及趋近于零的常数σ且σ＞0，选定权值R和Q，一种针对批次长度可变的优化ILC轨迹跟踪算法设计如下：

步骤1：根据边界条件K(N)＝0，通过式(31)计算时刻t∈[0,N-1]的状态反馈矩阵K(t)；

步骤2：对于批次k+1，根据边界条件ξ_k+1(N)＝0，通过式(32)计算时刻t∈[0,N-1]的前馈项ξ_k+1(t)；

步骤3：对于批次k+1的时刻t，根据边界条件p_k+1(N)＝0，通过式(28)计算共态向量p_k+1(t)；

步骤4：根据控制律(26)计算时刻t的输入信号u_k+1(t)后，令t←t+1；若t≤N_k-1，则重新回到步骤3，否则进入步骤5；

步骤5：如果N_k+1＜N-1，则将时刻t∈[N_k+1,N-1]的输入信号设置为u_k+1(t)＝u_k(t)；

步骤6：若||e_k+1||/N_k+1＞σ，则重新回到步骤2，否则算法结束；

第六步、分析批次长度变化下的优化ILC轨迹跟踪算法的收敛性：

为了分析第五步设计的优化ILC轨迹跟踪算法的收敛性，首先应分析式(18)所设计的连续投影顺序下点序列{z_k}_k≥0的收敛性；

由于式(12)和(13)定义的闭子空间是相互包含的，即：

则对于常数b≥a≥1，存在如下等式成立：

且对于序列段n＞m≥1，一定存在

其中b∈{1,2,...,J+1}是在投影次序段[m,n]上出现的最大数；由于投影算子具有自伴随性和幂等属性，则下式等式成立：

因而存在：

根据式(34)和(36)，得到：

由于式(35)同样等价于

根据点序列{z_k}_k≥0的定义，则有：

z_k-z_k+1⊥z_k+1 (38)

由式(38)可知，z_k-z_k+1和z_k+1的内积为0，于是：

||z_k||²＝||z_k+1||²+||z_k-z_k+1||² (39)

对式(39)做递归处理，得到：

将式(40)变形后代入式(37)，得到：

||z_n-z_m||²≤||z_m||²-||z_n||² (41)

根据式(39)可知，序列{||z_k||}_k≥0是单调递减的，且下界为0，所以一定存在一个常数α≥0，使得

若给定一个趋近于零的常数ε＞0，一定存在

使得对于任意n≥k，满足0≤||z_n||-α＜ε/2；因此结合式(41)，得到：

||z_n-z_m||²≤||z_m||²-α+α-||z_n||²＜ε/2+ε/2＝ε (42)

由于点序列{z_k}_k≥0是Hilbert空间中的柯西序列，得到点序列{z_k}_k≥0依范数收敛；

继续分析其收敛边界，根据式(18)所设计的连续投影顺序可知，闭子空间M_J+1每隔一次就会作为一次被投影的集合，所以一定存在一个子序列{z_2k}_k≥0使得z_2k∈M_J+1，则对于任意的

都有〈z_2k,z′>＝0；根据内积的连续性及Hilbert空间中的柯西序列的性质，则有：

即存在一个属于闭子空间M_J+1的序列收敛点z_∞与该闭子空间中的每一个点都正交；根据式(33)可知，对任意的j∈{1,2,...,J}都有

则z_∞存在于每一个闭子空间M_j,j∈{1,2,...,J+1}，即

因此得到点序列{z_k}_k≥0依范数收敛到式(12)和(13)定义的闭子空间交集中的一个点；

令点序列{z_k}_k≥0依范数收敛到的点为z_∞＝(0,u_∞)，由于式(12)和(13)所定义的闭子空间均为离散情形下的有限维Hilbert空间，根据点序列{z_k}_k≥0的收敛性可知，依据式(18)所设计的连续投影顺序，两个闭子空间之间的投影距离最终收敛至0，即性能指标(19)收敛至0，即：

得到：

根据式(45)可知，本申请设计的算法能使跟踪误差范数最终收敛至0；

第七步、根据优化ILC轨迹跟踪算法确定足下垂FES康复系统每一迭代批次中每一个时间点的控制输入量，将得到的控制输入量输入批次长度变化的足下垂FES康复系统进行轨迹跟踪控制，使得康复系统受到控制输入量的控制作用追踪对应的期望输出。

本发明的有益技术效果是：

本申请公开了针对足下垂功能性电刺激康复系统此类具有重复运动特征的线性系统，将足下垂功能性电刺激康复系统作为被控对象，针对被控对象出现的批次长度变化情况，提出优化迭代学习控制设计框架，并在该框架下设计可实现的优化迭代学习控制轨迹跟踪算法，进而在保证跟踪需求的同时获得更快的收敛速度，提升了系统对于批次长度可变情况的跟踪性能表现。同时基于提出的多集合逐次投影框架，对设计的优化迭代学习控制轨迹跟踪算法进行收敛性分析，保证了跟踪误差的收敛性。

附图说明

图1是本申请公开的生物阻抗测量系统原理框图。

图2是本申请公开的生物阻抗测量和功能性电刺激的电极位置放置图。

图3是本申请公开的足下垂功能性电刺激康复系统的模型框图。

图4是本申请公开的忽略外部扰动时足下垂功能性电刺激康复系统踝关节背屈角度跟踪曲线图。

图5是本申请公开的忽略外部扰动时踝关节背屈角度沿批次轴跟踪误差向量2-范数图。

图6是本申请公开的非重复扰动下足下垂功能性电刺激康复系统踝关节背屈角度跟踪曲线图。

图7是本申请公开的非重复扰动下踝关节背屈角度沿批次轴跟踪误差向量2-范数图。

具体实施方式

下面结合附图对本发明的具体实施方式做进一步说明。

请参考图1，其示出了本申请使用的生物阻抗测量系统原理框图。该系统采用四个电极，包括两个电流激励电极与两个电压检测电极，四个电极放置位置如图2所示。在该生物阻抗测量系统中，首先利用电压控制型电流源和可编程函数发生器，产生一个频率为50kHz、幅值为0.25mA的正弦交流电，并通过电流激励电极进行输入。然后，利用两个电压检测电极对人体生物阻抗受电流激励而产生的变化进行测量，得到的电压信号经过定制的具有高度共模抑制的仪表放大器进行放大。为了实现激励过程中的实时测量，该仪表放大器要求静默周期后1ms以下的快速工件功能恢复，并且在激励脉冲期间由一个型号为AQV248的光继电器进行输入短路。经放大后的信号再由一个25kHz的高通滤波器进行滤波，滤除由肌电和运动伪影产生的低频干扰信号。最后，经由整流器和低通频率为2kHz的低通滤波器组成的幅度解调电路得到生物阻抗信号的绝对值，再经过一个12-bit的A/D转换器转换，通过光学隔离的串行接口输出至电脑端。整个生物阻抗测量系统由STM32F103RCT6单片机操作，在采样前控制生物阻抗信号的偏移和增益，以便对信号质量进行调谐。该生物阻抗测量系统的测量频率设定为50Hz，即采样时间为0.02s，以便对FES的实时控制。

图2示出了本申请用于生物阻抗测量和功能性电刺激的电极位置放置图。如图所示，在本申请的足下垂FES康复系统中，共需在足下垂患者腿部放置六个电极。其中，生物阻抗测量系统需在髌骨下方和小腿横韧带上的小腿前表面放置两个电流激励电极，在小腿正面胫骨上方源头处和小腿背面跟腱上侧放置两个电压检测电极。另外，用于功能性电刺激的两个电极则放置在胫骨前肌上。

图3示出了本申请公开的足下垂功能性电刺激康复控制系统模型框图。第k批次t时刻的FES输入信号为u_k(t)，将其作用于足下垂患者腿部可以得到第k批次t时刻的踝关节背屈角度输出信号y_k(t)。针对批次长度变化问题，设计输出信号修正器，将该批次足部与地面接触之后的所有输出信号都设置为期望信号。之后，将输出信号y_k′(t)与储存在期望轨迹存储器的设定期望值y_d(t)进行比较，判断N_k时刻之前误差精度是否达到要求。若误差精度没有达到所设定的精度，则将误差e_k(t)与当前控制器输入u_k(t)传递到优化迭代学习控制器生成下一批次的控制器输入u_k+1(t)，如此循环运行直至康复系统实际输出与期望值之间的误差达到精度要求，则停止运行，此时的控制器输入即为最优控制输入。

结合图4-7所示，其示出了一次足下垂功能性电刺激康复系统的数值仿真。根据足下垂患者行走时的实际情况，设定康复系统一个批次的批次时间为1s，采样时间为0.02s，则N＝50。针对式(1)所示的足下垂功能性电刺激康复系统的实际物理模型，通过图1所示的生物阻抗测量系统，得到一组实际参数为：

φ₁＝-0.8097002,

φ₂＝-0.0777289,

根据该组实际参数，得到离散状态空间模型的系统参数为：

C＝[0.6634 0]

在该康复系统实际运行过程中，需要设定一个期望的踝关节背屈角度曲线。定义水平方向为0度，上抬脚背为正方向，则本实施方式设定的期望轨迹为：

y_d(t)＝16sin[(t-1)/20]-14.5

其中，输出信号的单位为度。同时令初始状态满足x_k(0)＝[0 14.5/0.6634]^T，即每一批次均开始于某一个设定的踝关节背屈角度。由于批次开始时间的不确定性，以及足下垂患者体力、安全和环境等因素导致的批次结束时间的不一致，实际批次长度会随着批次的变化而变化。作为一个简单的示例，令实际批次长度N_k满足离散均匀分布，且在40～50之间随机变化。不失一般性，令初始输入u₀(t)＝0,0≤t≤N-1。

本申请采用外部可控刺激装置RehaStim3生成FES刺激信号来刺激胫骨前肌使得踝关节产生背屈角度。该装置通过光学隔离的USB接口与PC机相连，为生物阻抗测量系统在刺激脉冲产生时对仪表放大器进行静置提供触发信号，刺激频率选择50Hz。选取权值Q＝100，R＝80，当权值与康复系统参数确定时，优化迭代学习控制跟踪轨迹算法的参数也随之确定。

首先，忽略外部扰动d_k(t)的影响，足下垂功能性电刺激康复系统的动态模型(1)运行时，请参考图4，其示出了足下垂功能性电刺激康复系统应用本申请提出的优化迭代学习控制算法的轨迹跟踪效果图，经过一定批次后，系统的输出值能在可接受的精度范围内跟踪上期望轨迹。请注意，前几个批次的实际运行长度均不到1s，这符合系统批次长度变化的情景假设。图5表明系统经过一定的迭代批次后其跟踪误差向量的2-范数能够收敛。当设置外部扰动为一个非重复扰动时，即令d_k(t)＝0.2Δ，其中Δ是一个在0到1范围内随机变化的随机数。图6和图7分别示出了在非重复扰动下优化迭代学习控制算法的轨迹跟踪效果图和跟踪误差向量的2-范数沿批次变化图。由于本申请设计的优化迭代学习控制算法引入了反馈机制，因而尽管在非重复扰动下算法的跟踪效果变差，但仍能经过一定批次后在可接受的精度范围内跟踪上期望值。

本申请提出了一种针对批次长度变化问题的优化迭代学习控制跟踪轨迹算法，该算法基于一个多集合逐次投影框架，通过确定对于不同集合的投影顺序来导出实际算法步骤。针对迭代学习控制用于批次长度变化情形时收敛速度较差的现象，本申请提出的算法结合了优化思想，在保证误差收敛的同时提高了收敛速度，并引入反馈机制，降低了批次过程中外部扰动的影响。

以上所述的仅是本申请的优选实施方式，本发明不限于以上实施例。可以理解，本领域技术人员在不脱离本发明的精神和构思的前提下直接导出或联想到的其他改进和变化，均应认为包含在本发明的保护范围之内。

Claims (1)

1.足下垂功能性电刺激康复系统变长度迭代学习控制方法，其特征在于，所述方法包括：

第一步、确定足下垂FES康复系统的输入输出模型：

由于足下垂患者病症的多样性，一种生物阻抗测量系统被用于实时测量踝关节背屈角度，以确定康复系统的输入输出模型；

所述足下垂FES康复系统的输出信号为踝关节背屈角度，记为y，由所述生物阻抗测量系统测量得到；输入信号为FES信号，记为u，具有固定幅值和频率的调制脉冲宽度，由外部可控刺激装置产生，作用于足下垂患者的胫骨前肌；由于踝关节背屈在腿部摆动阶段的运动幅度在20度以内，因而腿部摆动时所述输入信号与所述输出信号之间的关系由如下二阶渐近稳定离散时间传递函数来描述：

其中，t表示时间，q^-1定义为向后位移一步算子，即q^-1u(t)＝u(t-1)，q^-a表示有a个采样间隔的死区时间；y(t)、u(t)和d(t)分别表示t时刻的输出信号、输入信号和扰动信号；Φ(·)和Γ(·)是关于q^-1的实值多项式，根据式(1)，令Φ(q^-1)＝1+φ₁q^-1+φ₂q^-2，

其中

φ₁、φ₂、

是数据驱动模型，用二阶离散传递函数来近似；此外，所述u(t)应限制在幅值约束

内，其中u表示所述FES信号作用于目标肌肉使得目标肌肉单位开始运动的最小脉冲宽度，

表示患者所能承受的最大刺激强度；

第二步、构建所述足下垂FES康复系统的离散状态空间方程：

针对第一步确定的所述足下垂FES康复系统，将其从离散传递函数形式改写为离散状态空间方程，以适用于之后设计的基于模型的ILC方法；忽略死区时间，即令a＝0，则离散时间传递函数(1)改写为：

定义

则经反变换得到θ(t)与输入信号u(t)的关系为：

θ(t)+φ₁θ(t-1)+φ₂θ(t-2)＝u(t) (3)

定义状态向量为x(t)＝[θ(t-2) θ(t-1)]^T，输入输出变量分别为FES信号u(t)和踝关节背屈角度y(t)；用下标k表示批次，令每次摆动期间从足部与地面脱离到再次接触为一个批次，每个批次的实际运行时间长度为N_k，则建立如下足下垂FES康复重复过程的离散状态空间方程：

其中：

对于该重复过程，存在如下两点假设：一是该重复过程各批次的初始状态保持严格一致，即x_k(0)＝x₀，其中x₀是一个常数向量；二是CB≠0；

第三步、建立足下垂功能性电刺激康复系统变批次长度的输入输出提升矩阵模型：

考虑系统沿批次轴的性能，针对式(4)形式的所述足下垂FES康复重复过程，将其状态空间方程表达式转换为时间序列的输入输出提升矩阵模型，所述输入输出提升矩阵模型要求输出的时间长度一致；针对足下垂FES康复系统各批次的实际运行时间长度不一致的问题，依据实际情况事先设定好最大时间长度N及最小时间长度N_-，并令实际批次运行时间长度N_k在该范围内变化，即N_k∈{N_-,N_-+1,...,N}；由于是离散情形，一定存在有限个批次长度，定义出现的批次长度的个数为J，满足J＝N-N_-+1；于是，具有相同输出长度的输入输出提升矩阵模型为：

y_k＝Gu_k+d_k (5)

其中：

d_k＝[d_k(1),d_k(2),K,d_k(N)]^T

u_k＝[u_k(0),u_k(1),K,u_k(N-1)]^T

y_k＝[y_k(1),y_k(2),K,y_k(N)]^T

上述变量中，G是时间序列上的输入输出传递矩阵，d_k是外部干扰对输出的影响；输入Hilbert空间l₂[0,N-1]及输出Hilbert空间l₂[1,N]由如下内积及相关的诱导范数定义：

其中，u,v∈l₂[0,N-1]，y,e∈l₂[1,N]，权值R和Q为正常数；定义期望输出y_d∈l₂[1,N]为：

y_d＝[y_d(1) y_d(2) … y_d(N)]^T (8)

尽管该提升矩阵模型的输出长度均为N，但各批次实际输出信号仍然有所缺失，并且若仍采用足部落地后的踝关节背屈角度测量信号值作为输出，则会影响ILC控制效果，甚至使系统不可控；因此，将每次足部与地面接触后的踝关节角度直接设置为期望信号，即

y_k(t)＝y_d(t),t∈[N_k+1,N] (9)

在这种情况下，当N_k＜N时，e_k≠y_d-y_k，其中e_k为跟踪误差序列，于是引入如下随机矩阵F_k来消除该不等关系：

即使在[N_k+1,N]时刻测量到的踝关节背屈角度信号值对于ILC学习过程是无益的，将所述跟踪误差序列利用如下期望输出序列与实际输出序列来表示：

e_k＝F_k(y_d-y_k)∈l₂[0,N] (11)

第四步、针对批次长度变化问题提出优化ILC设计框架：

批次长度变化下优化ILC轨迹跟踪问题的设计目标是迭代地找到一个最优控制输入u_∞，使得所述跟踪误差序列收敛到零，等价于迭代地在属于Hilbert空间的多个闭子空间M_j和M_J+1的交集中寻找点(0,u_∞)；闭子空间M_j和M_J+1的定义如下：

M_j＝{(e,u)∈H:e＝F_j(y_d-Gu-d)}∈{M₁,M₂,...,M_J} (12)

M_J+1＝{(e,u)∈H:e＝0} (13)

其中，闭子空间M_j表示批次长度变化的系统动态，M_j包含的闭子空间代表了离散情况下每个批次能够到达的实际运行时间长度，闭子空间M_J+1表示跟踪需求；d∈l₂[1,N]表示一个扰动量，随机矩阵F_j的定义如下：

根据康复系统的能控性要求，假设多个所述闭子空间M_j和M_J+1在Hilbert空间中存在交集，即

Hilbert空间H的定义如下：

H＝l₂[0,N-1]×l₂[1,N] (15)

其中，×表示两个空间的笛卡尔积；定义的Hilbert空间由跟踪误差序列及输入信号的Hilbert空间组成，其内积和相关的诱导范数由(6)和(7)导出，如下所示：

定义投影算子P_M(z)表示z在集合M上的投影，定义{j_k}_k≥0是一个由批次k确定、在集合{1,2,...,J}中取值的整数序列；通过选择一个初始点z₀∈H，定义Hilbert空间中点序列{z_k}_k≥0满足如下关系：

当使用上述定义的Hilbert空间中的点z来表示式(12)及(13)所代表的闭子空间中的点时，即z＝(e,u)时，设计不同的取值序列{j_k}_k≥0，以不同的顺序对式(12)及(13)所代表的多个闭子空间进行投影，利用逐次投影思想，形成针对批次长度变化问题的优化ILC设计框架；

第五步、设计针对批次长度可变的优化ILC轨迹跟踪算法：

根据第四步提出的所述优化ILC设计框架，针对批次长度变化问题，只需设计对于多个闭子空间合理的投影顺序，即导出优化ILC轨迹跟踪算法；一种投影顺序的设计如下：

即选定所述初始点z₀∈H后，先对{M₁,M₂,...,M_J}中的一个闭子空间进行投影，具体投影空间的选定取决于当前实际运行时间长度，然后再对M_J+1进行投影，依次进行；在此过程中，序列{z_k}_k≥0中任意相邻两点的距离均由式(16)的内积进行定义；根据式(18)下Hilbert空间H中距离的变化，设计如下性能指标：

所述性能指标满足实际工业问题不同的实际要求，此处考虑了跟踪误差序列及控制信号批次间变化；对于具有性能指标(19)的ILC问题，采用如下范数优化ILC律解决：

u_k+1＝u_k+G^*(I+GG^*)^-1e_k (20)

其中，I表示适当维数的单位矩阵，G^*表示输入输出传递矩阵G在Hilbert空间中的伴随矩阵；

由于伴随矩阵的性质，式(20)所示的范数优化ILC律实际上是非因果的，即现实中无法直接实现，因此采用一种因果前馈加反馈框架进行实现；式(20)中存在如下等价形式：

e_k+1＝(I+GG^*)^-1e_k (21)

根据伴随算子的定义，存在如下等式：

<e,Gu>_Q＝e^TQGR^-1Ru＝<R^-1G^TQe,u>_R＝<G^*e,u>_R (22)

其中：

结合式(21)和(22)，控制律(20)改写为：

u_k+1＝u_k+R^-1G^TQe_k+1 (23)

即：

定义共态向量p_k+1(t)为：

则式(24)表示为：

u_k+1(t)＝u_k(t)+R^-1B^Tp_k+1(t) (26)

若确定边界条件p_k+1(N)＝0，则所述共态向量p_k+1(t)由如下递推关系计算：

p_k+1(t)＝A^Tp_k+1(t+1)+C^TQe_k+1(t+1),t∈[0,N-1] (27)

假设康复系统是可观测的，则式(27)的递推关系存在如下一种实现方式，即所述共态向量p_k+1(t)还表示为：

p_k+1(t)＝-K(t)[I_n+BR^-1B^TK(t)]^-1A[x_k+1(t)-x_k(t)]+ξ_k+1(t) (28)

其中，I_n表示n维单位矩阵，K(t)为状态反馈矩阵，ξ_k+1(t)为一个前馈项；为了得到状态反馈矩阵K(t)和前馈项ξ_k+1(t)的计算方式，首先结合式(4)、(27)和(28)得到如下等式：

然后将式(28)和(29)代入式(27)，得到：

其中，为了简化表述，使用f₁(·)和f₂(·)代表括号内存在项的一种多项式；若令f₁(·)和f₂(·)均等于0，则式(30)成立且与康复系统的状态无关，采取此步骤，得到所述状态反馈矩阵K(t)和前馈项ξ_k+1(t)的计算形式如下：

K(t)＝A^TK(t+1)[I_n+BR^-1B^TK(t+1)]^-1A+C^TQC (31)

ξ_k+1(t)＝[I_n+K(t)BR^-1B^T]^-1[A^Tξ_k+1(t+1)+C^TQe_k(t+1)] (32)

并且，将所述边界条件p_k+1(N)＝0代入式(28)，当t＝0时，若式(28)成立且与康复系统状态无关，则得到状态反馈矩阵和前馈项的边界条件为K(N)＝0和ξ_k+1(N)＝0；

因此给定批次长度变化系统初始输入u₀及趋近于零的常数σ且σ＞0，选定权值R和Q，一种针对批次长度可变的优化ILC轨迹跟踪算法设计如下：

步骤1：根据边界条件K(N)＝0，通过式(31)计算时刻t∈[0,N-1]的状态反馈矩阵K(t)；

步骤2：对于批次k+1，根据边界条件ξ_k+1(N)＝0，通过式(32)计算时刻t∈[0,N-1]的前馈项ξ_k+1(t)；

步骤3：对于批次k+1的时刻t，根据边界条件p_k+1(N)＝0，通过式(28)计算所述共态向量p_k+1(t)；

步骤4：根据控制律(26)计算时刻t的输入信号u_k+1(t)后，令t←t+1；若t≤N_k-1，则重新回到步骤3，否则进入步骤5；

步骤5：如果N_k+1＜N-1，则将时刻t∈[N_k+1,N-1]的输入信号设置为u_k+1(t)＝u_k(t)；

步骤6：若||e_k+1||/N_k+1＞σ，则重新回到步骤2，否则算法结束；

第六步、分析批次长度变化下的所述优化ILC轨迹跟踪算法的收敛性：

为了分析第五步设计的所述优化ILC轨迹跟踪算法的收敛性，首先应分析式(18)所设计的连续投影顺序下点序列{z_k}_k≥0的收敛性；

由于式(12)和(13)定义的闭子空间是相互包含的，即：

则对于常数b≥a≥1，存在如下等式成立：

且对于序列段n＞m≥1，一定存在

其中b∈{1,2,...,J+1}是在投影次序段[m,n]上出现的最大数；由于投影算子具有自伴随性和幂等属性，则下式等式成立：

因而存在：

根据式(34)和(36)，得到：

由于式(35)同样等价于

根据点序列{z_k}_k≥0的定义，则有：

由式(38)可知，z_k-z_k+1和z_k+1的内积为0，于是：

||z_k||²＝||z_k+1||²+||z_k-z_k+1||² (39)

对式(39)做递归处理，得到：

将式(40)变形后代入式(37)，得到：

||z_n-z_m||²≤||z_m||²-||z_n||² (41)

根据式(39)可知，序列{||z_k||}_k≥0是单调递减的，且下界为0，所以一定存在一个常数α≥0，使得

若给定一个趋近于零的常数ε＞0，一定存在

使得对于任意n≥k，满足0≤||z_n||-α＜ε/2；因此结合式(41)，得到：

||z_n-z_m||²≤||z_m||²-α+α-||z_n||²＜ε/2+ε/2＝ε (42)

由于点序列{z_k}_k≥0是Hilbert空间中的柯西序列，得到点序列{z_k}_k≥0依范数收敛；

继续分析其收敛边界，根据式(18)所设计的连续投影顺序可知，闭子空间M_J+1每隔一次就会作为一次被投影的集合，所以一定存在一个子序列{z_2k}_k≥0使得z_2k∈M_J+1，则对于任意的

都有z_2k,z′＝0；根据内积的连续性及Hilbert空间中的柯西序列的性质，则有：

即存在一个属于闭子空间M_J+1的序列收敛点z_∞与该闭子空间中的每一个点都正交；根据式(33)可知，对任意的j∈{1,2,...,J}都有

则z_∞存在于每一个闭子空间M_j,j∈{1,2,...,J+1}，即

因此得到点序列{z_k}_k≥0依范数收敛到式(12)和(13)定义的闭子空间交集中的一个点；

令点序列{z_k}_k≥0依范数收敛到的点为z_∞＝(0,u_∞)，由于式(12)和(13)所定义的闭子空间均为离散情形下的有限维Hilbert空间，根据点序列{z_k}_k≥0的收敛性可知，依据式(18)所设计的连续投影顺序，两个闭子空间之间的投影距离最终收敛至0，即所述性能指标(19)收敛至0，即：

得到：

根据式(45)可知，本申请设计的算法能使跟踪误差范数最终收敛至0；

第七步、根据所述优化ILC轨迹跟踪算法确定所述足下垂FES康复系统每一迭代批次中每一个时间点的控制输入量，将得到的控制输入量输入批次长度变化的足下垂FES康复系统进行轨迹跟踪控制，使得康复系统受到控制输入量的控制作用追踪对应的期望输出。

Images