工业机械臂建模
这个例子展示了一个工业机器人手臂动力学的灰盒建模。机械臂采用非线性三质量柔性模型描述，如图1所示。这个模型是理想化的，在某种意义上，运动被假设是围绕一个不受重力影响的轴。为了简单起见，建模也是在齿轮传动比r = 1的情况下进行的，然后通过与真实齿轮传动比直接缩放来获得真实的物理参数。下面详细介绍的建模和识别实验是基于发表在
E. Wernholt和S. Gunnarsson。物理参数化机器人模型的非线性辨识。第14届国际会计师联合会系统识别研讨会预印本，143-148页，澳大利亚纽卡斯尔，2006年3月。
图1:工业机械臂示意图。
机器人手臂建模
机器人的输入是电机产生的施加力矩u(t)=tau(t)，电机由此产生的角速度y(t) = d/dt q_m(t)为测量输出。齿轮箱后和臂结构末端的质量角位置q_g(t)和q_a(t)是不可测量的。齿轮箱内的柔性由一个非线性弹簧建模，由弹簧扭矩tau_s(t)描述，该弹簧位于电机和第二个质量之间，而位于最后两个质量之间的“线性”弹簧建模臂结构中的柔性。系统的摩擦主要作用于第一质量，这里用非线性摩擦力矩tau_f(t)来表示。
状态介绍:
(x1 (t)) (q_m (t) - q_g (t)) (x2 (t)) (q_g (t) - q_a (t) x (t) = (x3 (t)) = (d / dt q_m (t)) (x4 (t)) (d / dt q_g (t)) (x5 (t)) (d / dt q_a (t))
对三个质量施加力矩平衡，得到非线性状态空间模型结构如下:
d / dt x1 (t) = x3 (t) - x4 (t) d / dt x2 (t) = x4 (t) - x5 (t) d / dt x3 (t) = 1 / J_m * (-tau_s (t) - d_g * (x3 (t) x4 (t))——tau_f (t) + u (t)) d / dt x4 (t) = 1 / J_g * (tau_s (t) + d_g * (x3 (t) x4 (t) - k_a * x2 (t) - d_a * (x4 (t) x5 (t))) d / dt x5 (t) = 1 / J_a * (k_a * x2 (t) + d_a * (x4 (t) x5 (t)))
Y (t) = x3(t)
式中，J_m、J_g、J_a分别为电机、变速箱和臂结构的惯性矩，d_g、d_a为阻尼参数，k_a为臂结构的刚度。
齿轮箱摩擦力矩tau_f(t)的建模包括了实践中遇到的许多摩擦现象，其中包括所谓的库仑摩擦和斯特里贝克效应:
tau_f (t) =阵线* x3 (t) + (Fc + Fcs *双曲正割(α* x3 (t))) *双曲正切(β* x3 (t))
其中Fv和Fc是粘性系数和库仑摩擦系数，Fcs和alpha是反映Stribeck效应的系数，beta是用于获得x3(t)从负速度平滑过渡到正速度的参数。(类似的方法，但基于稍微不同的模型结构，用于描述速度和摩擦力矩/力之间的静态关系，将在名为idnlgreydemo5的教程中进一步讨论:“摩擦的静态建模”。)
弹簧的扭矩tau_s(t)，假设用x1(t)中没有平方项的三次多项式来描述:
Tau_s (t) = k_g1*x1(t) + k_g3*x1(t)^3
式中k_g1和k_g3为变速箱弹簧的两个刚度参数。
在Wernholt和Gunnarsson论文中讨论的其他类型的辨识实验中，有可能辨识出总的转动惯量J = J_m+J_g+J_a。由此，我们可以引入未知的比例因子a_m和a_g，并执行以下重新参数化:
J_m = J*a_m J_g = J*a_g J_a = J*(1-a_m-a_g)
其中(如果J已知)只需要估计a_m和a_g。
总之，这给出了以下状态空间结构，涉及13个不同的参数:Fv, Fc, Fcs, alpha, beta, J, a_m, a_g, k_g1, k_g3, d_g, k_a和d_a。(根据定义，我们还使用了sech(x) = 1/cosh(x)的事实。)
tau_f (t) =阵线* x3 (t) + (Fc + Fcs / cosh(α* x3 (t))) *双曲正切(β* x3 (t)) tau_s (t) = k_g1 * x1 (t) + k_g3 * x1 (t) ^ 3
d / dt x1 (t) = x3 (t) - x4 (t) d / dt x2 (t) = x4 (t) - x5 (t) d / dt x3 (t) = 1 / (J * a_m) * (-tau_s (t) - d_g * (x3 (t) x4 (t))——tau_f (t) + u (t)) d / dt x4 (t) = 1 / (J * a_g) * (tau_s (t) + d_g * (x3 (t) x4 (t) - k_a * x2 (t) - d_a * (x4 (t) x5 (t))) d / dt x5 (t) = 1 / (J (1-a_m-a_g)) * (k_a * x2 (t) + d_a * (x4 (t) x5 (t)))
Y (t) = x3(t)
IDNLGREY机器人手臂模型对象
上面的模型结构被输入到一个名为robotarm_c.c的C mex文件中，具有如下的状态和输出更新功能(整个文件可以通过命令“type robotarm_c.c”查看)。在状态更新函数中，注意我们在这里使用了两个中间双变量，一方面增强了方程的可读性，另一方面提高了执行速度(taus在方程中出现了两次，但只计算了一次)。
/*状态方程。*/ void compute_dx(double *dx, double *x, double *u, double **p){/*模型参数和中间变量的声明。* /双*阵线,* Fc, * Fcs, *α,β*,*,*,*,* kg1, * kg3, * dg, * ka, *达;双tauf, taus;/*中间变量。* /
/*检索模型参数。*/ Fv = p[0];/*粘性摩擦系数。*/ Fc = p[1];/*库仑摩擦系数。*/ Fcs = p[2];/* Striebeck摩擦系数。*/ alpha = p[3];/* Striebeck平滑系数。*/ beta = p[4]; /* Friction smoothness coefficient. */ J = p[5]; /* Total moment of inertia. */ am = p[6]; /* Motor moment of inertia scale factor. */ ag = p[7]; /* Gear-box moment of inertia scale factor. */ kg1 = p[8]; /* Gear-box stiffness parameter 1. */ kg3 = p[9]; /* Gear-box stiffness parameter 3. */ dg = p[10]; /* Gear-box damping parameter. */ ka = p[11]; /* Arm structure stiffness parameter. */ da = p[12]; /* Arm structure damping parameter. */
/*确定中间变量。*/ /* tauf:齿轮摩擦力矩。(sech(x) = 1/cosh(x)!*/ /* taus:弹簧扭矩。* / tauf =阵线[0]* x [2] + (Fc [0] + Fcs [0] / (cosh(α[0]* x[2]))) *双曲正切(β[0]* x [2]);Taus = kg1[0]*x[0]+kg3[0]*pow(x[0]，3);
/* x[0]:电机与变速箱的转速差。*/ /* x[1]:变速箱与手臂的转速差。*/ /* x[2]:电机转速。*/ /* x[3]:变速箱后转速。*/ /* x[4]:机械臂转速。*/ dx[0] = x[2]-x[3];Dx [1] = x[3]-x[4];dx [2] = 1 / (J [0] * [0]) * (-taus-dg [0] * (x [2] - x [3]) -tauf + u [0]);dx [3] = 1 / (J [0] * ag)[0]) *(τ+ dg [0] * (x [2] - x [3]) - ka [0] * x [1] - da [0] * (x [3] - x [4]));dx [4] = 1 / (J[0] *(1.0是[0]ag) [0])) * (ka [0] * x [1] + da [0] * (x [3] - x [4])); }
/*输出方程。*/ void compute_y(double y[]， double x[]) {/* y[0]:电机转速。*/ y[0] = x[2];｝
下一步是创建一个反映建模情况的IDNLGREY对象。这里应该指出，为机械臂找到合适的初始参数值需要一些额外的努力。在Wernholt和Gunnarsson的论文中，这是在前面的两个步骤中进行的，其中使用了其他模型结构和识别技术。下面使用的初始参数值是这些识别实验的结果。
文件名=“robotarm_c”；描述模型结构的文件。Order = [1 1 5];%模型订单[ny nu nx]。参数= [0.00986346744839 0.74302635727901 ..．.3.98628540790595 - 3.24015074090438.．.0.79943497008153 - 0.03291699877416.．.0.17910964111956 - 0.61206166914114.．.20.59269827430799 - 0.00000000000000.．.0.06241814047290 - 20.23072060978318.．.0.00987527995798];初始参数向量。InitialStates = 0 (5,1);初始状态。Ts = 0;%时间连续系统。nlgr = idnlgrey(文件名，顺序，参数，InitialStates, Ts，.．.“名字”，机器人手臂的，.．.“InputName”，“应用电机扭矩”，.．.“InputUnit”，“纳米”，.．.“OutputName”，“电机角速度”，.．.“OutputUnit”，“rad / s”，.．.“TimeUnit”，“年代”）;
为了便于记帐，提供了各州的名称和单位:
NLGR = setinit(NLGR，“名字”, {“电机和变速箱之间的角位置差”.．.变速箱和手臂之间的角度位置差.．.“电机角速度”.．.“变速箱角速度”.．.“机械臂角速度”} ');NLGR = setinit(NLGR，“单位”, {rad的rad的“rad / s”“rad / s”“rad / s”});
参数名也被详细指定。此外，建模是这样做的，所有参数都应该是正的，即每个参数的最小值应该设置为0(因此约束估计将稍后执行)。在Wernholt和Gunnarsson的论文中，我们也考虑了前6个参数，即Fv, Fc, Fcs, alpha, beta和J，它们非常好，不需要估计。
NLGR = setpar(NLGR，“名字”, {Fv:粘性摩擦系数.．.% 1。“Fc:库仑摩擦系数”.．.% 2。“Fcs: Striebeck摩擦系数”.．.% 3。'alpha: Striebeck平滑系数'.．.% 4。“beta:摩擦平滑系数”.．.% 5。J:总转动惯量.．.% 6。'a_m:电机转动惯量比例因子'.．.% 7。'a_g:变速箱转动惯量比例因子'.．.% 8。'k_g1:变速箱刚度参数1'.．.% 9。'k_g3:变速箱刚度参数3'.．.% 10。d_g:变速箱阻尼参数.．.% 11。“k_a:臂结构刚度参数”.．.% 12。d_a:臂结构阻尼参数.．.% 13。});NLGR = setpar(NLGR，“最低”num2cell(0(大小(nlgr“np”), 1)));%所有参数>= 0!为帕诺= 1:6修复前六个参数。nlgr.Parameters (parno)。Fixed = true;结束
目前进行的建模步骤给我们留下了一个初始的机械臂模型，其属性如下:
礼物(nlgr);
nlgr =由'robotarm_c' (MEX-file)定义的连续时间非线性灰盒模型:dx/dt = F(t, u(t)， x(t)， p1，…， p13) y(t) = H(t, u(t)， x(t)， p1，…，p13) + e(t) with 1 input(s), 5 state(s), 1 output(s), and 7 free parameter(s) (out of 13). Inputs: u(1) Applied motor torque(t) [Nm] States: Initial value x(1) Angular position difference between the motor and the gear-box(t) [rad] xinit@exp1 0 (fixed) in [-Inf, Inf] x(2) Angular position difference between the gear-box and the arm(t) [rad] xinit@exp1 0 (fixed) in [-Inf, Inf] x(3) Angular velocity of motor(t) [rad/s] xinit@exp1 0 (fixed) in [-Inf, Inf] x(4) Angular velocity of gear-box(t) [rad/s] xinit@exp1 0 (fixed) in [-Inf, Inf] x(5) Angular velocity of robot arm(t) [rad/s] xinit@exp1 0 (fixed) in [-Inf, Inf] Outputs: y(1) Angular velocity of motor(t) [rad/s] Parameters: Value p1 Fv : Viscous friction coefficient 0.00986347 (fixed) in [0, Inf] p2 Fc : Coulomb friction coefficient 0.743026 (fixed) in [0, Inf] p3 Fcs : Striebeck friction coefficient 3.98629 (fixed) in [0, Inf] p4 alpha: Striebeck smoothness coefficient 3.24015 (fixed) in [0, Inf] p5 beta : Friction smoothness coefficient 0.799435 (fixed) in [0, Inf] p6 J : Total moment of inertia 0.032917 (fixed) in [0, Inf] p7 a_m : Motor moment of inertia scale factor 0.17911 (estimated) in [0, Inf] p8 a_g : Gear-box moment of inertia scale factor 0.612062 (estimated) in [0, Inf] p9 k_g1 : Gear-box stiffness parameter 1 20.5927 (estimated) in [0, Inf] p10 k_g3 : Gear-box stiffness parameter 3 0 (estimated) in [0, Inf] p11 d_g : Gear-box damping parameter 0.0624181 (estimated) in [0, Inf] p12 k_a : Arm structure stiffness parameter 20.2307 (estimated) in [0, Inf] p13 d_a : Arm structure damping parameter 0.00987528 (estimated) in [0, Inf] Name: Robot arm Status: Created by direct construction or transformation. Not estimated. More information in model's "Report" property. Model Properties
输入输出数据
实验机器人收集了大量真实世界的数据集。为了使机器人保持在其工作点附近，但也出于安全原因，数据是使用实验反馈控制安排收集的，随后允许脱机计算关节控制器的参考信号。
在这个案例研究中，我们将继续讨论四个不同的数据集，一个用于估计，其余的用于验证。在每种情况下，使用大约10秒持续时间的周期性激励信号为控制器生成参考速度。所选采样频率为2 kHz(采样时间Ts = 0.0005秒)。对于数据集，使用了三种不同类型的输入信号:(ue:估计数据集的输入信号;Uv1, uv2, uv3:三个验证数据集的输入信号)
ue, uv1:多正弦信号，频率区间为1-40 Hz，振幅平坦，峰值为16 rad/s。多正弦信号叠加在振幅为20 rad/s、截止频率为1hz的滤波方波上。
uv2:类似于ue和uv1，但没有方波。
uv3:多正弦信号(正弦波的和)，频率为0.1、0.3和0.5 Hz，峰值为40 rad/s。
让我们加载可用数据，并将所有四个数据集放入一个IDDATA对象z中:
负载(fullfile (matlabroot“工具箱”，“识别”，“iddemos”，“数据”，“robotarmdata”));Z = iddata({ye yv1 yv2 yv3}， {ue uv1 uv2 uv3}， 0.5e-3，“名字”，机器人手臂的）;z.InputName =“应用电机扭矩”；z.InputUnit =“纳米”；z.OutputName =“电机角速度”；z.OutputUnit =“rad / s”；z.ExperimentName = {“估计”“验证1”《验证2》'验证3 '};z.Tstart = 0;z.TimeUnit =“年代”；礼物(z);
th =包含4个实验的时域数据集。实验样本样本时间估计19838 0.0005验证1 19838 0.0005验证2 19838 0.0005验证3 19838 0.0005名称:机械臂输出单位(如指定)电机角速度rad/s输入单位(如指定)电机施加扭矩Nm
下图是四个实验中使用的输入输出数据。
图(“名字”[z。的名字':输入输出数据'),.．.“DefaultAxesTitleFontSizeMultiplier”, 1.．.“DefaultAxesTitleFontWeight”，“正常”，.．.“位置”，[100 100 900 600]);为I = 1:z。Ne zi = getexp(z, i);次要情节(z。Ne, 2,2 *i-1);%的输入。情节(zi.u);标题([z。ExperimentName{我}“:”子。InputName {1}),“FontWeight”，“正常”）;如果(i < z.Ne) xlabel(”）;其他的包含([z。域“(”子。TimeUnit“)”]);结束次要情节(z。Ne, 2,2 *i);%输出。情节(zi.y);标题([z。ExperimentName{我}“:”子。OutputName {1}),“FontWeight”，“正常”）;如果(i < z.Ne) xlabel(”）;其他的包含([z。域“(”子。TimeUnit“)”]);结束结束
图2:实验机械臂的输入输出测量数据。
机械臂初始模型的性能
最初的机械臂模型有多好?让我们使用COMPARE来模拟模型输出(对于所有四个实验)，并将结果与相应的测量输出进行比较。对于所有四个实验，我们都知道前两个状态的值为0(固定)，而其余三个状态的值初始设置为开始时间的测量输出(不固定)。然而，在默认情况下，COMPARE估计所有的初始状态，如果z持有4个不同的实验，这意味着4*5 = 20个初始状态需要估计。即使在确定了前两个状态后，4*3 = 12个初始状态仍有待估计(如果遵循内部模型初始状态策略)。因为数据集相当大，这将导致冗长的计算，为了避免这种情况，我们使用PREDICT(如果初始状态作为初始状态结构传递)估计初始状态的4*3个自由成分，但将估计限制在可用数据的前10:th。然后，我们指示COMPARE使用得到的5 × 4初始状态矩阵X0init，而不执行任何初始状态估计。
zred = z(1:round(zi.N/10));NLGR = setinit(NLGR，“固定”，{真真假假假});X0 = nlgr.InitialStates;(X0。值]=交易(0(1、4),0(1、4),(你们(1)yv1 (1) yv2 (1) yv3(1)]。.．.[你们(1)yv1 (1) yv2 (1) yv3(1)]、[你们(1)yv1 (1) yv2 (1) yv3 (1)]);[~， X0init] = predict(zred, nlgr， []， X0);NLGR = setinit(NLGR，“价值”， num2cell(X0init(:， 1)));clf compare(z, nlgr， []， compareOptions(“InitialCondition”X0init));
图3:机械臂初始模型的测量输出与仿真输出的比较。
可以看出，初始机器人手臂模型的性能还算不错或相当不错。三种数据集的拟合率分别为ye和yv1约79%，yv2约37%，yv3约95%。请注意，与yv2相比，ye/yv1的高拟合在很大程度上是由于初始模型捕捉方波的能力，而多正弦部分捕捉得不是同样好。我们还可以看看四个实验的预测误差:
pe(z, nlgr, peOptions(“InitialCondition”X0init));
图4:机械臂初始模型的预测误差。
参数估计
现在让我们尝试通过估计z(估计数据集)的第一次实验的7个自由模型参数和3个自由初始状态来改善初始机器人手臂模型的性能。这个估计将花费一些时间(通常是几分钟)。
nlgr = nlgreyest(nlgr, getexp(z, 1)， nlgreyestOptions(“显示”，“上”));
估计机械臂模型的性能
再次使用COMPARE来评估估计的机械臂模型的性能。我们还在这里指示COMPARE不执行任何初始状态估计。对于第一个实验，我们用NLGREYEST估计的初始状态替换猜测的初始状态，对于其余三个实验，我们使用PREDICT来估计基于减少的IDDATA对象zred的初始状态。
X0init(:， 1) = cell2mat(getinit(nlgr，“价值”));X0 = nlgr.InitialStates;(X0。值]=交易(0 (1,3),0 (1,3),[yv1 (1) yv2 (1) yv3(1)]。.．.[yv1(1) yv2(1) yv3(1)]， [yv1(1) yv2(1) yv3(1)]);[yp X0init(: 2:4)] =预测(getexp (zr, 2:4) nlgr, [], X0);clf compare(z, nlgr， []， compareOptions(“InitialCondition”X0init));
图5:测量输出与估计机器人手臂模型的模拟输出的比较。
比较图显示了更好的拟合方面的改进。对于ye和yv1，现在的拟合度约为85%(以前:79%)，yv2约为63%(以前:37%)，yv3略低于95.5%(以前:也略低于95.5%)。第二个验证数据集的改进最为明显，其中应用了没有任何方波的多正弦信号作为输入。然而，估计模型跟踪ye和yv1的多正弦部分的能力也有了很大的提高(但这没有反映在拟合图中，因为这些更大程度上受到对方波拟合的影响)。预测误差图还显示，残差现在总体上小于初始机械臂模型:
图;pe(z, nlgr, peOptions(“InitialCondition”X0init));
图6:估计的机械臂模型的预测误差。
我们通过文本总结估计的机械臂模型的各种特性来结束案例研究。
礼物(nlgr);
nlgr =由'robotarm_c' (MEX-file)定义的连续时间非线性灰盒模型:dx/dt = F(t, u(t)， x(t)， p1，…， p13) y(t) = H(t, u(t)， x(t)， p1，…，p13) + e(t) with 1 input(s), 5 state(s), 1 output(s), and 7 free parameter(s) (out of 13). Inputs: u(1) Applied motor torque(t) [Nm] States: Initial value x(1) Angular position difference between the motor and the gear-box(t) [rad] xinit@exp1 0 (fixed) in [-Inf, Inf] x(2) Angular position difference between the gear-box and the arm(t) [rad] xinit@exp1 0 (fixed) in [-Inf, Inf] x(3) Angular velocity of motor(t) [rad/s] xinit@exp1 -19.0684 (estimated) in [-Inf, Inf] x(4) Angular velocity of gear-box(t) [rad/s] xinit@exp1 -22.2346 (estimated) in [-Inf, Inf] x(5) Angular velocity of robot arm(t) [rad/s] xinit@exp1 -23.4238 (estimated) in [-Inf, Inf] Outputs: y(1) Angular velocity of motor(t) [rad/s] Parameters: Value Standard Deviation p1 Fv : Viscous friction coefficient 0.00986347 0 (fixed) in [0, Inf] p2 Fc : Coulomb friction coefficient 0.743026 0 (fixed) in [0, Inf] p3 Fcs : Striebeck friction coefficient 3.98629 0 (fixed) in [0, Inf] p4 alpha: Striebeck smoothness coefficient 3.24015 0 (fixed) in [0, Inf] p5 beta : Friction smoothness coefficient 0.799435 0 (fixed) in [0, Inf] p6 J : Total moment of inertia 0.032917 0 (fixed) in [0, Inf] p7 a_m : Motor moment of inertia scale factor 0.26648 0.000286925 (estimated) in [0, Inf] p8 a_g : Gear-box moment of inertia scale factor 0.647531 0.00019097 (estimated) in [0, Inf] p9 k_g1 : Gear-box stiffness parameter 1 20.0776 0.0263505 (estimated) in [0, Inf] p10 k_g3 : Gear-box stiffness parameter 3 24.1775 0.332324 (estimated) in [0, Inf] p11 d_g : Gear-box damping parameter 0.0305126 0.000282693 (estimated) in [0, Inf] p12 k_a : Arm structure stiffness parameter 11.7553 0.031113 (estimated) in [0, Inf] p13 d_a : Arm structure damping parameter 0.00283482 8.0518e-05 (estimated) in [0, Inf] Name: Robot arm Status: Termination condition: Change in cost was less than the specified tolerance.. Number of iterations: 10, Number of function evaluations: 11 Estimated using Solver: ode45; Search: lsqnonlin on time domain data "Robot arm". Fit to estimation data: 85.21% FPE: 9.503, MSE: 9.494 More information in model's "Report" property. Model Properties
结束语
系统识别技术广泛应用于机器人领域。“好的”机器人模型对于现代机器人控制概念至关重要，通常被认为是满足不断增长的速度和精度需求的必要条件。这些模型也是各种机器人诊断应用中的关键组件，其中模型用于预测与磨损相关的问题，并用于检测机器人故障的实际原因。
您也可以从以下列表中选择一个网站: