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

这个例子展示了估计算法的选择如何影响非线性灰盒模型估计的结果。我们使用由非线性摆系统产生的数据，如图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);初始初始状态。t = 0;%的时间连续系统。nlgr = idnlgrey(FileName, Order, Parameters, 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) eps(0)});%所有参数> 0。

三种初始摆模型的性能

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

A.用一阶欧拉正向ODE求解器计算的模型。nlgref = nlgr;nlgref.SimulationOptions.Solver =“ode1”；%欧拉向前。nlgref.SimulationOptions.FixedStep = z.Ts * 0.1;%步长。B.用自适应龙格-库塔23 ODE求解器计算的模型。nlgrrk23 = nlgr;nlgrrk23.SimulationOptions。解算器=“ode23”；%龙格-库塔23。C.用自适应龙格-库塔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。我们从基于欧拉正向的模型结构开始。

选择= nlgreyestOptions (“显示”，“上”）;微软目前=时钟;Nlgref = nlgreyest(z, Nlgref, opt);执行参数估计。Tef = etime(clock, 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'微软,nlgref.Parameters (3) value);流(' ode23: %3.1f %1.3f\n'、trk23 nlgrrk23.Parameters (3) value);流(' ode45: %3.1f %1.3f\n'、trk45 nlgrrk45.Parameters (3) value);

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

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

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

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

消防工程(nlgref nlgrrk23 nlgrrk45);

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

结论

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

额外的信息

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