将其作为单独发布，因为我让您的代码正常工作。好吧，运行并产生输出：P
实际上，一个大问题是一些我没有注意到的隐形广播，但我一路上改变了很多其他事情。
首先，隐形广播是，如果将一维函数与一个参数集成，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 and 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()
