软体机器人作为新一代的机器人，在很多领域有着得天独厚的优势^[1]. 软体机器人具有大范围连续变形的能力，对环境有着极强的适应能力，近年来已经成为了机器人研究的热点^[2]. 柔性关节驱动器作为软体机器人的核心技术^[3-5]，已经被开发用于各种应用，包括软夹具、人工肌肉、可穿戴设备、触觉设备和医疗设备等^[6-7]，具有广泛的应用前景. 柔性关节驱动器通常采用柔性材料制造，由于柔性材料内部摩擦力、周期性松弛和强非线性等方面的原因，柔性关节驱动器具有非常复杂的输出特性^[8]. 目前国内外对柔性驱动器的迟滞特性研究还处于初步阶段，研究人员将压电陶瓷、形状记忆合金中的迟滞理论引入到柔性驱动器中，将其迟滞看作是一种扰动量进行研究，尚未有针对柔性驱动器迟滞特性的统一建模方法^[9].

根据超弹性材料理论和驱动器结构的分析，柔性关节驱动器的迟滞原因主要如下：1）硅橡胶的内摩擦（内黏性）；2）诱导结晶；3）柔性关节驱动器外界摩擦；4）纤维增强或约束情况下，纤维编织网与硅胶基体间的干摩擦^[6-7]. 目前的迟滞建模方法主要有唯象模型和物理模型^[10]2大类. 唯象模型是仅根据输入和输出之间的关系而建立的抽象模型，常用的唯象模型主要有Preisach模型、KP模型、PI模型等^[11]. 物理模型是从不同的迟滞物理机制出发，来具体描述迟滞过程. 物理模型的参数具有物理意义，因此物理模型更加直观精确，但通用性不足，往往只能描述特定的对象，且模型的参数往往需要与实际系统的各个参数相对应，在建模过程中须根据实际情况对模型参数进行改进和修正.

目前常用的物理迟滞模型主要有Jile-Atherton模型、Maxwell模型、Duhem模型、Bouc-Wen模型等^[12]. Jile-Atherton模型是最常见的物理迟滞模型之一，通常用于描述磁化过程中的磁滞特性，无法用于动态驱动力的描述^[13]. Maxwell模型是针对黏弹性材料提出的微积分模型，由于有记忆项，参数过多且不便被辨识，其复杂程度远高于不具备记忆功能的模型^[14]. Duhem模型是由Duhem提出的一种基于微分方程的动态模型，其输出与输入信号频率及方向相关，具有明确的函数表达式，由于其模型参数获取方式复杂，不易表征复杂的迟滞非线性^[15]. Bouc-Wen模型是一种微分模型，所描述的迟滞呈双曲线对称，已被广泛应用于迟滞建模. 此外，Bouc-Wen模型具有方向性，能够满足关节驱动器中关节伸展和收缩的物理描述需求^[16].

本研究选择Bouc-Wen作为柔性关节驱动器动态输出力的迟滞建模本体模型，其包含A、β、γ、n共4个参数，Bouc-Wen模型参数辨识易收敛，可以通过最小二乘法、遗传算法、粒子群优化等实现. 在迟滞实验的基础上，提出Bouc-Wen修正模型，并通过龙格-库塔-费尔贝格法进行参数辨别，能够将柔性关节驱动器的转角φ与输出力F之间的迟滞拟合精度显著提高. 对柔性关节驱动器的力迟滞进行准确建模，可以为力闭环控制提供基础，也为其他超弹性材料柔性驱动器的迟滞建模提供参考方法.

蜘蛛腿的运动主要是由腿节-膝节关节和胫节-后跗节关节（统称为液压关节）的驱动来实现，这2个关节处只有屈肌而没有伸肌^[17]. 液压关节由坚硬的外骨骼和柔软的折叠关节膜组成，外骨骼中的多个通道中充满蜘蛛的血淋巴，血淋巴即蜘蛛的血液. 血淋巴通过蜘蛛前体的肌肉加压，进而驱动柔软的关节膜产生膨胀，从而推动腿部伸展，如图1所示. 腿部屈曲则主要依赖于外骨骼内的屈肌，以及关节膜本身的弹性恢复^[19].

蜘蛛腿关节的驱动特征与现有的流体弹性驱动器相类似^[20]. 目前，关节式气动驱动器研究还相对较少，尤其是具有微小型结构的气动驱动器. 且大多数气动柔性驱动器为了简化分析和控制的难度，结构均设计为规则的形状，例如圆柱状、矩体状、尺蠖型、象鼻型等^[21]. 不规则形状的驱动器由于理论方面较复杂，研究相对较少. 本研究基于蜘蛛腿液压关节的结构和驱动原理，设计了一种新型的仿生柔性关节驱动器，呈不规则的扇形结构，运动原理与蜘蛛腿关节膜相似，并采用气动驱动方法，满足了质量轻、结构简单、性能优良等要求，其结构及原理如图2所示.

