使用integrate.odeint()时出现ValueError和odepack.error

Question

我正在尝试编写一个方程来建模和绘制积分控制系统(特别是关于巡航控制)。然而，每当我运行它时，我都会收到两个错误：

ValueError: 对象太深，无法得到所需的数组 odepack.error: 函数调用结果不是浮点数数组。

我已经阅读了以下问题：

- scipy curve_fit error: Result from function call is not a proper array of floats - 如何使用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)

请让我知道您需要我提供的其他任何信息！

- Dan Oberlam

---

你能在这里发布一段相关的代码片段吗？那可能会有所帮助。 - aIKid

对不起，你是什么意思？我正在使用的确切代码在帖子中，并且在这里http://pastebin.com/5pubqMg1。 - Dan Oberlam

我不确定你的结构是否可行，因为每个集成取决于其他集成的方式。此外，您似乎希望 controller 取一个 t 的值，但是 controller 只是 odeint 的包装器，它期望 t 是要集成的所有时间值的数组，但是您在 vderivs 中使用标量来调用它。您可能只需将 vderivs 向量化即可解决此问题，但我不确定。 - askewchan

是的，当我开始出现错误时，我意识到了这一点。积分控制系统应该采取自己的输出，将其与预期/期望的输出进行比较，然后将该误差输入到自身中。我尽力处理这个问题，但正如你所看到的，效果不佳。我还在考虑使用延迟（这更加真实）来实现它可能会解决这个问题，但是引入延迟会带来新的问题...当你说矢量化vderivs时，你指的是什么？ - Dan Oberlam

是的，但一个积分控制系统集成到当前时间，而你的代码要求在所有时间上进行积分。 - askewchan

Answer 1

这不是一个答案，而是我注意到你代码中的一些评论。

正如@qmorgan所指出的那样，您在controller函数中命名了一个参数与另一个函数相同：velocity。尽量保持变量名称的一致性，以避免此类情况发生。例如，您所有的常量都是简短的缩写词。因此，在函数中命名参数时，可以使用单词，这样您就会习惯于知道您在哪里使用了什么。

您在函数中有一些类似的错误，其中您在函数中为某些内容定义了参数，但在函数体中忽略了它，并使用了常量。例如，您已经定义了velocity以接受(mass, desired, theta, t)，但您传递了(m, v, theta, time)，其中v是您的起始速度。要小心！您没有注意到这个错误，因为实际上，velocity 忽略了 desired，而是只使用全局常量vr。相反，整个代码块应该像这样：

def velocity(mass, desired, theta, t):
    return integrate.odeint(vderivs, desired, t)

plt.plot(time, velocity(m, vr, theta, time), 'k-')

---

要将列表转换为numpy数组，您只需使用np.array（vel，float）而不是np.array（[x for x in vel]，float），因为[x for x in vel]与vel本身相同。

---

在vderivs中，您循环遍历t可能完全可以被消除，但至少我认为它应该更像这样：

ti = 0
while ti < t:
    ...
    ti += 1

不会忽略输入的 t。

---

目前，uderivs 接收当前速度和全局定义的期望速度，并返回它们之间的差异：

def uderivs(v, t):  
    return vr - v   

但是你总是将其传递给controller（请参见controller的定义），因此每次集成它时，它只会返回一个长度为t.size且值为vr的平坦数组。

---

你可能想使用 theta = np.radians(4)，而不是 theta = 4

---

在numpy中已经存在函数np.sign，无需自行实现。

- askewchan

---

好的，我已经按照你建议所做的更改除了关于传递 uderivs vr 部分。我不确定如何传递实际当前速度，因为当前速度取决于控制器。我应该尝试实现一些延迟吗？如果是这样，应该放在哪里？我在想，如果我可以在 vderivs 函数中生成某种初始条件，那么这样会更好，只是我不确定如何着手。此外，为什么 vderivs 中的循环是不必要的？目标是创建一个“上坡前”和“上坡后”的情况。 - Dan Oberlam

Answer 2

我将这个作为独立的帖子发布，因为我已经让你的代码运行并产生了输出 :P

实际上一个大问题是我没有注意到的一些隐秘的广播，但是我沿途改变了很多其他的东西。

首先，隐秘广播是如果你用一个参数来集成一个1d函数，odeint返回一个列向量，然后当你对该结果进行一些行向量操作时，你会得到一个2d数组（矩阵）。例如：

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()

- askewchan

---

我不再收到错误信息了，有时它似乎在做某些事情，但它一直被卡在uderivs()的循环中。你知道原因吗？我是直接从你在这里写的内容复制/粘贴的。 - Dan Oberlam

@Dannnno，我不确定是什么原因导致这种情况，这段代码在我的电脑上运行良好。 - askewchan

我找到问题了，我的电脑似乎不喜欢ts。移除它后，一切都正常了，谢谢！ - Dan Oberlam

我从未听说过那个。你知道它不喜欢什么吗？ - askewchan

我不是很确定。我能想到的唯一可能就是在使用早期的time = np.linspace(0,20,50)来形成ts = np.arange(0,t,dt)时制作了一个矩阵或类似的东西。 - Dan Oberlam

Answer 3

我看到的一个错误出现在函数controller中; 您正在尝试从整数(vr)中减去一个函数(velocity)，这可能是一些问题的根源。

"对象对于所需数组太深"的错误可能来自另一个问题，即controller和velocity函数的返回值形状错误；它们需要是1-d numpy数组。您可以使用flatten()方法来解决此问题：

In [82]: z.shape
Out[82]: (50, 1)

In [83]: z_flat=z.flatten()

In [84]: z_flat.shape
Out[84]: (50,)

但是为了完全调试，您需要在控制器函数中修复上述错误。

- qmorgan

---

实际上，velocity是函数controller的参数名称。 这是不良习惯，但并非错误。 - askewchan

就像askewchan所说的那样，那不是使用函数velocity而是一个名为velocity的变量。它已经被重命名为currentV。我还添加了.flatten()但是仍然出现相同的错误。问题已经被编辑以反映这些内容。 - Dan Oberlam