线缆是工业机器人的重要组成部分，具有传输电动力和控制信号等作用^[1]。除需满足电气性能及与设备相匹配的基本要求之外，伴随着机器人关节长时间、大角度的运动及频繁启停等工况，挂载于机器人关节两端的线缆呈现大角度、高频率的弯曲运动，导致其局部应力疲劳，从而严重影响其服役寿命。因此，明确工业机器人线缆的动态变形特性，对其关节处线缆的空间布局进行优化，有利于改善线缆的机械性能，并提升工业机器人的整体可靠性^[2-3]。

工业机器人线缆为通体柔性，会随着边界条件的改变发生局部弯曲、扭转等非线性变形，同时其受空间约束力的影响，受力状况复杂，而目前尚无成熟的方法能够全面分析其在服役过程中的动态特性。

目前，国内外许多学者对线缆等柔性体进行了建模和仿真研究。如：Terzopoulos等^[4]提出了一种基于Lagrange方程的线缆物理模型，考虑了外载荷，并对几何约束进行了处理，分析了包含线缆重力、抗拉伸和抗弯曲三部分的物理变形能；Grégoire等^[5]基于Cosserat理论对机电产品的线缆进行了仿真分析，在仿真模型中考虑了电缆的重力和接触摩擦力，但没考虑线缆的拉伸变形；Loock等^[6]将质点-弹簧系统应用于线缆装配仿真中，用扭转弹簧表征线缆的动态弯曲特性，但是由于模型的参数较多，数值求解过程较复杂；金望韬等^[7]基于弹性细杆静力学理论，提出了一种光滑平面约束下的活动线缆物理建模与仿真方法。以上研究考虑了线缆的物理属性，并对其静态特性进行了建模和分析，使柔性线缆模型在仿真和布局过程中更具真实性，但没考虑线缆活动时所产生的动态应力问题。杨炜烽等^[8]基于位置动力学，对线缆离散Cosserat杆模型的位姿进行连续的预估和修正，完成了线缆的运动过程仿真，仿真结果受位姿修正迭代次数的影响。此外，还有不少利用多刚体链法^[9-11]、有限元法^[12-13]等进行线缆类柔性体建模的研究。然而，在上述研究中，建模过程都较复杂，数值求解效率较低。因此，需要对工业机器人线缆的动态特性进行建模优化。

在线缆布局设计方面，国内外学者也进行了研究。如：Park等^[14]提出了利用多智能体系统进行协作规划的方法，研发了一种用于导弹线缆布线的线缆并行设计系统；Conru^[15]利用人工智能方法研究了在复杂空间环境中线缆路径的搜索问题；Zhu等^[16]利用基于知识的工程（knowledge based engineering, KBE）技术对飞机线缆的自动布局设计进行了研究，利用离散优化等技术分析了线缆长度约束和布线区域约束等，进行了线缆路径搜索；刘潇等^[17-18]利用改进的运动规划算法，快速获得了单根线缆和分支线满足一定约束的线缆路径；杨啸东等^[19]针对复杂机电产品的线缆装配序列规划问题，提出了装配优先关系矩阵（assembly priority matrix, APM）的概念，对标准粒子群算法的初始化加以改进，并引入了混沌算法，提高了算法的搜索效率；王发麟等^[20]提出了基于数字孪生的动态建模与仿真方法，将基于随机有限集（random finite set, RFS）的同步定位与地图构建（simultaneous localization and mapping, SLAM）方法应用于线缆装配路径与位姿的估计，来构建以动态数据驱动的线缆装配工艺模型。上述线缆布局规划研究主要集中在复杂空间中线缆的可行路径搜索上，同样未考虑线缆本身处于实时运动状态的工况。特别对于工业机器人线缆而言，其在机器人第2、第3关节处通常仅有2个固定点，大部分线缆的服役特性为动态反复弯曲导致的应力动态变化和疲劳损伤，因此，上述研究成果并不适用于工业机器人关节处运动线缆的空间布局与优化。