所设计的柔性关节驱动器是一种双连通式气囊弹性驱动器，考虑到气囊弹性驱动器体积较小、基体制备难度大，将驱动器设计为一体式气腔结构，具有不易破损、结构简单、性能稳定、制备相对容易等优势. 柔性关节驱动器气腔内部形状为一个具有2个腔室的扇形，2个腔室紧密连通，2个腔室之间的间隙侧壁厚度比柔性关节驱动器的其他侧壁厚度要小. 这种非对称式几何结构对驱动器产生局部约束、定向膨胀的作用，当驱动器空腔内的压力增加时，2个间隙侧壁膨胀从而相互挤压，实现驱动器转动. 其驱动原理如图2（b）所示.

柔性关节驱动器结构的主要参数如表1、图2（c）所示. 表中，b为两侧面扇形层的厚度；w为驱动器宽度；w₁为气囊腔室厚度，w₁=w−2b；r为驱动器内腔的曲率半径.

柔性关节驱动器制造材料为Dragon Skin^TM30硅橡胶. 为了提高柔性驱动器的驱动性能和抗压强度，采用一体化失蜡铸造方法进行制造，如图3所示，避免了由于局部受力不均匀而导致驱动器的破损，大大提高了驱动器使用寿命. 其主要步骤如下：1）3D打印软关节失蜡铸造模具；2）制作气腔蜡芯模具；3）硅胶模具浇筑出气腔蜡芯；4）将气腔蜡芯装入柔性关节驱动器外模具；5）浇筑柔性关节驱动器半成品；6）取下柔性关节驱动器半成品，放入水浴中加热融化蜡芯，得到柔性关节驱动器.

影响柔性关节驱动器迟滞的主要因素有材料特性^[6-7,10]、驱动器结构、充气压力等. 制造柔性关节驱动器的均为Dragon Skin30硅橡胶材料，其特性对迟滞的影响在此不进行讨论. 主要讨论的是充气压力对柔性关节驱动器迟滞的影响.

构建如图4所示的迟滞实验平台，柔性关节驱动器的充气由气泵完成，并通过精密调节阀（IR2000-02BG）进行压力调节. 采用3D打印方法制造的刚性支架将柔性关节驱动器底面固定，支架侧面同关节驱动器侧面相约束，随关节转动而转动. 在转动中心处，同轴安装有一个42步进电机，用以控制关节驱动器转动角度. 支架侧面连接有高精度力传感器，测量驱动器在不同转动角度时的输出力大小. 控制器为Arduino，控制步进电机转动，并对力传感器数据进行采集.

在进行迟滞实验之前，须获得驱动器在不同充气压力下的最大转动角度. 即将关节驱动器底面固定，在侧面无支架约束情况下，对于不同恒定充气压力分别测定驱动器的转动角度. 当驱动器充气压力小于60 kPa时，驱动器转动角度过小，测试误差较大，而当气压大于120 kPa时，由于材料特性及气密性不足，驱动器容易爆裂或漏气，因此将充气压力限定在60~120 kPa.

柔性关节驱动器的迟滞实验原理如图5所示. 首先，向驱动器通入恒定压力压缩空气，并控制步进电机以5˚/s逆时针恒速转动，如图5（b）所示，直至驱动器伸展至该气压下的最大转角处. 由于关节驱动器在步进电机带动下被动匀速伸展，将该过程称为回程. 然后，控制步进电机以5˚/s顺时针恒速转动，如图5（c）所示，压缩驱动器直至其初始状态（底面和侧面呈90˚），该过程称为去程. 为了保证实验结果的准确性，制造了5个柔性驱动器进行测试，取平均测试数据.

分别通入60、80、100、120 kPa的恒定气压进行测试，测得周期下转角与驱动力的关系，结果如图6（a）~（d）所示. 此外，为了更进一步分析驱动器的力特性，控制步进电机每5˚转动一次，每次停顿5 s，测得驱动器的静态输出力. 并将静态输出力与动态输出力进行比较，如图6（e）、（f）所示.

