我正在尝试编写一个方程式来建模，然后绘制一个整体控制系统(特别是关于巡航控制)。但是，每当我运行它时，我都会收到两个错误:
ValueError:对象对于所需数组来说太深 odepack.error:函数调用的结果不是正确的 float 组。
我读过这些问题:
- > scipy curve_fit error: Result from function call is not a proper array of floats
- > How to solve this differential equation using scipy odeint?
- > Object Too Deep for Desired Array - scipy.integrate.odeint
它们似乎应该有所帮助，但我不确定如何将它们应用到我的问题中。我是 python 的新手，所以如果我错过了一些明显的事情或做了一些非常愚蠢的事情，请多多包涵。我对绘制它没有任何问题，所以一旦我弄清楚如何真正让它发挥作用，我想我已经准备好了。
import numpy as np
import scipy.integrate as integrate
##Parameters
kp=.5 #proportional gain
ki=.1 #integral gain
vr=30 #desired velocity in m/s
Tm=190 #Max Torque in Nm
wm=420 #engine speed
B=0.4 #Beta
an=12 #at gear 4
p=1.3 #air density
Cd=0.32 #Drag coefficient
Cr=.01 #Coefficient of rolling friction
A=2.4 #frontal area
##Variables
m=18000 #weight
v=20 #starting velocity
time=np.linspace(0,10,50) #time
theta=np.radians(4) #Theta
def vderivs(state,t):
v = state
vel=[]
ti=0
while ti < t:
v1 = an*controller(ti,vr,v)*torque(v)
v2 = m*Cr*np.sign(v)
v3 = 0.5*p*Cd*A*v**2
v4 = m*np.sin(theta)
if t < 10:
vtot = v1+v2+v3
vfin = np.divide(vtot,m)
else:
vtot = v1+v2+v3+v4
vfin = np.divide(vtot,m)
vel.append(vfin)
ti+=1
trueVel = np.array(vel, float)
return trueVel
def uderivs(state,t):
v = state
deltax = vr - v
return deltax
def controller(time,desired,currentV):
z = integrate.odeint(uderivs, currentV, time)
u = kp*(vr-currentV)+ki*z
return u.flatten()
def torque(v):
return Tm*(1-B*(np.divide(an*v,wm)-1)**2)
def velocity(mass,desired,theta,t):
v = integrate.odeint(vderivs, desired, t)
return v.flatten()
test = velocity(m,vr,theta,time)
print(test)
如果您还需要我，请告诉我!
最佳答案
将其单独发布，因为我已让您的代码正常工作。好吧，运行并产生输出:P
实际上一个大问题是一些我没有注意到的 secret 广播，但我一路上改变了很多其他东西。
首先，隐形广播是，如果您将一个一维函数与一个参数集成，odeint 返回一个列向量，然后当您对该结果进行处理时，该结果是一个行向量，然后您得到二维数组(矩阵)。例如:
In [704]: a
Out[704]: array([0, 1, 2, 3, 4])
In [705]: b
Out[705]:
array([[0],
[1],
[2]])
In [706]: a+b
Out[706]:
array([[0, 1, 2, 3, 4],
[1, 2, 3, 4, 5],
[2, 3, 4, 5, 6]])
您获得的速度输出是一个像 b 这样的列向量，并将其添加到其他时间函数，并获得一个矩阵。
关于递归，我想我解决了那个问题。两个导数函数应该采用单个标量 t，在该点计算导数。为此，vderivs 需要进行积分，它应该一直进行直到 t。所以我重新定义了它们:
dt = .1 # another global constant parameter
def vderivs(v, t):
ts = np.arange(0, t, dt)
v1 = an * controller(v, ts) * torque(v)
v2 = m*Cr*np.sign(v)
v3 = 0.5*p*Cd*A*v**2
v4 = m*np.sin(theta)
vtot = v1+v2+v3+v4*(ts>=10) # a vector of times includes incline only after ts = 10
return vtot/m
当然 uderivs 也可以，但可以写得更简单:
def uderivs(v, t):
return vr - v
然后，确保 velocity 和 controller 传递正确的值(使用 v0 而不是 v起始速度):
def controller(currentV, time):
z = integrate.odeint(uderivs, currentV, time)
return kp*(vr-currentV) + ki*z.squeeze()
def velocity(desired, theta, time):
return integrate.odeint(vderivs, desired, time)
谁知道物理学是否正确，但这给出了:
请注意，它没有达到所需的速度，所以我增加了解决它的时间
time = np.linspace(0,50,50) #time
这是我运行的所有代码:
import matplotlib.pylab as plt
import numpy as np
import scipy.integrate as integrate
##Parameters
kp = .5 #proportional gain
ki = .1 #integral gain
vr = 30 #desired velocity in m/s
Tm = 190 #Max Torque in Nm
wm = 420 #engine speed
B = 0.4 #Beta
an = 12 #at gear 4
p = 1.3 #air density
Cd = 0.32 #Drag coefficient
Cr = .01 #Coefficient of rolling friction
A = 2.4 #frontal area
##Variables
m = 18000.0 #weight
v0 = 20. #starting velocity
t = np.linspace(0, 20, 50) #time
dt = .1
theta = np.radians(4) #Theta
def torque(v):
return Tm * (1 - B*(an*v/wm - 1)**2)
def vderivs(v, t):
ts = np.arange(0, t, dt)
v1 = an * controller(v, ts) * torque(v)
v2 = m*Cr*np.sign(v)
v3 = 0.5*p*Cd*A*v**2
v4 = m*np.sin(theta)
vtot = v1+v2+v3+v4*(ts>=10)
return vtot/m
def uderivs(v, t):
return vr - v
def controller(currentV, time):
z = integrate.odeint(uderivs, currentV, time)
return kp*(vr-currentV) + ki*z.squeeze()
def velocity(desired, theta, time):
return integrate.odeint(vderivs, desired, time)
plt.plot(t, velocity(v0, theta, t), 'k-', lw=2, label='velocity')
plt.plot(t, controller(v0, t), 'r', lw=2, label='controller')
plt.legend()
plt.show()
关于python - ValueError 和 odepack.error 使用 integrate.odeint()，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/20435432/