两物体间的干摩擦建模

使用牛顿第三运动定律，运动物体的运动由以下公式描述：

F_tot（t）=m*a（t）=m*d/dtv（t）=m*d^2/dt^2s（t）

式中，F_tot（t）等于外力F（t）减去两个物体之间接触产生的摩擦力。摩擦力假定为滑动摩擦力F_s（t）和干摩擦力F_d（t）之和。前者通常被建模为速度的线性函数，即F_s（t）=k*v（t），其中k是未知的滑动摩擦参数。另一方面，干摩擦是一种相当复杂的现象。论文中：

A.克拉维尔、M.索林和Q.张。“钢板弹簧系统的建模和识别”。在第三届IFAC汽车控制进展研讨会上，2001年。

它由一个常微分方程建模：

d/dt F_d(t) = -1/e*|v(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运动物体的质量[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.c。其状态和输出更新方程compute_dx和compute_y如下:

/ *状态方程。void compute_dx(double *dx, double *x, double *u, double **p, const mxArray *auxvar){/*检索模型参数。*/ *m， *k， *e， *f;m = p [0];/*运动物体的质量。* k = p;/*滑动摩擦力系数。*/ e = p;/*与距离相关的干摩擦参数。*/ f = p[3]; /* Force-related dry friction parameter. */

/*x[0]：位置.//*x[1]：速度.//*x[2]：干摩擦力.*/dx[0]=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（double*y，double t，double*x，double*u，double**p，const mxArray*auxvar）{/*y[0]：位置。*/y[0]=x[0]；}

在编写了描述模型结构的文件之后，下一步是创建反映建模情况的IDNLGREY对象。我们还添加了关于模型结构的输入、输出、状态和模型参数的名称和单位的信息。注意，这里的Parameters和InitialStates被指定为向量，这意味着在调用NLGREYEST时，所有模型参数和初始状态向量都将被估计。

文件名=“两个人”;%描述模型结构的文件。顺序=[13]；%型号订单[ny nu nx]。参数=[380；2200；0.00012；1900]；%初始参数向量。初始状态=[0；0；0]；%初始状态。Ts=0；%时间连续系统。nlgr=IDNLGRY（文件名、顺序、参数、初始状态、Ts、，...“姓名”,“两体系统”,...“InputName”,“外生力量”,...“输入单元”,“不”,...“OutputName”,“动体位置”,...“输出单元”,“我是,...“时间单位”,'s'); nlgr=setinit（nlgr，“姓名”, {“动体位置”...“运动物体的速度”...“物体之间的干摩擦力”})；nlgr=setinit（nlgr，“单位”, {“我是“m/s”“不”}); nlgr=setpar（nlgr，“姓名”, {“运动物体的质量”...“滑动摩擦力系数”,...“与距离相关的干摩擦参数”...“力相关干摩擦参数”}); nlgr=setpar（nlgr，“单位”, {“我是公斤/年代”“我是“不”});

输入输出数据

此时，我们装入可用的(模拟的)输入-输出数据。该文件包含来自三个不同(模拟)测试运行的数据，每个测试运行都包含1000个噪声损坏的输入-输出样本，使用0.005秒的采样率(Ts)生成。输入u(t)是作用于运动物体的外生力[N]。在实验中，输入是一个对称的锯齿形信号，其中波形重复频率在所有实验中是相同的，但最大信号幅度在测试运行之间是不同的。输出y(t)为移动体(相对于固定体)的位置[m]。出于建模的目的，我们创建一个IDDATA对象，其中包含三个不同的实验:

加载（完整文件）（matlabroot，“工具箱”,“识别”,“iddemos”,“数据”,“双体数据”))；z=合并（iddata（y1，u1，0.005）、iddata（y2，u2，0.005）、iddata（y3，u3，0.005））；z.名称=“两体系统”; z、 InputName=nlgr.InputName；z、 InputUnit=nlgr.InputUnit；z、 OutputName=nlgr.OutputName；z、 OutputUnit=nlgr.OutputUnit；z、 Tstart=0；z、 TimeUnit=nlgr.TimeUnit；

用于向前识别实验的输入-输出数据显示在绘图窗口中。