如图6（a）~（d）所示，柔性关节驱动器的输出力和转动角度之间存在明显的迟滞现象，具体表现在去程与回程时的力-位特征不同. 当驱动器为去程时，随着转角的减小，输出力随之增大，去程的各转角对应输出力均大于静态输出力；当驱动器为回程时，随着转角的增大，输出力随之减小，回程的各转角对应输出力在大部分情况下要小于静态输出力. 同时，不同气压下的迟滞曲线均不闭合，回程起点的输出力略微小于去程终点的输出力，这是硅橡胶松弛和惯性力的影响所导致的.

由图6（e）可以看出，静态测量的关节驱动器输出力F与转角φ之间大体呈现线性关系. 当驱动器角度φ一定时，关节的输出力F与充气压力p近似成正比关系，充入气体的压力越大，输出力也越大. 输出力与转角关系曲线的斜率可以看作是转动刚度，由图6（e）可知，在静态输出时，柔性关节驱动器的转动刚度与充气压力及转角均无关.

由图6（f）可以发现，在动态输出时，转动刚度随着驱动器转角变化而呈非线性变化. 在回程过程中，转动刚度减小速度会随着转角增大先快后慢，且驱动器整体刚度较静态时偏小；在去程过程中，转动刚度增大的速度会随着转角的减小而先快后慢，且与静态时相比整体刚度偏大.

由图6综合可得，在静态情况下，柔性关节驱动器的输出性能比较稳定，且可以看作是线性变化，而动态下的力输出与转动角度之间，则具有明显的迟滞效应. 驱动器内部气压越大，其输出力的迟滞效应越明显，在60、80、100 、120 kPa这4种充气压力情况下，只有60 kPa时的迟滞曲线才近乎对称.

柔性关节驱动器采用的是超弹性非线性硅橡胶材料，在运动过程中具有较强的时变性. 将运动过程中驱动器的输出力F( $\varphi $, p)分解为静态力F_s( $\varphi $, p)和非线性项F_n( $\varphi $, p)^[22]，即

(1) $ F\left(\varphi ,p\right)={F}_{\mathrm{s}}\left(\varphi ,p\right)+{F}_{\mathrm{n}}\left(\varphi ,p\right) .$

式中：F_s( $\varphi $, p)为驱动器静态输出力，F_n( $\varphi $, p)为驱动器运动过程中的迟滞力， $\varphi $为驱动器转角，p为驱动器内部的气体压力.

如图6所示，通过静态实验详细测得驱动器的静态输出力. 为了进一步描述运动过程中的动态输出力，引入Bouc-Wen 模型，将动态输出力z =F_n( $\varphi $, p)表示为非线性微分方程的形式^[22]：

(2) $ \dot{z}=A\Delta \dot{\varphi }-\beta \left|\Delta \dot{\varphi }\right|{\left|z\right|}^{n-1}z-\gamma \Delta \dot{\varphi }{\left|z\right|}^{n} . $

式中：z为柔性关节驱动器的输出力， $ \Delta \dot{\varphi } $为柔性关节驱动器的角速度，A为比例系数，β、γ分别为滞回环宽度调节参数和高度调节参数，n为曲线圆滑系数，

Bouc-Wen 模型是一种光滑双曲线迟滞模型，被广泛应用于非线性结构单元的迟滞建模. 但是通过经典的Bouc-Wen 模型所建立的迟滞曲线属于对称滞回曲线^[23]，而柔性关节驱动器的输出力-转角迟滞并不是严格的对称曲线，因此为了建立更精确的迟滞模型，须对经典Bouc-Wen 模型进行修正.

根据柔性驱动器的超弹性材料构成可知，造成驱动器动态输出曲线在去程和回程时非对称的一个很大因素是硅橡胶分子链间的内摩擦，而内摩擦力一般具有明确的方向性. 另外，在柔性关节驱动器伸展/收缩时，转动速度较快，囤积在硅橡胶内的热能无法完全消散，从而导致驱动器能量损失的不确定性，造成应变能储存与释放的非均匀性，引起柔性关节驱动器不同的输出力.

在Bouc-Wen模型的基础上，添加一个具有方向性的修正项 $\phi \left({\rm{sign}}\;(\mathrm{\Delta }\dot{\varphi }\right)+{\rm{sign}}\;\left(z\right))$，能够极大地克服滞回的不确定性因素，且精度较高、计算简单高效，便于参数的拟合. 修正后的Bouc-Wen模型可以表示为

(3) $ \begin{split} \dot{z}=&A\Delta \dot{\varphi }-\beta \left|\Delta \dot{\varphi }\right|{\left|z\right|}^{n-1}z-\gamma \Delta \dot{\varphi }{\left|z\right|}^{n}+ \\&\phi \left(\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\;\left(\mathrm{\Delta }\dot{\varphi }\right)+\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\;\left(z\right)\right). \end{split} $

