Simulink 仿真的行为与使用 rk4 的相同仿真有很大不同

Question

我在 simulink 中模拟由以下状态 space 方程描述的对象时遇到困难：

状态 space 等式的右侧由下面的函数描述。

function dxdt = RHS( t, x, F)
% parameters
b = 1.22;          % cart friction coeffitient 
c = 0.0027;           %pendulum friction coeffitient
g = 9.81;           % gravity
M = 0.548+0.022*2; % cart weight
m = 0.031*2; %pendulum masses
I = 0.046;%0.02*0.025/12+0.02*0.12^2+0.011*0.42^2; % moment of inertia
l = 0.1313;
% x(1) = theta
% x(2) = theta_dot
% x(3) = x
% x(4) = x_dot


dxdt = [x(2);
    (-(M+m)*c*x(2)-(M+m)*g*l*sin(x(1))-m^2*l^2*x(2)^2*sin(x(1))*cos(x(1))+m*l*b*x(4)*cos(x(1))-m*l*cos(x(1))*F)/(I*(m+M)+m*M*l^2+m^2*l^2*sin(x(1))^2);
    x(4);
    (F - b*x(4) + l*m*x(2)^2*sin(x(1)) + (l*m*cos(x(1))*(c*x(2)*(M + m) + g*l*sin(x(1))*(M + m) + F*l*m*cos(x(1)) + l^2*m^2*x(2)^2*cos(x(1))*sin(x(1)) - b*l*m*x(4)*cos(x(1))))/(I*(M + m) + l^2*m^2*sin(x(1))^2 + M*l^2*m))/(M + m)];
end

下面显示了具有简单可视化的对应 rk4 函数。

function [wi, ti] = rk4 ( RHS, t0, x0, tf, N )

%RK4 approximate the solution of the initial value problem
%
% x'(t) = RHS( t, x ), x(t0) = x0
%
% using the classical fourth-order Runge-Kutta method - this 
% routine will work for a system of first-order equations as 
% well as for a single equation
%
% calling sequences:
% [wi, ti] = rk4 ( RHS, t0, x0, tf, N )
% rk4 ( RHS, t0, x0, tf, N )
%
% inputs:
% RHS string containing name of m-file defining the 
% right-hand side of the differential equation; the
% m-file must take two inputs - first, the value of
% the independent variable; second, the value of the
% dependent variable
% t0 initial value of the independent variable
% x0 initial value of the dependent variable(s)
% if solving a system of equations, this should be a 
% row vector containing all initial values
% tf final value of the independent variable
% N number of uniformly sized time steps to be taken to
% advance the solution from t = t0 to t = tf
%
% output:
% wi vector / matrix containing values of the approximate 
% solution to the differential equation
% ti vector containing the values of the independent 
% variable at which an approximate solution has been
% obtained
%

% x(1) = theta
% x(2) = theta_dot
% x(3) = x
% x(4) = x_dot
t0 = 0; tf = 5; x0 = [pi/2; 0; 0; 0]; N = 400;
neqn = length ( x0 );
ti = linspace ( t0, tf, N+1 );
wi = [ zeros( neqn, N+1 ) ];
wi(1:neqn, 1) = x0';

h = ( tf - t0 ) / N;
% force
u = 0.0;

%init visualisation
h_cart = plot(NaN, NaN, 'Marker', 'square', 'color', 'red', 'LineWidth', 6);
hold on
h_pend = plot(NaN, NaN, 'bo', 'LineWidth', 3);
axis([-5 5 -5 5]);
axis manual;
xlim([-5 5]);
ylim([-5 5]);

for i = 1:N
    k1 = h * feval ( 'RHS', t0, x0, u );
    k2 = h * feval ( 'RHS', t0 + (h/2), x0 + (k1/2), u);
    k3 = h * feval ( 'RHS', t0 + h/2, x0 + k2/2, u);
    k4 = h * feval ( 'RHS', t0 + h, x0 + k3, u);
    x0 = x0 + ( k1 + 2*k2 + 2*k3 + k4 ) / 6;
    t0 = t0 + h;
    % model output
    wi(1:neqn,i+1) = x0';

    % model visualisation
    %plotting cart
    l = 2;
    set(h_cart, 'XData', x0(3), 'YData', 0, 'LineWidth', 5);
    %plotting pendulum
    %hold on;
   set(h_pend, 'XData', sin(x0(1))*l+x0(3), 'YData', -cos(x0(1))*l, 'LineWidth', 2);
   %hold off;
    % regulator
    pause(0.02);
end;

figure;
plot(ti, wi);
legend('theta', 'theta`', 'x', 'x`');

这为手推车上的钟摆提供了逼真的效果。

现在进入正题。我想在 simulink 中重新创建完全相同的方程式。我认为它会像创建以下 simulink 模型一样简单。我用 RHS 文件中的第二个和第四个方程填充 fcn 块。像这样。

(-(M+m)*c*u(2)-(M+m)*g*l*sin(u(1))-m^2*l^2*u(2)^2*sin(u(1))*cos(u(1))+m*l*b*u(3)*cos(u(1))-m*l*cos(u(1))*u(4))/(I*(m+M)+m*M*l^2+m^2*l^2*sin(u(1))^2)

(u(5) - b*u(4) + l*m*u(2)^2*sin(u(1)) + (l*m*cos(u(1))*(c*u(2)*(M + m) + g*l*sin(u(1))*(M + m) + u(5)*l*m*cos(u(1)) + l^2*m^2*u(2)^2*cos(u(1))*sin(u(1)) - b*l*m*u(4)*cos(u(1))))/(I*(M + m) + l^2*m^2*sin(u(1))^2 + M*l^2*m))/(M + m)

问题是这并没有给出上面的正确结果，而是下面的

有人知道我做错了什么吗？

编辑：@am304 评论后我决定添加以下信息。我将 simulink 求解器的设置更改为使用固定步长 rk4 求解器，以获得相同的结果。上述模型中的第二个 integrator3 已初始化为 pi/2.

Edit2：如果有人想自己查看 simulink 模型，请单击 link 下载文件。

Edit3：正如您在下面的答案中看到的那样，这个问题很简单。您可以下载正确的模型here

Answer 1

我查看了您的 Simulink 模型，看来您可能混淆了您使用的两个函数。您在要放置 x_dd 的地方使用了 theta_dd 函数，反之亦然。

在您的模型中，您还强制将 x_d 设置为常数值 0。我假设您实际上打算将初始条件设置为 0，您可以看到这是通过 Integrator 模块完成的。 x_d（作为 f 的输入）应该是你的状态向量，它也是你的积分器的输出。这只是您定义 x_d 的结果，即 x_dd 的积分。这也是 RK4 的工作方式；您首先使用初始状态向量，然后使用预测状态向量来驱动下一个 RK4 步骤。

范围的结果输出（我在这里输出了你的整个状态向量）如下所示，看起来像你所期望的：

我认为我不应该 link 外部到 simulink 文件，所以如果您想要该文件的副本，您可以打开聊天并索要它。否则上图应该足以帮助您重现相同的结果。

Simulink 仿真的行为与使用 rk4 的相同仿真有很大不同

Simulink simulation behaving massively different than the same simulation using rk4

matlab

simulink

runge-kutta