Identification and compensation of friction for modular joints based on grey wolf optimizer

基于灰狼算法的模块化关节摩擦辨识和补偿 Identification and compensation of friction for modular joints based on grey wolf optimizer 1 2 first-author 崔靖凯（1997-），男，山东聊城人，博士研究生，2019年于中国海洋大学获得学士学位，主要从事超冗余机器人、机器人振动抑制等方面的研究。E-mail： cuijingkai19@mails.ucas.ac.cn http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008330&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008336&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008333&type= 22.01332855 32.00399399 崔靖凯 （1997-），男，山东聊城人，博士研究生，2019年于中国海洋大学获得学士学位，主要从事超冗余机器人、机器人振动抑制等方面的研究。E-mail： cuijingkai19@mails.ucas.ac.cn cuijingkai19@mails.ucas.ac.cn 1 2 1 1 1 3 corresp E-mail： xuzhenbang@ciomp.ac.cn E-mail： xuzhenbang@ciomp.ac.cn 徐振邦（1982-），男，内蒙古通辽人，研究员，2005年、2010年于中国科学技术大学分别获得学士、博士学位，主要从事空间智能机器人、空间遥感器设计、振动控制方面的研究。E-mail： xuzhenbang@ciomp.ac.cn http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008338&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008342&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008340&type= 22.01332855 32.00399399 徐振邦 （1982-），男，内蒙古通辽人，研究员，2005年、2010年于中国科学技术大学分别获得学士、博士学位，主要从事空间智能机器人、空间遥感器设计、振动控制方面的研究。E-mail： xuzhenbang@ciomp.ac.cn xuzhenbang@ciomp.ac.cn 中国科学院 长春光学精密机械与物理研究所，吉林 长春130033 中国科学院 长春光学精密机械与物理研究所，吉林 长春130033 Changchun Institute of Optics， Fine Mechanics and Physics，Chinese Academy of Sciences， Changchun 130033， China Changchun Institute of Optics， Fine Mechanics and Physics，Chinese Academy of Sciences， Changchun 130033， China 中国科学院大学，北京 100049 中国科学院大学，北京 100049 University of Chinese Academy of Sciences， Beijing 100049， China University of Chinese Academy of Sciences， Beijing 100049， China 中国科学院大学 材料与光电研究中心，北京 100049 中国科学院大学 材料与光电研究中心，北京 100049 Center of Materials Science and Optoelectronics Engineering， University of Chinese Academy of Sciences， Beijing 100049， China Center of Materials Science and Optoelectronics Engineering， University of Chinese Academy of Sciences， Beijing 100049， China 针对模块化关节，设计了一套离线辨识方法，对关节的摩擦模型进行了辨识，并对摩擦力进行了补偿。首先，介绍了该模块化关节的结构和控制系统，建立了关节的动力学模型。然后，基于LuGre摩擦模型对关节摩擦力进行了建模，分别用灰狼算法和分段最小二乘配合伪随机序列法辨识了模型的参数，并进行对比分析。最后，设计了基于LuGre摩擦模型的前馈补偿算法，并进行了实验验证。实验结果表明：相较于分段最小二乘法，灰狼算法的辨识精度提高了19.2%；给定幅值为1 （°）/s，频率为10 Hz的正弦速度信号，摩擦补偿使关节速度跟踪误差从0.295 （°）/s减少为0.183 （°）/s；关节速度环带宽从12 Hz提高到18 Hz。多次实验显示，辨识的数据具有较高的重复性，验证了辨识方法的有效性。设计的前馈摩擦补偿算法可以改善关节控制系统的动态性能。 To identify the friction model parameters of a modular joint， an off-line identification method that compensates the joint friction is proposed. First， the structure and control system of the modular joint are presented， and the dynamic model of the joint is established. Second， the LuGre friction model is developed. The grey wolf algorithm and piecewise least-square algorithm with a pseudo random sequence are then used to identify the respective model parameters. The results of two methods are compared and analyzed， and a feed-forward compensation algorithm based on the LuGre friction model is designed and verified experimentally. The experimental results indicate that compared with the piecewise least-square method， the identification accuracy of the grey wolf algorithm improved by 19.2%； the joint velocity tracking error decreased from 0.295 （°）/s to 0.183 （°）/s when the given velocity signal was a sine wave with an amplitude of 1 （°）/s and a frequency of 10 Hz； and the velocity loop bandwidth increased from 12 Hz to 18 Hz after friction compensation. Several experiments are repeated， and the identified data exhibit a high repeatability， which verifies the suitability of the proposed method. The proposed feed-forward friction compensation algorithm can be used to improve the dynamic performance of the joint control system. 模块化关节 参数辨识 LuGre摩擦模型 灰狼算法 摩擦补偿 modular joint parameter identification LuGre friction model grey wolf optimizer friction compensation 1 引　言 1 模块化关节能通过重构组合，使机械臂形成不同功能的构型，满足不同环境和任务的需求。具有模块化关节的机械臂制造成本低、设计周期短、互换性强并易于改进，被广泛应用于各个领域 ［ 1 1 ］ 。为了获得高负载自重比，模块化关节常具有谐波减速器等非线性环节，这增大了摩擦干扰，影响机械臂的定位和跟踪精度。 摩擦补偿是减少或消除关节非线性环节的不利影响、提高系统动态性能的重要途经。摩擦补偿包括基于摩擦模型的补偿和不基于摩擦模型的补偿。其中，基于摩擦模型的补偿由于具有可预见、响应速度快等特点 ［ 2 2 ］ ，获得了广泛应用。模块化关节的摩擦模型包括库伦-黏性摩擦模型 ［ 3 3 ］ 、Stribeck摩擦模型 ［ 4 4 ］ 、LuGre摩擦模型 ［ 5 5 ］ 等。库仑-黏性摩擦模型和Stribeck摩擦模型无法描绘电机速度接近于零时的预滑动区域的摩擦特性，导致电机速度换向时补偿的摩擦力不连续，引起力矩的跳变。为了充分考虑摩擦的非线性特性，保证补偿力矩的连续性，本文基于LuGre摩擦模型进行摩擦补偿。 摩擦模型参数的准确度直接影响了摩擦模型的精度，进而影响摩擦补偿的效果 ［ 6 6 ］ ，因此需要对模型参数进行辨识。参数辨识是指在满足激励的条件下，采集系统的状态信息对物理参数进行估计的过程 ［ 7 7 ］ 。参数辨识可以分为在线辨识和离线辨识。在线辨识适用于对环境变化敏感的参数，如温度的变化、材料的磨损等，虽然具有高精度和实时响应的优点，但可能存在参数收敛不完全的问题，且受到实验对象能否实时执行的限制。离线辨识相对简单易行，在对环境变化不敏感的参数如力矩系数、摩擦系数等上更具优势，因此，本文采用离线辨识的方法。 针对LuGre模型，学者们提出了多种辨识方法。肖前进 ［ 8 8 ］ 等用分段最小二乘法辨识了LuGre模型的静态参数，将电机预滑动区域的特征作线性化处理，等价为一个二阶系统，计算得到动态参数。该方法需要对实验数据分段辨识，预滑动区域的线性化处理会引入建模误差、计算误差等，影响辨识精度。为了提高辨识精度，元启发式算法成为研究热点。于伟 ［ 9 9 ］ 等用一般辨识方法辨识了静态参数，并以转速变化为目标函数，用遗传算法拟合出动态参数；曾德林 ［ 10 10 ］ 等用分步夹逼和分步搜索法初步辨识了摩擦参数，然后用遗传算法提高辨识精度。 灰狼算法是由SEYEDALI M ［ 11 11 ］ 提出的元启发式优化算法，模仿自然界中灰狼的领导阶层和捕猎机制。灰狼算法具有简单易用、灵活性强和避免局部最优等优点，适用于多参数辨识的情况。相较于工程常用的最小二乘法，无需对数据分段处理，在辨识非线性模型方面具有优势。 本文针对自主研发的模块化关节，选用能较为全面地描述摩擦特性的LuGre模型，分别用基于分段最小二乘法的一般辨识方法和基于灰狼算法的元启发式辨识方法对摩擦参数进行辨识，并将辨识结果对比分析，最后通过前馈补偿算法对摩擦力进行补偿。实验结果表明，相比于一般辨识方法，灰狼算法在摩擦辨识领域精度更高、适用性更强；本文提出的摩擦补偿方法能显著减小关节速度跟踪误差，提高控制品质。 2 模块化关节系统及建模 1 2.1　关节结构 2 图1 图1 为模块化关节的内部构造示意图。关节主要由外壳、无框力矩电机、谐波减速器、电磁制动器、增量式编码器、绝对式编码器和驱动电路板构成。左侧和下侧的电连接板用于给相邻关节供电通讯，方便关节的装卸和组合，体现模块化的特点。关节采用双编码器设计，增量式编码器获取的相对位置信息通过差分运算得到速度信息，绝对式编码器用于获取关节负载端的绝对位置信息。 2.2　控制系统 2 关节采用位置-速度-电流三闭环控制，控制框图如 图2 图2 所示。计算机使用倍福公司的TwinCAT软件对关节发送指令，关节的驱动电路板接收到位置或速度指令，通过控制 q 轴电流驱动电机转子和输入轴转动，谐波减速器的波发生器与输入轴固连，经减速后的输出轴与负载连接，实现力矩的输出。 2.3　动力学模型 2 建立模块化关节的简化动力学模型如 图3 图3 所示，其动力学方程为： <math display="block" id="M1"><mi>J</mi><mover accent="true"><mrow><mi>θ</mi></mrow><mi>¨</mi></mover><mo>=</mo><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi></mrow></msub><mo>-</mo><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi><mi>f</mi></mrow></msub><mo>-</mo><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi><mi>l</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008360&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008358&type= 28.70199966 4.48733330 ， <math display="block" id="M2"><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi></mrow></msub><mo>=</mo><msub><mrow><mi>K</mi></mrow><mrow><mi>t</mi></mrow></msub><msub><mrow><mi>i</mi></mrow><mrow><mi>q</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008364&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008362&type= 14.13933277 3.72533321 ， <math display="block" id="M3"><mi>τ</mi><mo>=</mo><mi>M</mi><mo stretchy="false">(</mo><mi>q</mi><mo stretchy="false">)</mo><mover accent="true"><mrow><mi>q</mi></mrow><mi>¨</mi></mover><mo>+</mo><mi>C</mi><mo stretchy="false">(</mo><mi>q</mi><mo>,</mo><mover accent="true"><mi>q</mi><mo>˙</mo></mover><mo stretchy="false">)</mo><mover accent="true"><mi>q</mi><mo>˙</mo></mover><mo>+</mo><mi>G</mi><mo stretchy="false">(</mo><mi>q</mi><mo stretchy="false">)</mo><mo>+</mo><msub><mrow><mi>τ</mi></mrow><mrow><mi>f</mi></mrow></msub><mo>+</mo><msub><mrow><mi>τ</mi></mrow><mrow><mi>l</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008368&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008366&type= 64.09266663 4.82600021 ， <math display="block" id="M4"><mi>τ</mi><mo>=</mo><mi>η</mi><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008372&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008370&type= 11.43000031 4.57200003 ， <math display="block" id="M5"><msub><mrow><mi>τ</mi></mrow><mrow><mi>f</mi></mrow></msub><mo>=</mo><mi>η</mi><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi><mi>f</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008375&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008374&type= 12.78466606 4.57200003 ， <math display="block" id="M6"><msub><mrow><mi>τ</mi></mrow><mrow><mi>l</mi></mrow></msub><mo>=</mo><mi>η</mi><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi><mi>l</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008379&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008377&type= 12.86933327 4.57200003 ， 其中： <math id="M7"><mi>J</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008383&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008381&type= 1.60866666 2.79399991 为电机转子及波发生器折算到电机端的转动惯量； <math id="M8"><mi>θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008601&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008599&type= 1.77800000 2.79399991 为电机转角； <math id="M9"><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008391&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008389&type= 3.21733332 3.72533321 为电机端输出力矩， <math id="M10"><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi><mi>f</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008395&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008393&type= 3.80999994 3.72533321 和 <math id="M11"><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi><mi>l</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008399&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008397&type= 3.80999994 3.72533321 分别为折算到电机端的摩擦力矩和负载力矩； <math id="M12"><msub><mrow><mi>K</mi></mrow><mrow><mi>t</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008403&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008401&type= 3.55599999 3.72533321 为力矩系数； <math id="M13"><msub><mrow><mi>i</mi></mrow><mrow><mi>q</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008407&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008405&type= 2.28600001 3.72533321 为 q 轴电流； <math id="M14"><mi>τ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008411&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008409&type= 1.60866666 2.79399991 为关节端输出力矩； <math id="M15"><mi>q</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008414&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008412&type= 1.86266661 3.64066648 ， <math id="M16"><mover accent="true"><mi>q</mi><mo>˙</mo></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008418&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008416&type= 1.86266661 3.64066648 ， <math id="M17"><mover accent="true"><mrow><mi>q</mi></mrow><mi>¨</mi></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008422&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008420&type= 1.86266661 3.64066648 分别为关节端角度、角速度、角加速度； <math id="M18"><mi>M</mi><mo stretchy="false">(</mo><mi>q</mi><mo stretchy="false">)</mo><mover accent="true"><mrow><mi>q</mi></mrow><mi>¨</mi></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008426&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008424&type= 11.59933376 3.97933316 为连杆惯性项； <math id="M19"><mi>C</mi><mo stretchy="false">(</mo><mi>q</mi><mo>,</mo><mover accent="true"><mi>q</mi><mo>˙</mo></mover><mo stretchy="false">)</mo><mover accent="true"><mi>q</mi><mo>˙</mo></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008429&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008428&type= 14.13933277 3.97933316 为离心力和科氏力项； <math id="M20"><mi>G</mi><mo stretchy="false">(</mo><mi>q</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008520&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008518&type= 8.38199997 3.97933316 为重力项； <math id="M21"><msub><mrow><mi>τ</mi></mrow><mrow><mi>f</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008437&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008435&type= 2.28600001 3.72533321 和 <math id="M22"><msub><mrow><mi>τ</mi></mrow><mrow><mi>l</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008441&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008439&type= 2.37066650 3.72533321 分别为关节端摩擦力矩和负载力矩； <math id="M23"><mi>η</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008445&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008443&type= 1.86266661 3.64066648 为谐波减速器的减速比。为了统一度量，本文均采用关节端的物理量作为实验和计算数据。 3 摩擦参数辨识 1 3.1　LuGre摩擦模型 2 摩擦来源于接触器件的相对运动。对于模块化关节，谐波减速器的波发生器和柔轮的相对运动、柔轮和钢轮的相对运动都会产生摩擦。此外，电机输入轴、输出轴的支撑轴承也存在摩擦。 LuGre摩擦模型为： <math display="block" id="M24"><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi><mi>f</mi></mrow></msub><mo>=</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mi>z</mi><mo>+</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mover accent="true"><mi>z</mi><mo>˙</mo></mover><mo>+</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mover accent="true"><mi>θ</mi><mo>˙</mo></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008455&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008453&type= 33.10466385 4.57200003 ， <math display="block" id="M25"><mover accent="true"><mi>z</mi><mo>˙</mo></mover><mo>=</mo><mover accent="true"><mi>θ</mi><mo>˙</mo></mover><mo>-</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mfrac><mrow><mfenced open="|" close="|" separators="|"><mrow><mover accent="true"><mi>θ</mi><mo>˙</mo></mover></mrow></mfenced></mrow><mrow><mi>g</mi><mo stretchy="false">(</mo><mover accent="true"><mi>θ</mi><mo>˙</mo></mover><mo stretchy="false">)</mo></mrow></mfrac><mi>z</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008459&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008457&type= 26.92399979 11.34533310 ， <math display="block" id="M26"><mi>g</mi><mo stretchy="false">(</mo><mover accent="true"><mi>θ</mi><mo>˙</mo></mover><mo stretchy="false">)</mo><mo>=</mo><msub><mrow><mi>F</mi></mrow><mrow><mi>c</mi></mrow></msub><mo>+</mo><mo stretchy="false">(</mo><msub><mrow><mi>F</mi></mrow><mrow><mi>s</mi></mrow></msub><mo>-</mo><msub><mrow><mi>F</mi></mrow><mrow><mi>c</mi></mrow></msub><mo stretchy="false">)</mo><msup><mrow><mi>e</mi></mrow><mrow><mo>-</mo><msup><mrow><mfenced separators="|"><mrow><mfenced open="|" close="|" separators="|"><mrow><mover accent="true"><mi>θ</mi><mo>˙</mo></mover></mrow></mfenced><mo>/</mo><msub><mrow><mi>ω</mi></mrow><mrow><mi>s</mi></mrow></msub></mrow></mfenced></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008463&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008461&type= 47.15933609 6.85799980 ， 其中： z 为鬃毛平均形变量， <math id="M27"><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008468&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008465&type= 2.87866688 3.72533321 为鬃毛刚度， <math id="M28"><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008471&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008469&type= 2.87866688 3.72533321 为鬃毛阻尼系数， <math id="M29"><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008475&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008473&type= 2.87866688 3.72533321 为黏滞摩擦系数， <math id="M30"><mi>g</mi><mo stretchy="false">(</mo><mover accent="true"><mi>θ</mi><mo>˙</mo></mover><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008479&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008477&type= 7.45066643 4.40266657 为Stribeck效应， <math id="M31"><msub><mrow><mi>F</mi></mrow><mrow><mi>c</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008494&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008492&type= 3.47133350 3.72533321 为库仑摩擦力， <math id="M32"><msub><mrow><mi>F</mi></mrow><mrow><mi>s</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008486&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008484&type= 3.38666677 3.72533321 为最大静摩擦力， <math id="M33"><msub><mrow><mi>ω</mi></mrow><mrow><mi>s</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008490&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008488&type= 3.38666677 3.72533321 为Stribeck速度。模型包括四个静态参数 <math id="M34"><msub><mrow><mi>F</mi></mrow><mrow><mi>c</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008494&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008492&type= 3.47133350 3.72533321 、 <math id="M35"><msub><mrow><mi>F</mi></mrow><mrow><mi>s</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008498&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008496&type= 3.38666677 3.72533321 、 <math id="M36"><msub><mrow><mi>ω</mi></mrow><mrow><mi>s</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008502&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008500&type= 3.38666677 3.72533321 和 <math id="M37"><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008506&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008504&type= 2.87866688 3.72533321 ，以及两个动态参数 <math id="M38"><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008624&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008622&type= 2.87866688 3.72533321 和 <math id="M39"><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008513&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008511&type= 2.87866688 3.72533321 。 3.2　基于分段最小二乘法的一般辨识方法 2 一般辨识方法分为两步，第一步将实验数据分段，用最小二乘法分别辨识电机转速为正和转速为负时的摩擦静态参数；第二步将转速为零附近的预滑动区域作线性化处理，使其近似为二阶系统，然后辨识该系统的幅频特性，计算得到动态参数。 3 . 2 . 1　静态参数的辨识 3 单个关节水平放置，可以忽略离心力和科氏力项 <math id="M40"><mi>C</mi><mo stretchy="false">(</mo><mi>q</mi><mo>,</mo><mover accent="true"><mi>q</mi><mo>˙</mo></mover><mo stretchy="false">)</mo><mover accent="true"><mi>q</mi><mo>˙</mo></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008516&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008515&type= 14.13933277 3.97933316 以及重力项 <math id="M41"><mi>G</mi><mo stretchy="false">(</mo><mi>q</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008520&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008518&type= 8.38199997 3.97933316 的影响，再令电机匀速转动，由于加速度为0，惯性项 <math id="M42"><mi>M</mi><mo stretchy="false">(</mo><mi>q</mi><mo stretchy="false">)</mo><mover accent="true"><mrow><mi>q</mi></mrow><mi>¨</mi></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008524&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008522&type= 11.59933376 3.97933316 置零。将关节空载，根据 式（3） 式（3） ，得到： <math display="block" id="M43"><mi>τ</mi><mo>=</mo><msub><mrow><mi>τ</mi></mrow><mrow><mi>f</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008528&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008526&type= 8.72066593 3.72533321 . 此时关节的输出力矩主要克服摩擦力矩。 电机匀速转动即 <math id="M44"><mover accent="true"><mi>z</mi><mo>˙</mo></mover><mo>=</mo><mn mathvariant="normal">0</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008531&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008530&type= 8.38199997 2.79399991 时，稳态摩擦力为： <math display="block" id="M45"><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi><mi>f</mi></mrow></msub><mo>=</mo><mi>g</mi><mo stretchy="false">(</mo><mover accent="true"><mi>θ</mi><mo>˙</mo></mover><mo stretchy="false">)</mo><mi mathvariant="normal">s</mi><mi mathvariant="normal">g</mi><mi mathvariant="normal">n</mi><mo stretchy="false">(</mo><mover accent="true"><mi>θ</mi><mo>˙</mo></mover><mo stretchy="false">)</mo><mo>+</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mover accent="true"><mi>θ</mi><mo>˙</mo></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008535&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008533&type= 37.16867065 5.33400011 . 结合 式（1） 式（1） 和（2），在电机空载匀速转动的条件下，有： <math display="block" id="M46"><mtable><mtr><mtd><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi></mrow></msub><mo>=</mo><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi><mi>f</mi></mrow></msub><mo>=</mo><msub><mrow><mi>K</mi></mrow><mrow><mi>t</mi></mrow></msub><msub><mrow><mi>i</mi></mrow><mrow><mi>q</mi></mrow></msub><mo>=</mo></mtd></mtr><mtr><mtd><mo stretchy="false">[</mo><msub><mrow><mi>F</mi></mrow><mrow><mi>c</mi></mrow></msub><mo>+</mo><mo stretchy="false">(</mo><msub><mrow><mi>F</mi></mrow><mrow><mi>s</mi></mrow></msub><mo>-</mo><msub><mrow><mi>F</mi></mrow><mrow><mi>c</mi></mrow></msub><mo stretchy="false">)</mo><msup><mrow><mi>e</mi></mrow><mrow><mo>-</mo><msup><mrow><mfenced separators="|"><mrow><mfenced open="|" close="|" separators="|"><mrow><mover accent="true"><mi>θ</mi><mo>˙</mo></mover></mrow></mfenced><mo>/</mo><msub><mrow><mi>ω</mi></mrow><mrow><mi>s</mi></mrow></msub></mrow></mfenced></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mrow></msup><mo stretchy="false">]</mo><mi mathvariant="normal">s</mi><mi mathvariant="normal">g</mi><mi mathvariant="normal">n</mi><mo stretchy="false">(</mo><mover accent="true"><mi>θ</mi><mo>˙</mo></mover><mo stretchy="false">)</mo><mo>+</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mover accent="true"><mi>θ</mi><mo>˙</mo></mover></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008539&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008537&type= 60.11333084 11.43000031 . 因此，可以测出电机在不同速度下匀速转动时的 q 轴电流，乘以力矩系数 K t = 0.21 N·m/A得到摩擦力矩，绘制摩擦力-速度曲线。 图4 图4 是摩擦参数辨识和补偿平台。关节空载且水平放置，计算机的TwinCAT软件用于控制电机匀速转动，MATLAB软件用于数据的处理和分析。计算机和关节驱动板用网线连接，通信协议为倍福公司的EtherCAT。实验时，采样周期设置为5 ms，关节端速度变化范围为 - 35 （°）/s到 - 0.5 （°）/s和0.5 （°）/s到35 （°）/s，每次改变0.5 （°）/s，每个速度维持15 s。为了消除电机在启动和停止阶段的电流波动，取中间5 s的 q 轴电流值的平均值作为该速度下的电流值。根据 式（12） 式（12） 计算出每个速度对应的摩擦力，用分段最小二乘法拟合摩擦力-速度曲线，结果如 图5 图5 所示。 1;2; 按照 式（12） 式（12） 辨识的静态摩擦参数如 表1 表1 所示。 模型参数 速度区间/（°/s） 辨识结果 平均值 2 F c /（N·m） ［-35，-0.5］ 14.300 6 2 14.409 8 ［0.5，35］ 14.519 0 2 F s /（N·m） ［-35，-0.5］ 7.174 1 2 7.200 9 ［0.5，35］ 7.227 7 2 ω s /（°/s） ［-35，-0.5］ 11.868 2 2 12.410 5 ［0.5，35］ 12.952 7 2 σ 2 /（N·m·s/°） ［-35，-0.5］ 0.400 2 2 0.398 8 ［0.5，35］ 0.397 3 3 . 2 . 2　动态参数的辨识 3 在电机的预滑动区域，即驱动力矩小于最大静摩擦力矩时，对方程（7）在零状态和零速度附近作线性化处理，即令 <math id="M47"><mi>z</mi><mo>=</mo><mi>θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008567&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008565&type= 8.29733276 2.79399991 ， <math id="M48"><mover accent="true"><mi>z</mi><mo>˙</mo></mover><mo>=</mo><mover accent="true"><mi>θ</mi><mo>˙</mo></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008571&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008569&type= 8.29733276 3.55599999 ，得到： <math display="block" id="M49"><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi><mi>f</mi></mrow></msub><mo>=</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mi>θ</mi><mo>+</mo><mo stretchy="false">(</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>+</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo stretchy="false">)</mo><mover accent="true"><mi>θ</mi><mo>˙</mo></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008575&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008573&type= 33.86666870 4.91066647 . 关节空载即 <math id="M50"><msub><mrow><mi>τ</mi></mrow><mrow><mi>m</mi><mi>l</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">0</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008579&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008577&type= 10.58333302 3.72533321 时，结合 式（1） 式（1） 和（2），得到： <math display="block" id="M51"><mi>J</mi><mover accent="true"><mrow><mi>θ</mi></mrow><mi>¨</mi></mover><mo>=</mo><msub><mrow><mi>K</mi></mrow><mrow><mi>t</mi></mrow></msub><msub><mrow><mi>i</mi></mrow><mrow><mi>q</mi></mrow></msub><mo>-</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mi>θ</mi><mo>-</mo><mo stretchy="false">(</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>+</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo stretchy="false">)</mo><mover accent="true"><mi>θ</mi><mo>˙</mo></mover></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008583&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008581&type= 44.19599915 4.91066647 . 式（14） 式（14） 在频域中的表达可写成： <math display="block" id="M52"><mfrac><mrow><mi>θ</mi><mo stretchy="false">(</mo><mi>j</mi><mi>ω</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mrow><mi>i</mi></mrow><mrow><mi>q</mi></mrow></msub><mo stretchy="false">(</mo><mi>j</mi><mi>ω</mi><mo stretchy="false">)</mo></mrow></mfrac><mo>=</mo><mi>H</mi><mo stretchy="false">(</mo><mi>j</mi><mi>ω</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msub><mrow><mi>K</mi></mrow><mrow><mi>t</mi></mrow></msub></mrow><mrow><mo>-</mo><mi>J</mi><msup><mrow><mi>ω</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><mo>+</mo><mo stretchy="false">(</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>+</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo stretchy="false">)</mo><mi>j</mi><mi>ω</mi><mo>+</mo><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008587&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008585&type= 70.86599731 9.48266602 . 由此可见，在预滑动区域，LuGre模型表现为带有阻尼的质量-弹簧二阶系统。可以通过伪随机序列测试该二阶系统的幅频特性，分析伯德图和拟合的频响函数获得动态参数 <math id="M53"><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008624&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008622&type= 2.87866688 3.72533321 和 <math id="M54"><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008628&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008626&type= 2.87866688 3.72533321 的估计值。 实验时，以幅值为0.19 A，长度为511位，周期为2 s的M序列作为 q 轴电流值输入关节系统。该序列满足准白噪声要求，其自相关函数可近似为 <math id="M55"><mi>δ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008598&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008596&type= 1.94733346 2.79399991 函数。序列幅值的均方根值对应的关节驱动力矩为6.65 N·m，小于最大静摩擦力7.200 9 N·m，可防止系统离开预滑动区域。测量的系统响应为电机的转角 <math id="M56"><mi>θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008601&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008599&type= 1.77800000 2.79399991 。为了排除非零初始条件的影响，输入4个周期的M序列，测出后3个周期的系统响应并做平均处理。 图6 图6 展示了其中一个周期的M序列和系统响应曲线。 用相关分析法将系统辨识为二阶系统，得到频响函数 <math id="M57"><mi>H</mi><mo stretchy="false">(</mo><mi>j</mi><mi>ω</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008612&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008609&type= 10.83733368 3.97933316 。在 <math id="M58"><mi>ω</mi><mo>=</mo><mn mathvariant="normal">0</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008616&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008614&type= 9.22866726 2.79399991 时，由 式（15） 式（15） 得到： <math display="block" id="M59"><mi>H</mi><mo stretchy="false">(</mo><mn mathvariant="normal">0</mn><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msub><mrow><mi>K</mi></mrow><mrow><mi>t</mi></mrow></msub></mrow><mrow><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008621&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008618&type= 17.52599907 8.80533314 . 因此可以得到 <math id="M60"><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008624&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008622&type= 2.87866688 3.72533321 的估计值。再根据 式（15） 式（15） 计算 <math id="M61"><msub><mrow><mi>σ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008628&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008626&type= 2.87866688 3.72533321 的估计值，辨识结果如 表2 表2 所示。 模型参数 辨识结果 σ 0 /（N·m/°） 1 160 σ 1 /（N·m·s/°） 0.44 3.3　基于灰狼算法的元启发式辨识方法 2 灰狼算法（GWO）是一种新的元启发式算法，灵感来自灰狼社会的狩猎行为。GWO的种群代表一群狼，根据灰狼的适应度将其分为4个阶级：最优解称为α狼，次优解称为β狼，再次优解称为γ狼，其余解称为ω狼。 α ， β 和 γ 是搜索猎物的领导者， ω 狼通过以下算子来更新它们的位置，以跟随三只头狼到达搜索空间的最佳区域： <math display="block" id="M62"><mtable columnalign="left"><mtr><mtd><msub><mrow><mi>D</mi></mrow><mrow><mi>α</mi></mrow></msub><mo>=</mo><mo stretchy="false">|</mo><msub><mrow><mi mathvariant="bold-italic">C</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>·</mo><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mi>α</mi></mrow></msub><mo>-</mo><mi mathvariant="bold-italic">X</mi><mo stretchy="false">|</mo></mtd></mtr><mtr><mtd><msub><mrow><mi>D</mi></mrow><mrow><mi>β</mi></mrow></msub><mo>=</mo><mo stretchy="false">|</mo><msub><mrow><mi mathvariant="bold-italic">C</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo>·</mo><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mi>β</mi></mrow></msub><mo>-</mo><mi mathvariant="bold-italic">X</mi><mo stretchy="false">|</mo></mtd></mtr><mtr><mtd><msub><mrow><mi>D</mi></mrow><mrow><mi>γ</mi></mrow></msub><mo>=</mo><mo stretchy="false">|</mo><msub><mrow><mi mathvariant="bold-italic">C</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub><mo>·</mo><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mi>γ</mi></mrow></msub><mo>-</mo><mi mathvariant="bold-italic">X</mi><mo stretchy="false">|</mo></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008632&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008631&type= 26.58533096 15.07066536 ， <math display="block" id="M63"><mtable columnalign="left"><mtr><mtd><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>=</mo><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mi>α</mi></mrow></msub><mo>-</mo><msub><mrow><mi mathvariant="bold-italic">A</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>·</mo><msub><mrow><mi>D</mi></mrow><mrow><mi>α</mi></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo>=</mo><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mi>β</mi></mrow></msub><mo>-</mo><msub><mrow><mi mathvariant="bold-italic">A</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo>·</mo><msub><mrow><mi>D</mi></mrow><mrow><mi>β</mi></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub><mo>=</mo><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mi>γ</mi></mrow></msub><mo>-</mo><msub><mrow><mi mathvariant="bold-italic">A</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub><mo>·</mo><msub><mrow><mi>D</mi></mrow><mrow><mi>γ</mi></mrow></msub></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008634&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008633&type= 26.67000008 12.95400047 ， <math display="block" id="M64"><mi mathvariant="bold-italic">X</mi><mo stretchy="false">(</mo><mi>t</mi><mo>+</mo><mn mathvariant="normal">1</mn><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>+</mo><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo>+</mo><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub></mrow><mrow><mn mathvariant="normal">3</mn></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008636&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008635&type= 40.72467041 7.95866632 ， 其中： D α ， D β 和 D γ 分别是狼群 ω 与头狼 α ， β 和 γ 的距离； X ， X α ， X β 和 X γ 分别是狼群 ω 、头狼 α ， β 和 γ 当前迭代次数下的位置向量； X （ t+ 1）是狼群 ω 下一次迭代的位置向量， t 是当前迭代次数； A i 和 C i （ i =1，2，3）是随机生成的系数向量，如下所示： <math display="block" id="M65"><msub><mrow><mi mathvariant="bold-italic">A</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">2</mn><mi mathvariant="bold-italic">a</mi><mo>·</mo><msub><mrow><mi mathvariant="bold-italic">r</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>-</mo><mi mathvariant="bold-italic">a</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008639&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008637&type= 22.77533340 3.72533321 ， <math display="block" id="M66"><msub><mrow><mi mathvariant="bold-italic">C</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">2</mn><mo>·</mo><msub><mrow><mi mathvariant="bold-italic">r</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008641&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008640&type= 13.80066681 3.72533321 ， 其中： r 1 和 r 2 是区间［0，1］中的随机向量； a 的值在迭代过程中按照下式从2线性减小到0： <math display="block" id="M67"><mi mathvariant="bold-italic">a</mi><mo>=</mo><mn mathvariant="normal">2</mn><mo>·</mo><mfenced separators="|"><mrow><mn mathvariant="normal">1</mn><mo>-</mo><mfrac><mrow><mi>t</mi></mrow><mrow><mi>M</mi><mi>a</mi><mi>x</mi><mi>I</mi><mi>t</mi><mi>e</mi><mi>r</mi></mrow></mfrac></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008645&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008643&type= 33.18933105 9.22866726 ， 其中 MaxIter 是最大迭代次数。 图7 图7 展示了狼群如何在二维搜索空间内更新其位置。可以观察到，最终位置将位于由搜索空间中头狼α、β和γ所定义的圆内的随机位置。换句话说，α、β和γ估计猎物的位置，其他狼随机更新至猎物周围的位置。随机向量 A i 和 C i 是该算法需要调节的两个参数。| A i | < 1时，灰狼会向猎物进攻，而| A i | > 1时，灰狼会被迫逃离猎物，这允许灰狼进行全局搜索。| C i |始终是区间［0，2］中的随机值，它将提供随机权重，以强化（| C i | > 1）或弱化（| C i | < 1）猎物在 式（17） 式（17） 中定义距离的效果。 C i 也可视为自然界中障碍物对灰狼接近猎物的影响，它有助于灰狼突破局部最优停滞。 在LuGre摩擦模型的辨识中，该算法所处理的猎物是 式（7） 式（7） 到（9）所包含的6个参数。定义目标函数为： <math display="block" id="M68"><mi>J</mi><mo>=</mo><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><mi>N</mi></mrow></mfrac><mstyle displaystyle="true"><munderover><mo largeop="true">∑</mo><mrow><mi>i</mi><mo>=</mo><mn mathvariant="normal">1</mn></mrow><mrow><mi>N</mi></mrow></munderover></mstyle><mrow><mo stretchy="false">(</mo><msub><mrow><mi>f</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>-</mo><msub><mrow><mi>f</mi></mrow><mrow><mi>d</mi><mi>i</mi></mrow></msub><msup><mrow><mo stretchy="false">)</mo></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mrow></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008656&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008654&type= 30.14133072 8.63599968 ， 其中： N 是样本数， f i 和 f di 分别是第 i 次采样的模拟摩擦力和实际摩擦力。该问题是最小化寻优问题，即目标函数越小，灰狼的适应度越高。 图8 图8 展示了灰狼算法的操作流程。实验时，参数约束设置为： F c ∈［0，20］， F s ∈［0，10］， ω s ∈［0，20］， σ 2 ∈［0，1］， σ 0 ∈［0，30000］， σ 1 ∈［0，1］，狼群规模设置为50。因为算法是随机搜索的，故用不同的初始种群进行了10次重复运算，并在200次迭代后停止运行。 图9 图9 展示了在不同的初始种群下，α狼（最优解）的目标函数随迭代次数的演变过程。 可以看出，十次重复运行的结果都成功地收敛到最优参数附近，最佳运行的第38次迭代和最差运行的第158次迭代的目标函数都达到了10 -0.84 。一些运行在搜索过程中陷入了局部最优，但在最后的几次迭代中都成功地突破停滞，走向全局最优。 基于灰狼算法的辨识结果如 表3 表3 所示。 模型参数 辨识结果 F c /（N·m） 14.548 1 F s /（N·m） 7.286 8 ω s /（°/s） 12.577 7 σ 2 /（N·m·s/°） 0.394 2 σ 0 /（N·m/°） 29 931 σ 1 /（N·m·s/°） 0.68 3.4　两种算法的对比 2 选择 式（23） 式（23） 作为辨识精度的衡量指标。 J 代表辨识结果和实验数据之间的误差， J 越小，辨识精度越高。基于分段最小二乘法的一般辨识方法和基于灰狼算法的元启发式辨识方法的对比如 表4 表4 所示。由此可见，灰狼算法的辨识精度高于一般辨识方法，且无需对实验数据分段辨识，无需对预滑动区域作线性化处理，避免了建模误差和计算误差，在摩擦辨识领域具有更强的竞争力。 辨识算法 辨识误差 使用条件 最小二乘法+预滑动区域线性化 0.154 53 最小二乘法要求函数连续可微，只适应于速度同向的辨识，需要分段辨识；将预滑动区域线性化来辨识动态参数会引入建模误差和计算误差 灰狼算法 0.124 71 对辨识的数学模型无特殊要求，使用范围广，适用性强 4 摩擦补偿及实验 1 为验证摩擦辨识的准确性，改善模块化关节的动态性能，设计了 图10 图10 所示的基于摩擦模型的前馈补偿算法。关节仍采用位置-速度-电流三闭环控制，取期望的电机端速度值为摩擦补偿器的输入，按照 式（7） 式（7） 到（9）计算相应的摩擦补偿值 <math id="M69"><msub><mrow><mi>τ</mi></mrow><mrow><mi>c</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008689&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008685&type= 2.53999996 3.72533321 ，除以力矩系数得到补偿电流值 <math id="M70"><msub><mrow><mi>I</mi></mrow><mrow><mi>q</mi><mi>c</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008696&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008692&type= 3.13266683 3.72533321 ，作为前馈量叠加到电流环的输入端，以消除摩擦扰动。 关节输出角度的跟踪精度很大程度上取决于速度的跟踪精度，因此取关节的角速度作为观测量。实验采用 图4 图4 所示的设备，摩擦补偿前，测得模块化关节的速度环带宽为12 Hz，实验时，给定幅值为1 （°）/s，频率为10 Hz的正弦关节端速度输入信号sin（20 <math id="M71"><mi mathvariant="normal">π</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008716&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=24008712&type= 1.94733346 2.79399991 t ），加入摩擦补偿前后的关节端速度曲线、摩擦补偿值曲线、速度跟踪误差曲线如 图11 图11 所示。 图11 图11 中，在箭头所标注的时刻开启摩擦补偿器，可以看到摩擦补偿后，实线代表的速度反馈值能更好地跟踪虚线代表的速度期望值，跟踪误差显著减小，由补偿前的峰值0.295 （°）/s减小到补偿后的峰值0.183 （°）/s，降低了38%。改变输入正弦波的频率，测得摩擦补偿后，关节速度环带宽为18 Hz，比补偿前提高了50%。由此可见，建立的摩擦补偿算法效果显著，能提高速度跟踪精度和速度环控制带宽，同时验证了摩擦模型辨识的准确性。 5 结　论 1 本文针对模块化关节，设计了摩擦模型的离线辨识方法和摩擦补偿算法。分别使用基于分段最小二乘法的一般辨识方法和基于灰狼算法的元启发式辨识方法对LuGre摩擦模型进行了参数辨识，并对辨识结果进行对比分析。设计前馈摩擦补偿算法，以期望速度作为补偿器的输入，以补偿的电流值作为输出前馈给电流闭环。在模块化关节上进行实验，结果表明，相比于一般辨识方法，灰狼算法的辨识精度提高了19.2%；给定幅值为1 （°）/s，频率为10 Hz的速度正弦输入信号，摩擦补偿后，关节速度跟踪误差减小了38%；摩擦补偿后，速度环带宽提高了50%。所提出的参数辨识方法准确有效，摩擦补偿算法可以显著提高关节速度跟踪精度和速度环带宽。 参考文献 DOYLE M J ， MARQUES J V A ， VANDERMEULEN I ， et al . Modular fluidic propulsion robots ［J］. IEEE Transactions on Robotics ， 2020 ， 37 （ 2 ）： 532 - 549 . doi: 10.1109/tro.2020.3031880 http://dx.doi.org/10.1109/tro.2020.3031880 Modular fluidic propulsion robots CHU Z Y ， CHEN G ， CUI J ， et al . Classifier-based approximator for friction compensation in high accelerated positioning system ［J］. IEEE Transactions on Industrial Electronics ， 2021 ， 68 （ 5 ）： 4090 - 4098 . doi: 10.1109/tie.2020.2987268 http://dx.doi.org/10.1109/tie.2020.2987268 Classifier-based approximator for friction compensation in high accelerated positioning system ZHANG D W ， OU H G . Relationship between friction parameters in a Coulomb-Tresca friction model for bulk metal forming ［J］. Tribology International ， 2016 ， 95 ： 13 - 18 . doi: 10.1016/j.triboint.2015.10.030 http://dx.doi.org/10.1016/j.triboint.2015.10.030 Relationship between friction parameters in a Coulomb-Tresca friction model for bulk metal forming CONG S ， DENG K ， SHANG W W ， et al . Isolation control for inertially stabilized platform based on nonlinear friction compensation ［J］. Nonlinear Dynamics ， 2016 ， 84 （ 3 ）： 1123 - 1133 . doi: 10.1007/s11071-015-2557-4 http://dx.doi.org/10.1007/s11071-015-2557-4 Isolation control for inertially stabilized platform based on nonlinear friction compensation LU Y J ， ZHANG J N ， YANG S P ， et al . Study on improvement of LuGre dynamical model and its application in vehicle handling dynamics ［J］. Journal of Mechanical Science and Technology ， 2019 ， 33 （ 2 ）： 545 - 558 . doi: 10.1007/s12206-019-0108-5 http://dx.doi.org/10.1007/s12206-019-0108-5 Study on improvement of LuGre dynamical model and its application in vehicle handling dynamics YEH S S ， SU H C . Development of friction identification methods for feed drives of CNC machine tools ［J］. The International Journal of Advanced Manufacturing Technology ， 2011 ， 52 （ 1/2/3/4 ）： 263 - 278 . doi: 10.1007/s00170-010-2720-5 http://dx.doi.org/10.1007/s00170-010-2720-5 Development of friction identification methods for feed drives of CNC machine tools BAI J ， FAN L ， ZHANG S Y ， et al . The parameter identification model considering both geometric parameters and joint stiffness ［J］. Industrial Robot： the International Journal of Robotics Research and Application ， 2019 ， 47 （ 1 ）： 76 - 81 . doi: 10.1108/ir-11-2018-0223 http://dx.doi.org/10.1108/ir-11-2018-0223 The parameter identification model considering both geometric parameters and joint stiffness 肖前进 ， 贾宏光 ， 章家保 ， 等 . 电动舵机伺服系统非线性辨识及补偿 ［J］. 光学 精密工程 ， 2013 ， 21 （ 8 ）： 2038 - 2047 . doi: 10.3788/OPE.20132108. 2038 http://dx.doi.org/10.3788/OPE.20132108. 2038 电动舵机伺服系统非线性辨识及补偿 XIAO Q J ， JIA H G ， ZHANG J B ， et al . Identification and compensation of nonlinearity for electromechanical actuator servo system ［J］. Opt. Precision Eng. ， 2013 ， 21 （ 8 ）： 2038 - 2047 . （in Chinese） . doi: 10.3788/OPE.20132108. 2038 http://dx.doi.org/10.3788/OPE.20132108. 2038 Identification and compensation of nonlinearity for electromechanical actuator servo system 于伟 ， 马佳光 ， 李锦英 ， 等 . 基于LuGre模型实现精密伺服转台摩擦参数辨识及补偿 ［J］. 光学 精密工程 ， 2011 ， 19 （ 11 ）： 2736 - 2743 . doi: 10.3788/ope.20111911.2736 http://dx.doi.org/10.3788/ope.20111911.2736 基于LuGre模型实现精密伺服转台摩擦参数辨识及补偿 YU W ， MA J G ， LI J Y ， et al . Friction parameter identification and friction compensation for precision servo turning table ［J］. Opt. Precision Eng. ， 2011 ， 19 （ 11 ）： 2736 - 2743 . （in Chinese） . doi: 10.3788/ope.20111911.2736 http://dx.doi.org/10.3788/ope.20111911.2736 Friction parameter identification and friction compensation for precision servo turning table 曾德林 ， 肖凯 ， 林竹翀 ， 等 . 航空遥感惯性稳定平台LuGre摩擦参数的分步辨识 ［J］. 光学 精密工程 ， 2016 ， 24 （ 5 ）： 1148 - 1158 . doi: 10.3788/ope.20162405.1148 http://dx.doi.org/10.3788/ope.20162405.1148 航空遥感惯性稳定平台LuGre摩擦参数的分步辨识 ZENG D L ， XIAO K ， LIN ZH CH ， et al . Sub-step identification of LuGre friction parameters of inertially stabilized platform for airborne remote sensing ［J］. Opt. Precision Eng. ， 2016 ， 24 （ 5 ）： 1148 - 1158 . （in Chinese） . doi: 10.3788/ope.20162405.1148 http://dx.doi.org/10.3788/ope.20162405.1148 Sub-step identification of LuGre friction parameters of inertially stabilized platform for airborne remote sensing MIRJALILI S ， MIRJALILI S M ， LEWIS A . Grey wolf optimizer ［J］. Advances in Engineering Software ， 2014 ， 69 ： 46 - 61 . doi: 10.1016/j.advengsoft.2013.12.007 http://dx.doi.org/10.1016/j.advengsoft.2013.12.007 Grey wolf optimizer 10.37188/OPE.20212911.2683 TH117.1 A 1004-924X（2021）11-2683-09 崔靖凯,赛华阳,张恩阳等.基于灰狼算法的模块化关节摩擦辨识和补偿[J].光学精密工程,2021,29(11):2683-2691. 崔靖凯,赛华阳,张恩阳等.基于灰狼算法的模块化关节摩擦辨识和补偿[J].光学精密工程,2021,29(11):2683-2691. CUI Jing-kai,SAI Hua-yang,ZHANG En-yang,et al.Identification and compensation of friction for modular joints based on grey wolf optimizer[J].Optics and Precision Engineering,2021,29(11):2683-2691. CUI Jing-kai,SAI Hua-yang,ZHANG En-yang,et al.Identification and compensation of friction for modular joints based on grey wolf optimizer[J].Optics and Precision Engineering,2021,29(11):2683-2691. 国家自然科学基金项目（No.11672290）；国家重点研发计划项目（No.2016YFC0803000）；吉林省科技发展计划重点研发项目（No.2018020102GX）；中国科学院与吉林省高技术合作与产业化项目（No.2018SYHZ0004） 2020-12-02 2021-02-01 2021-11-15 2022-01-25 光学精密工程 11 29