1 引言（Introduction）

至2030年，我国将有超过3 000万脑卒中患者，各类残疾患者的总数将超过8 000万人^[1-2]。下肢康复外骨骼^[3-6]作为服务于下肢伤残及运动障碍患者的辅助康复设备，能够完成更高精度的重复性康复训练，相比于固定康复设备能显著提高患者的能动性。

近十年来国内外多个研究团队进行了下肢康复外骨骼的研究。Kwa等^[7]设计了IHMC-MAE（IHMC mobility assist exoskeleton）系统，该系统具有多个主被动自由度以适应下肢关节运动，以健康人群步态数据为参考轨迹，牵引穿戴者复现步行运动。Sankai ^[8]设计了HAL（hybrid assistive limb）外骨骼康复系统，利用肌电信号对关节力矩进行估计，并结合阻抗模型计算人机交互关节的输出力矩，从而帮助偏瘫患者完成康复运动训练。以色列ReWalk公司研制了ReWalk系列外骨骼^[9]，该外骨骼能够检测穿戴者的运动重心变化，进而模拟自然步态运动并防止在运动过程中发生意外跌倒。中国科学院合肥智能机械研究所、哈尔滨工业大学、中国科学院深圳先进技术研究院、中国科学院沈阳自动化研究所^[10]等科研院所也开展了相关研究，并设计了一系列实验样机。尽管围绕下肢康复外骨骼的研究很多，但其应用的可行性仍不显著，原因在于刚性外骨骼在结构设计和人机交互方面均存在技术难点，具体表现在以下两方面：

(1) 在结构设计方面，刚性外骨骼以人体下肢骨骼为构形，当刚性连杆与关节位置未对齐时将产生非期望交互力矩，从而严重影响系统控制性能。

(2) 在人机交互方面，外骨骼的刚性远大于人体下肢，系统柔顺性决定了穿戴者的舒适性和康复治疗效率。除硬件系统的柔顺设计外，也要求控制算法能实现步态的准确识别和参考轨迹的准确追踪。同时应综合考虑外骨骼驱动力矩的约束或输出约束，确保交互过程的安全性与舒适性。

模块化、柔顺性的外骨骼关节可有效解决上述问题，原因在于：1) 模块化关节能够优化连杆—关节配置方案，改善结构设计的灵活性，提高系统功率密度，降低总重量与末端惯量；2) 对刚性系统引入柔顺性结构从而降低系统机械阻抗，缓冲瞬时过载冲击，提高交互柔顺性与安全性。

串联弹性驱动器（SEA）是柔顺性模块化关节的重要设计内容。该概念最早由Pratt等^[11]在1995年提出，即在刚性传动机构中引入弹性环节。相比于传统刚性驱动器，SEA能够改善系统柔顺性，具有低阻抗、抗冲击、力控稳定等特点，但对控制算法设计提出了更高的要求。Pratt等率先提出应用于SEA的比例—积分—微分（PID）力矩控制算法，Tagliamonte等^[12]在此基础上引入了级联PD以优化系统柔顺性。为提升控制算法的鲁棒性，Kong等^[13]提出了基于线性扰动观测器的控制方法，Calanca等^[14]设计了基于参考模型的自适应控制方法。前述研究对象均为SEA线性系统，而实际系统存在非线性环节。为此，朱秋国等^[15]利用反向传播（BP）神经网络补偿负载重力，王萌等^[16]先后利用有限时间输出反馈控制、鲁棒误差积分最优控制、有限冲击响应滤波等方法对非线性SEA进行力矩控制和抖动抑制。

在康复机器人领域中，不仅需要对SEA模块进行控制，还需要考虑系统集成问题和人机交互设计^[17]。在系统集成问题上，SEA应兼顾临床应用需求与结构的小型化、轻量化、低惯量；在交互过程中应确保安全性和舒适性^[18]。Gerlach-Hahn等^[19]设计了集成SEA的膝关节康复外骨骼，并应用了级联PD反馈控制，但其SEA体积较大，需要额外的弹性元件固定器件。Kim等^[20]针对下肢外骨骼的非线性SEA系统采用了一种基于逆向模型的控制方法，但未处理输入饱和导致的控制抖动。

综上所述，上述研究中设计的SEA弹性元件难以实现小型化与轻量化，在控制算法上缺少控制量的约束设计，仅以输入饱和限制控制量可能引发系统抖动甚至系统发散。因此，本文设计了一种应用于下肢康复外骨骼的串联弹性关节，并提出了迭代线性化模型预测控制方法，主要创新工作在于：

(1) 设计了面向康复外骨骼的串联弹性关节，该关节具有结构紧凑、质量较轻、承载能力大的特点，在优化系统柔顺性的同时，能够根据外部载荷实现线性/非线性刚度的切换。

(2) 针对串联弹性关节研究了带约束非线性优化问题，提出基于模型预测控制的迭代线性化方法，结合Laguerre网络优化求解非凸优化控制问题。通过对系统控制量的约束，避免了控制量出现尖峰和抖动，改善了人机交互过程的舒适性。

2 串联弹性关节设计（Design of the series elastic joint）

康复外骨骼需要与穿戴者进行深度的人机交互，因此对关节轻量化与紧凑化要求较高。现有的康复外骨骼倾向于在关节处集中布置驱动元件，导致关节轴向尺寸较大而连杆空间未能有效利用。基于此设计了一种具有多分支弹性单元的串联弹性关节。

2.1 多分支弹性单元结构设计

通常情况下，标准的SEA需要包含弹性元件以串联驱动电机实现运动传递，因此弹性元件是SEA设计的主要对象。本文基于多分支平面涡卷扭簧提出了一种新型弹性单元，该弹性单元由多个不接触的涡卷弹簧并联组成，如图 1所示。

基于并联的结构配置，多分支涡卷扭簧能够实现结构的紧凑化与轻量化，从而能在关节有限的空间内作为SEA的弹性元件。扭簧在预紧装配后呈现阿基米德螺旋线（涡状线）形状，两端固定于SEA金属结构的内外轮。加载外力后，内外轮发生相对转动并将载荷传递至涡卷扭簧，导致涡卷扭簧发生弹性形变：

$ \rho =\alpha +\beta \cdot \theta $

(1)

其中，$\theta $为涡状线旋角，$\rho $为涡状线半径，$\alpha $、$\beta $为涡状线设计参数。

扭簧的弯曲变形主要由扭簧横截面所受弯矩引起，弹性体受力分析如图 2所示。其中，$B$点为弹簧形变部分外端点，$\phi $为扭簧总旋角，$F_{\rm t}$、$F_{\rm f}$分别为作用在$B$点的等效切向力与等效法向力，另取$\mu $为$B$点位置约束的修正系数。

任取弹簧形变部分内任一点$A$进行弹性体受力分析，则扭簧外轮在外力矩$M$作用下，$A$点受到的等效切向力矩$M_{\rm t}$和等效法向力矩$M_{\rm f}$分别为

$ M_{\rm t} =[(\alpha +\beta \theta)\cos (\phi -\pi -\theta)+\alpha +\beta \theta]F_{\rm t} \\ M_{\rm f} =\mu (\alpha +\beta \theta)\sin (\phi -\pi -\theta)F_{\rm t} $

弹簧外端点$B$在$M_{\rm t}$与$M_{\rm f}$作用下发生相对位移，并由位置约束计算法向力$F_{\rm f}$沿$B$点切线方向的位移，最终结合所受力矩计算弹簧刚度。

为实现刚度计算，基于Ansys-WorkBench有限元分析软件对多分支涡卷扭簧进行了仿真试验^[21]。仿真模型包括弹簧内外轮及弹簧弹性体，接触部分添加绑定接触。仿真过程中在对内轮施加固定约束的同时对外轮施加等效力矩。假设弹簧受力矩$M_{{\text{ext}}}$作用，扭簧内轮不动，扭簧弹性体发生形变，使得扭簧外轮转动角度$\theta_{\rm s}$。此时扭簧外轮外径$R$处的某点由扭簧外轮转动前的$A$点移动到转动后的$B$点。$A$、$B$两点之间的距离$ d $即为扭簧在有限元分析中的最大位移，距离$ d $从仿真结果中获得。则弹簧刚度$ k $为

$ k =\frac{M_{{\text{ext}}}} {\theta_{\rm s}}, \quad \theta_{\rm s} =\arcsin \frac{d}{2R} $