图(“姓名”[z。的名字“：输入输出数据”]);对于i = 1: z。nezi = getexp(z, i); / /输入z次要情节(z。不2张2 *);%的输入。绘图（zi.samplingstants，zi.InputData）；标题（[z.ExperimentName{i}）': '输入名称{1}]）；如果（i''）;其他的包含([z。域' ('时间单位')']);终止轴(“紧”)分段（z.Ne，2，2*i）；%输出。绘图（zi.samplingstants，zi.OutputData）；标题（[z.ExperimentName{i}）': 'zi.OutputName{1}]）；如果（i''）;其他的包含([z。域' ('时间单位')']);终止轴(“紧”）;终止

图2：两体系统的输入输出数据。

初始两体模型的性能

在估计四个未知参数之前，我们使用初始参数值模拟系统。我们使用默认微分方程求解器（ode45）使用默认解算器选项。当仅使用两个输入参数调用时，COMPARE将估计完整的初始状态向量，在这种情况下，每个实验一个，即3个实验，每个实验具有3×1的状态向量，总共意味着9个估计的初始状态。模拟和真实输出显示在绘图窗口中，可以看出拟合为还算不错，但不如预期的好。

clf比较(z, nlgr);

图3：比较初始两体模型的真实输出和模拟输出。

参数估计

为了提高拟合度，接下来对四个参数进行估计。我们在估计阶段使用第一个和最后一个实验的数据，并保留第二个实验的数据，以进行纯粹的验证。

opt=nlgreyestOptions(“显示”,“开”); nlgr=nlgreyest（getexp（z[13]），nlgr，opt）；

估计两体模型的性能

为了研究估计模型的性能，最后对其进行了仿真。通过裁剪nlgr的初始状态结构阵列，可以完全指定在每个实验中估计哪些状态，例如COMPARE。让我们在这里定义并使用一个结构，其中实验1的初始状态x1(0)和x2(0)是估计的，实验2的初始状态x2(0)是估计的，实验3的初始状态x3(0)是估计的。经过此修改后，测量值和模型输出之间的比较将显示在绘图窗口中。

nlgr.InitialStates=struct(“姓名”，getinit（nlgr，“姓名”），...“单位”，getinit（nlgr，“单位”），...“价值”，零（1,3），“最低限度”无穷(1、3)...“最大值”，Inf（1,3），“固定的”,...{(假假真);(真的假的真的);[真真正假]});比较(z, nlgr compareOptions (“初始条件”,“模型”））;

图4：比较估计的两体模型的真实输出和模拟输出。

特别令人感兴趣的是第二个实验的数据结果，这些数据没有用于参数估计。所有实验都清楚地模拟了真实系统的动力学。估计的参数也非常接近用于生成实验数据的参数：

disp(“真实估计参数向量”）;

真估计参数向量

ptrue=[400；2e3；0.0001；1700]；fprintf(“% 10.5 f % 10.5 f \ n”，[ptrue'；getpvec（nlgr'）；

400.00000 402.16624 2000.00000 1987.96042 0.00010 0.00010 1700.00000 1704.51209

最后使用PRESENT命令，我们可以得到关于估计模型的附加信息:

出席（nlgr）；

nlgr=由“两个物体”定义的连续时间非线性灰箱模型（MEX文件）：dx/dt=F（t，u（t），x（t），p1，…，p4）y（t）=H（t，u（t），x（t），p1，…，p4）+e（t），具有1个输入，3个状态，1个输出和4个自由参数（共4个）。输入：u（1）外力（t）[N]状态：初始值x（1）运动物体的位置（t）[m] xinit@exp10（估计）英寸[-Inf，Inf]xinit@exp20（估计）英寸[-Inf，Inf]xinit@exp30（固定）in[-Inf，Inf]x（2）运动物体的速度（t）[m/s]xinit@exp1[-Inf，Inf]中的0（固定）xinit@exp20（估计）英寸[-Inf，Inf]xinit@exp30（固定）in[-Inf，Inf]x（3）物体之间的干摩擦力（t）[N]xinit@exp10（固定）英寸[-Inf，Inf]xinit@exp2[-Inf，Inf]中的0（固定）xinit@exp3[-Inf，Inf]中的0（估计值）输出：y（1）移动体的位置（t）[m]参数：值p1移动体的质量[m]402.166（估计值）in[-Inf，Inf]p2滑动摩擦力系数[kg/s]1987.96（估计值）in[-Inf，Inf]p3距离相关干摩擦参数[m][-Inf，Inf]中的0.000104534（估计值）p4力相关干摩擦参数[N]1704.51（估计值）in[-Inf，Inf]名称：两体系统状态：模型在估计后修改。更多信息请参见模型的“报告”属性。

结论

这个例子描述了在进行IDNLGREY建模时如何使用多个实验数据。可以使用任意数量的实验，对于每一个这样的实验，都可以完全指定要在NLGREYEST、COMPARE、PREDICT等中估计哪个或哪个初始状态。

补充资料

有关使用系统识别工具箱识别动态系统的详细信息™ 参观系统识别工具箱产品信息页面。