来自 scipy.integrate 的 Odeint 函数给出了错误的结果

Odeint function from scipy.integrate gives wrong result

我使用 odeint 函数求解耦合微分方程组,并在系统求解后绘制其中一个变量 (theta_i)。我的变量 (theta_i) 来自等式:

theta_i = np.arctan2(g1,g2)

其中 g1 和 g2 是在同一函数中计算的变量。结果必须在 -pi 和 pi 之间,它们应该看起来像这样(来自 matlab 模拟的图):

但是,当我在 odeint 完成后尝试绘制 theta_i 时,我得到了这个(从我的 python 代码绘制):

这真的很奇怪。当我在计算后立即打印 theta_i 的值(仍在函数内部)时,它们看起来是正确的(在 -0.2 和 0.5 之间),因此它必须与结果的存储和我的 odeint 实现有关。来自 odeint 解决方案的所有其他变量都是正确的。我搜索了类似的帖子,但没有人对我有同样的问题。这可能是什么问题?我是 python 的新手,我使用 python 2.7.12。先感谢您。

import numpy as np 
from scipy.integrate import odeint 
import matplotlib.pyplot as plt

added_mass_x = 0.03 # kg 
added_mass_y = 0.04 
mb = 0.3 # kg 
m1 = mb-added_mass_x 
m2 = mb-added_mass_y 
l1 = 0.07 # m 
l2 = 0.05 # m 
J = 0.00050797 # kgm^2 
Sa = 0.0110 # m^2 
Cd = 2.44 
Cl = 3.41 
Kd = 0.000655 # kgm^2 
r = 1000 # kg/m^3

c1 = 0.5*r*Sa*Cd 
c2 = 0.5*r*Sa*Cl 
c3 = 0.5*mb*(l1**2) 
c4 = Kd/J 
c5 = (1/(2*J))*(l1**2)*mb*l2 
c6 = (1/(3*J))*(l1**3)*mb

theta_0 = 10*(np.pi/180) # rad 
theta_A = 20*(np.pi/180) # rad 
f = 2 # Hz

t = np.linspace(0,100,8000) # s