Ansys-WorkBench软件中扭簧的位移结果如图 3所示。由前述计算得到的弹簧刚度为60 N$\cdot $m/rad。

2.2 弹性单元刚度标定

为验证四分支平面涡卷弹簧的性能，设计了针对该多分支弹性单元的刚度标定实验。刚度标定平台如图 4所示，弹簧外轮与平台固定安装，弹簧内轮由伺服电机加载，其中伺服电机基于比例—微分（PI）控制器实现速度闭环控制。扭矩传感器用于测量电机对弹簧内轮施加的扭矩。扭簧材料选择高强度弹簧钢60Si2CrA，其弹性模量206 GPa、屈服极限1 568 Mpa，扭簧加工方式采用中走丝线切割，切割后喷丸处理。考虑弹性元件在外骨骼关节运动的实际工作场景，对刚度标定测试设计如下工作范式：

(1) 以扭矩传感器为力矩反馈，将标定弹簧零力矩状态下的电机位置信息作为单次标定零点。

(2) 选取速度作为单次标定实验的加载速度。

(3) 电机执行匀速旋转运动，在形变测定区间$[-0.35, $ $0.35]$ rad内对扭簧进行双向加载。

(4) 双向加载结束后电机牵引扭簧至零力矩位置，计算当前电机位置与标定零点间的位置偏差。若位置偏差超过1$^{\circ}$，则扭簧发生塑性形变，更换弹簧并缩小形变测定区间，回到步骤(1) 进行重新标定；否则，继续以当前加载速度重复执行标定测试3次。

(5) 调整加载速度，重复标定步骤(1)~(4)。

刚度标定实验中，在$[10, 50]$ r/min区间内均匀选取弹簧的加载速度，其中传感器测量力矩均进行了低通滤波与时间平均处理。为简化刚度计算，对循环加载标定过程中出现的弹性体迟滞现象进行数值平均，即忽略刚度在正、负向加载过程中由于迟滞现象产生的数值波动，以平均刚度作为弹簧均匀刚度。受力—形变关系如图 5、6所示。

2.3 弹性单元刚度切换特性

由图 5可知，弹性单元刚度存在显著分段现象，即刚度特性随弹性形变发生切换：当扭簧形变量小于0.224 rad时，扭簧受力—形变呈线性关系，刚度为57 N$\cdot $m。此后形变速度随扭矩增大而逐渐减缓，该区间内具有较大刚度变化范围，最大刚度为2 500 N$\cdot $m/rad。

定义刚度切换平面的额定形变区间为$[-0.22, 0.22]$ rad，弹性单元等效形变量为$|\Delta d |=0.22$。额定形变区间内弹性体具有线性刚度，而额定形变区间外呈现非线性刚度。采用多项式模型对非线性刚度进行数值拟合，则该弹性单元受力—形变模型为

$ \tau_{{\text{spr}}} = \begin{cases} k\cdot \Delta d, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; | \Delta d |<0.22 \\ p_{1} \Delta d^{4}+p_{2} \Delta d^{3}+p_{3} \Delta d^{2}+p_{4} \Delta d+p_{5}, \;\; {\text{其他}} \end{cases} $

其中$ p_{i}$（$i=1, \cdots, 5$）为多项式模型各阶系数。

事实上，变形量小于0.22 rad时各层扭簧间存在间隙，各层扭簧具有独立的线性刚度；当变形量大于0.22 rad后，相邻扭簧发生接触，且其接触面积随所受扭矩增大而增大，进而导致弹性元件的自由形变区间越来越小。得益于其结构特点，受外力影响，弹性单元在不同的形变区间将改变自身的受力状态，故其在较大形变下具有非线性刚度特性。相较于具有固定刚度的弹性元件，所提出的设计方案在提高承载性能的同时，具有线性/非线性刚度多态特性。

所设计的弹性单元根据所受载荷情况切换线性/非线性刚度，其具体特性体现为：(1) 弹性单元具有非线性刚度当且仅当系统处于过载状态。受非线性刚度影响，弹性单元对于持续的过载力矩能快速响应，形变量在较短时间内收敛至稳态值，实现过载情况下的限位保护。对于瞬时过载情形，弹性元件的刚度特性由线性过渡至非线性以实现力矩平衡，并在载荷卸除后快速收敛至额定形变区间，因此弹性单元在冲击载荷下具有较大阻尼响应。(2) 弹性单元在额定形变区间内具有线性刚度，该区间内非线性特性不显著，等效刚度为常值，弹性元件对外部载荷呈现被动柔顺。

2.4 下肢外骨骼串联弹性关节

现有的外骨骼机器人通常将驱动元件集中布置在关节处，使得关节轴向尺寸较大，同时连杆空间难以有效利用。故本文提出了一种分体式串联弹性关节，以优化关节紧凑性、提高连杆空间利用率。

所设计的串联弹性关节如图 7、8所示，其额定功率430 W，额定力矩51.8 N$\cdot $m，厚度65 mm，关节宽度85 mm，减速器模块质量1 010 g，电机模块质量860 g，具有结构扁平、模块化、轻量化特点。

串联弹性关节包含电机模块和减速器模块两部分，采用分体式设计。减速器输出轴与串联弹性关节旋转轴对齐，电机配置在连杆方向，模块间由同步齿形带传动。所设计的串联弹性关节搭载TQ-borodrive ILM50x08盘式电机与Harmonic CSD-20-100-GR-BB谐波减速器，其中新型弹性单元与谐波减速器中的波发生器配合构成串联弹性驱动器。串联弹性关节在外骨骼矢状面具有主动自由度，为防止人机交互过程中其他方向的力矩对关节造成损坏，在关节前后连杆间安装交叉滚子轴承从而保护弹性单元与谐波减速器。串联弹性关节通过双编码器分别测量电机角度（采用Posital公司的KCD-BC00B传感器测量）与弹簧形变量（采用RLS公司RM08SC传感器测量），电机及各传感器数据由关节内的电机驱动器（Synapticon公司，型号NODE E.1000）通过EtherCAT通信协议传输至主控制器。

串联弹性关节中，运动首先通过电机转子—转轴传递至电机端同步带轮，再由同步带传递至减速器端同步带轮，进而由减速器输出端传递至扭簧，最终驱动末端连杆。模块化关节等效传动效率$\eta =0.7$，实际电机力矩计算满足：

$ \tau_{{\text{mot}}} = \begin{cases} \tau_{{\text{current}}} \cdot \eta, \;\;\; {\rm sgn }(\dot{\theta} \cdot \tau_{{\text{current}}})>0 \\ \tau_{{\text{current}}} /\eta, \;\;\; {\rm sgn }(\dot{\theta} \cdot \tau_{{\text{current}}})\leqslant 0 \end{cases} $

其中，$\tau_{{\text{current}}}$为电流环力矩，$\dot{\theta}$为电机角速度，$\tau_{{\text{current}}}$和$\dot{\theta}$的方向与规定的正方向一致时为正数，相反为负数，sgn$(\cdot)$为符号函数。

3 动力学建模与参数辨识（Dynamics modeling and parameter identification）

所设计的串联弹性关节兼具刚性连杆结构与柔性元件的动力学特性，构成刚柔耦合动力学模型。该关节的动力学模型基于以下基本假设：

假设 1：在刚度切换特性下，关节发生较大形变时受其结构特性影响会进行收敛运动，而关节形变较小时对关节柔性的影响近似为线性。

假设 2：驱动电机与负载连杆具有独立自由度。

假设 3：串联弹性关节的驱动元件被建模为质心位于关节旋转轴上的均匀刚体。

单关节串联弹性关节传动模型如图 9所示，其中连杆位置记为$ \boldsymbol q $，减速器（减速比$N$为特定数值）输出端位置记为${\boldsymbol\theta}_{\rm M} =N{\boldsymbol\theta}_{\rm m}$。关节电机直接作用于弹簧，并使弹簧产生力矩，弹簧作为中间环节再进一步为连杆提供运动所需的动力。刚性连杆与柔性元件间通过弹簧串联。

3.1 串联弹性关节动力学模型

为实现基于模型的控制算法设计，首先对串联弹性关节的动力学模型特性进行分析。对广义物理系统，由拉格朗日法建立系统动力学模型：

$ F_{i} =\frac{\rm d}{{\rm d} t}\frac{\partial L}{\partial \dot{q}_{i}} -\frac{\partial L}{\partial q_{i}}, \quad i=1, 2, \cdots, n \\ L =T-P $

(2)

其中，$L$为拉格朗日函数，$T$为系统总动能，$P$为系统总势能，$ q_{i}$为关节$ i $处连杆变量。串联弹性关节的弹性单元将弹簧形变产生的势能引入到系统总势能：

$ P=P_{{\text{grav}}} +P_{{\text{elas}}} $

(3)

由假设2可知，重力引入的势能$P_{{\text{grav}}}$由连杆质量和电机质量共同提供，而弹性势能$P_{{\text{elas}}}$由弹簧形变量${\boldsymbol{\alpha}} ={\boldsymbol{q}}-{\boldsymbol{\theta}}_{\rm M}$产生。又由假设3可知，$P_{{\text{grav}}}$由连杆位置决定，电机质量产生的重力势能部分与电机旋转位置无关，仅由连杆位置决定。定义$P_{{\text{grav, link}}} ({\boldsymbol{q}})$为连杆重力势能，$P_{{\text{grav, motor}}} ({\boldsymbol{q}})$为电机的重力势能，则：

$ P_{{\text{grav}}} =P_{{\text{grav, link}}} ({\boldsymbol{q}})+P_{{\text{grav, motor}}} ({\boldsymbol{q}}) $

(4)

对于弹性势能$P_{{\text{elas}}}$，令${\boldsymbol{K}}$为弹簧刚度对角矩阵，则由假设1可得：

$ P_{{\text{elas}}} =\frac{1}{2}({\boldsymbol{q}}-{\boldsymbol{\theta}}_{\rm M})^{\rm T} {\boldsymbol{K}}( {\boldsymbol{q}}-{\boldsymbol{\theta}}_{\rm M}) $

(5)

类似的，由假设2，系统总动能$T$由连杆和电机的动能共同构成：

$ T =T_{{\text{link}}} +T_{{\text{motor}}} \\ T_{{\text{link}}} =\frac{1}{2}{\boldsymbol{\dot{q}}}^{\rm T}{\boldsymbol{M}}_{\rm L} ({\boldsymbol{q}}){\boldsymbol{\dot{q}}} \\ T_{{\text{motor}}} =\sum\limits _{i=1}^n \frac{1}{2}\left(m_{{\rm r}i} {\boldsymbol{v}}_{{\rm r}i}^{\rm T} {\boldsymbol{v}}_{{\rm r}i} +{\boldsymbol{\omega}}_{{\rm r}i}^{\rm T} {\boldsymbol{I}}_{{\rm r}i} {\boldsymbol{\omega}}_{{\rm r}i}\right) $

(6)

其中，${\boldsymbol{M}}_{\rm L} ({\boldsymbol{q}})$为正定对称的连杆惯量矩阵，${\boldsymbol{v}}_{{\rm r}i}$与${\boldsymbol{\omega}}_{{\rm r}i}$分别为第$ i $个电机的质心线速度与角速度，${\boldsymbol{I}}_{{\rm r}i}$为电机惯量对角矩阵，$ m_{{\rm r}i}$为第$ i $个电机的质量。

串联弹性关节动力学模型由刚性连杆子模型、电机—弹簧子模型构成。为建立完整的动力学模型，对串联弹性关节进行分步解耦建模。

(1) 刚性连杆子模型

利用投影牛顿－欧拉方法（projected Newton-Euler）建立刚性连杆的动力学模型并设置模型参数。连杆动力学基本模型为

$ {\boldsymbol{M}}_{\rm L} ({\boldsymbol{q}}){\boldsymbol{\ddot{q}}}+{\boldsymbol{C}}({\boldsymbol{q}}, \dot{\boldsymbol q}){\boldsymbol{\dot{q}}}+ {\boldsymbol{G}} ({\boldsymbol{q}}) =\boldsymbol \tau $

(7)

其中，${\boldsymbol{C}}({\boldsymbol{q}}, \dot{\boldsymbol q})$为科氏力及离心力矩阵，${\boldsymbol{G}}({\boldsymbol{q}})$为重力矩阵，${\boldsymbol{\tau}}$为关节力矩向量。

模型中各矩阵参数由下式计算：

$ {\boldsymbol{M}}_{\rm L} ({\boldsymbol{q}}) =\sum\limits _{i=1}^N \left( {\boldsymbol{J}}_{{\rm S}i}^{\rm T} m_{i} {\boldsymbol{J}}_{{\rm S}i} +{\boldsymbol{J}}_{{\rm R}i}^{\rm T} {\boldsymbol{I}}_{i} {\boldsymbol{J}}_{{\rm R}i}\right) \\ {\boldsymbol{C}}({\boldsymbol{q}}, \dot{\boldsymbol q}){\boldsymbol{\dot{q}}} =\sum\limits _{i=1}^N \left({{\boldsymbol{J}}_{{\rm S}i}^{\rm T} m_{i} {\boldsymbol{\dot{J}}}_{{\rm S}i} {\boldsymbol{\dot{q}}}+{\boldsymbol{J}}_{{\rm R}i}^{\rm T} ({\boldsymbol I_{i} {\boldsymbol{\dot{J}}}_{{\rm R}i} {\boldsymbol{\dot{q}}}+{\boldsymbol{\varOmega}}_{i} \times {\boldsymbol{I}}_{i} {\boldsymbol{\varOmega}}_{i}})}\right) \\ {\boldsymbol{G}}({\boldsymbol{q}}) =\sum\limits _{i=1}^N \left({-{\boldsymbol{J}}_{{\rm S}i}^{\rm T} m_{i} {{\boldsymbol{g}}}}\right) $

其中，${\boldsymbol{\varOmega}}_{i}$为惯性坐标系下的绝对角速度，${\boldsymbol{J}}_{{\rm S}i}$与${\boldsymbol{J}}_{{\rm R}i}$分别为对应连杆$ i $的绝对线速度与角速度的雅可比矩阵，${\boldsymbol{I}}_{i}$为连杆$ i $的惯量矩阵。

(2) 电机—弹簧子模型

由拉格朗日动力学模型式(2)、系统势能模型式(3)与动能模型式(6)，进一步地，电机引入的动能项描述为

$ T_{{\text{motor}}} =\frac{1}{2}{\boldsymbol{\dot{q}}}^{\rm T} \left({\boldsymbol{M}}_{\rm R} ({\boldsymbol{q}})+{{\boldsymbol{S}}}({\boldsymbol{q}}){{\boldsymbol{B}}}^{-1}{{\boldsymbol{S}}}^{\rm T}({\boldsymbol{q}})\right){\boldsymbol{\dot{q}}}+\notag \\ {\boldsymbol{\dot{q}}}^{\rm T} {\boldsymbol{S}} ({\boldsymbol{q}}){\boldsymbol{\dot{\theta}}}_{\rm M} +\frac{1}{2}{\boldsymbol{\dot{\theta}}}_{\rm M}^{\rm T} {\boldsymbol{B}}\dot{\boldsymbol \theta}_{\rm M} $

(8)

其中，${\boldsymbol{M}}_{\rm R} ({\boldsymbol{q}})$为电机引入的连杆附加惯量矩阵，${\boldsymbol{B}}$为考虑减速比后的电机惯量矩阵，${\boldsymbol{S}}({\boldsymbol{q}})$为电机－连杆耦合矩阵。

结合弹簧柔性引入的弹性势能式(5) 与拉格朗日方程(2)，导出串联弹性关节动力学模型为

$ \begin{pmatrix} {\boldsymbol{H}}({\boldsymbol{q}}) & {\boldsymbol{S}}({\boldsymbol{q}}) \\ {\boldsymbol{S}}^{\rm T}({\boldsymbol{q}}) & {\boldsymbol{B}} \end{pmatrix} \begin{pmatrix} {\boldsymbol{\ddot{q}}} \\ {\boldsymbol{\ddot{\theta}}}_{\rm M} \end{pmatrix}+ \begin{pmatrix} {\boldsymbol{C}}({\boldsymbol{q}}, \dot{\boldsymbol q})+{\boldsymbol{c}}_{1} ({\boldsymbol{q}}, \dot{\boldsymbol q}, \dot{\boldsymbol \theta}) \\ {\boldsymbol{c}}_{2} ({\boldsymbol{q}}, \dot{\boldsymbol q}) \end{pmatrix}+ \notag \\ \quad \begin{pmatrix} {\boldsymbol{G}}({\boldsymbol{q}})+{\boldsymbol{K}}({\boldsymbol{q}}-{\boldsymbol{\theta}}_{\rm M}) \\ {\boldsymbol{K}}({\boldsymbol{\theta}}_{\rm M} -{\boldsymbol{q}}) \end{pmatrix}= \begin{pmatrix} 0 \\ {\boldsymbol{\tau}}_{{\text{mot}}} \end{pmatrix} $

(9)

其中，${\boldsymbol{H}}({\boldsymbol{q}})$为总惯量矩阵：

$ {\boldsymbol{H}}({\boldsymbol{q}})={\boldsymbol{M}}_{\rm L} ({\boldsymbol{q}})+{\boldsymbol{M}}_{\rm R} ({\boldsymbol{q}})+{\boldsymbol{S}}({\boldsymbol{q}}){\boldsymbol{B}}^{-1}{\boldsymbol{S}}^{\rm T}({\boldsymbol{q}}) $

3.2 串联弹性关节动力学退化模型

假设 4：对于较大的关节减速比，电机—连杆间惯量耦合引入的能量相对于系统总能量不显著，电机角速度仅由自身旋转决定。

在外骨骼康复应用场景中，串联弹性关节运动被约束在矢状面内，且连杆与电机连接处存在摩擦力。基于假设4对动力学模型式(9) 进行简化与修正：

(1) 科氏力项退化为平面形式：

$ {\boldsymbol{C}}({\boldsymbol{q}}, {\boldsymbol{\dot{q}}}){\boldsymbol{\dot{q}}}=\sum\limits _{i=1}^n \left({\boldsymbol{J}}_{{\rm S}i}^{\rm T} m_{i} {\boldsymbol{\dot{J}}}_{{\rm S}i} {\boldsymbol{\dot{q}}}+{\boldsymbol{J}}_{{\rm R}i}^{\rm T} {\boldsymbol{I}}_{i} {\boldsymbol{\dot{J}}}_{{\rm R}i} {\boldsymbol{\dot{q}}}\right) $

此外，科氏力项$ {\boldsymbol{c}}_{1} ({\boldsymbol{q}}, {\boldsymbol{\dot{q}}}, {\boldsymbol{\dot{\theta}}})={\boldsymbol{c}}_{2} ({\boldsymbol{q}}, {\boldsymbol{\dot{q}}})\equiv {\bf 0}$。

(2) 由假设4，电机—连杆耦合矩阵${\boldsymbol{S}}({\boldsymbol{q}}$) 退化为常数矩阵${\boldsymbol{S}}$，且${\boldsymbol{S}}\equiv {\bf 0}$，即电机总角动能为$\dfrac{1}{2}{\boldsymbol{\dot{\theta}}}_{\rm M}^{\rm T} {\boldsymbol{B}}\dot{{\boldsymbol\theta}}_{\rm M}$。

(3) 当关节进行反向驱动时，摩擦力对外力的阻碍作用将被放大。考虑库仑力与黏滞摩擦力模型，引入库仑摩擦系数$\tau_{\rm k}$以及黏滞摩擦系数$D_{\rm v}$，则关节$ i $处摩擦力模型为

$ \tau_{\rm f}^{i} =\tau_{\rm k}^{i} {\rm{sgn}}({\dot{\theta}}_{{\rm M}, i})+D_{\rm v}^{i} {\dot{\theta}}_{{\rm M}, i} $

(4) 串联弹性关节中弹簧运动具有阻尼$D_{\rm s}$。

基于上述讨论，面向控制系统设计的退化动力学模型为

$ \begin{gathered} {\boldsymbol{H}}({\boldsymbol{q}}){\boldsymbol{\ddot{q}}}+{\boldsymbol{C}}({\boldsymbol{q}}, {\boldsymbol{\dot{q}}}){\boldsymbol{\dot{q}}}+{\boldsymbol{G}}({\boldsymbol{q}})+{\boldsymbol{\tau}}_{\rm f} -{\boldsymbol{K}}({\boldsymbol{\theta}}_{\rm M} -{\boldsymbol{q}}) ={\bf 0} \\ {\boldsymbol{B}}\ddot{{\boldsymbol\theta}}_{\rm M} +{\boldsymbol{K}}( {\boldsymbol{\theta}}_{\rm M} - {\boldsymbol{q}})+{\boldsymbol{D}}_{\rm s} ({\boldsymbol{\dot{\theta}}}_{\rm M} -{\boldsymbol{\dot{q}}}) ={\boldsymbol{\tau}}_{{\text{mot}}} \end{gathered} $

(10)

其中，${\boldsymbol{\tau}}_{\rm f}$为关节摩擦力向量。

3.3 动力学参数辨识 3.3.1 激励轨迹设计

一般地，机器人系统辨识问题需要设计激励轨迹以充分激励系统并获得充分的辨识数据^[22]。由退化动力学模型式(10) 可知，待辨识模型包含连杆变量$({\boldsymbol{q}}, {\boldsymbol{\dot{q}}}, {\boldsymbol{\ddot{q}}})$与电机变量$({\boldsymbol{\theta}}_{\rm M}, {\dot{{\boldsymbol\theta}}}_{\rm M}, {\ddot{{\boldsymbol\theta}}}_{\rm M})$，为处理变量间的耦合关系，仅考虑连杆变量建立次优问题设计激励轨迹。此时动力学线性化模型为

$ {\boldsymbol{Y}}_{{\text{red}}} ({\boldsymbol{q}}, {\boldsymbol{\dot{q}}}, {\boldsymbol{\ddot{q}}}) {\boldsymbol{\varPhi}}_{{\text{red}}} ={\boldsymbol{H}}({\boldsymbol{q}}){\boldsymbol{\ddot{q}}}+{\boldsymbol{C}}({\boldsymbol{q}}, {\boldsymbol{\dot{q}}}){\boldsymbol{\dot{q}}}+{\boldsymbol{G}}({\boldsymbol{q}}) $

其中，${\boldsymbol{Y}}_{{\text{red}}} ({\boldsymbol{q}}, {\boldsymbol{\dot{q}}}, {\boldsymbol{\ddot{q}}})$为回归矩阵，${\boldsymbol{\varPhi}}_{{\text{red}}}$为辨识参数向量。

激励轨迹模型选择高阶傅里叶级数式(11)，以获得微分连续且具有周期性的激励信号。

$ q_{i} (t)=\sum\limits _{l=1}^{n_{\rm f}} a_{l}^{i} \sin (w_{\rm f} \textit{ltTs})-b_{l}^{i} \cos (w_{\rm f} \textit{ltTs}) +a_{i0} $

(11)

选取D-optimality准则^[23]建立最优化问题，目标函数式(12) 为回归矩阵行列式的最大值。最优化问题中决策变量为傅里叶级数参数组$(a_{l}^{i}, $ $b_{l}^{i}, $ $a_{i0})$，约束条件为傅里叶模型相关的等式约束以及关节角度、速度和加速度相关的不等式约束。

$ \min\limits _{a_{\rm l}^{i}, b_{\rm l}^{i}, a_{i0}} -\log \left(\det ({\boldsymbol{Y}}_{{\text{red}}}^{\rm T} {\boldsymbol{Y}}_{{\text{red}}})\right) $

(12)

辨识激励实验中，通过参数整定后的PD控制器进行参考轨迹的位置追踪。为有效激励弹簧特性，对连杆输出端配置等效负重，激励轨迹频率选取200 Hz，并持续进行20次的轨迹激励。对数据低通滤波和时间平均处理后，激励实验结果如图 10所示。

3.3.2 解耦动力学辨识

动力学模型式(10) 中，电机侧动力学与连杆侧动力学通过弹簧力矩${\boldsymbol{K}}({\boldsymbol{\theta}}_{\rm M} -{\boldsymbol{q}})$耦合。进行参数辨识时，需要首先进行动力学模型的解耦。

(1) 将弹簧力矩作为系统内力建立等效动力学模型：

$ {\boldsymbol{H}}({\boldsymbol{q}}){\boldsymbol{\ddot{q}}}+{\boldsymbol{B}}\ddot{{\boldsymbol\theta}}_{\rm M} +{\boldsymbol{C}}({\boldsymbol{q}}, {\boldsymbol{\dot{q}}}){\boldsymbol{\dot{q}}}+{\boldsymbol{G}}({\boldsymbol{q}})+{\boldsymbol{\tau}}_{\rm f} ={\boldsymbol{\tau}}_{{\text{mot}}} $

(13)

进而对解耦模型式(13) 作线性分离：

