0 引言
相比于传统的驱动器，气动人工肌肉关节通过气囊的收缩与膨胀执行指定动作，除具有气压传动技术低成本，清洁、安装简便的优点外，还具有高功率/质量比、自然柔顺性和与生物肌肉类似的力学特性，因而在康复理疗装置、仿生机器人等多个领域中得到了广泛使用，其中应用最广泛的是MicKibben型气动人工肌肉[1-2].为了让气动人工肌肉关节能够更好地在工程中得到应用，近年来国内外学者针对气动人工肌肉关节的静/动态特性建模、迟滞特性建模和控制策略等方面进行了大量研究[3-4].
气动人工肌肉关节具有高度非线性和时变性，从而难以获得其精确的动态模型.基于模型的方法通常对建模的精确度都具有较高的敏感性，相比之下，数据驱动控制方法直接使用在线数据，这些数据不包括任何被控对象的模型参数，因此能够有效解决控制不确定性问题.迭代反馈整定(IFT)作为一种无模型的数据驱动学习方法，由Hjalmarsson率先提出[5]，该方法预先给定控制器结构，设计LQG型性能准则函数，通过在闭环系统上迭代采集输入输出数据，计算性能准则函数对控制器参数的梯度，然后采用Gauss-Newton法自动搜索最小化性能准则函数对应的控制器参数. IFT算法的优点在于其给出了一种计算性能准则函数对控制器参数无偏梯度估计的无模型方法，克服了往常需要通过对被控对象和扰动动态模型的估计来求取这些微分信号的弊端.文[6]将工业中广泛使用的4种经典PID调优方法，包括Ziegler-Nichols规则、内模控制以及积分平方误差等方法应用于4个线性系统并进行了仿真测试，证明了IFT算法可以达到甚至超过经典PID控制器参数整定方法的性能，而且可以使系统具有更快的收敛速度和更小的超调.此外，经典PID控制器参数整定方法多适用与无噪声环境，与之相比IFT算法可以对输入信号进行滤波，从而有效减少噪声带来的干扰.基于以上这些优点，IFT算法已从线性系统的单回路控制逐渐发展到MIMO系统、级联控制和非线性系统等[7-8]，并从提高系统的暂态性能，加快算法收敛速度以及增强鲁棒性等方面不断发展完善.
针对气动人工肌肉充气变形过程中存在的强非线性、时变性和控制参数难以确定的问题，本文将迭代反馈整定(IFT)算法推广应用于非线性系统的PID控制器参数的整定和控制性能优化，目的是提高气动人工肌肉关节的跟踪控制性能.本文首先建立气动人工肌肉关节驱动的单自由度运动数学模型，然后在PID控制策略框架下，基于IFT算法的基本原理、参数最优化理论和实验整定方法，根据闭环系统的性能准则函数以及输入输出信号自动整定PID控制器参数，运用Gauss-Newton算法得到PID控制器参数的最优值，并引入辅助因子不断迭代性能准则函数的权重因子达到改善系统鲁棒性的目的，最后通过实验仿真验证了算法的有效性.
1 单自由度气动人工肌肉关节建模 1.1 气动人工肌肉数学模型MicKibben型气动人工肌肉由编织网和橡胶管组成，通过开关阀调节进气与排气流量来控制气动人工肌肉的伸缩，为此Tondu等利用虚功原理推导出气动人工肌肉关节的静态理论模型[9]：
|
(1) |
其中，p是人工肌肉内部的相对压力，F是人工肌肉的收缩力，l0和l分别是人工肌肉的初始长度和当前长度，r0、α0和ε分别表示人工肌肉运动的初始半径、编制网初始夹角和收缩率.可以看出，式(1)所示气动人工肌肉的收缩力与气压成正比，与收缩率成非线性关系.假设气动人工肌肉在运动过程中的静摩擦力f保持不变，则最终气动人工肌肉在收缩伸展时的输出力F(p)与内部气压p的关系为[10]
|
(2) |
同时根据理想气体的多变过程原理得出人工肌肉内部气体的状态方程为[11]
|
(3) |
其中，Rg为气体常数，T为绝对温度，V为人工肌肉的工作容积，mq为腔室内气体质量，其微分信号为进入腔室的气体质量流量qm.
假设气动人工肌肉的充放气为恒温过程，则工作容积V可表示为
|
(4) |
由式(4)的微分方程可以得到气压变化率
|
(5) |
在建立单个气动人工肌肉数学模型的基础上，气动人工肌肉关节的设计原理如图 1所示[12].该关节使用一对对抗性人工肌肉控制一个同步带轮以驱动关节旋转，使用两个PWM信号u1、u2驱动高速开关阀，进而调节进出一对人工肌肉关节内部的气体质量流量qm1和qm2，以及内部压强p1和p2，从而改变人工肌肉的长度，进而产生不同的牵引力来驱动同步带轮带动力臂旋转[13].信号u1控制气动人工肌肉1的进气与气动人工肌肉2的出气，u2相反，因此两个气动人工肌肉伸缩完全相反.因为关节正负运动对称，为减少控制量，不妨以控制信号u0为基准点，令u1=u0+Δu、u2=u0-Δu，则最终输入可以表示为PWM信号Δu，输出为关节旋转角度θ的数学模型.可以看出，式(3)建立了输出力F(p)与内部气压P的关系，式(6)则建立了气压变化率
|
| 图 1 单自由度气动人工肌肉关节设计原理图 Fig.1 Design principle of 1-DOF PMAs joint |
气体通过高速开关阀的过程通常是基于Sanville公式来描述的[14]，则两个气动人工肌肉关节的进气质量流量可以表示为
|
(6) |
|
(7) |
其中，kv是正常数，ρ是气体密度，ps是气源的绝对压力.定义p0是大气压力，进一步得到进气流量qin(p)为
|
(8) |
出气流量为qout(p)：
|
(9) |
令g1(p)=kvρpsqin(p)-kvρpqout(p)，g2(p)=kvρpsqin(p)+kvρpqout(p)，并代入u1=u0+Δu、u2=u0-Δu后，两个气动人工肌肉气体质量流量qm1和qm2与控制信号Δu的关系分别为
|
(10) |
|
(11) |
假设lb为关节杆的长度，lc为关节杆质心到转轴的长度，负载固定于关节杆末端，mb和M分别为关节杆和负载的质量，F(p1)、F(p2)为根据式(3)得到的两个气动人工肌肉在内部压强p1、p2下产生的牵引力，Cf为粘滞摩擦系数.
因而，基于Newton-Euler法可建立气动人工肌肉的动力学方程为
|
(12) |
|
(13) |
|
(14) |
联立式(1)~式(14)进一步可推导出气动人工肌肉关节的整体数学模型为
|
(15) |
将推导出的气动人工肌肉关节的数学模型(17)作为控制系统的未知不确定性被控对象G，通过获取实际角度与给定角度的跟踪误差，设计PID控制器用以输出PWM信号Δu，其中PID控制器参数通过IFT算法自动获取.根据IFT算法理论，假定PID控制器是一个具有通用参数的二自由度控制器[15]，控制系统的结构如图 2所示.图中C(ρ)=[Cr(ρ) Cy(ρ)]，Cr(ρ)、Cy(ρ)是线性时不变传递函数，r(t)是参考输入，u(t)是控制器输出，y(t)是系统输出，v(t)是外部随机扰动，PID控制器参数ρ=[kP kI kD].
|
| 图 2 气动肌肉关节的控制系统 Fig.2 Control system of PMAs joint |
在图 2所示反馈控制方案作用下的系统响应为
|
(16) |
为便于描述，下文将T0(ρ)和S0(ρ)简写为T0和S0.
接下来，定义yd是系统的期望轨迹信号，则实际输出与期望输出之间的误差为
|
(17) |
对于受参数ρ直接影响的固定结构PID控制器，通过
|
(18) |
式中，
|
(19) |
|
(20) |
接下来，在闭环控制系统上进行3次响应实验求出
|
(21) |
式中，ri(1)、ri(2)和ri(3)为第i次迭代下第一次、第二次和第三次的实验输入，vi(1)、vi(2)和vi(3)是三次实验过程中相互独立的外部干扰，因而
|
(22) |
|
(23) |
进而得到基于实验数据的第i次迭代的性能指标J(ρi)的估计梯度为
|
(24) |
根据J(ρi)的估计梯度以及上一次迭代的控制参数ρi，使用Gauss-Newton算法计算下一次迭代更新值：
|
(25) |
其中，γi>0表示步长，Ri为正定Hessian矩阵表示优化搜索方向，表示为
|
(26) |
为IFT算法选择的步长需要注意把握优化算法的收敛性和系统稳定性之间的平衡，步长越小且迭代次数足够多将更有利于确保算法的收敛性，在实践中步长通常根据经验定义约束步长为0 < γi≤1[18].
另需指出的是，性能指标J(ρi)的定义决定了IFT算法的性能，通常情况下J(ρi)中每个元素的权重因子都是由经验得出的.文[19]表明λ在不同区间下IFT算法有着不同的收敛效果，为此考虑引入辅助因子Li使得λ在不同系统中都具有可比性和随机性，该方法将λ整定到一个合适的区间，从而使得IFT算法能够加速收敛.此时将式(18)和式(26)修改为
|
(27) |
|
(28) |
其中：
|
(29) |
式中，yd，max和yd，min为期望输出的最大值与最小值，ut，max和ut，min表示第i次迭代过程中控制信号在所有N个采样点中的最大值与最小值.由于这些值都是在每次迭代结束时给出，所有的采样点都会被考虑在内，因而引入辅助因子使得λ在不同系统中都代表了当前迭代优化后
步1：选择控制器最初的参数ρ1并根据式(27)设计性能指标J(ρi).
步2：根据式(21)进行三次实验，实验的三次输入分别为：ri(1)=r、ri(2)=y(1)(ρi)和ri(3)=r；将第一次实验所得y(1)(ρi)作为第二次实验的参考输入，第二次和第三次实验获取控制器输入值u(2)(ρi)和u(3)(ρi)，以及系统的输出值y(2)(ρi)和y(3)(ρi).
步3：运用第二、三次实验的结果根据式(22)和式(23)计算估计梯度
步4：利用Hessian阵Ri并选取步长，根据式(25)获取新的控制器参数ρi+1.
步5：判断系统性能指标J(ρi)是否达到要求，则转步6结束，否则重新转步2.
步6：结束.
2.3 迭代反馈整定算法的收敛性为保证算法的收敛性，一方面应尽量保证目标函数的估计梯度是无偏的，另一方面步长序列γi要求能够收敛到零[20].进一步由式(19)~式(21)可得估计梯度为
|
(30) |
|
(31) |
在IFT实验中设定三次实验的vi(m)是同一系统相互独立的零均值有界随机噪声，即|vi(m)| < C，m=1，2，3，其中假设三次实验中界限值C和随机噪声的均方值均保持不变，则可以保证式(22)和式(23)的估计是无偏的.
另一方面需要保证算法收敛的条件通常要求序列γi的所有元素满足[21]：
|
(32) |
该收敛条件的基本要求是在整个优化迭代过程中输入信号r(t)始终保持有界.虽然确定更新方向的矩阵Ri不会影响IFT的收敛能力，但理想选择是通过Gauss-Newton方向来加速收敛速度.因此正如Huusom等人指出，使用Gauss-Newton优化算法可以保证算法的收敛性，使得设计的IFT算法能够很快覆盖到一个固定的优化点[22].该结论除了时不变性条件外，对系统的性质没有任何假设，因此该结论适用于简单的PID控制器或者更复杂的控制器.
3 仿真结果为验证IFT优化方法的有效性，本文在Matlab环境下针对单自由度气动人工肌肉关节模型进行控制仿真，模型参数如表 1中所示.假定关节臂摆动的期望轨迹是振幅为0.2 rad、频率为0.15 Hz的正弦曲线，它在工作过程中受到均值为零的随机白噪声vi影响.
| 参数 | 数值 |
| r0/mm | 10 |
| α0/(°) | 68 |
| l0/m | 0.2 |
| T/K | 298 |
| Cf/(kg·m2/s) | 0.025 |
| Rg /(J·kg-1·K-1) | 287.15 |
| p0/MPa | 0.101 |
| kv/(m4·s·kg-1) | 5.1×10-9 |
| g/(m·s-2) | 9.8 |
| ps/MPa | 0.7 |
| u0 | 0.5 |
| lb/m | 0.5 |
| lc/m | 0.25 |
| mb/kg | 0.5 |
| M/kg | 2.45 |
| Rj/mm | 25 |
| f/N | 40 |
为与本文控制优化算法进行比对，首先使用常用的Ziegler-Nichols参数整定规则获得PID增益为ρ=[3.35 5.5×10-2 2.25×10-3]，将其作为ρ1，在该PID控制器作用下气动人工肌肉关节的角跟踪效果及其误差采样点频数分布直方图如图 3(a)所示.为进一步提升输出跟踪控制效果，在此基础上采用IFT算法进行迭代优化，其中参数Kp、KI以及KD的整定步长分别为0.2、1×10-3和1×10-4，选取λ=1×10-2，经过20次PID控制器参数自动迭代优化后气动人工肌肉的关节角跟踪控制效果和误差采样点频数分布直方图如图 3(b)所示.
|
| 图 3 Ziegler-Nichols法与IFT法下跟踪效果曲线和误差频数直方图 Fig.3 The tracking curve and histogram of control error between Ziegler-Nichols and IFT |
由图 3可以得到，当关节臂驱动2.45 kg的负载，以振幅为0.2 rad、频率为0.15 Hz的轨迹摆动时，经过20次迭代后，累计跟踪误差的平方和从5.29×10-2 rad2减小到了1.78×10-2 rad2，跟踪效果得到了明显改善.控制过程中相对应的性能准则函数J(ρ)与PID控制器参数KP、KI以及KD的变化趋势如图 4所示.
|
| 图 4 迭代过程中性能指标J(ρ)与PID控制器参数的变化趋势 Fig.4 The trend of performance criterion function J(ρ) and controller parameter |
为了进一步验证结论的适用性，改变仿真实验条件以求逼近于工程实际要求，本文接下来通过改变负载大小并驱动关节臂在不同的轨迹下运动，从而得到多组对照实验结果，对比不同条件下系统在15 s内的累计跟踪误差平方，部分仿真实验结果如表 2所示.
| 负载/kg & 摆动轨迹 | 累计误差平方 | 减小百分比/% | |
| Z-N /×10-2 rad2 | IFT /×10-2 rad2 | ||
| 2.45 & 0.2 rad、频率为0.15 Hz的正弦曲线 | 7.41 | 2.57 | 65.32 |
| 4.45 & 0.2 rad、频率为0.15 Hz的正弦曲线 | 8.23 | 5.05 | 33.89 |
| 2.45 & 0.2 rad、频率为0.20 Hz的正弦曲线 | 17.15 | 4.27 | 75.10 |
| 2.45 & 0.2 rad、频率为0.25 Hz的正弦曲线 | 27.86 | 5.57 | 80.01 |
由四组对比实验可以看出，随着关节臂负载的重量增加，IFT算法的优化效果有所降低，这是由于随着气动肌肉关节输出牵引力的增加，相对应的高速开关阀门PWM调节信号Δu的幅度由±0.21变为±0.15，并出现了部分震荡，从而降低了控制器输出ut(ρ)无偏梯度的计算精度.同时也可以看出，在摆动频率变大，系统跟踪变化速度更快的轨迹时，IFT算法相比较Ziegler-Nichols法拥有更强的适应性.
另外需指出的是，在性能准则函数J(ρ)中，合适的
|
| 图 5 对比实验中性能指标J(ρ)与PID控制器参数的变化趋势 Fig.5 The trend of performance criterion function J(ρ) and controller parameter in comparison experiments |
由图 5可见，性能准则函数J(ρ)随着迭代过程的进行稳步下降，并且在引入辅助因子Li后，J(ρ)进一步减小，可以很好地适应第10次迭代时关节臂摆动频率的突变情况，表明系统的整体控制性能和鲁棒性得到了进一步提升.与经典IFT算法相比，引入辅助因子Li后可以得到一个更佳的参数λ，便于获取更好的优化方向，最终通过不断调节控制器参数补偿关节臂扭矩在摆动过程中的变化.
4 结论本文在建立非线性气动人工肌肉数学模型的基础上，直接基于系统的输入输出数据，提出了一种通过迭代反馈整定算法自动调节PID控制器参数的方法，然后定义包含跟踪误差及加权控制信号的性能准则函数，采用Gauss-Newton算法提高迭代反馈整定算法的收敛速度，同时引入辅助因子使得权重因子在不同系统中都具有适应性，进一步改善了系统的跟踪控制精度和鲁棒性.仿真实验表明，通过IFT算法自动参数整定后的PID控制器能够很好得完成单自由度人工肌肉关节的轨迹跟踪任务.但也需指出的是，IFT算法虽然相比较于Ziegler-Nichols等传统PID参数整定方法简便，可以实现自动调参，但也增加了算法的计算复杂度，这无疑限制了IFT算法的工程适用性.另外由于实验条件限制，本文尚未在实际人工肌肉装置上进行控制实验，开展实物仿真验证进一步验证本文算法的有效性将是今后需努力开展的一项工作.
| [1] |
谢胜龙, 梅江平, 刘海涛.
McKibben型气动人工肌肉研究进展与趋势[J]. 计算机集成制造系统, 2018, 24(5): 1065–1080.
Xie S L, Mei J P, Liu H T. Achievements and trends of research on McKibben pneumatic artificial muscles[J]. Computer Integrated Manufacturing Systems, 2018, 24(5): 1065–1080. |
| [2] |
朱雪峰, 王建辉, 王晓峰.
基于人工免疫算法的滚动优化在康复机器人中的应用研究[J]. 信息与控制, 2017, 46(3): 335–341.
Zhu X F, Wang J H, Wang X F. Application of the rolling optimization based on artificial immune algorithm to rehabilitation robot[J]. Information and Control, 2017, 46(3): 335–341. |
| [3] | Tondu B. Modelling of the McKibben artificial muscle:A review[J]. Journal of Intelligent Material Systems and Structures, 2012, 23(3): 225–253. DOI:10.1177/1045389X11435435 |
| [4] | Kelasidi E, Andrikopoulos G, Nikolakopoulos G, et al. A survey on pneumatic muscle actuators modeling[C]//Proceedings of 2011 IEEE International Symposium on Industrial Electronics. Piscataway, NJ, USA: IEEE, 2011: 1263-1269. https://www.researchgate.net/publication/259310435_A_survey_on_pneumatic_muscle_actuators_modeling |
| [5] | Hjalmarsson H, Gevers M, Gunnarsson S, et al. Iterative feedback tuning:theory and applications[J]. IEEE Control Systems Magazine, 1998, 18(4): 26–41. DOI:10.1109/37.710876 |
| [6] | Lequin O, Gevers M, Mossberg M, et al. Iterative feedback tuning of PID parameters:Comparison with classical tuning rules[J]. Control Engineering Practice, 2003, 11(9): 1023–1033. DOI:10.1016/S0967-0661(02)00303-9 |
| [7] | Radac M B, Precup R E, Preitl S, et al. Iterative feedback tuning in MIMO systems. Signal processing and application[C]//5th International Symposium on Applied Computational Intelligence and Informatics. Piscataway, NJ, USA: IEEE, 2009: 77-82. https://www.researchgate.net/publication/221172814_Iterative_Feedback_Tuning_in_MIMO_Systems_Signal_Processing_and_Application |
| [8] | Tesch D, Eckhard D, Bazanella A S. Iterative feedback tuning for cascade systems[C]//2016 European Control Conference (ECC). Piscataway, NJ, USA: IEEE, 2016: 495-500. https://www.researchgate.net/publication/312326497_Iterative_feedback_tuning_for_cascade_systems |
| [9] | Tondu B, Ippolito S, Guiochet J, et al. A seven-degrees-of-freedom robot-arm driven by pneumatic artificial muscles for humanoid robots[J]. The International Journal of Robotics Research, 2005, 24(4): 257–274. |
| [10] |
张宏立, 申珉珉, 彭光正.
气动人工肌肉静态数学模型与实验研究[J]. 液压与气动, 2009(4): 17–19.
Zhang H L, Shen M M, Peng G Z. Static mathematical model and experimental study of pneumatic muscle actuator[J]. Chinese Hydraulics & Pneumatics, 2009(4): 17–19. DOI:10.3969/j.issn.1000-4858.2009.04.008 |
| [11] | Minh T V, Tjahjowidodo T, Ramon H, et al. Cascade position control of a single pneumatic artificial muscle-mass system with hysteresis compensation[J]. Mechatronics, 2010, 20(3): 402–414. DOI:10.1016/j.mechatronics.2010.03.001 |
| [12] |
谢建蔚, 陶国良, 周洪.
气动人工肌肉关节的建模与仿真[J]. 浙江大学学报(工学版), 2007(3): 450–455.
Xie J W, Tao G L, Zhou H. Modelling and simulation of joint actuated by pneumatic muscles[J]. Journal of Zhejiang University (Engineering Science), 2007(3): 450–455. DOI:10.3785/j.issn.1008-973X.2007.03.017 |
| [13] | Najjari B, Barakati S M, Mohammadi A, et al. Modelling and controller design of electro-pneumatic actuator based on PWM[J]. IAES International Journal of Robotics and Automation, 2012, 1(3): 125. |
| [14] | Erkaya S, Uzmay I. Experimental investigation of joint clearance effects on the dynamics of a slider-crank mechanism[J]. Multibody System Dynamics, 2010, 24(1): 81–102. |
| [15] | Pan H, Zhang Y, Wang L. Iterative feedback tuning for two-degree-of-freedom system[M]. Berlin, Germany: Springer, 2018: 365-379. |
| [16] | Kora K, Lu C Z, McDaid A J. Automatic tuning with feedforward compensation of the HuREx rehabilitation system[C]//2014 IEEE/ASME International Conference. Piscataway, NJ, USA: IEEE, 2014: 1504-1509. https://www.researchgate.net/publication/271484085_Automatic_tuning_with_feedforward_compensation_of_the_HuREx_rehabilitation_system |
| [17] | Meng W, Xie S Q, Liu Q, et al. Robust iterative feedback tuning control of a compliant rehabilitation robot for repetitive ankle training[J]. IEEE-ASME Transactions on Mechatronics, 2017, 22(22): 173–184. |
| [18] | Heertjes M F, Van der Velden B, Oomen T. Constrained iterative feedback tuning for robust control of a wafer stage system[J]. IEEE Transactions on Control Systems Technology, 2016, 24(1): 56–66. DOI:10.1109/TCST.2015.2418311 |
| [19] | Lu C Z, Xie S Q, Deng C. Optimal normalized weighting factor in iterative feedback tuning of step input responses[J]. IFAC Proceedings Volumes, 2017, 22(22): 173–184. |
| [20] | Rădac M B, Precup R E, Petriu E M, et al. Data-driven reference trajectory tracking algorithm and experimental validation[J]. IEEE Transactions on Industrial Informatics, 2013, 9(4): 2327–2336. DOI:10.1109/TII.2012.2220973 |
| [21] | Huusom J K, Poulsen N K, Jørgensen S B. Improving convergence of iterative feedback tuning[J]. Journal of Process Control, 2009, 19(4): 570–578. DOI:10.1016/j.jprocont.2008.09.004 |
| [22] | Hjalmarsson H. Iterative feedback tuning-an overview[J]. International Journal of Adaptive Control and Signal Processing, 2002, 16(5): 373–395. DOI:10.1002/acs.714 |