def direct(u,t):
    vcx = u[0]
    vcy = u[1]
    wz = u[2]
    psi = u[3]
    x = u[4]
    y = u[5]
    vcx_i = u[6]
    vcy_i = u[7]
    psi_i = u[8]
    wz_i = u[9]
    theta_i = u[10]
    theta_deg_i = u[11]

    # Subsystem 1
    omega = 2*np.pi*f # rad/s
    theta = theta_0 + theta_A*np.sin(omega*t) # rad
    theta_deg = (theta*180)/np.pi # deg
    thetadotdot = -(omega**2)*theta_A*np.sin(omega*t) # rad/s^2

    # Subsystem 2
    vcxdot = (m2/m1)*vcy*wz-(c1/m1)*vcx*np.sqrt((vcx**2)+(vcy**2))+(c2/m1)*vcy*np.sqrt((vcx**2)+(vcy**2))*np.arctan2(vcy,vcx)-(c3/m1)*thetadotdot*np.sin(theta)
    vcydot = -(m1/m2)*vcx*wz-(c1/m2)*vcy*np.sqrt((vcx**2)+(vcy**2))-(c2/m2)*vcx*np.sqrt((vcx**2)+(vcy**2))*np.arctan2(vcy,vcx)+(c3/m2)*thetadotdot*np.cos(theta)
    wzdot = ((m1-m2)/J)*vcx*vcy-c4*wz*wz*np.sign(wz)-c5*thetadotdot*np.cos(theta)-c6*thetadotdot
    psidot = wz

    # Subsystem 3
    xdotdot = vcxdot*np.cos(psi)-vcx*np.sin(psi)*wz+vcydot*np.sin(psi)+vcy*np.cos(psi)*wz # m/s^2
    ydotdot = -vcxdot*np.sin(psi)-vcx*np.cos(psi)*wz+vcydot*np.cos(psi)-vcy*np.sin(psi)*wz # m/s^2
    xdot = vcx*np.cos(psi)+vcy*np.sin(psi) # m/s
    ydot = -vcx*np.sin(psi)+vcy*np.cos(psi) # m/s

    # Subsystem 4
    vcx_i = xdot*np.cos(psi_i)-ydot*np.sin(psi_i)
    vcy_i = ydot*np.cos(psi_i)+xdot*np.sin(psi_i)
    psidot_i = wz_i
    vcxdot_i = xdotdot*np.cos(psi_i)-xdot*np.sin(psi_i)*psidot_i-ydotdot*np.sin(psi_i)-ydot*np.cos(psi_i)*psidot_i
    vcydot_i = ydotdot*np.cos(psi_i)-ydot*np.sin(psi_i)*psidot_i+xdotdot*np.sin(psi_i)+xdot*np.cos(psi_i)*psidot_i

    g1 = -(m1/c3)*vcxdot_i+(m2/c3)*vcy_i*wz_i-(c1/c3)*vcx_i*np.sqrt((vcx_i**2)+(vcy_i**2))+(c2/c3)*vcy_i*np.sqrt((vcx_i**2)+(vcy_i**2))*np.arctan2(vcy_i,vcx_i)
    g2 = (m2/c3)*vcydot_i+(m1/c3)*vcx_i*wz_i+(c1/c3)*vcy_i*np.sqrt((vcx_i**2)+(vcy_i**2))+(c2/c3)*vcx_i*np.sqrt((vcx_i**2)+(vcy_i**2))*np.arctan2(vcy_i,vcx_i)

    A = 12*np.sin(2*np.pi*f*t+np.pi) # eksiswsi tail_frequency apo simulink
    if A>=0.1:
        wzdot_i = ((m1-m2)/J)*vcx_i*vcy_i-c4*wz_i**2*np.sign(wz_i)-c5*g2-c6*np.sqrt((g1**2)+(g2**2))
    elif A<-0.1:
        wzdot_i = ((m1-m2)/J)*vcx_i*vcy_i-c4*wz_i**2*np.sign(wz_i)-c5*g2+c6*np.sqrt((g1**2)+(g2**2))
    else:
        wzdot_i = ((m1-m2)/J)*vcx_i*vcy_i-c4*wz_i**2*np.sign(wz)-c5*g2


    if g2>0:
        theta_i = np.arctan2(g1,g2)
    elif g2<0 and g1>=0:
        theta_i = np.arctan2(g1,g2)-np.pi
    elif g2<0 and g1<0:
        theta_i = np.arctan2(g1,g2)+np.pi
    elif g2==0 and g1>0:
        theta_i = -np.pi/2
    elif g2==0 and g1<0:
        theta_i = np.pi/2
    elif g1==0 and g2==0:
        theta_i = 0
    theta_deg_i = (theta_i*180)/np.pi
    #print theta_deg_i


    return [vcxdot, vcydot, wzdot, psidot, xdot, ydot, vcxdot_i, vcydot_i, psidot_i, wzdot_i, theta_i, theta_deg_i]

# arxikes synthikes 
vcx_0 = 0.1257 
vcy_0 = 0 
wz_0 = 0 
psi_0 = 0 
x_0 = 0 
y_0 = 0 
vcx_i_0 = 0.1257 
vcy_i_0 = 0 
psi_i_0 = 0 
wz_i_0 = 0 
theta_i_0 = 0.1745 
theta_deg_i_0 = 9.866 
u0 = [vcx_0, vcy_0, wz_0, psi_0, x_0, y_0, vcx_i_0, vcy_i_0, psi_i_0, wz_i_0, theta_i_0, theta_deg_i_0]

u = odeint(direct, u0, t, tfirst=False)

vcx = u[:,0] 
vcy = u[:,1] 
wz = u[:,2] 
psi = u[:,3] 
x = u[:,4] 
y = u[:,5] 
vcx_i = u[:,6] 
vcy_i = u[:,7] 
psi_i = u[:,8] 
wz_i = u[:,9] 
theta_i = u[:,10] 
theta_deg_i = u[:,11] 
print theta_i

plt.figure(17) 
plt.plot(t,theta_i,'r-',linewidth=1,label='theta_i') 
plt.xlabel('t [s]') 
plt.title('theta_i [rad] (Main body CF)') 
plt.legend() 
plt.show()

您所说的问题是 theta_i 不是渐变的一部分。当你制定你的 direct 时,它应该是这样的形式:

def direct(vector, t):
    return vector_dot

最快和最脏的解决方案(不清理代码)是使用您已经定义的函数:

theta_i = [direct(u_i, t_i)[10] for t_i, u_i in zip(t, u)]

我使用了一个较短的间隔:t = np.linspace(0,10,8000)。它产生了这个:

编辑:如何从积分器中删除 theta:

def direct(u, t):
    # your original function as it is

def direct2(u,t):
    return direct(u,t)[:9]

#now integrate the second function
u = odeint(direct2, u0, t)