$ {\boldsymbol{Y}}_{{\text{dec}}}{\boldsymbol{\varPhi}}_{{\text{dec}}} = {\boldsymbol{H}}({\boldsymbol{q}}){\boldsymbol{\ddot{q}}}+{\boldsymbol{B}}\ddot{{\boldsymbol\theta}}_{\rm M}+ {\boldsymbol{C}}({\boldsymbol{q}}, {\boldsymbol{\dot{q}}}){\boldsymbol{\dot{q}}}+{\boldsymbol{G}}({\boldsymbol{q}})+{\boldsymbol{\tau}}_{\rm f} $

(14)

(2) 结合退化动力学模型式(10) 与线性分离模型式(14)，对弹簧刚度与弹簧阻尼系数进行参数辨识。2.2节中，对串联弹性关节的弹性单元进行了离线刚度标定，而本节中对弹簧刚度进行系统辨识的另一目的是补偿关节装配误差。此时，待辨识动力学模型为

$ {\boldsymbol{B}}\ddot{{\boldsymbol\theta}}_{\rm M} +{\boldsymbol{K}} ({\boldsymbol{\theta}}_{\rm M} -{\boldsymbol{q}})+{\boldsymbol{D}}_{\rm s} ({\boldsymbol{\dot{\theta}}}_{\rm M} -{\boldsymbol{\dot{q}}})={\boldsymbol{\tau}}_{{\text{mot}}} $

对应的线性分离参数化形式为

$ {\boldsymbol{Y}}_{{\text{spr}}} {\boldsymbol{\varPhi}}_{{\text{spr}}} ={\boldsymbol{B}}\ddot{{\boldsymbol\theta}}_{\rm M} +{\boldsymbol{K}}({\boldsymbol{\theta}}_{\rm M} -{\boldsymbol{q}})+{\boldsymbol{D}}_{\rm s} ({\boldsymbol{\dot{\theta}}}_{\rm M} -{\boldsymbol{\dot{q}}}) $

对动力学模型式(10) 解耦后，采用带约束的最小二乘法进行参数辨识，将串联弹性关节在CAD设计软件中的参考计算值作为约束条件中的下界。为验证参数辨识结果的有效性，基于辨识结果进行逆动力学计算以获取电机力矩估计值，并与辨识过程中电机实际力矩进行对比，其中对电流环测量值进行零相位滤波以处理测量噪声。本文以单自由度串联弹性关节为对象进行了参数辨识实验，其逆动力学验证如图 11所示，参数辨识结果如表 1所示。

4 迭代线性化模型预测控制设计（Design of iterative linearized model predictive control）

串联弹性关节适用于在下肢康复外骨骼矢状面运动，受重力项与摩擦力作用，系统具有非线性。同时具有自主运动能力的外骨骼需要考虑系统控制功耗，并有效抑制交互过程中的外部环境扰动。控制器应满足如下设计需求：

(1) 控制器能够实现控制约束，避免控制量出现尖峰与抖动以提高交互舒适性；

(2) 控制器能够有效抑制交互过程中的输入输出扰动，具有鲁棒性。

基于退化动力学模型式(10)，以单自由度串联弹性关节为被控对象设计控制算法，从而建立非线性最优化问题，并在模型预测控制（MPC）框架下进行数值求解。针对非线性模型预测控制（NMPC）的数值计算难点，利用线性化方法沿系统状态轨迹对原系统进行迭代变换和周期性更新，最终结合Laguerre网络对控制律进行滚动优化。

4.1 非线性优化问题设计

将前述动力学模型式(10) 等效为单输入单输出（SISO）非线性系统：

$ \begin{cases} {\boldsymbol{\dot{x}}}={\boldsymbol{f}}({\boldsymbol{x}})+{\boldsymbol{g}}({\boldsymbol{x}})u \\ y=x_{1} \end{cases} $

(15)

其中，

$ {\boldsymbol{x}} =[x_{1} \; x_{2} \; x_{3} \; x_{4} ]^{\rm T} =[q \; \dot{q} \; \theta_{\rm M} \; \dot{\theta}_{\rm M} ]^{\rm T}\\ {\boldsymbol{f}}({\boldsymbol{x}}) = \begin{pmatrix} x_{2} \\ \dfrac{K(x_{3} -x_{1})-Cx_{2} -mgd\sin x_{1}-\tau_{\rm f} (x_{4})}{H} \\ x_{4} \\ -\dfrac{K(x_{3} -x_{1})+D_{\rm s} (x_{4} -x_{2})}{B} \end{pmatrix} \\ {\boldsymbol{g}}({\boldsymbol{x}}) = \left(0 \; 0 \; 0 \; \frac{1}{B} \right)^{\rm T}, \qquad u=\tau_{\rm M} $

系统在$ k $时刻，带约束的非线性优化问题为$\min \limits_{{u}} J(k) $，其中：

$ J(k)=\| {\boldsymbol{x}}^{ k}-{\boldsymbol{x}}_{\rm d}^{ k}\|^{2}+| u^{ k}-u_{\rm d}^{ k} |^{2} \\ {\text{s.t. }} \begin{cases} u_{\min}^{ k} \leqslant u^{ k} \leqslant u_{\max}^{ k} \notag\\ \dot{{\boldsymbol{x}}}^{ k}={\boldsymbol{f}}({\boldsymbol{x}}^{ k})+{\boldsymbol{g}}({\boldsymbol{x}}^{ k})u^{ k} \end{cases} $

(16)

其中，上标$ k $表示对应时刻，$ {\boldsymbol{x}}_{\rm d}^{ k}$为$ k $时刻系统期望状态，$ u_{\rm d}^{ k}$为$ k $时刻系统期望控制量，$ u_{\max}^{ k}$与$ u_{\min}^{ k}$分别为$ k $时刻系统控制量的上下界。

系统状态方程(15) 存在重力项和摩擦力项而引入非线性因素，因而无法确保优化问题(16) 为凸优化问题，更进一步地，对优化问题(16) 的直接求解不能确保全局最优解的存在性。NMPC的求解首先需要对非凸问题进行凸化，并利用数值方法进行在线计算^[24]。然而目前主流的NMPC求解速度与问题规模相关，尤其对预测时域$N$较为敏感。下肢外骨骼属于实时人机交互系统，对控制频率具有较高要求（不低于500 Hz），对NMPC的直接求解难以满足计算时间要求。为解决上述问题，设计了一种迭代线性化模型预测控制算法，在每一控制周期内对当前系统的工作状态进行线性化，将非线性系统沿实际状态轨迹等效为线性系统，并采用线性MPC算法逐次求解原非线性问题。

4.2 迭代线性化模型预测控制律

任意时刻，名义连续系统一致满足非线性关系式(15)，同时期望轨迹$({\boldsymbol{x}}_{\rm d}, {\boldsymbol{\dot{x}}}_{\rm d})$应满足

$ {\boldsymbol{\dot{x}}}_{\rm d} ={\boldsymbol{f}}({\boldsymbol{x}}_{\rm d})+{\boldsymbol{g}}({\boldsymbol{x}}_{\rm d})u $

(17)

考虑系统工作点${\boldsymbol{x}}={\boldsymbol{x}}_{\rm d}$，则对于该平衡点处状态微分方程具有1阶泰勒展开线性逼近：

$ {\boldsymbol{\dot{x}}} \approx {\boldsymbol{f}}({\boldsymbol{x}}_{\rm d})+{\boldsymbol{g}}({\boldsymbol{x}}_{\rm d})u+\frac{\partial {\boldsymbol{f}}}{\partial {\boldsymbol{x}}} \bigg| _{{\boldsymbol{x}}={\boldsymbol{x}}_{\rm d}} ({\boldsymbol{x}}-{\boldsymbol{x}}_{\rm d}) \notag \\ ={\boldsymbol{f}}({\boldsymbol{x}}_{\rm d})+{\boldsymbol{g}}({\boldsymbol{x}}_{\rm d})u+{\boldsymbol{J}}_{{\boldsymbol{f}}} ({\boldsymbol{x}})({\boldsymbol{x}}-{\boldsymbol{x}}_{\rm d}) $

(18)

其中，${\boldsymbol{J}}_{{\boldsymbol{f}}} ({\boldsymbol{x}})$为$ {\boldsymbol{f}}({\boldsymbol{x}}_{\rm d})$相对于${\boldsymbol{x}}$的雅可比矩阵。

