3

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

ValueError: object too deep for required array odepack.error: 函数调用的结果不是正确的浮点数组。

我读过这些问题：

看起来它们应该会有所帮助，但是我不确定如何将它们应用于我的问题。我对 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)

请让我知道您是否还有其他需要我的东西！

4

3 回答 3

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

这不会忽略 input t。

就目前而言，uderivs采用当前速度和全局定义的期望速度并返回它们的差：

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

但是你总是传递它vr（见 的定义controller），所以每次你集成它时，它都会简单地返回一个长度t.size和值的平面数组vr。

而不是theta = 4你可能想要theta = np.radians(4)

numpy 中已经存在该功能np.sign，无需自己实现。

于 2013-12-07T00:37:53.047 回答
1

将此作为单独发布，因为我让您的代码正常工作。好吧，运行并产生输出：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()
于 2013-12-07T01:14:17.433 回答
-1

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

“对象对于所需数组来说太深”错误可能来自一个单独的问题，即函数返回的controller形状velocity错误；它们需要是一维 numpy 数组。您可以使用以下flatten()方法解决此问题：

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

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

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

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

于 2013-12-07T00:06:39.010 回答