综上，尽管目前对于线缆动态建模和布局设计的研究取得了一定的成果，但仍没有一种适合工业机器人线缆建模和布局设计的优化算法。因此，本文提出了基于弹簧阻尼链式等效的柔性线缆建模方法，并考虑了线缆动力学模型参数的取值对最终仿真精度的影响，避免了根据经验确定线缆动力学模型的参数值而造成仿真误差。通过线缆仿真参数识别技术来确定线缆力学模型的最优参数，使模型具有所需参数少而精度高的优势，实现了工业机器人线缆的运动仿真；以线缆最大应力最小为优化目标，对工业机器人线缆布局进行优化设计；通过实验，将线缆动力学模型仿真形态与工业机器人线缆的真实运动形态进行对比，来验证线缆模型的准确性。

工业机器人本体线缆如图1所示。本文主要针对工业机器人本体关节处柔性线缆的变形特性进行研究，并对其进行合理布局，以提高线缆寿命。

柔性线缆通过多个线性系统与扭转系统的结合进行模拟，如图2所示。其中，线性系统包括线性弹簧和线性阻尼，扭转系统包括扭转弹簧和扭转阻尼^[21]。线缆被均匀地分为n段，每一段都用一根连杆来表示，线性系统置于连杆之中，使得连杆可以沿着其轴向伸长或压缩。区别于作者先前的研究^[21]，本文从柔性线缆的基座开始，连杆之间的连接处及连杆首末两端都设有一个扭转系统。柔性线缆的总质量m被均匀分散到n个质点上，每个质点位于连杆的末端。质点i (i=1, 2, …, n)与质点i-1之间连杆的长度用Li表示，连杆i与连杆i-1之间的夹角用θ_i 表示。线缆无主动运动能力，其动态变形源于首端和尾端两固定点随机械臂的相对运动。

采用牛顿法对线缆进行受力分析。以任意非特定质点i (i≠n)为例进行分析，如图3所示，以确保分析的广泛适用性。在已知线缆上各质点当前位置和速度的前提下，通过受力分析推导出各质点所受的合力，进而计算出相应的加速度；随后，经过微小时间段的迭代，来更新各质点的位置和速度信息。

针对图3所示的质点i，设其线性系统的拉力为ηi，扭转系统的扭矩为τi，则：

1）质点i受到其下方相邻质点i-1处扭转系统的作用力F₁，F1=τi/Li，方向垂直于两质点之间连杆的长度方向。

2）质点i受到其上方相邻质点i+1处扭转系统的作用力F₂，F2=τi+2/Li+1，方向垂直于两质点之间连杆的长度方向。

3）质点i受到自身扭转系统施加于下方相邻质点i-1力的反作用力，即质点i-1受到质点i处扭转系统的作用力F31'，F31'=τi+1/Li，方向垂直于两质点之间连杆的长度方向。那么，质点i会受到一个与F31'大小相等、方向相反的反作用力F31。

4）质点i受到自身扭转系统施加于上方相邻质点i+1力的反作用力，即质点i+1受到质点i处扭转系统的作用力F32'，F32'=τi+1/Li+1，方向垂直于两质点之间连杆的长度方向。那么，质点i会受到一个与F32'大小相等、方向相反的反作用力F32。

5）质点i受到其下方线性系统所产生的力F₄₁，F41=ηi，方向为沿着连杆i由质点i指向质点i-1。

6）质点i受到其上方线性系统所产生的力F₄₂，F42=ηi+1，方向为沿着连杆i+1由质点i指向质点i+1。

因此，当i≠n时，线缆上的质点会受到上述6个作用力；当i=n时，质点仅受到1）、3）和5）所述的3个作用力。其中，扭转系统的扭矩τi由扭转弹簧的回复扭矩τ_{tor, i} 和阻尼扭矩τ_{dam, i} 构成，线性系统力η_i 由线性弹簧的回复力η_{lin, i} 和线性阻尼力η_{dam, i} 构成。

τ_{tor, i} 仅与其当前连杆转角θ_i 有关，即：

