1

我正在使用带有 ODE45 的 Octave 来模拟 ODE 方程组。但问题是 ODE 模拟给出了错误的值。看看这个 Octave 代码：

function dx = dynamik(t, x)
b1 = 1000;
b2 = 2000;
m1 = 10;
m2 = 7;
M = 2000;
g = 9.82;
mu = 0.3;
L = 0.1;
Ap = 0.004;
Am = 0.002;
Pp = 2*10^6;
Pm = 2.1*10^6;
k1 = 1.78e+4;
k2 = 4.04e+4;
k3 = 8.86e+3;

dx= [ x(2);
    -k1/m1*x(1) + k1/m1*x(3) - b1/m1*x(4) + b1/m1*x(4) + Ap/m1*Pp - Pm*Am/m1*x(2);
    x(4);
    k1/M*x(1) - k1/M*x(3) + b1/M*x(2) - b1/M*x(4) - g*mu*x(4) - k2/M*x(3) + k3*L/M*sin(x(5));
    x(6);
    3*k2/(m2*L)*x(3) - 3*k2/m2*sin(x(5)) - 3*k3/(m2*L^2)*x(5) - 3*b2/(m2*L^2)*x(6) + 3*g/(2*L)*sin(x(5))];
endfunction

tspan = 0:0.5:10;
init = [0, 0, 0, 0, 0, 0];
[t, y] = ode45(@dynamik, tspan, init);

这给出了：

>> y
y =

0.0000e+00    0.0000e+00    0.0000e+00    0.0000e+00    0.0000e+00    0.0000e+00
1.6659e+08   -6.9253e+10   -1.9336e+05    8.0380e+07   -4.4787e+07    3.8388e+12
5.9331e+18   -2.4665e+21   -6.8864e+15    2.8628e+18   -1.3333e+32    1.1428e+37
2.1131e+29   -8.7843e+31   -2.4526e+26    1.0196e+29   -3.9691e+56    3.4019e+61
7.5258e+39   -3.1286e+42   -8.7350e+36    3.6313e+39   -1.1816e+81    1.0127e+86
2.6803e+50   -1.1142e+53   -3.1110e+47    1.2933e+50  -3.5174e+105   3.0148e+110
9.5460e+60   -3.9684e+63   -1.1080e+58    4.6060e+60  -1.0471e+130   8.9747e+134
3.3998e+71   -1.4134e+74   -3.9461e+68    1.6404e+71  -3.1171e+154   2.6717e+159
1.2109e+82   -5.0337e+84   -1.4054e+79    5.8425e+81  -9.2794e+178   7.9533e+183
4.3125e+92   -1.7928e+95   -5.0054e+89    2.0808e+92  -2.7624e+203   2.3676e+208
1.5359e+103  -6.3849e+105  -1.7827e+100   7.4109e+102  -8.2234e+227   7.0482e+232
5.4701e+113  -2.2740e+116  -6.3491e+110   2.6394e+113  -2.4480e+252   2.0982e+257
1.9482e+124  -8.0989e+126  -2.2612e+121   9.4003e+123  -7.2875e+276   6.2461e+281
6.9386e+134  -2.8844e+137  -8.0534e+131   3.3479e+134  -2.1694e+301   1.8594e+306
       NaN           NaN           NaN           NaN           NaN           NaN
       NaN           NaN           NaN           NaN           NaN           NaN
       NaN           NaN           NaN           NaN           NaN           NaN
       NaN           NaN           NaN           NaN           NaN           NaN
       NaN           NaN           NaN           NaN           NaN           NaN
       NaN           NaN           NaN           NaN           NaN           NaN
       NaN           NaN           NaN           NaN           NaN           NaN

但在 OpenModelica 中，我有以下代码：

model Hydraulik_System_1
// Types for variables
type PositionCylinder = Real(unit="m");
type Position = Real(unit="m");
type Velocity = Real(unit="m/s");
type DegreesPosition = Real(unit="rad");
type DegreesVelocity = Real(unit="rad/s");
type pressure = Real(unit="Pa");
type flows = Real(unit="l/min", min=0.0);

// Types for parameters
type spring = Real(unit="N/m");
type damper = Real(unit="Ns/m");
type mass = Real(unit="kg");
type friction = Real(unit="%");
type length = Real(unit="m");
type gravitation = Real(unit="m/s^2");
type area = Real(unit="cm^2");

// Parameters
parameter spring k1 = 1.78*10^4;
parameter spring k2 = 4.04*10^4; 
parameter spring k3 = 8.86*10^3;
parameter mass m1 = 10;
parameter mass m2 = 7;
parameter mass M = 2000;
parameter damper b1 = 1000;
parameter damper b2 = 2000; 
parameter gravitation g = 9.82;
parameter friction mu = 0.3;
parameter area Am = 0.002;
parameter area Ap = 0.004;
parameter length L = 0.1;
parameter pressure Pp = 2*10^6;
parameter pressure Pm = 2.1*10^6;

// Variables
PositionCylinder x1;
Position x3;
Velocity x2 , x4;
DegreesPosition x5;
DegreesVelocity x6;
initial equation
x1 = 0;
x2 = 0;
x3 = 0;
x4 = 0;
x5 = 0;
x6 = 0;
equation

der(x1) = x2;
der(x2) = - k1/m1*x1 + k1/m1*x3 - b1/m1*x4 + b1/m1*x4 + Ap/m1*Pp - Pm*Am/m1*x2;
der(x3) = x4;
der(x4) = k1/M*x1 - k1/M*x3 + b1/M*x2 - b1/M*x4 - g*mu*x4 - k2/M*x3 + k3*L/M*sin(x5);
der(x5) = x6;
der(x6) = 3*k2/(m2*L)*x3 - 3*k2/m2*sin(x5) - 3*k3/(m2*L^2)*x5 - 3*b2/(m2*L^2)*x6 + 3*g/(2*L)*sin(x5);

end Hydraulik_System_1;

结果如下所示：

在此处输入图像描述

你能告诉我为什么我的模拟会发生这种情况吗？OpenModelica 仿真结果和 Octave 仿真结果之间存在巨大差异。为什么？我是否更改 ODE 求解器并不重要。结果将几乎相同。

4

1 回答 1

1

我使用了lsode，也得到了正确的答案，但是调用参数必须切换，并且tspan更小。

首先，交换函数中的参数：

function dx = dynamik(x, t)

设置时间跨度：

tspan = 0:0.0625:2;

然后lsode调用：

[y,t] = lsode(@dynamik, init, tspan );

更新：已安装odepkg，并且能够重现您的错误。也看到了错误ode23，但没有看到ode23s。这表明您的 ODE 是“僵硬的”，因此 Runge-Kutta 4/5 并不是真正合适的算法。

于 2017-06-09T22:36:28.440 回答