令矩阵${\boldsymbol{A}}_{\rm L} (t)$与矩阵${\boldsymbol{B}}_{\rm L} (t)$为$ t $时刻的线性化矩阵，$\tilde{{\boldsymbol{x}}}={\boldsymbol{x}}-{\boldsymbol{x}}_{\rm d}$为状态偏差，则由式(17)(18) 可得系统状态偏差模型：

$ {\boldsymbol{\dot{\tilde{x}}}}={\boldsymbol{J}}_{\boldsymbol f} ({\boldsymbol{x}})\tilde{{\boldsymbol{x}}}+{\boldsymbol{g}}({\boldsymbol{x}})\tilde{u} ={\boldsymbol{A}}_{\rm L} (t)\tilde{{\boldsymbol{x}}}+{\boldsymbol{B}}_{\rm L} (t)\tilde{u} $

(19)

进一步，对连续系统式(19) 进行欧拉离散化，得到离散线性时变（LTV）系统：

$ \begin{split} \tilde{{\boldsymbol{x}}}(k+1)& ={\boldsymbol{A}}_{k} \tilde{{\boldsymbol{x}}}(k)+{\boldsymbol{B}}_{k} \tilde{{u}}(k) \\ \tilde{y}(k)& ={\boldsymbol{C}}_{k} \tilde{{\boldsymbol{x}}}(k) \end{split} $

(20)

考虑$ k $时刻至$ k+N$时刻内系统状态量与控制量，分别对状态偏差$\tilde{{\boldsymbol{x}}}$与控制量偏差$\tilde{u}$配置状态权重矩阵${\boldsymbol{Q}}$和控制权重因子${R}$，控制可行域记为${ U}$，则非线性优化问题式(16) 能够通过线性凸优化问题求解：

$ \min\limits _{\tilde{u}_{k} \in { U}} J(k) \\ {\text{s.t. }} \begin{cases} \tilde{{\boldsymbol{x}}}(k+1)={\boldsymbol{A}}_{k} \tilde{{\boldsymbol{x}}}(k)+{\boldsymbol{B}}_{k} \tilde{u}(k) \\ \tilde{{\boldsymbol{x}}}(k)={\boldsymbol{x}}(k)-{\boldsymbol{x}}_{\rm d} (k) \\ \tilde{u}(k)=u(k)-u_{\rm d} (k) \end{cases}\notag $

(21)

其中，$ k $时刻目标函数为

$ \tilde{{\boldsymbol{x}}}^{\rm T}(k+N |k){\boldsymbol{P}} \tilde{{\boldsymbol{x}}}(k+N |k)+\sum\limits _{i=1}^{N-1} \tilde{{\boldsymbol{x}}}^{\rm T}(k+i| k){\boldsymbol{Q}} \tilde{{\boldsymbol{x}}}(k+i| k)+ \\ \quad \sum\limits _{i=1}^{N-1} R\tilde{u}^{2}( k+i |k)=J(k) $

为使系统在$ k+N$步收敛，令${\boldsymbol{P}}$为终值权重矩阵且${\boldsymbol{P}}$应满足李雅普诺夫方程：

$ {\boldsymbol{P}}-{\boldsymbol{A}}^{\rm T}{\boldsymbol{PA}} = {\boldsymbol{Q}} $

(22)

一般地，在$ k $时刻离散MPC关于最优化问题(21) 的解为最优控制向量$\tilde{{\boldsymbol{U}}}_{k}$，并保留向量$\tilde{{\boldsymbol{U}}}_{k}$中首项的控制量作为实际控制输入。为提高优化求解效率、降低在线求解复杂度，结合Laguerre网络^[25]对最优化问题(21) 引入一组离散正交基函数，从而构造系统最优解。

基于Laguerre网络，$ k $时刻后的任意采样时刻$ k+m $处系统控制量可描述为

$ \Delta u(k+m) ={\boldsymbol{L}}^{\rm T}(m){\boldsymbol{\eta}}_{k} \\ {\boldsymbol{L}}(m) =[l_{1} (m), l_{2} (m), \cdots, l_{N_{\rm L}} (m)]^{\rm T} \\ {\boldsymbol{\eta}}_{k} =[c_{1}, c_{1}, \cdots, c_{N_{\rm L}}]^{\rm T} $

(23)

其中$m=1, 2, \cdots, N$，$c_{i}$为Laguerre网络参数，$ l_{i} (m)$为Laguerre函数，$N_{\rm L}$为Laguerre网络阶数，$ m $为预测步长。

则当前时刻的最优控制输入可描述为

$ \Delta \tilde{u}^{*} (k)={\boldsymbol{L}}^{\rm T}(0){\boldsymbol{\eta}}_{k} $

假设 5：存在$ u\in { U}$使LTV系统式(20) 镇定。

假设 6：预测时域$N$的选取应在满足追踪效果的条件下取最小值以满足实时计算需求。

假设 7：在较高的系统控制频率下，系统状态量${\boldsymbol{x}}$的变化率较小。

定理 1：在假设5~7成立的条件下，对LTV系统式(20) 求1阶差分，系统可等效为带有积分器的增广模型：

$ {\boldsymbol{X}}(k+1) ={\boldsymbol{AX}}(k)+{\boldsymbol{B}}\Delta u(k) \\ \tilde{y}(k) ={\boldsymbol{CX}}(k) $

(24)

其中，

$ \begin{align*} {\boldsymbol{X}}(k) = \begin{pmatrix} \Delta {\boldsymbol{x}}(k) \\ \tilde{y}(k) \end{pmatrix}\\ {\boldsymbol{A}} =\begin{pmatrix} {\boldsymbol{A_{k}}} & {\bf 0} \\ {\boldsymbol{CA}}_{k} &{\boldsymbol{I}} \end{pmatrix}, \;\;{\boldsymbol{B}}=\begin{pmatrix} {\boldsymbol{B_{k}}} \\ {\boldsymbol{CB}}_{k} \end{pmatrix}, \;\;{\boldsymbol{C}}=\begin{pmatrix} {\bf 0} \\ {\boldsymbol{I}} \end{pmatrix}^{\rm T} \end{align*} $

若满足条件：

1) 最优化问题式(21) 在$ k+N$时刻满足等式约束$\tilde{{\boldsymbol{x}}}(k+N)=0$

2) 算法1对最优化问题具有可行解

那么，增广系统式(24) 关于原点是渐近稳定的。

算法1迭代线性化MPC算法
输入: 给定权重矩阵$({\boldsymbol{Q}}, {R})$，控制量可行域$U$以及初值条件$ {\boldsymbol{x}}(k)$。
1: 计算线性化轨迹${\boldsymbol{x}}_{\rm d}, {u}_{\rm d}$；
2: 根据初值条件$ {\boldsymbol{x}}(k)$初始化$ k $时刻离散LTV系统式(20)，由李雅普诺夫方程(22) 更新终值权重矩阵；
3: 结合Laguerre网络对凸优化问题式(21) 求解最优控制增量$\Delta \tilde{u}^{*} (k)$；
4: $\tilde{u}^{*} (k)\leftarrow \tilde{u}(k-1)+\Delta \tilde{u}^{*} (k)$；
5: $ u^{*} (k)\leftarrow \tilde{u}^{*} (k)+u_{\rm d} (k)$；
6: 更新非线性系统式(15) 控制量$ u(k)=u^{*} (k)$；
7: $ k\leftarrow k+1$，回到步骤1。

证明：首先考察LTV系统式(20)，易知$({\boldsymbol{B}}_{k}, {\boldsymbol{C}}_{k})$为常向量，${\boldsymbol{A}}_{k}$为时变矩阵。对${\boldsymbol{A}}_{k}$求1阶差分后矩阵时变参数项为$\cos x(k) -\cos x(k-1)$，由假设5及拉格朗日中值定理可知${\boldsymbol{A}}_{k}$为慢时变矩阵，即存在小常数$\varepsilon $使不等式成立：

$ | \cos x(k) -\cos x(k-1) |\leqslant \Delta x(k) <\varepsilon $

对系统1阶作差分计算，则

$ \Delta \tilde{{\boldsymbol{x}}}(k+1) ={\boldsymbol{A}}_{k} \tilde{{\boldsymbol{x}}}(k) -{\boldsymbol{A}}_{k-1} \tilde{{\boldsymbol{x}}}(k-1) +{\boldsymbol{B}}_{k} \tilde{u}_{k} -{\boldsymbol{B}}_{k-1} \tilde{u}(k-1) \\ \approx {\boldsymbol{A}}_{k} (\tilde{{\boldsymbol{x}}}(k) -\tilde{{\boldsymbol{x}}}(k-1))+{\boldsymbol{B}}_{k} (\tilde{u}(k) -\tilde{u}(k-1)) \\ ={\boldsymbol{A}}_{k} \Delta \tilde{{\boldsymbol{x}}}(k) +{\boldsymbol{B}}_{k} \Delta \tilde{u}(k) $