式中：k_tor为扭转弹簧的弹性系数。

η_{lin, i} 仅与其当前连杆长度L_i 有关，即：

式中：k_lin为线性弹簧的弹性系数，L_i0为连杆原长。

求解τ_{dam, i} 和η_{dam, i}，首先需要得到质点i的速度。质点i的速度分解如图4所示。设某时刻质点i及其下方两质点i-1和i-2的速度分别为v_i 、v_i-1和v_i-2。

以连杆i下方的质点i-1作为参照，则连杆i的相对速度vrlt,i=vi-vi-1。将该相对速度分解为2个分量：一个为沿连杆长度方向的v_{A, i}，另一个为垂直于连杆长度方向的v_{B, i}，即vrlt,i=vA,i+vB,i。因此，可求得连杆i的线速度为v_{A, i}，角速度为v_{B, i} /Li，则：

式中：C_lin为线性阻尼系数。

同理，以连杆i-1的下方质点i-2为参照，则连杆i-1的相对速度vrlt,i-1=vi-1-vi-2。同样，将该相对速度分解为沿连杆长度方向和垂直于连杆长度方向的2个分量，即vrlt,i-1=vA,i-1+vB,i-1，进而可求得连杆i-1的线速度为v_{A, i-1}，角速度为v_{B, i-1}/Li_-1。扭转系统i处的角速度ω_i 为连杆i与连杆i-1的角速度之差，即：

则：

式中：C_tor为扭转阻尼系数。

基于上述质点分力的计算，可求得每个质点所受的合力。随后，利用这些合力，可进一步求得当前时刻线缆上所有质点的加速度。

在传统基于质点-弹簧的线缆建模与仿真中，由于弹簧数量较多，仿真时求解速度慢且容易出现失真的情况。本文提出的基于弹簧阻尼链式等效的柔性线缆建模方法，可以用较少的弹簧来模拟线缆，并提高求解效率。同时，该模型还能描述线缆两端固定点运动引起的各质点作用力的变化，从而适用于描述工业机器人线缆的被动变形特征。

工业机器人线缆可分为2种形态：机器人处于稳固静止状态时的静态形态和机器人在执行工作任务时的动态形态。在不同的形态下，选取的线缆模型的弹性系数和阻尼系数不同。因此，利用粒子群算法对上述参数进行识别。通过实验测得线缆的真实形态，并与线缆仿真模型的形态进行对比，对仿真结果与实验结果之间的误差不断迭代来调整模型参数，直至误差值足够小，从而精确识别出最优的模型参数。对线缆弹性系数和阻尼系数的精确识别可以有效提升线缆模型的准确性，对于后续工业机器人线缆仿真与优化具有重要意义。

运用最小势能原理求解线缆模型在约束条件下的静态形态。线缆处于稳定平衡状态时，其势能为最小值，线缆静态形态即线缆模型的重力势能、线性弹簧和扭转弹簧的弹性势能之和为最小值时的线缆形态。

线缆模型的重力势能E_p为：

式中：M为线缆重力，hi为质点i与基座的高度差。

线缆模型线性弹簧的弹性势能E_lin为：

式中：k_lin为线性弹簧的弹性系数。

线缆模型扭转弹簧的弹性势能E_tor为：

式中：k_tor为扭转弹簧的弹性系数。

通过线缆的静态实验识别线缆的弹性系数。取一个线束进行静态实验，测量在多种工况下线缆的空间变形。将机器人线缆的下端固定，上端选择4个不同的位置进行4组静态实验，并用DIC（digital image correlation，数字图像相关）相机测量线缆静态时中间4个标记点的位置坐标。线缆实验装置如图5所示。以固定端为坐标原点O，实验中机器人线缆上端位置信息如表1所示，线缆各标记点的坐标如表2所示。

采用2.1.1节中介绍的线缆静态形态仿真方法对线缆模型进行相应场景的仿真，得到与标记点相同位置的坐标。以仿真结果与实验结果的差值最小为优化目标，目标函数F为：

