两体间干摩擦:基于多实验数据的参数估计
这个例子展示了如何利用多个实验数据估计非线性灰盒模型的参数。在两个固体之间表现干摩擦的系统将被用作讨论的基础。如图1所示，在这个系统中，一个物体是固定的，而另一个物体在外力的作用下在固定的物体上前后移动。
图1:二体系统的示意图。
模拟两体间的干摩擦
根据牛顿第三运动定律，运动物体的运动可以用:
f_t (t) = m*a(t) = m*d/dt v(t) = m*d²/dt²s(t)
式中F_tot(t)等于外生力F(t)减去两个物体接触所产生的摩擦力。假定摩擦力为滑动摩擦力F_s(t)和干摩擦力F_d(t)的总和。前者通常建模为速度的线性函数，即F_s(t) = k*v(t)，其中k为未知滑动摩擦参数。干摩擦，另一方面，是一个相当复杂的现象。在论文中:
A. Clavel, M. Sorine和Q. Zhang。钢板弹簧系统的建模与辨识。2001年，第三届IFAC汽车控制进展研讨会。
用常微分方程表示:
d/dt F_d(t) = -1/e*|v(t)|*F_d(t) + f/e*v(t)
式中，e和f分别为两个未知参数，尺寸分别为距离和力。
令输入信号u(t) = F(t) [N]，引入状态为:
x1(t) = s(t)运动体的位置[m]。x2(t) = v(t)运动物体的速度[m/s]。x3(t) = F_d(t)体间干摩擦力[N]。
模型参数为:
运动物体的质量[m]。k滑动摩擦力系数[kg/s]。e距离相关干摩擦参数[m]。f力相关干摩擦参数[N]。
我们得到如下状态空间模型结构:
d / dt x1 (t) = x2 (t) d / dt x2 (t) = 1 / m * (u (t) - k * x2 (t) - x3 (t)) d / dt x3 (t) = 1 / e * (- | x2 (t) | * x3 (t) + f * x2 (t))
Y (t) = x1(t)
这些方程被输入到一个C-MEX模型文件中，twobodies_c。其状态和输出更新方程compute_dx和compute_y如下:
/*状态方程。*/ void compute_dx(double *dx, double t, double *x, double *u, double **p, const mxArray *auxvar){/*检索模型参数。*/ double *m， *k， *e， *f;M = p[0];/*运动物体的质量。*/ k = p[1];/*滑动摩擦力系数。*/ e = p[2];/*距离相关干摩擦参数。*/ f = p[3]; /* Force-related dry friction parameter. */
/* x[0]:位置。*/ /* x[1]:速度。*/ /* x[2]:干摩擦力。*/ dx[1] = x[1];Dx [1] = (u[0]-k[0]*x[1]-x[2])/m[0];Dx [2] = (-fabs(x[1])*x[2]+f[0]*x[1])/e[0];}
/*输出方程。*/ void compute_y(双*y，双t，双*x，双*u，双**p, const mxArray *auxvar) {/* y[0]:位置。*/ y b[0] = x b[0];}
在编写了描述模型结构的文件之后，下一步是创建反映建模情况的IDNLGREY对象。我们还添加了关于模型结构的输入、输出、状态和模型参数的名称和单位的信息。注意，这里的Parameters和InitialStates被指定为向量，默认情况下，这意味着当调用NLGREYEST时，将估计所有模型参数而不估计初始状态向量。
文件名=“twobodies_c”;描述模型结构的文件。Order = [1 1 3];%模型订单[ny nu nx]。参数= [380;2200;0.00012;1900);%初始参数向量。InitialStates = [0;0;0);%初始状态。Ts = 0;%时间连续系统。nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts，…“名字”，“双体系统”，…“InputName”，“外生力量”，…“InputUnit”，“N”，…“OutputName”，“运动体位置”，…“OutputUnit”，“米”，…“TimeUnit”，“年代”);NLGR = setinit(NLGR);“名字”， {“运动体位置”…“运动物体的速度”…"物体间的干摩擦力"}）;NLGR = setinit(NLGR);“单位”， {“米”“米/秒”“N”}）;NLGR = setpar(NLGR，“名字”， {运动物体的质量…滑动摩擦力系数，…“与距离相关的干摩擦参数”…力相关干摩擦参数}）;NLGR = setpar(NLGR，“单位”， {“米”公斤/年代”“米”“N”}）;
输入输出数据
此时，我们加载可用的(模拟的)输入-输出数据。该文件包含来自三个不同(模拟)测试运行的数据，每个测试运行使用0.005秒的采样率(Ts)生成1000个噪声损坏的输入-输出样本。输入u(t)是作用在运动物体上的外生力[N]。在实验中，输入是对称锯齿形信号，所有实验的波形重复频率相同，但在不同的测试运行中，最大信号幅度不同。输出y(t)是运动物体的位置[m](相对于固定物体)。出于建模目的，我们创建一个IDDATA对象，其中包含三个不同的实验:
负载(fullfile (matlabroot“工具箱”，“识别”，“iddemos”，“数据”，“twobodiesdata”));Z = merge(iddata(y1, u1, 0.005)， iddata(y2, u2, 0.005)， iddata(y3, u3, 0.005));z.Name =“双体系统”;z.InputName = nlgr.InputName;z. inpuunit = nlgr. inpuunit;z.OutputName = nlgr.OutputName;z. outpuunit = nlgr. outpuunit;z.Tstart = 0;z.TimeUnit = nlgr.TimeUnit;
用于后续识别实验的输入输出数据显示在绘图窗口中。
图(“名字”[z。名字':输入输出数据']）;为I = 1:z。nezi = getexp(z, i);次要情节(z。Ne, 2,2 *i-1);%的输入。情节(子。SamplingInstants zi.InputData);标题([z。ExperimentName{我}“:”zi.InputName {1}));如果(i < z.Ne) xlabel(＇＇);其他的包含([z。域“(”子。TimeUnit“)”]）;结束轴(“紧”);次要情节(z。Ne, 2,2 *i);%输出。情节(子。SamplingInstants zi.OutputData);标题([z。ExperimentName{我}“:”zi.OutputName {1}));如果(i < z.Ne) xlabel(＇＇);其他的包含([z。域“(”子。TimeUnit“)”]）;结束轴(“紧”);结束
图2:二体系统的输入输出数据。
初始二体模型的性能
在估计四个未知参数之前，我们使用初始参数值对系统进行了仿真。我们使用默认的微分方程求解器(ode45)和默认的求解器选项。当只有两个输入参数调用时，COMPARE将估计完整的初始状态向量，在这种情况下，每个实验一个，即3个实验，每个3 × 1的状态向量意味着总共9个估计的初始状态。模拟和真实的输出显示在绘图窗口中，可以看出，拟合程度不错，但不如预期的那么好。
CLF比较(z, nlgr);
图3:初始两体模型真实输出与仿真输出的比较。
参数估计
为了改善拟合，接下来对四个参数进行估计。我们在估计阶段使用来自第一次和最后一次实验的数据，并保留来自第二次实验的数据用于纯粹的验证目的。
opt = nlgreyestOptions(“显示”，“上”);NLGR = nlgreyest(getexp(z， [13])， NLGR, opt);
估计二体模型的性能
为了研究估计模型的性能，最后对其进行了仿真。通过剪裁nlgr的初始状态结构数组，可以完全指定每个实验中估计的状态，例如COMPARE。让我们在这里定义并使用一个结构，其中初始状态x1(0)和x2(0)在实验1中估计，x2(0)在实验2中估计，x3(0)在实验3中估计。通过这种修改，测量输出和模型输出之间的比较显示在绘图窗口中。
nlgr。InitialStates = struct(“名字”getinit (nlgr“名字”），…“单位”getinit (nlgr“单位”），…“价值”， 0 (1,3)，“最低”， -Inf(1,3)，…“最大”， Inf(1,3)，“固定”，…{[假假真];[真假真];[真真假]});compare(z, nlgr, compareOptions)“InitialCondition”，“模型”));
图4:估计的两体模型的真实输出与模拟输出的比较。
特别令人感兴趣的是第二次实验数据的结果，这些数据没有用于参数估计。真实系统的动力学在所有实验中都得到了很好的模拟。估计的参数也与生成实验数据的参数相当接近:
disp (“真实估计参数向量”);
估计参数向量
Ptrue = [400;2 e3;0.0001;1700);流(' %10.5f %10.5f\n'[ptrue ';getpvec (nlgr) '));
400.00000 402.16624 2000.00000 1987.96042 0.00010 0.00010 1700.00000 1704.51209
最后使用PRESENT命令，我们可以得到关于估计模型的额外信息:
礼物(nlgr);
nlgr =由'twobodies_c' (MEX-file)定义的连续时间非线性灰盒模型:dx/dt = F(t, u(t)， x(t)， p1，…， p4) y(t) = H(t, u(t)， x(t)， p1，…，p4) + e(t) with 1 input(s), 3 state(s), 1 output(s), and 4 free parameter(s) (out of 4). Inputs: u(1) Exogenous force(t) [N] States: Initial value x(1) Position of moving body(t) [m] xinit@exp1 0 (estimated) in [-Inf, Inf] xinit@exp2 0 (estimated) in [-Inf, Inf] xinit@exp3 0 (fixed) in [-Inf, Inf] x(2) Velocity of moving body(t) [m/s] xinit@exp1 0 (fixed) in [-Inf, Inf] xinit@exp2 0 (estimated) in [-Inf, Inf] xinit@exp3 0 (fixed) in [-Inf, Inf] x(3) Dry friction force between the bodies(t) [N] xinit@exp1 0 (fixed) in [-Inf, Inf] xinit@exp2 0 (fixed) in [-Inf, Inf] xinit@exp3 0 (estimated) in [-Inf, Inf] Outputs: y(1) Position of moving body(t) [m] Parameters: Value p1 Mass of the moving body [m] 402.166 (estimated) in [-Inf, Inf] p2 Sliding friction force coefficient [kg/s] 1987.96 (estimated) in [-Inf, Inf] p3 Distance-related dry friction parameter [m] 0.000104534 (estimated) in [-Inf, Inf] p4 Force-related dry friction parameter [N] 1704.51 (estimated) in [-Inf, Inf] Name: Two body system Status: Model modified after estimation. More information in model's "Report" property. Model Properties
结论
本例描述了在执行IDNLGREY建模时如何使用多个实验数据。可以使用任意数量的实验，并且对于每个这样的实验，可以完全指定要在NLGREYEST、COMPARE、PREDICT等中估计哪些初始状态或状态。
额外的信息
有关使用System identification Toolbox™识别动态系统的更多信息，请访问系统识别工具箱产品信息页面。