根据假设4，当预测时域$N$较小时存在：

$ \prod\limits _k^N {\boldsymbol{A}}_{k} \approx {\boldsymbol{A}}_{k}^{N} $

在预测时域为$N$、控制时域为$N_{\rm c}$条件下对系统式(20) 进行状态预测，则系统输出序列的向量形式为

$ {\boldsymbol{Y}}={\boldsymbol{F}}{\boldsymbol x}(k)+{\boldsymbol{\varPhi}} \Delta {\boldsymbol{U}} $

(25)

其中，

$ \begin{align*} {\boldsymbol{\varPhi}} & = \begin{pmatrix} {\boldsymbol{C}}_{k} {\boldsymbol{B}}_{k} & 0 & \cdots & 0 \\ {\boldsymbol{C}}_{k} {\boldsymbol{A}}_{k} {\boldsymbol{B}}_{k} & {\boldsymbol{C}}_{k} {\boldsymbol{B}}_{k} & \cdots & 0 \\ \vdots & \vdots & \cdots & 0 \\ {\boldsymbol{C}}_{k} \displaystyle \prod _k^{N-1} {\boldsymbol{A}}_{k} {\boldsymbol{B}}_{k} & {\boldsymbol{C}}_{k} \displaystyle \prod _k^{N-2} {\boldsymbol{A}}_{k} {\boldsymbol{B}}_{k} & \cdots & {\boldsymbol{C}}_{k} \displaystyle \prod _k^{N-N_{\rm c}} {\boldsymbol{A}}_{k} {\boldsymbol{B}}_{k} \end{pmatrix}_{N\times N_{\rm c}} \\ & \approx \begin{pmatrix} {\boldsymbol{C}}_{k} {\boldsymbol{B}}_{k} & 0 & \cdots & 0 \\ {\boldsymbol{C}}_{k} {\boldsymbol{A}}_{k} {\boldsymbol{B}}_{k} & {\boldsymbol{C}}_{k} {\boldsymbol{B}}_{k} & \cdots & 0 \\ \vdots & \vdots & \cdots & 0 \\ {\boldsymbol{C}}_{k} {\boldsymbol{A}}_{k}^{N-1}{\boldsymbol{B}}_{k} & {\boldsymbol{C}}_{k} {\boldsymbol{A}}_{k}^{N-{2}}{\boldsymbol{B}}_{k} & \cdots & {\boldsymbol{C}}_{k} {\boldsymbol{A}}_{k}^{N-N_{\rm c}}{\boldsymbol{B}}_{k} \end{pmatrix}_{N\times N_{\rm c}} \\ {\boldsymbol{F}}& = \begin{pmatrix} {\boldsymbol{C}}_{k} {\boldsymbol{A}}_{k} \\ {\boldsymbol{C}}_{k} {\boldsymbol{A}}_{k}^{2} \\ \vdots \\ {\boldsymbol{C}}_{k} {\boldsymbol{A}}_{k}^{N} \end{pmatrix}_{N\times 1}, \quad \Delta {\boldsymbol{U}}= \begin{pmatrix} \Delta u(k) \\ \Delta u(k+1) \\ \vdots \\ \Delta u(k+N_{\rm c} -1) \end{pmatrix}_{N_{\rm c} \times 1} \\ {\boldsymbol{Y}}& =( y(k+1), y(k+2), \cdots, y(k+N) )^{\rm T} \end{align*} $

即LTV系统式(20) 等效为增广系统式(24)。

现证明增广系统式(24) 关于原点是渐近稳定的。在$ k $时刻，对系统取李雅普诺夫函数：

$ V(\tilde{{\boldsymbol{x}}}(k), k) =\bar{J}(k) \notag \\ =\tilde{{\boldsymbol{x}}}^{\rm T}(k+N |k){\boldsymbol{P}}\tilde{\boldsymbol x}(k+N |k)+ \notag \\ \quad \sum\limits _{i=1}^{N-1} \tilde{{\boldsymbol{x}}}^{\rm T}(k+i| k ){\boldsymbol{Q}}\tilde{\boldsymbol x}(k+i| k) +\notag \\ \quad \sum\limits _{i=0}^{N-1} {R}(\Delta \tilde{u}^{*} (k+i|k))^{2} $

(26)

易知李雅普诺夫函数$V(\tilde{\boldsymbol x}(k), k)$具有正定性且随$\tilde{{\boldsymbol{x}}}(k)$增大趋于无穷。此时系统沿$ k $时刻的最优解向量为

$ \Delta \tilde{{\boldsymbol{U}}}_{k}^{*} =[\Delta \tilde{u}^{*} (k |k), \Delta \tilde{u}^{*} (k+1 |k), \cdots, \Delta \tilde{u}^{*} (k+N-1 |k)] \\ =[{\boldsymbol{L}}^{\rm T}(0)\boldsymbol \eta_{k}, {\boldsymbol{L}}^{\rm T}(1)\boldsymbol \eta_{k}, \cdots, {\boldsymbol{L}}^{\rm T}(N-1)\boldsymbol \eta_{k}] $

将$\Delta \tilde{{\boldsymbol{U}}}_{k}^{*}$沿采样时刻由$ k $平移至$ k+1$，并取$ k+N+1$时刻的控制量为0，则：

$ \Delta {\boldsymbol{\bar{U}}}_{k+1}^{*} =[{\boldsymbol{L}}^{\rm T}(1)\boldsymbol \eta_{k}, {\boldsymbol{L}}^{\rm T}(2)\boldsymbol\eta_{k}, \cdots, {\boldsymbol{L}}^{\rm T}(N-1)\boldsymbol\eta_{k}, 0] $

为$ k+1$时刻的一组可行解。将可行解$\Delta {\boldsymbol{\bar{U}}}_{k+1}^{*}$代入式(26)，得到类李雅普诺夫函数$\bar{V}(\tilde{{\boldsymbol{x}}}(k+1), k+1)$，且下式成立：

$ \begin{align*} V(\tilde{{\boldsymbol{x}}}(k+1), k+1)-V(\tilde{{\boldsymbol{x}}}(k), k) \\ \leqslant \bar{V}(\tilde{{\boldsymbol{x}}}(k+1), k+1)-V(\tilde{{\boldsymbol{x}}}(k), k) \end{align*} $

注意到$\bar{V}(\tilde{{\boldsymbol{x}}}(k+1), k+1)$与$V(\tilde{{\boldsymbol{x}}}(k), k)$在$ k+1$时刻后具有相同的控制输入，因此

$ \begin{align*} \bar{V}(\tilde{{\boldsymbol{x}}}(k+1), k+1)-V(\tilde{{\boldsymbol{x}}}(k), k) \\ =-\tilde{{\boldsymbol{x}}}^{\rm T}(k+1){\boldsymbol{Q}}\tilde{{\boldsymbol{x}}}(k+1)-\Delta \tilde{u}^{2}(k){R}<0 \end{align*} $

由李雅普诺夫定理可知，

$ \tilde{{\boldsymbol{x}}}(k+1)\to 0, \quad k\to \infty $

即闭环系统关于原点$(\tilde{x}=0, \tilde{y}=0)$渐近稳定。

5 实验结果（The experimental result）

集成串联弹性关节的下肢外骨骼机器人如图 12所示。外骨骼机器人具备单侧三关节主动自由度，以实现人体矢状面内髋—膝—踝关节运动。主控制器对下肢运动轨迹进行解耦，将关节间的耦合运动划分为单关节子任务，最后由驱动器实现底层控制。本文以外骨骼机器人的髋关节模块为对象，基于所提出的算法设计了控制实验，并对实验结果进行分析。实验中，各控制参数如表 2所示。以PD前馈控制律^[26]作为对比算法：

$ \begin{split} u =K_{\rm P, m} (\theta_{\rm d} -\theta_{\rm M})-K_{\rm D, m} (\dot{\theta}_{\rm d} -\dot{\theta}_{\rm M})+\tau_{\rm ff} \\ \tau_{\rm ff} =B \ddot{\theta}_{\rm d} +K(\theta_{\rm d} -q_{\rm d}) \end{split} $

(27)