式中：Et(j)和St(j)分别为第t组实验中第j个标记点的实验结果和仿真结果。

将所有数据和目标函数值代入粒子群优化算法，对线缆的弹簧参数进行识别。根据识别结果调整线性弹簧和扭转弹簧的弹性系数，得到k_lin＝218 485 N/m，k_tor＝28 N/m。种群的适应度变化如图6所示。由图可知，第50代种群的适应度为202，表示上述参数与实验结果的误差为202 mm，各点的平均误差为12.6 mm。

给定线缆模型某一质点的运动位置和速度，根据1.2节中线缆质点受力分析，可得到当前时刻所有质点的加速度。将线缆运动分为若干个微小时间单元，由牛顿法可得到下一时间单元的质点位置、姿态、速度并进行受力分析，得到新的质点位置，按此方法依次迭代便可实现线缆的动态形态仿真。

通过线缆的动态实验识别线缆的阻尼系数。取一个线束进行动态实验，测量线缆在运动过程中的空间变形。将线缆下端固定，线缆上端运动，并用DIC相机测量线缆运动时中间4个标记点相对初始位置的位移。在线缆运动过程中每隔0.1 s选取一组数据，共取6个时刻，得到6组数据，如表3所示。

对线缆模型进行相应场景的动态仿真，得到与实验相同时刻的标记点位移。以仿真结果与实验结果的差值最小为优化目标，将所有数据和目标函数值代入粒子群优化算法，对线缆的阻尼参数进行识别。根据识别结果调整线性阻尼系数和扭转阻尼系数，得到C_lin＝45 N/(m·s^-1)，C_tor＝28 N/(m·s^-1)。种群的适应度变化如图7所示。由图可知，第50代种群的适应度为57.1，表示上述参数与实验结果的误差为57.1 mm，各点的平均误差为14.3 mm，其占线缆总长度的4.8%。

为了验证线缆物理简化模型与运动仿真的准确性，进行实际工况下工业机器人关节处线缆运动的仿真和实验。实验装置如图5所示。将线缆下端在台架上用卡箍固定，上端夹持在六轴机器人的末端执行器中。通过调整机器人的末端轨迹，可模拟实际工况下线缆在机械臂上安装位置不同导致的运动特性差异。以六轴机器人末端执行器为线缆运动端，通过编程即可实现模拟工业机器人某一关节处线缆的运动。在实验过程中，用DIC相机对线缆表面的标记点进行追踪拍摄，从而得到标记点的位移。

以埃斯顿工业机器人ER12B-1510第2关节处的线缆为例进行分析。机器人第2关节处的线缆如图8所示。调节线缆实验装置中线缆固定端与运动端的相对距离，使其与图8所示一致。以第2关节与线缆运动端的距离为半径，通过编程控制实验装置中六轴机器人的末端轨迹，使其模拟第2关节逆时针旋转90°，再顺时针旋转90°。安装的可运动线缆的长度为0.3 m，在线缆中间均匀粘贴4个标记点，用DIC相机拍摄线缆运动过程，结果如图9所示。

根据所建立的线缆动态模型及所识别的仿真参数，在MATLAB软件中进行线缆运动仿真。仿真中将线缆模型的质点数量设定为10个，线缆长度为0.3 m，直径为0.028 m，质量为0.9 kg，与真实线缆保持一致，并且线缆两端的相对位置及其与工业机器人第2关节的相对位置均与实验一致。线缆运动仿真结果如图10所示。

线缆运动实验与仿真结果的对比如图11所示。由图可知，实验与仿真结果存在一定差值，最大误差为10.5 mm，占线缆总长度的3.5%，符合预期要求，验证了线缆模型的准确性。

在工业机器人工作时，其关节处线缆的应力随线缆两端的固定位置、运动端出线方向及活动线缆长度的变化而变化。为了优化线缆布局，提高线缆寿命，可通过优化线缆两端的固定位置、运动端出线方向及活动线缆的长度来减小线缆应力。基于线缆动态模型求得各个质点处的应力，建立使线缆最大应力最小为优化目标的线缆布局优化模型。目标函数F₁为：

