一种基于刚度性能优化的移动机械臂站位规划方法

1.本发明属于移动机械臂技术领域，具体涉及一种基于刚度性能优化的移动机械臂站位规划方法。


背景技术：

2.现如今，随着技术的发展，越来越多实际作业场景中将六自由度工业机器人的操作能力和移动平台的移动能力相结合，形成了一种冗余自由度移动机械臂系统；该系统实现了末端工作空间大、工作效率和灵活性高、系统低功耗和低成本的有机结合；在进行作业任务时，由于系统的冗余性，面对同一作业任务，存在多个移动平台相对于该作业任务的停靠点均能完成任务，当其运动到提前规划好的停靠点静止后，此时在该停靠点对应移动平台相对于墙体的位姿被称为站位；再由工业机器人完成该作业任务，如此循环交替工作，直至完成全部作业，但该系统同时也将面临如下问题：与固定基座机械臂相比，该系统为融合了六自由度工业机器人和移动平台的冗余自由度系统，附加了移动平台相对于任务的站位移动环节；现如今关于移动机械臂站位规划还是依靠人工根据现场情况进行移动平台的站位布置，难以规划至最优的站位，而产生了移动平台站位移动次数多、整体作业精度低和效率低等问题；因此，必须针对各种作业任务开发一种普适性的算法来对移动平台相对于具体任务的站位进行优化，从而以最高的工作效率、作业精度实现作业任务。


技术实现要素：

3.本发明的目的在于提供一种移动机械臂系在统执行多种实际作业场景中为提高整体作业精度，而提出的基于刚度性能优化的移动机械臂站位规划方法。
4.基于上述目的，本发明采用如下技术方案：一种基于刚度性能优化的移动机械臂站位规划方法，包括如下步骤：
5.步骤(1)建立坐标系：确定移动平台坐标系、机械臂基坐标系、机械臂各连杆坐标系、法兰坐标系、打印单元坐标系和墙体坐标系；
6.步骤(2)根据步骤(1)建立的坐标系，建立整个系统的运动学模型，根据整个系统的运动学模型及实际作业任务确定移动平台相对所需打印墙体的可达站位空间；
7.步骤(3)根据步骤(1)建立的坐标系，建立机械臂的力学模型：已知机械臂末端所期望的运动轨迹，计算连杆构型、连杆质量、末端载荷和各关节角速度的因素作用下的各关节驱动力矩；
8.步骤(4)根据步骤(3)建立的机械臂的力学模型，确定机械臂末端变形量作为机械臂刚度评价指标；
9.步骤(5)根据步骤(4)得到的机械臂刚度评价指标做为优化目标，以移动平台相对实际工作任务在可达站位空间内的站位做为设计变量，建立单次作业的优化函数，利用遗传算法进行合理的站位优化。
10.优选的，机械臂为六连杆机械臂，根据改进的dh参数法分别确定机械臂各连杆坐
标系；机械臂包括第一连杆、第二连杆、第三连杆、第四连杆、第五连杆和第六连杆；机械臂各连杆坐标系包括第一连杆坐标系、第二连杆坐标系、第三连杆坐标系、第四连杆坐标系、第五连杆坐标系和第六连杆坐标系；第一连杆坐标系、第二连杆坐标系、第三连杆坐标系、第四连杆坐标系、第五连杆坐标系和第六连杆坐标系分别对应机械臂的第一连杆、第二连杆、第三连杆、第四连杆、第五连杆和第六连杆。
11.根据权利要求2所述的方法，其特征在于，步骤(2)中，建立整个系统的运动学模型的步骤具体为：首先根据机械臂各连杆的坐标系和机械臂各连杆的结构尺寸确定机械臂运动学模型；然后根据移动平台的尺寸、机械臂安装在移动平台上的位置和作业任务确定整个系统运动学模型，即
[0012][0013]
式中：为移动平台坐标系相对于机械臂基坐标系的齐次变换矩阵；为打印工具坐标系相对于法兰坐标系的齐次变换矩阵；为打印工具坐标系相对墙体坐标系的齐次变换矩阵；为移动平台坐标系{m}相对墙体坐标系{w}的齐次变换矩阵；为法兰坐标相对于第六连杆坐标系的齐次变换矩阵；为第六连的坐标系相对机械臂基坐标系的齐次变换矩阵；
[0014]
其中，第六连杆的坐标系相对机械臂基坐标系的齐次变换矩阵由整个系统的运动学模型求得，为
[0015][0016]
根据机械臂运动学模型列出机械臂的逆运动学方程，根据机械臂的逆运动学方程校验移动平台相对实际作业任务在可达站位空间内规划的某一站位下，判断机械臂的末端位姿序列中是否各点均存在运动学逆解且满足各轴关节角度限制。
[0017]
优选的，步骤(2)中，根据整个系统的运动学模型及实际作业任务确定移动平台相对实际作业任务的可达站位空间的具体过程为：首先根据机械臂运动学模型绘制机械臂的可达工作空间云图；然后根据实际作业任务的尺寸、机械臂安装在移动平台上的位置和机械臂末端打印工具的尺寸确定移动平台相对作业任务的可达站位空间；移动平台相对作业任务的可达站位空间转化为移动平台坐标系在墙体坐标系下的位姿：
[0018][0019]
式中：x,y为移动平台坐标系的原点在墙体坐标系下的位置，β为移动平台坐标系相对于墙体坐标系竖直方向z轴的姿态旋转角度，r表示移动平台相对于作业任务的站位空间的最远距离。
[0020]
优选的，根据机械臂运动学模型列出机械臂的逆运动学方程，具体步骤为：首先根据机械臂运动学模型列出机械臂的正运动学方程，然后根据代数法进行机械臂逆运动学求解，得机械臂逆运动学方程如公式(2)：
[0021][0022]
式中：θ1、θ2、θ3、θ4、θ5、θ6为各个关节的角度值，为第一连杆坐标系相对机械臂基坐标系的齐次变换矩阵，为第二连杆坐标系相对第一连杆的坐标系的齐次变换矩阵，为第三连杆坐标系相对第二连杆坐标系的齐次变换矩阵，为第四连杆坐标系相对第三连杆坐标系的齐次变换矩阵，为第五连杆坐标系相对第六连杆坐标系的齐次变换矩阵，为第六连杆坐标系相对第五连杆坐标系的齐次变换矩阵，为第六连杆坐标系相对机械臂基坐标系的齐次变换矩阵：
[0023]
优选的，步骤(3)中所述机械臂力学模型的建立的步骤具体为：首先由外递推计算各连杆的速度与加速度，再由牛顿-欧拉公式计算出各连杆的惯性力和力矩；向内递推计算各连杆相互作用的力和力矩，以及各关节驱动力矩；
[0024]
机械臂各连杆的惯性力矩：
[0025]fj
fa＝m[(a
j-g) α
×bc
ω
×
(ω
×bc
)]
[0026]
式中：m表示连杆质量；α表示为角加速度；ω表示角速度；fj表达连杆i的作用力；fa表示连杆i 1的作用力；aj表示连杆j的线加速度；bc表示关节与连杆质心的偏距；
[0027]
机械臂各连杆的力矩：
[0028][0029]
式中：式中：式中：是bc的反对称矩阵，表示的转置矩阵。
[0030]
求得各关节绕转轴的驱动力矩：
[0031][0032]
式中：为连杆坐标系{i}在机械臂基坐标系下的方向，i＝1、2、3、4、5、6，mi为关节i上的约束力矩矢量，mi＝[m
ix m
iy m
iz
]，m
ix
，m
iy
和m
iz
表示沿x,y,z三个方向上的约束力矩，zi＝[0 0 1]
t
。
[0033]
优选的，步骤(4)中，确定机械臂末端变形量作为机械臂刚度评价指标的具体过程为：以机械臂末端沿z轴方向的变形量为机械臂刚度评价指标，其中，械臂末端沿z轴方向的变形量记作μz[0034]
μz＝z
·
δx＝z
·
j(θ)
·
(k
θ
)-1
·
τ
[0035]
式中：z＝[0 0 1]，δx为机械臂末端产生的广义变形量，j(θ)为机械臂的雅可比矩阵，τ为τi的矩阵式，k
θ
为机械臂关节空间的刚度矩阵。
[0036]
优选的，步骤(5)建立的优化函数，以步骤(2)确定的移动平台相对实际作业任务的可达站位空间为界限约束，以移动平台在移动平台相对实际作业任务的可达站位空间内相对实际作业任务的位姿满足整个系统的逆运动学方程为性态约束。
[0037]
优选的，步骤1中移动平台坐标系、机械臂基坐标系、机械臂各连杆坐标系、法兰坐标系、打印单元坐标系和墙体坐标系均为右手坐标系。
[0038]
与现有技术相比，本发明的有益效果如下：
[0039]
1.本发明提供一种基于刚度性能优化的移动机械臂站位规划方法，该方法可应用于多种实际作业场景中；本发明根据建立的坐标系，建立整个系统的运动学模型，根据整个系统的运动学模型计算移动平台相对作业任务的理论可达站位空间，在根据机械臂的逆运动学方程对理论可达站位空间进行验证；然后根据整个系统的力学模型提出以机械臂末端沿z轴方向的变形量为机械臂刚度评价指标，求出移动平台坐标系在墙体坐标系不同位姿下的机械臂末端沿z轴方向的变形量，以移动平台坐标系在墙体坐标系位姿下的机械臂末端沿z轴方向的最小变形量为最优站位；从而极大减少了移动机械臂的站位移动次数，使得机器人的作业效率得到显著提升；
[0040]
2.本发明应用牛顿-欧拉方法建立机械臂动力学方程过程中，除了考虑连杆构型与末端载荷对各关节驱动力矩影响外，还引入了机械臂各连杆质量、各关节角速度和角加速度这些影响因素，这使得后续根据该动力学方程建立的机械臂刚度模型以及求解得到的机械臂末端变形量更加精确；
[0041]
3.由于本发明提供的一种移动机械臂站位规划方法可应用于多种实际场景中，将针对具体实际作业场景提出机械臂刚度评价指标，即机械臂末端沿某个轴方向变形量来进行移动机械臂的站位规划，可以使得机械臂末端沿该轴方向的最大变形量在移动平台理论上的可达站位空间内达到最小，整体的机械臂末端作业精度达到最优；
[0042]
4.本发明使用遗传算法在整个站位空间中搜寻移动平台的最优站位，能够跳出局部最优而找到全局最优点，而且遗传算法允许使用非常的目标函数，并对变量的变化范围可以加以限制。基于遗传算法在整个站位空间中搜寻移动平台的最优站位的过程是利用matlab遗传算法工具箱，其扩展了matlab在处理优化问题方面的能力，使结果更具有准确性与高效性，解决了依靠人工根据现场情况进行移动平台的站位布置而使得移动平台站位移动次数多、整体打印精度低和效率低等问题。
附图说明
[0043]
图1为本发明的流程图；
[0044]
图2为本械臂系统结构示意图以及系统各坐标系关系图；
[0045]
图3为kuka kr90 r3100型机械臂各连杆坐标系；
[0046]
图4为kuka kr90 r3100型机械臂可达工作空间云朵图；
[0047]
图5为移动平台可达性站位空间分析示意图；
[0048]
图6为机械臂简化示意图；
[0049]
图7为单连杆力平衡示意图；
[0050]
图8为相邻连杆间力关系示意图；
[0051]
图9为adams仿真虚拟样机；
[0052]
图10为在选取机械臂不同位形下机械臂末端沿x、y、z方向的变形量大小；
[0053]
图11为民用建筑实例的平面结构及尺寸示意图；
[0054]
图12为单次作业任务四种类型墙体站位示意图；
[0055]
图13为针对一字墙体3d打印任务移动机械臂系统的最优站位示意图；
[0056]
图14为针对直角墙体3d打印任务移动机械臂系统的最优站位示意图；
[0057]
图15为针对圆弧墙体3d打印任务移动机械臂系统的最优站位示意图；
[0058]
图16为针对t形墙体3d打印任务移动机械臂系统的最优站位示意图。
[0059]
图中：1-移动平台；2-kuka kr90 r3100型的机械臂；3-混凝土3d打印单元；4-直角墙体。
具体实施方式
[0060]
下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的建筑3d打印场景下的实施例是本发明一部分实施例，而不是全部的实施例。
[0061]
如图2所示该实施例以直角墙体作为示例，为便于对本实施例进行理解，首先对本发明实施例所公开的一种移动机械臂整车系统进行简单介绍；移动机械臂系统如图2所示，移动机械臂系统包括移动平台、机械臂和混凝土3d打印单元；机械臂为kuka kr90 r3100机械臂，包括机械臂基座、第一连杆、第二连杆、第三连杆、第四连杆、第五连杆、第六连杆和法兰；机械臂基座固定安装在移动平台上，混凝土3d打印单元固定安装在机械臂末端法兰上；作业任务为所需打印的墙体。
[0062]
基于上述移动机械臂系统，本发明提供一种如图1所示的基于刚度性能优化的移动机械臂站位规划方法，具体步骤如下：
[0063]
步骤(1)建立坐标系：确认移动平台的坐标系为移动平台坐标系，机械臂基座的坐标系为机械臂基坐标系，法兰的坐标系为法兰坐标系，混凝土3d打印单元的坐标系为工具坐标系，所需打印墙体的坐标系为墙体坐标系；机械臂各连杆的坐标系为机械臂各连杆坐标系；移动平台坐标系、机械臂基坐标系、机械臂各连杆坐标系、法兰坐标系、混凝土3d打印单元坐标系和墙体坐标系均为右手坐标系；
[0064]
选取移动平台中心为移动平台坐标系的原点，机械臂基座中心为机械臂基坐标系的原点，机械臂基座中心与移动平台中心有一距离偏置l；机械臂基坐标系与移动平台坐标系的x轴方向均平行于移动平台长边方向，移动平台坐标系的x轴正方向由移动平台坐标系的原点指向机械臂基坐标系的原点，即移动平台坐标系与机械臂基坐标系x轴方向共线；移动平台坐标系和机械臂基坐标系的z轴正方向平行且竖直向上；移动平台坐标系与机械臂基坐标系的y轴方向符合右手法则，根据右手法则可对机械臂基坐标系和移动平台坐标系的y轴方向进行判断。
[0065]
根据改进的dh参数法确定机械臂各连杆坐标系；机械臂包括第一连杆、第二连杆、第三连杆、第四连杆、第五连杆和第六连杆；机械臂各连杆对应的坐标系如图3，分别为第一连杆坐标系{1}、第二连杆坐标系{2}、第三连杆坐标系{3}、第四连杆坐标系{4}、第五连杆坐标系{5}和第六连杆坐标系{6}；移动平台坐标系记作{m}，机械臂基坐标系记作{0}，法兰坐标系记作{f}，打印工具坐标系记作{t}，墙体坐标系记作{w}。
[0066]
dh参数法是由denavit和hartenberg于1955年提出的一种描述串联式链路上连杆
和关节的系统方法；1986年khalil，kleinfinger提出了一种改进的dh参数，在建立关节坐标系时将坐标系固结于该连杆的近端，而非远端，更符合直观理解；为了便于区分，这种方法成为改进dh(modified dh)，而将之前的方法称为标准dh(standard dh)。
[0067]
步骤(2)，建立整个系统的运动学模型：,根据机械臂各连杆对应的坐标系和机械臂各连杆的结构尺寸确定机械臂运动学模型；然后根据移动平台的尺寸、机械臂安装在移动平台上的位置和作业任务确定整个系统运动学模型：
[0068]
根据建立的机械臂各连杆坐标系和机械臂各连杆的结构尺寸得到该型号机械臂的dh参数表，如表1所示，其中θ1,θ2,θ3,θ4,θ5,θ6为未知数；各连杆间位姿的关系就能够通过坐标系间的变换得到，即可明确表示出各连杆之间的运动关系，相邻两坐标系间的位移关系由连杆扭角α、连杆长度a、关节转角θ和关节偏置d这4个参数描述。
[0069]
表1：dh参数表：
[0070]
关节ia
i-1
(mm)α
i-1
(
°
)di(mm)θi(
°
)10180-675θ1(0)2350900θ2(0)3135000θ3( 90)441-90-1400θ4(0)50900θ5( 180)60900θ6(0)
[0071]
连杆变换可定义为从连杆坐标系{i-1}到连杆坐标系{i}的变换，这个其次变换可以当作坐标系{i}经过四个子变换得到，因这些子变换是相对动坐标系描述，按照坐标变换的原则，得到：
[0072][0073]
式中：表示连杆坐标系{i}相对连杆坐标系{i-1}的齐次变换矩阵；rot表示旋转变换；tran表示平移变换；x表示方向x；z表示方向z；di表示关节i偏置；θi表示关节i转角；a
i-1
表示连杆i长度；
[0074]
由式(1)可以得到连杆变换的通式：
[0075][0076]
将所有连杆变换依次相乘，得末端齐次变换矩阵即机械臂正运动学方程
[0077][0078]
式中：为第一连杆坐标系{1}相对机械臂基坐标系{0}的齐次变换矩阵，为第二连杆坐标系{2}相对第一连杆{1}的坐标系的齐次变换矩阵，为第三连杆坐标系{3}相对第二连杆坐标系{2}的齐次变换矩阵，为第四连杆坐标系{4}相对第三连杆坐标系{3}的齐次变换矩阵，为第五连杆坐标系{5}相对第六连杆坐标系{4}的齐次变换矩阵，为
第六连杆坐标系{6}相对第五连杆坐标系{5}的齐次变换矩阵，为第六连杆坐标系{6}相对机械臂基坐标系{0}的齐次变换矩阵。
[0079]
机械臂的逆运动学问题是在已知工作时最后一个连杆坐标系即连杆坐标系{6}相对于基坐标系{0}的齐次变换矩阵求解满足要求作业任务时机械臂各关节变量值(θ1,θ2,θ3,θ4,θ5,θ6)的问题，dh参数表中其他的参数均为已知。
[0080]
采用代数法(也称反变换法)进行机械臂逆运动学求解，机械臂逆运动学方程为。
[0081][0082]
式中：θ1、θ2、θ3、θ4、θ5、θ6为各个关节的角度值，式中：p
x
、py、pz表示平移数据：r
11
、r
12
、r
13
、r
21
、r
22
、r
23
、r
31
、r
32
、r
33
表示旋转角度数据；
[0083]
在公式(4)中，已知等号左侧的位置与姿态组成的矩阵中16个元素，需要求解等号右侧的6个关节角度变量；因此，逆运动学可以看作依次将等号右侧含有待求角度变量的矩阵求逆，将关节角度变量分离，且求得各个关节角度的过程。
[0084]
根据建系结果可以得到法兰坐标系{f}相对于连杆坐标系{6}的连杆参数，故可得到法兰坐标系{f}相对于连杆坐标系{6}的齐次变换矩阵
[0085]
根据移动平台坐标系{m}在墙体坐标系{w}下的位姿[x,y,θ]，其中：x,y为移动平台坐标系{m}的原点在墙体坐标系{w}下的位置，θ为移动平台坐标系{m}相对于墙体坐标系{w}竖直方向z轴的姿态旋转角度；按照坐标变换的原则，得到为移动平台坐标系{m}相对墙体坐标系{w}的齐次变换矩阵
[0086]
因为在进行墙体3d打印过程中，打印工具坐标系{t}走过的坐标即为墙体离散出的任务点，该任务点可以从墙体坐标系中获取，并且可以得到工具坐标系{t}相对墙体坐标系{w}的齐次变换矩阵
[0087]
由于混凝土3d打印单元相对于机械臂末端法兰是固定安装，在不考虑其工作中撞击等原因造成变形的情况下，其结构尺寸固定不变，所以可以得出打印工具坐标系{t}相对于法兰坐标系{f}的齐次变换矩阵
[0088]
根据移动平台坐标系{m}与机械臂基坐标系{0}之间的偏置距离和连杆变换可以得出移动平台坐标系{m}相对于机械臂基坐标系{0}的齐次变换矩阵
[0089]
综上所述，列出整个系统的运动学方程，即
[0090][0091]
根据公式(5)得，
[0092][0093]
公式(6)中变量为移动平台坐标系{m}在墙体坐标系{w}下的位姿[x,y,θ]，根据移动平台坐标系{m}在墙体坐标系{w}下不同位姿[x,y,θ]求出再根据已建立的机械臂的逆运动学方程(4)，校验在规划的某一站位下，末端打印位姿序列中是否各点均存在运动学逆解且满足各轴关节角度限制。
[0094]
根据系统运动学模型及实际作业任务确定移动平台的站位空间：如图4所示，采用几何法可绘制出kuka kr90 r3100机械臂的可达工作空间云朵图，包络线内的区域表示机械臂末端可达空间；在3d打印墙体时，为保证储料桶不发生倾漏，混凝土3d打印单元始终竖直向下；当混凝土3d打印单元出料口恰好接触地面打印最低层墙体时，应保证混凝土3d打印单元的打印头不与机械臂本体发生干涉(碰撞)；根据机械臂各连杆参数、移动平台高度以及打印头尺寸参数，通过几何作图可以得出料口接触地面时，出料口距机械臂基座中心的最小距离为1328mm；当打印最高层墙体时，出料口距地面达到2800mm，且要求墙体厚度200mm在机械臂工作范围内，通过几何作图可以得到出料口距机械臂基座中心最大距离为2679mm，分析简图如图5(a)所示；如图5(b)为本项目建筑3d打印移动机械臂系统的目标墙体打印范围俯视图。
[0095]
依据上述分析可知，为满足目标墙体的打印，墙体在机械臂基坐标系下的范围1328-2679mm的圆环区域，因为本文目标是对移动平台相对于墙体的站位进行规划，根据建立的系统各坐标系中机械臂基座中心距移动平台中心的距离548mm，得到目标墙体在移动平台坐标系下的最远可达距离为3227mm；综上移动平台相对于墙体的站位空间为半径为3227mm的圆形区域；将移动平台相对于墙体的理论站位空间转化为移动平台坐标系{m}在墙体坐标系{w}下的表达为：
[0096][0097]
式中：x,y为移动平台坐标系{m}的原点在墙体坐标系{w}下的位置，θ为移动平台坐标系{m}相对于墙体坐标系{w}竖直方向z轴的姿态旋转角度；r为移动平台相对于墙体的站位空间为半径，计算得r＝3227mm。
[0098]
受移动平台相对于墙体的姿态以及墙体的类型、尺寸影响，该区域仅为移动平台理论上的可达站位空间，在该空间内并不是所有的移动平台站位都能保证目标墙体所有点均可达，在下文中选取移动平台站位进行站位优化时，首先会利用机械臂逆运动学进行可达性校验：根据公式(6)得出移动平台坐标系{w}相对墙体坐标系{w}不同位姿[x,y,θ]下的然后通过公式(4)机械臂逆运动学方程进行验证，存在逆解即表示在此站位下，打印工具可达作业单元的所有点；反之，若逆解不存在，则此站位不符合可达性要求，应不予考虑。
[0099]
步骤(3)，建立机械臂动力学方程，本文主要采用牛顿-欧拉方法进行机械臂的动力学建模：主要研究已知机械臂末端所期望的运行轨迹，计算连杆构型、连杆质量、末端载荷、各关节角速度多种因素作用下的各关节驱动力矩；整个方法主要由两个步骤组成：首先由外递推计算各连杆的速度与加速度，由牛顿-欧拉方法计算出各连杆的惯性力和力矩，第二步向内递推计算各连杆相互作用的力和力矩，以及各关节驱动力矩。
[0100]
由于在机械臂末端额定负载情况下，机械臂前三关节比后三关节对于末端精度的影响更大，所以在动力学计算部分，将仅考虑机械臂前三关节；该型号机械臂前三个连杆相关参数如表1，其中连杆3质量参数合并了机械臂的肘部、前臂、手腕机构和末端法兰。
[0101]
图6为机械臂简化示意图，关节与连杆的运动学建模包含了静态部分和动态部分；静态部分代表了每个连杆的初始配置；动态部分代表了每个关节或连杆的运动效用。
[0102]
静态部分：即每个连杆的初始配置，包含两个部分，静态旋转矩阵r
si
和静态体向量b
′
si
，表达如下：
[0103]rsi
＝r
xryrz
ꢀꢀꢀꢀꢀꢀꢀ
(8)
[0104]b′
si
＝xix
′i yiy
′i ziz
′iꢀꢀꢀꢀꢀꢀꢀ
(9)
[0105]
式中：r
si
为静态旋转矩阵；b
′
si
为静态体向量；r
x
、ry与rz是局部坐标系{i}相对局部坐标系{i-1}关于x、y、z三个轴的旋转矩阵；x
′i、y
′i和z
′i是局部坐标系{i}的三个单位矢量；b
′
si
是在局部坐标系{i}下连杆i初始构型的矢量；“'”表示局部坐标系；xi、yi和zi表示局部坐标系{i}的三个方向上的分量。
[0106]
动态部分，同样也包含两个部分：动态旋转矩阵r
mi
与动态体向量b
′
mi
。
[0107]
连杆i的总平移和旋转表达如下：
[0108][0109]b′i＝b
′
si
b
′
mi
ꢀꢀꢀꢀ
(11)
[0110]
式中：r
si
为静态旋转矩阵；b
′
si
为静态体向量；r
mi
为动态旋转矩阵；b
′
mi
为动态体向量；表达相邻连杆之间的旋转矩阵；b
′i表示相邻连杆之间的平移变换；
[0111]
实际上，动态旋转矩阵r
mi
可视为机械臂关节的旋转电机所驱动的旋转副，r
mi
与rz(θi)等价，θi为关节i的旋转角；动态体向量b
′
mi
可视为直线电机驱动的移动副；rz(θi)表示绕z轴旋转θi的旋转矩阵。
[0112]
连杆坐标系{i}的位置pi和方向在基坐标系下可表达如下：
[0113][0114][0115]
式中：pi为连杆坐标系{i}在机械臂基坐标系{0}下的方向位置；为连杆坐标系{i}在机械臂基坐标{0}系下的方向；b
‘
i-1
表示表示相邻连杆之间的平移变换；表示相邻连杆之间的旋转矩阵。
[0116]
对公式(12)和公式(13)分别求导可得到线速度与角速度的向前迭代方程：
[0117][0118][0119]
式中：vi为连杆坐标系{i}原点的线速度矢量；ωi为连杆i的角速度矢量；表示相邻连杆之间的平移变换的导数，表示相邻连杆之间的平移变换的导数，其中表示动态体向量的导数；表示连杆i相对于连杆i-1的角速度，
其中为关节i的角速度，即为关节i的角速度，即表示连杆i相对于连杆i-1的静态旋转矩阵。
[0120]
进一步，对公式(14)与公式(15)求导，可得到线加速度和角加速度的向前迭代方程：
[0121][0122][0123]
式中：ai为连杆坐标系{i}原点的线加速度矢量；αi为连杆i的角加速度矢量；表示连杆i相对于连杆i-1的角加速度矢量其中为关节i的角加速度，即度，即表示为连杆i相对于连杆i-1的角加速度矢量；表示相邻连杆之间的平移变换的二阶导数；。
[0124]
将机械臂的连杆看成刚体，它的质心加速度ac、总质量m与产生这一加速度的作用力f之间的关系满足牛顿第二运动定律：
[0125]
f＝macꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0126]
当刚体绕过质心的轴线旋转时，角速度ω、角加速度α、惯性张量ic与作用力矩m之间满足欧拉方程：
[0127]
m＝icα ω
×
(icω)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(19)
[0128]
如图7所示，根据公式(18)和公式(19)推导出单连杆质心处的力平衡方程和力矩平衡方程为：
[0129]
mac＝mg fj faꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
[0130]
icα ω
×
(icω)＝ma mj fj×bc
fa×bf
ꢀꢀꢀꢀꢀꢀꢀ
(21)
[0131]
式中：fj表示连杆i的作用力；fa表示连杆i 1的作用力；bf表示关节i相对于连杆i的偏距；bc表示关节i-1相对于连杆i的偏距；ma表示作用在关节i上的作用力矩；mj表示作用在关节i-1上的作用力矩；g表示重力加速度。
[0132]
因为连杆质心处线速度为：
[0133]
vc＝vj ω
×bc
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(22)
[0134]
式中：vc表示为连杆质心处线速度；vj表示为关节处的线速度。
[0135]
对公式(22)求导得质心处线加速度为：
[0136]ac
＝aj α
×bc
ω
×
(ω
×bc
)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(23)
[0137]
式中：ac表示质心处的加速度；aj表示关节处的线加速度；
[0138]
将公式(23)带入公式(20)，得：
[0139]fj
fa＝m[(a
j-g) α
×bc
ω
×
(ω
×bc
)]
ꢀꢀꢀꢀꢀꢀꢀ
(24)
[0140]
通过公式(24)求出fj，将其带入公式(21)，得：
[0141]
ma mj fa×
(b
f-bc)＝icα m[(a
j-g) α
×bc
ω
×
(ω
×bc
)]
×
(-bc) ω
×
(icω) (25)
[0142]
综合公式(24)与(25)，得：
[0143]
[0144]
式中：式中：式中：是bc的反对称矩阵，表示的转置矩阵。
[0145]
所以根据以上推导与计算，如图8得出机械臂各关节处约束力与约束力矩向后迭代方程如下：
[0146][0147]
式中：wi包含关节i上的关节约束力矢量fi和约束力矩矢量mi，且wi是在基坐标系{0}下的表示，即：
[0148][0149]
式中：f
it
表示作用力矩阵，m
it
表示作用力矩矩阵，f
ix
，f
iy
，f
iz
表示作用力三个方向的分量，m
ix
，m
iy
，m
iz
表示作用力矩三个方向的分量；
[0150]
mi是连杆i的广义质量矩阵：
[0151][0152]
式中：mi为连杆i的质量；b
ic
为连杆i的质心相对于关节i的矢量；是b
ic
的反对称矩阵；1为3
×
3的单位矩阵；ii为连杆i相对关节i旋转轴线的惯量矩阵在基坐标系{0}下的表示。
[0153]
是关节i的线加速度矢量与角加速度矢量：
[0154][0155]bi
是代表包括向心力和陀螺力矩的交叉项矩阵：
[0156][0157]
代表关节i 1相对于关节i的力与力矩的转换矩阵，表达如下：
[0158][0159]
式中：是bi的反对称矩阵；转换矩阵是对机械臂构型的一种描述，即：在给定机械臂的情况下，该矩阵仅受到各杆件尺寸的影响。
[0160]
由于，编号为n 1的力和力矩为作用在机械臂末端点上的力和力矩，即：
[0161]
[0162]
式中：表示末端的作用力矩阵，表示末端的作用力矩矩阵，f
tipx
，f
tipy
，f
tipz
表示三个方向的作用力，m
tipx
，m
tipy
，m
tipz
表示三个方向的作用力矩；
[0163]wn 1
由混凝土3d打印单元以及打印单元储料桶内混凝土材料在内负载的重力与惯性力和力矩组成，因为由表2可知机械臂末端额定有效负载为90kg，所以将机械臂末端所受负载假设为恒定质量90kg。在建筑3d打印任务中，混凝土3d打印单元出料口挤出混凝土材料时相互作用力非常小，所以此项不予考虑。
[0164]
表2
[0165]
参数质量(kg)长度(m)相对关节的转动惯量(kg
·
m2)连杆15120.350diag(15.4 30.9 30.2)连杆22491.350diag(80 150.6 150.8)连杆32231.400diag(125.9 45.9 125.9)额定负载90——
[0166]
通过公式(27)反向迭代计算各关节约束力/力矩，该方法首先从编号为n 1的机械臂末端执行器计算到机械臂最后一个关节n，直至计算到机械臂关节1。
[0167]
为进一步求解机械臂各关节驱动力矩，因为上一节中各关节约束力/约束力矩wi是机械臂基坐标系{0}下的表示，得将wi中的关节约束力矩mi利用坐标系间的齐次变换转化为局部坐标系{i}下的表示，再取z轴方向分量得到各关节绕关节转轴的驱动力矩τi，即：
[0168][0169]
式中：mi＝[m
ix m
iy m
iz
]，m
ix
，m
iy
和m
iz
表示沿x,y,z三个方向上的约束力矩；zi＝[0 0 1]
t
，为连杆坐标系{i}在机械臂基坐标系下的方向。
[0170]
如图9所示为在adams软件中建立adams虚拟样机，可在adams软件中验证所建立的动力学模型的正确性。
[0171]
步骤(4)将基于步骤(3)中动力学理论，确定机械臂末端变形量作为机械臂刚度评价指标。
[0172]
建筑3d打印移动机械臂执行任务过程中机械臂末端受到混凝土3d打印单元以及储料桶内混凝土材料在内较大负载的作用，机械臂各个关节驱动系统需要输出相应的力矩与该末端载荷平衡，同时，驱动系统本身由于反馈系统的增益有限而具有的刚度是有限的，即各个关节刚度有限；系统执行建筑3d打印任务过程，末端会产生相应的变形，机械臂末端抵御负载变形的能力称为机械臂笛卡尔空间刚度。
[0173]
为了精确计算出工业机器人的刚度特性，建立合理的刚度模型十分重要，本发明在进行机械臂刚度建模时，做出如下假设：
[0174]
(1)混凝土3d打印单元刚体假设：建筑3d打印移动机械臂系统在执行打印任务过程中，所有的变形均来自工业机器人本体，末端执行器变形忽略不计；
[0175]
(2)工业机器人关节柔性假设：机械臂受力后产生的变形全部是由于关节发生弹性扭转引起的，机械臂刚度模型旨在建立关节空间刚度矩阵与笛卡尔空间刚度矩阵的映射关系，此时认为机械臂连杆是刚性的，末端受力后连杆不发生形变；
[0176]
(3)弹性变形假设：影响机械臂关节刚度的因素主要来源于关节处的伺服驱动系
统和传动结构刚度，将两者的影响视作关节处发生的线弹性扭转变形。
[0177]
机械臂n个关节的驱动力矩组成n维矢量τ＝[τ1,τ2,
…
,τn]
t
，n个关节对应的角度值记为θ＝[θ1,θ2,
…
,θn]
t
，在传递驱动力矩时，每个关节都会产生变形，δθ是各关节角度变形矢量，记为δθ＝[δθ1,δθ2,
…
,δθn]
t
；每个关节的刚度用弹簧系数k
θi
表示，则关节所受力矩可表示为：
[0178]
τi＝k
θi
·
δθi(i＝1,2,
…
,n)
ꢀꢀꢀ
(36)
[0179]
公式(36)可写成矩阵形式：
[0180]
τ＝k
θ
·
δθ
ꢀꢀ
(37)
[0181]
式中：k
θ
＝diag(k
θ1
,k
θ2
,
…kθn
)，是机械臂关节空间的刚度矩阵，反映了每个关节的刚度系数。
[0182]
再由微分运动关系可得：
[0183]
δx＝j(θ)
·
δθ
ꢀꢀꢀꢀꢀꢀ
(38)
[0184]
式中：j(θ)是机械臂的雅可比矩阵。
[0185]
联立公式(37)和公式(38)求解得出机械臂末端产生的广义变形量：
[0186]
δx＝j(θ)
·
(k
θ
)-1
·
τ
ꢀꢀꢀꢀꢀ
(39)
[0187]
式中：δx为机械臂末端产生的广义变形量。
[0188]
由于机械臂的前三个关节决定机械臂末端的空间位置，后三个关节决定机械臂末端的空间姿态，在机械臂执行任务过程中，末端载荷引起的末端角位移变形量很小，通常情况下可以忽略其作用，且机械臂第四、五、六关节对末端刚度的影响相比于前三关节要小的多，即可以将机械臂第四、五、六关节看作是刚体；综上所述，影响建筑3d打印过程中机械臂末端定位精度和打印质量的主要因素是机械臂前三关节产生的末端线位移变形。
[0189]
由于实验条件的限制，为了继续展开进一步的研究工作，本文将参考hoai nam huynh等学者在2020年对本项目中同型号工业机器人进行的关节刚度辨识实验结果，将其作为本文中kuka kr90 r3100 extra型工业机器人的各个关节刚度值；kuka kr90 r3100 extra型工业机器人前三关节刚度值：k
θ
＝diag(5.0
×
106,6.9
×
106,2.8
×
106)，单位为：n
·
m/rad。
[0190]
对于建筑3d打印工况中，机械臂末端受到的负载包含混凝土3d打印单元以及打印单元储料桶内混凝土材料，总体负载较大，在本文中将机械臂末端受到的总负载设为该型号机械臂的额定负载，即90kg，且始终竖直向下；下面将考虑机械臂末端x、y、z方向的变形情况与末端所受负载方向之间的关系，为了定量分析两者之间的关系，对四种类型的墙体分别随机选取移动平台在某一站位下三个建筑3d打印的任务点，编号为1～12，其任务点对应的各关节角度见表3，并根据墙体轨迹计算出以下任务点对应的末端笛卡尔空间速度，通过运动微分关系式计算出对应的各关节角速度；
[0191]
表3
[0192][0193]
基于各任务点对应的各关节角度、各关节角速度运用上述推到的机械臂动力学模型求解机械臂各关节力矩值。
[0194]
其中根据墙体轨迹计算出以下任务点对应的末端笛卡尔空间速度步骤如下：
[0195]
机械臂末端执行器混凝土3d打印单元在操作空间中的运动速度是一个矢量，在进行墙体3d打印过程中，混凝土3d打印单元出料口始终竖直向下，运行过程仅包含在操作空间中的平移运动，并不包含旋转运动，故混凝土3d打印单元末端的角速度大小为0；每一个任务点对应的末端线速度单位方向向量可通过运行轨迹前后两点在墙体坐标系{w}下坐标值计算得出，并将其转换为在机械臂基坐标系{0}下的表达如下：
[0196]
(1)当运行轨迹前后两点在同一高度时，
[0197][0198]
式中：0pv表示在机械臂基坐标系{0}下末端的线速度；表示墙体坐标系{w}相对于机械臂基坐标系{0}的旋转矩阵；wpv表示在墙体坐标系{w}下的末端线速度；表示机械臂基坐标系{0}相对于墙体坐标系{w}的旋转矩阵；y
i 1
表示轨迹点i 1的y轴坐标分量；yi表示轨迹点i的y轴坐标分量；x
i 1
表示轨迹点i 1的x轴坐标分量；xi表示轨迹点i的x轴坐标分量；
[0199]
(2)当运行轨迹前后两点不在同一高度时，此时恰好打印结束某一层，混凝土3d打印单元竖直方向抬升，继续下一层的打印作业。
[0200]0pv＝[0 0 1]
t
ꢀꢀꢀꢀꢀꢀ
(41)
[0201]
混凝土3d打印单元出料口始终竖直向下，运行过程仅包含在操作空间中的平移运动，并不包含旋转运动，故混凝土3d打印单元末端的角速度大小为0。考虑建筑3d打印材料的粘稠度、工作环境以及天气等因素，其末端在操作空间中的线速度大小设置为0.1m/s。根据上述求出的运动方向，以及线速度大小即可以得到任务点对应的末端笛卡尔空间速度；
[0202]
图10所示为上表各任务点利用公式(39)模型计算得出的机械臂末端产生的广义变形量，其包含x、y、z方向的变形量大小与方向；比较发现，机械臂末端x、y、z方向都有不同
程度的变形，但具体比较各个方向变形量的大小，可以得出机械臂末端沿x、y轴方向的变形量大小均在0～0.5mm之间，该变形量极小，对3d打印的实现影响不大，可忽略不计；因此，将以执行建筑3d打印过程中机械臂末端沿z轴方向的变形量作为刚度性能评价指标，记为μz。
[0203]
μz＝z
·
δx＝z
·
j(θ)
·
(k
θ
)-1
·
τ
ꢀꢀꢀꢀꢀꢀꢀ
(42)
[0204]
式中：z＝[0 0 1]，μz单位为：mm。
[0205]
步骤(5)根据步骤(4)得到的机械臂刚度评价指标做为优化目标，以移动平台相对所需打印墙体在可达站位空间内的站位做为设计变量，建立单次作业的优化函数，利用遗传算法进行合理的站位优化。
[0206]
针对民用建筑中常见的四种类型墙体基于遗传算法在整个可达站位空间中搜寻移动平台的最优站位；如图11所示为本发明提供的民用建筑实例的平面结构及尺寸示意图，将其大致划分为如图12所示的四种常见墙体类型，并且移动平台相对于墙体站位如图12所示。下面将针对每一类型墙体求解出一个移动平台的最优站位。
[0207]
基于遗传算法在整个站位空间中搜寻移动平台的最优站位，最优站位求解过程是利用matlab遗传算法工具箱，matlab遗传算法工具箱扩展了matlab在处理优化问题方面的能力，可以用于处理传统的优化技术难以解决的问题，还可以用于解决目标函数较复杂的问题。利用matlab遗传算法工具箱在系统执行3d打印任务移动平台的整个站位空间内搜寻移动平台的最优站位，该最优站位能使得系统在执行单次3d打印任务时所有任务点中机械臂末端沿z轴方向的最大变形量最小，建筑3d打印系统执行任务过程中整体的末端作业精度得到显著提高。
[0208]
下面是优化模型的建立过程。
[0209]
设计变量：将移动平台坐标系{m}在墙体坐标系{w}下的位姿[x,y,θ]作为最优站位的参数，也就是该优化模型的设计变量；
[0210]
优化目标：优化目的选择合理的站位使得单次建筑3d打印工作所有任务点中机械臂末端沿z轴方向的最大变形量在移动平台理论上的可达站位空间内达到最小；因为机械臂末端沿z轴方向的变形始终竖直向下，即μz在机械臂基坐标系下的表示总是负值，故取其绝对值表示变形量大小。
[0211]
约束条件：该问题的约束条件分为界限约束和性态约束。界限约束上述得出的移动平台理论上可达站位空间，即公式(7)；性态约束即满足系统执行建筑3d打印过程中混凝土3d打印单元的位置和姿态约束，即保证混凝土3d打印单元坐标系{t}的原点与墙体轨迹重合，且其出料口姿态始终竖直向下直至完成任务。也就是需要保证墙体所有的任务点均能通过公式(6)系统运动学方程与机械臂逆运动学方程求解出至少一组解。
[0212]
综上所述表明该问题是一个多变量、单目标、有约束的优化问题，遗传算法是一种全局搜索的方法，更适合工程优化问题的求解。使用遗传算法在整个站位空间中搜寻移动平台的最优站位，能够跳出局部最优而找到全局最优点，而且遗传算法允许使用非常复杂的目标函数，并对变量的变化范围可以加以限制。
[0213]
基于遗传算法在整个站位空间中搜寻移动平台的最优站位的过程是利用matlab遗传算法工具箱，其扩展了matlab在处理优化问题方面的能力，使结果更具有准确性与高效性，解决了依靠人工根据现场情况进行移动平台的站位布置，而使得移动平台站位移动次数多、整体打印精度低和效率低等问题。
[0214]
根据以上所述求解出了四种类型墙体的移动平台坐标系{m}在墙体坐标系{w}下的最优位姿[x,y,θ]。计算出的结果如下：
[0215]
一字墙体：[1878.9mm,-1252.4mm,162.79
°
]；
[0216]
直角墙体：[2135.8mm,-74.87mm,-170.9
°
]；
[0217]
圆弧墙体：[2171.8mm,679.1mm,-163.11
°
]；
[0218]
t型墙体：[891.4mm,-1875.8mm,90.8
°
]。
[0219]
如图13、图14、图15、图16为分别针对一字墙体、直角墙体、圆弧墙体、t形墙体作业任务建筑3d打印移动机械臂基于刚度性能提升的最优站位示意图。图中坐标值代表任务点在笛卡尔空间中的位置，点的颜色代表机械臂末端沿z轴方向的变形量大小，右侧竖直图例表明点的颜色对应的数值大小。且表4可以看出，四种墙体优化后机械臂末端沿z轴方向最大的变形量均很小。由此表明移动平台相对于墙体的站位经过遗传算法优化后，系统在执行墙体打印任务时，整体的机械臂末端作业精度得到了显著的提高，由此优化了建筑3d打印出的墙体质量。
[0220]
表4
[0221]