式中： $\phi $为迟滞模型的回环系数.

采用龙格-库塔-费尔贝格法(Runge-Kutta-Fehlberg method，RKF45)对修正迟滞模型的参数进行辨识^[24]，主要过程是确定合适的步长h，其原理如下：在每一次计算中，比较Runge-Kutta法的4阶近似解和5阶较优解. 如果上述2个解的绝对差符合规定的精度，则采用4阶近似解，并按变步长方法增大步长h，加快计算速度；反之，则减小步长h. 每一次计算都须使用以下6个值^[25]：

(4) $ \left. \begin{aligned} &{k}_{1}=hf\left({t}_{k},{y}_{k}\right) , {k}_{2}=hf\left({t}_{k}+\frac{1}{4}h,\;{y}_{k}+\frac{1}{4}{k}_{1}\right),\\ &{k}_{3}=hf\left({t}_{k}+\frac{3}{8}h,{y}_{k}+\frac{3}{32}{k}_{1}+\frac{9}{32}{k}_{2}\right), \\ &{k}_{4}=hf\left({t}_{k}+\frac{12}{13}h,{y}_{k}+\frac{1\;932}{2\;197}{k}_{1}-\frac{7\;200}{2\;197}{k}_{2}+\frac{7\;296}{2\;197}{k}_{3}\right),\\ &{k}_{5}=hf\left({t}_{k}+h,{y}_{k}+\frac{439}{216}{k}_{1}-8{k}_{2}+\frac{3\;680}{513}{k}_{3}-\frac{845}{4\;104}{k}_{4}\right),\\ &{k}_{6}=hf\left({t}_{k}+\frac{1}{2}h,{y}_{k}-\frac{8}{27}{k}_{1}+2{k}_{2}-\right. \\ & \left.\frac{3\;544}{2\;565}{k}_{3}+\frac{1\;859}{4\;104}{k}_{4}-\frac{11}{40}{k}_{5}\right). \end{aligned}\right\}$

式中：k₁~k₆为系数，h为步长，t_k为柔性关节驱动器转角，y_k为柔性关节驱动器输出力.

使用4 阶Runge-Kutta法对非线性微分方程的解进行近似^[25]：

(5) $ {y}_{k+1}={y}_{k}+\frac{25}{216}{k}_{1}+\frac{1\;408}{2\;565}{k}_{3}+\frac{2\;197}{4\;104}{k}_{4}-\frac{1}{5}{k}_{5}. $

使用5阶Runge-Kutta法确定较优解^[25]：

(6) $ {z}_{k+1}={y}_{k}+\frac{16}{135}{k}_{1}+\frac{6\;656}{12\;825}{k}_{3}+\frac{28\;561}{56\;430}{k}_{4}-\frac{9}{50}{k}_{5}+\frac{2}{55}{k}_{6}. $

Runge-Kutta-Fehlberg 方法的误差^[25]，即4阶近似解与5阶较优解的绝对差为

(7) $ |{z}_{k+1}-{y}_{k+1}|=\left|\frac{1}{360}{k}_{1}-\frac{128}{4\;275}{k}_{3}-\frac{2\;197}{75\;240}{k}_{4}+\frac{1}{50}{k}_{5}\frac{2}{55}{k}_{6}\right|.$

下一次计算过程中的步长通过将标量s乘以当前步长h来确定^[25]. s的表达式为

(8) $ \mathrm{s}\le {\left(\frac{\varepsilon h}{2\left|{z}_{k+1}-{y}_{k+1}\right|}\right)}^{{1}/{4}}\approx 0.84{\left(\frac{\varepsilon h}{\left|{z}_{k+1}-{y}_{k+1}\right|}\right)}^{{1}/{4}}. $

式中：ε为指定的误差控制值，y_k+1及z_k+1分别为使用4阶、5阶Rung-Kutta方法所确定的近似解和较优解.

模型参数的辨识结果如表2所示.

进一步，采用多项式数值拟合方法，建立了各参数与气压之间的关系：

(9) $ \left.\begin{aligned} &A=8\;229{p}^{3}-1\;501{p}^{2}+17p+5.71 ,\\& \beta =-106\;833{p}^{3}+25\;950{p}^{2}-1\;757p+23 ,\\& n=-301\;812{p}^{3}+77\;048{p}^{2}-6\;307p+166, \\ &\gamma =-167\;166{p}^{3}+41\;525{p}^{2}-4\;114p+63 ,\\& \varphi =-49\;895{p}^{3}+13\;430{p}^{2}-1\;206p+38. \end{aligned}\right\}$