式中：p₀为线缆固定端沿X向的位移，p_l为线缆运动端沿Z向的位移，q_l为线缆运动端的出线方向与机械臂在XOZ平面的夹角，L为活动线缆长度，σi为线缆质点i的应力。

p₀、p_l、q_l和L为线缆布局优化模型的主要优化参数，如图12所示。

根据第1.2节中线缆的受力分析，对线缆各质点所受的应力进行计算。以质点i为例，分析线缆运动中的应力。为了确保分析的普适性，取i≠n。

质点i受到前4个扭转系统产生的应力，为：

式中：M1、M2、M31和M32分别为F1、F2、F31和F32所产生的扭矩，b为质点到线缆轴心的距离，I为极惯性矩，d为线缆直径。

质点i受到后2个线性系统产生的应力，为：

式中：A为线缆的横截面面积。

当i=n时，质点仅受σ1，σ31和σ42三个应力。

对3.2节所示的工业机器人第2关节处的线缆进行布局优化设计，对线缆两端的固定位置、运动端出线方向及活动线缆长度进行调节。设定p₀=-5~5 cm，p_l=-5~5 cm，q_l=-45°~45°。根据企业调研总结，线缆长度的设计要求为：活动线缆长度应大于线缆两端距离的1.1倍，且不宜过长，故取L=1.1~1.4L_c，其中L_c为线缆两端之间的距离，如图12所示。

改变其中一个线缆布局优化参数，保持其余参数不变，则线缆最大应力的变化趋势如图13所示。

由图13可得：p₀越大，线缆最大应力越大；p_l对线缆最大应力的影响不大；q_l≈-20°时，线缆最大应力最小；线缆的最大应力随L的增大而减小。因此，线缆布局优化结果为：p₀=-0.05 m，q_l=-20°，L=1.4L_c。

通过线缆各质点应力的仿真计算，得到工业机器人第2关节处线缆在布局优化前后最大应力的变化，如图14所示。由图可知，相较于布局优化前，布局优化后线缆的最大应力约减小了83%。

为了验证线缆布局优化方案的可行性，通过实验对布局优化前后线缆运动过程中的最大应变进行测量。首先对线缆表面进行喷涂散斑处理，用DIC相机识别线缆表面散斑，并进行线缆表面曲面构建，如图15所示，从而对线缆运动时的变形表面进行图像采集与应变计算。线缆布局优化前后最大应变的变化趋势如图16所示。由图可知，相较于布局优化前，布局优化后线缆的最大应变约减小75%，故线缆的最大应力也相应减小，验证了线缆布局优化方案的可行性。

1）本文提出了基于弹簧阻尼链式等效的柔性线缆建模方法，将线缆通过线性弹簧、线性阻尼、扭转弹簧和扭转阻尼进行分割，实现了其质点动力学描述。

2）基于建立的线缆模型，通过最小势能法实现了线缆的静态形态仿真，并通过对线缆各质点的受力分析，由牛顿法解得此时质点的动力学参数，经过一个微小时间内的运动后即可得到下一时刻各质点的位置，迭代上述步骤即可实现线缆动态形态仿真。而后，通过线缆静态和动态实验，利用粒子群算法对仿真结果与实验结果的误差进行迭代，直至误差值足够小，从而得到线缆动力学模型的最优参数。搭建了模拟工业机器人第2关节处线缆运动的实验平台，用DIC相机测得线缆运动形态，将实验结果与仿真结果进行对比，验证了线缆动力学模型的准确性。

3）提出了工业机器人线缆布局优化策略，计算了线缆动力学模型各质点处的应力。通过改变线缆两端的安装位置、运动端出线方向及活动线缆的长度，使得线缆的最大应力最小，从而实现线缆的布局优化。通过对布局优化前后线缆运动过程中线缆最大应力、最大应变的对比，验证了所提出的工业机器人线缆布局优化方案的可行性。