经典钟摆:一些算法相关问题
这个例子显示了估计算法的选择如何影响非线性灰盒模型估计的结果。我们使用由非线性摆系统产生的数据，如图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(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)});%所有参数> 0。
三种初始摆模型的性能
在尝试估计任何参数之前，我们用猜测的参数值模拟系统。我们对三个可用的求解器这样做，固定步长的欧拉正向(ode1)，自适应步长的龙格-库塔23 (ode23)和自适应步长的龙格-库塔45 (ode45)。使用这三个解算器获得的输出显示在绘图窗口中。
用一阶欧拉正演ODE求解器计算的模型。Nlgref = nlgr;nlgref.SimulationOptions.Solver =“ode1”;%欧拉远期。nlgref.SimulationOptions.FixedStep = z.Ts*0.1;%步长。% B.用自适应龙格-库塔23 ODE求解器计算模型。Nlgrrk23 = nlgr;nlgrrk23.SimulationOptions。解算器=“ode23”;%龙格-库塔用自适应龙格-库塔45 ODE求解器计算模型。Nlgrrk45 = nlgr;nlgrrk45.SimulationOptions。解算器=“数值”;%龙格-库塔比较(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(clock, Tef);
然后，我们继续使用基于Runge-Kutta 23求解器(ode23)的模型。
Trk23 =时钟;Nlgrrk23 = nlgreyest(z, Nlgrrk23, opt);%执行参数估计。Trk23 = etime(clock, Trk23);
我们最后使用龙格-库塔45求解器(ode45)。
Trk45 =时钟;Nlgrrk45 = nlgreyest(z, Nlgrrk45, opt);%执行参数估计。Trk45 = etime(clock, 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'， trk23, nlgrrk23.Parameters(3).Value);流(' ode45: %3.1f %1.3f\n'， trk45, nlgrrk45.Parameters(3).Value);
然而，这并不意味着模拟的输出与实际输出差别很大，因为由不同的微分方程求解器产生的误差通常在估计过程中被考虑在内。通过模拟三种估计的钟摆模型，很容易看出这一点。
比较(z, nlgref, nlgrrk23, nlgrrk45, 1，…compareOptions (“InitialCondition”，“模型”));
图4:三种估算钟摆模型的真实输出与模拟输出的比较。
从图中可以看出，这些模型的最终预测误差(Final Prediction Error, FPE)判据值也非常接近:
Fpe (nlgref, nlgrrk23, nlgrrk45);
基于此，我们得出结论，所有三种模型都能够捕捉钟摆动力学，但基于欧拉的正演模型反映了底层物理相当差。提高其“物理”性能的唯一方法是减少求解器的步长，但这也意味着求解时间显著增加。实验表明，当考虑精度和计算速度时，Runge-Kutta 45解算器是求解非刚性问题的最佳解算器。
结论
龙格-库塔45 (ode45)求解器通常相对较快地返回高质量的结果，因此被选为IDNLGREY中的默认微分方程求解器。输入“help idnlgrey”。SimulationOptions”提供了有关可用求解器的更多详细信息，输入“help nlgreyestOptions”提供了可用于IDNLGREY建模的各种估计选项的详细信息。
额外的信息
有关使用System identification Toolbox™识别动态系统的更多信息，请访问系统识别工具箱产品信息页面。
