经典摆:一些与算法相关的问题

打开实时脚本

这个例子显示了估计算法的选择如何影响非线性灰盒模型估计的结果。我们使用非线性摆系统产生的数据，如图1所示。我们特别说明了微分方程求解器的选择如何影响结果。

图1:经典摆的示意图。

输出数据

我们通过加载输出数据(时间序列数据)开始建模之旅。数据包含一个输出，y，这是摆的角度位置[rad]。钟摆向下时，角度为零，逆时针方向增大。数据点有1001个(模拟的)样本，样本时间为0.1秒。钟摆受到恒定重力的影响，但没有其他外生力影响钟摆的运动。为了研究这种情况，我们创建了一个IDDATA对象:

负载(fullfile (matlabroot“工具箱”,“识别”,“iddemos”,“数据”,“pendulumdata”));Z = iddata(y， []， 0.1，“名字”,“摆”）;z.OutputName =“摆的位置”；z.OutputUnit =rad的；z.Tstart = 0;z.TimeUnit =“年代”；

摆的角度位置(输出)显示在一个绘图窗口中。

图(“名字”[z。的名字:输出数据]);情节(z);

图2:摆(输出)的角位置。

摆造型

下一步是为摆系统指定模型结构。它的动力学已经在许多书籍和文章中进行了研究，并得到了很好的理解。如果我们选择x1(t)作为摆的角位置[rad]， x2(t)作为摆的角速度[rad/s]，那么建立如下的状态空间结构就相当简单了:

d / dt x1 (t) = x2 (t) d / dt x2 (t) = - (g / l) * sin (x1 (t))——(b / l ^ 2)) (m * * x2 (t)

Y (t) = x1(t)

具有参数(或常量)

g -重力常数[m/s^2] l -摆杆长度[m] b -粘性摩擦系数[Nms/rad] m -摆锤的质量[kg]

我们将这些信息输入到MATLAB文件pendulum_m中。M，内容如下:

function [dx, y] = pendulum_m(t, x, u, g, l, b, m, varargin) % pendulum_m一个摆系统。

输出方程。Y = x(1);%角位置。

%状态方程。Dx = [x(2);．.．%角位置。- (g / l) * sin (x (1)) - b / (m * l ^ 2) * x(2)……%角速度。];

下一步是创建一个通用的IDNLGREY对象来描述摆。我们还输入有关输入、状态、输出和参数的名称和单位的信息。由于物理现实，所有模型参数必须为正，这是通过将所有参数的“最小”属性设置为MATLAB®中最小可识别的正值，eps(0)来实现的。

文件名=“pendulum_m”；描述模型结构的文件。Order = [1 0 2];%模型订单[ny nu nx]。参数= [9.81;1;0.2;1);初始参数。InitialStates = [1;0);初始状态。Ts = 0;%时间连续系统。nlgr = idnlgrey(文件名，顺序，参数，InitialStates, Ts);nlgr。OutputName =“摆的位置”；nlgr。OutputUnit =rad的；nlgr。TimeUnit =“年代”；NLGR = setinit(NLGR，“名字”, {“摆的位置”“摆速度”});NLGR = setinit(NLGR，“单位”, {rad的“rad / s”});NLGR = setpar(NLGR，“名字”, {“引力常数”“长度”.．.的摩擦系数“质量”});NLGR = setpar(NLGR，“单位”, {“米/秒^ 2”“米”“Nms / rad”“公斤”});NLGR = setpar(NLGR，“最低”， {eps(0) eps(0) eps(0) eps(0)};%所有参数> 0。

三种初始摆模型的性能

在尝试估计任何参数之前，我们用猜测的参数值模拟系统。我们对三个可用的求解器这样做，欧拉前进固定步长(ode1)，龙格-库塔23自适应步长(ode23)，龙格-库塔45自适应步长(ode45)。使用这三个求解器时获得的输出显示在一个绘图窗口中。

A.用一阶欧拉正向ODE求解器计算的模型。Nlgref = nlgr;nlgref.SimulationOptions.Solver =“ode1”；%欧拉正向。fixedstep = z.Ts*0.1;步长。B.自适应龙格-库塔23 ODE求解器计算的模型。Nlgrrk23 = nlgr;nlgrrk23.SimulationOptions。解算器=“ode23”；龙格-库塔23。用自适应龙格-库塔45 ODE求解器计算模型。Nlgrrk45 = nlgr;nlgrrk45.SimulationOptions。解算器=“数值”；龙格-库塔45。比较(z, nlgref, nlgrrk23, nlgrrk45, 1，.．.compareOptions (“InitialCondition”,“模型”));

图3:三种初始摆模型的真实输出与模拟输出的比较。

可以看到，欧拉正向法的结果(网格比默认使用的网格更细)与龙格-库塔求解器得到的结果有显著差异。在这种情况下，欧拉前向求解器产生了一个相当好的结果(就模型拟合而言)，而龙格-库塔求解器获得的输出有点有限。然而，这有点具有欺骗性，后面会很明显地看到，因为粘性摩擦系数b的初值是实际值的两倍。

参数估计

重力常数g，杆子的长度l，和杆子的质量m，可以很容易地测量或从表格中得到，而无需估计。然而，这通常不是摩擦系数的情况下，通常必须估计。因此，我们将前三个参数固定为g = 9.81, l = 1, m = 1，并且只考虑b作为自由参数:

Nlgref = setpar(Nlgref，“固定”，{真真假真});Nlgrrk23 = setpar(Nlgrrk23，“固定”，{真真假真});Nlgrrk45 = setpar(Nlgrrk45，“固定”，{真真假真});

接下来我们用三个微分方程求解器估计b。我们从基于欧拉正演的模型结构开始。

opt = nlgreyestOptions(“显示”,“上”）;Tef =时钟;Nlgref = nlgreyest(z, Nlgref, opt);进行参数估计。Tef = etime(时钟，Tef);

然后我们继续基于龙格-库塔23求解器(ode23)的模型。

Trk23 =时钟;Nlgrrk23 = nlgreyest(z, Nlgrrk23, opt);进行参数估计。Trk23 = etime(时钟，Trk23);

最后我们使用龙格-库塔45求解器(ode45)。

Trk45 =时钟;Nlgrrk45 = nlgreyest(z, Nlgrrk45, opt);进行参数估计。Trk45 = etime(时钟，Trk45);

三种估计摆模型的性能

下面总结了使用这三个求解器时的结果。b的真值是0.1，由ode45获得。Ode23返回一个非常接近0.1的值，而ode1返回一个与0.1相差很大的值。

disp ('估计时间估计b值'）;流(' ode1: %3.1f %1.3f\n'， tef, nlgref.Parameters(3).Value);流(' ode23: %3.1f %1.3f\n'参数(3).Value);流(' ode45: %3.1f %1.3f\n'， trk45, nlgrrk45.Parameters(3).Value);

然而，这并不意味着模拟输出与实际输出差异很大，因为不同的微分方程求解器产生的误差通常在估计过程中被考虑在内。这是一个可以通过模拟三个估计的摆模型很容易看到的事实。

比较(z, nlgref, nlgrrk23, nlgrrk45, 1，.．.compareOptions (“InitialCondition”,“模型”));

图4:三种估算摆模型的真实输出与模拟输出的比较。

如图所示，这些模型的最终预测误差(FPE)准则值也相当接近:

Fpe (nlgref, nlgrrk23, nlgrrk45);

基于此，我们得出结论，所有三个模型都能够捕捉到钟摆动力学，但基于欧拉正演的模型反映底层物理相当差。提高其“物理”性能的唯一方法是减少求解器的步长，但这也意味着求解时间显著增加。我们的实验表明，龙格-库塔45求解器在精度和计算速度方面都是非刚性问题的最佳求解器。

结论

龙格-库塔45 (ode45)求解器通常相对快速地返回高质量的结果，因此被选为IDNLGREY中的默认微分方程求解器。输入“help idnlgrey”。SimulationOptions”提供了有关可用解算器的更多详细信息，键入“help nlgreyestOptions”提供了关于可用于IDNLGREY建模的各种估计选项的详细信息。

额外的信息

有关使用系统识别工具箱™识别动态系统的更多信息，请访问系统识别工具箱产品信息页面。