其中$K_{\rm P, m}$、$K_{\rm D, m}$分别为PD前馈控制律中的比例调节增益与微分调节增益，$\tau_{\rm ff}$为前馈力矩项，$\theta_{\rm d}$、$\dot{\theta}_{\rm d}$分别为电机期望角度与期望速度。

为量化控制性能，定义如下评价指标：

(1) 平均绝对值误差为所有实际输出与期望输出的绝对值误差算术平均：

$ e_{\text{MAE}} =\frac{1}{N}\sum\limits _{t=1}^N | y_{{\text{actual}}, {t}} -y_{{\text{desired}}, {t}} | $

(28)

(2) 均方根误差为实际输出与期望输出的偏差平方相比于统计次数的算术平方根：

$ e_{\text{RMSE}} =\sqrt{\frac{1}{N}\sum\limits _{t=1}^N (y_{{\text{actual}}, {t}} -y_{{\text{desired}}, {t}})^{2}} $

(29)

(3) 控制能耗采用串联弹性关节中电机控制力矩的绝对平均值表征：

$ J_{\text{EST}} =\frac{1}{N}\sum\limits _{k=1}^N | u(k) | $

(30)

5.1 串联弹性关节姿态定位实验

外骨骼运动控制中对关节姿态具有定位需求，例如关节初始姿态复位、姿态保持等场景任务。以阶跃参考轨迹为实例，设计3组不同幅值的参考轨迹以模拟关节姿态定位任务。其中，幅值分别选取为0.1 rad、0.2 rad、0.3 rad，实验结果如图 13及表 3所示。

对于不同的信号幅值，本文方法都能够较快地（0.2 s内）收敛至目标值附近（稳态误差$<$ 5%），而对比算法存在更大稳态误差并具有更大的超调量。另一方面，本文方法对控制量实现了有效约束，而对比算法在位置跟踪误差较大时产生了远超约束的控制量（约3 N$\cdot $m）。一般地，外骨骼设备在姿态定位任务中由于复杂的人机交互环境存在过载或临界过载的情况。从控制算法设计角度对控制量进行约束可以主动降低系统负荷以及潜在的故障风险。通过评价指标式(3) 对控制过程中的能耗进行估计，本文方法相较于对比算法降低了约60% 能耗，从而提升了自主运动场景（内部电池驱动）下的续航能力。

5.2 串联弹性关节运动跟踪实验

如前述，本文算法验证中以独立的髋关节为研究对象。下肢步态运动中，不同步速下髋关节位置轨迹均呈现类正弦曲线。为验证串联弹性关节运动跟踪性能，以不同周期的正弦轨迹为参考位置，分别取正弦频率为0.5 Hz、0.75 Hz、1 Hz以对应不同步态速度，信号幅值为0.2 rad。实验结果如图 14及表 4所示，本文方法与对比算法相比具有更低的相位延迟，这是因为基于模型的控制方法有效提升了系统的控制带宽，并显著优于对比方法中前馈项对模型动力学的补偿。

2种方法在控制过程中均满足控制约束（电机输出力矩处于$[-1, 1]$ N$\cdot$m范围内），然而值得注意的是，PD前馈控制在位置方向改变时出现了控制量的尖峰，导致速度换向时关节运动存在抖动。人机交互中速度换向较为频繁，关节抖动现象将随关节数量增加而加剧（例如髋—膝双关节配置、髋—膝—踝三关节配置），在对控制器引入扰动的同时降低了交互过程中的舒适性与协调性。

5.3 外部扰动抑制实验

外骨骼设备在人机交互中存在外部环境及穿戴者的扰动，对系统提出了鲁棒性要求。在5.1节的基础上，设计了系统稳态下发生外部冲击扰动的实验。所设计的扰动发生时刻位于稳定的位置跟踪阶段，以模拟串联弹性关节进行位置定位和姿态保持时突然受到外部扰动的情景。扰动实验中，阶跃参考信号的幅值为0.2 rad，扰动发生时刻为0.6 s，设计扰动大小分别为参考信号幅值的15%、20%、30%，实验结果如图 15及表 5所示。

对比算法中缺乏对外部扰动的主动调节，但串联弹性关节具有被动柔顺性，在扰动结束后连杆位置仍能回到参考轨迹，因此在本实验中，对比算法退化为被动扰动抑制情形。本文方法在扰动发生后的0.1 s内回到参考轨迹，同时扰动发生前后控制量一致满足控制约束。在未进行主动扰动响应的对比算法中，连杆位置在扰动发生后的0.3 s内回到参考轨迹邻域，伴随小幅度的稳态误差，收敛时间随扰动增加而显著延长。实验结果表明本文方法对外部扰动具有抑制作用，验证了算法的鲁棒性。

对比算法中缺乏对外部扰动的主动调节，但串联弹性关节具有被动柔顺性，在扰动结束后连杆位置仍能回到参考轨迹，因此在本实验中，对比算法退化为被动扰动抑制情形。本文方法在扰动发生后的0.1 s内回到参考轨迹，同时扰动发生前后控制量一致满足控制约束。在未进行主动扰动响应的对比算法中，连杆位置在扰动发生后的0.3 s回到参考轨迹邻域，伴随小幅度的稳态误差，收敛时间随扰动增加而显著延长。实验结果表明本文方法对外部扰动具有抑制作用，验证了算法的鲁棒性。

5.4 系统带宽实验

所设计的串联弹性关节对系统引入弹性环节，进而对系统控制带宽产生影响。为计算系统在所提出算法下的实际控制带宽，设计了带宽实验。在迭代线性化模型预测控制中，控制内环为电机力矩控制，最外环为连杆位置控制。外骨骼设备在康复应用场景中的主要任务是对患者下肢的运动功能进行复现，即利用外骨骼引导患者下肢运动并实现对期望轨迹的跟踪，因此本文围绕系统位置带宽进行实验设计。

实验中以5 kg均匀质量块等效患者下肢大腿质量^[27]，并与串联弹性关节的连杆侧固定。设计正弦线性变频连杆位置信号$ r(t)$：

$ r(t)=A\sin \left(\phi_{0} +2\pi \left(f_{0} t+\frac{k}{2}t^{2}\right)\right) $

其中，$A$为幅值比例因子，$\phi_{0} =0$为初始相位，$ f_{0} =0.1$ Hz为初始频率，$ k=0.5$为正弦频率变化率。

实验对串联弹性关节进行时间长度为200 s的位置跟踪控制，并分别对期望信号和实际输出信号进行快速傅里叶变换以得到频域下的增益与相位。其中，采用零相位滤波器对所需数据进行了低通滤波以降低测量噪声。所得系统闭环Bode图如图 16所示。其中，幅值直流分量为$-$0.52 dB，系统控制带宽为11.12 Hz，此时对应相角为$-$75.34$^{\circ}$。

由于具有较高柔顺性，所提出的串联弹性关节难以获得刚性系统的高带宽，与SEA相关文献^[28-29]相比，本文方法实现了较好的带宽控制，能够覆盖人体下肢运动的频率范围^[30]。在康复外骨骼的临床应用中，下肢运动功能障碍患者相比于健康人群具有更低的运动频率。目前商用康复外骨骼频率为0.5~0.8 Hz^[31]，本文提出的控制方法能够满足康复外骨骼频率需求。

6 结论（Conclusion）

针对下肢康复外骨骼系统，提出了一种紧凑、轻量化的串联弹性关节，并设计了一种迭代线性化模型预测控制方法。该关节具有独有的刚度切换特性，在利用非线性刚度对抗过载形变以提供结构保护的同时，能够保证额定工作区间内具有一致的线性刚度。为实现控制器设计，本文首先对串联弹性关节进行了系统能量分析与动力学建模，并实现分步解耦动力学参数辨识以及对应的激励和验证实验；其次建立了最优化问题框架，设计了迭代线性化模型预测控制以求解带有控制量约束的非线性非凸优化问题。针对外骨骼任务场景，对姿态定位任务、髋关节步态跟踪任务设计了对比实验。实验结果表明，本文方法能够实现串联弹性关节的运动跟踪并对控制量进行有效约束，从而避免控制量出现尖峰与抖动，降低了系统控制能耗并抑制了外部扰动。此外，基于模型的控制方法有效提升了柔性关节的控制带宽以适应复杂人机交互场景。本文提出的串联弹性关节能够作为柔顺交互的关节模组，被拓展应用于不同自由度的外骨骼与相关柔顺交互控制任务中。