由此得到了动态输出力迟滞模型中的各参数.

为了验证Bouc-Wen修正模型的准确性，将驱动器在恒定充气压力为70、90、110 kPa下的实验数据与2种模型仿真数据进行对比，包括经典Bouc-Wen模型以及Bouc-Wen修正模型，结果如图7所示.

定义模型的最大相对误差：

(10) $ {\delta }_{\mathrm{M}\mathrm{A}\mathrm{X}}=\frac{\mathrm{\Delta }{F}_{\mathrm{M}\mathrm{A}\mathrm{X}}}{{F}_{\mathrm{M}\mathrm{A}\mathrm{X}}}\times 100{\text{%}} . $

式中：F_MAX为实验测量的输出力最大值， $ \Delta {F_{{\text{MAX}}}} $为模型与实验相比的最大输出力误差.

定义R²统计量来检验修正模型的拟合优度：

(11) $ R ^2 = \frac{{{\rm{SSR}}}}{{{\rm{SST}}}} = \frac{{\displaystyle\sum\limits_{i = 1}^n {{{\left( {{{\hat F}_i} - \bar F} \right)}^2}} }}{{\displaystyle\sum\limits_{i = 1}^n {{{\left( {{F_i} - \bar F} \right)}^2}} }} . $

式中：F_i为实验测得柔性关节驱动器的输出力， $ \bar{F} $为实验测得柔性关节驱动器输出力的平均值， ${{\hat F}_{i}}$为修正模型得到的柔性关节驱动器输出力. R²用来描述仿真模型与实验值的拟合程度，范围为[0, 1.0]，越接近1.0说明拟合优度越好.

如图7所示，相较于经典Bouc-Wen模型，所提出的Bouc-Wen修正模型，能够更为准确地描述柔性关节驱动器的迟滞特性. 由表3的3种充气压力情况下的最大误差∆F_MAX、最大相对误差δ_MAX、平均偏差A.D.、均方差MSE以及拟合优度R²都可以看出，与经典Bouc-Wen模型相比较，Bouc-Wen修正模型对于柔性关节驱动器迟滞实验的拟合结果有了明显提高. 在70、90、110 kPa的3种充气压力情况下，Bouc-Wen修正模型的迟滞拟合结果和实验数据相比，柔性关节驱动器去程输出力的最大偏差为0.847 N，去程输出力最大相对误差为7.75%；回程输出力最大偏差0.637 N，回程输出力最大相对误差为7.57%. 去程平均偏差均小于0.45 N，而回程平均偏差则小于0.30 N. 在均方差方面，Bouc-Wen修正模型对于迟滞曲线的拟合结果均小于0.42 N，较经典Bouc-Wen模型也有明显提升. 在模型的拟合优度方面，经典Bouc-Wen模型的R²均小于0.90，而Bouc-Wen修正模型的R²均大于0.99，显示了精确拟合效果.

此外，随着充气压力的增大，修正模型的迟滞偏差有增大的趋势，但均在较小误差范围内. 误差的来源除了实验误差外，最主要的是模型的参数个数和复杂程度，理论上来说，模型的参数越多，模型越复杂，其精度也越高，但过多的参数会导致模型复杂，计算可能不收敛. Bouc-Wen模型本身具有广泛适用性，其修正模型仅针对柔性关节驱动器迟滞曲线的非对称性而提出，可以为类似结构或材料的柔性驱动器迟滞建模研究提供借鉴.

设计了一款仿蜘蛛液压关节的柔性关节驱动器，针对其输出力的迟滞影响，基于Bouc-Wen修正方法对力迟滞进行建模. 分析柔性关节驱动器产生迟滞的原因，得到不同充气压力下柔性关节驱动器的迟滞特性. 采用实验建模与理论建模相结合的方法，提出Bouc-Wen修正模型，引入一个具有方向性的修正项 $\phi (\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\;\left(\mathrm{\Delta }\dot{\varphi }\right)+\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\;\left(z\right))$，克服滞回的非对称及不确定性，并采用龙格-库塔-费尔贝格方法进行模型参数辨识. 从最大误差、相对误差、平均偏差、均方差以及拟合优度这5个方面对Bouc-Wen修正模型的准确性进行验证，结果表明相较于经典Bouc-Wen模型，柔性关节驱动器的力迟滞建模精度有了极大提升.

下一步将开展柔性关节驱动器的动态迟滞特性研究，如阶跃气压作用下的输出力动态过程实验与建模分析. 并且建立基于力迟滞建模及力反馈的闭环控制系统，研究柔性关节之类的超弹性材料驱动器的控制特性.