ValueError and odepack.error using integrate.odeint()

Question

I'm trying to write an equation to model and then plot an integral control system (specifically regarding cruise control). However I'm receiving two errors whenever I run it:

ValueError: object too deep for desired array odepack.error: Result from function call is not a proper array of floats.

I've read these questions:

• 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

Which seem like they should be helpful, however I'm unsure how to apply those to my problem. I'm fairly new to python so please bear with me if I've missed something obvious or done something exceptionally silly. I have no problems with plotting it, so once I figure out how to actually get this working I think I'm set.

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)

Please let me know if there is anything else you need from me!

Answer 1

Posting this as separate, because I got your code to work. Well, to run and produce output :P

Actually one big problem is some stealth broadcasting that I didn't notice, but I changed a lot of other things along the way.

First the stealth broadcasting is that if you integrate a 1d function with one parameter, odeint returns a column vector, and then when you do stuff with that result that is a row vector, then you get a 2d array (matrix). For example:

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

You were getting output for velocity that was a column vector like b and adding it to some other function of time, and getting a matrix.

---

With regards to the recursion, I think I sovled that issue. The two derivative functions should take a single scalar t at which point they calculate the derivative. To do that, vderivs needs to do the integral, which it should do over all time up to t. So I redefined them as such:

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

And of course uderivs is fine as is but can be written more simply:

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

Then, make sure that velocity and controller pass the right values (using v0 instead of v for the starting velocity):

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)

Who knows if the physics is correct, but this gives:

Note that it hasn't reached the desired velocity, so I increased the time over which it was to be solved

time = np.linspace(0,50,50) #time

Here is all the code that I ran:

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

Answer 2

One error that I see is in the function controller; you are trying to subtract a function (velocity) from an integer (vr) which is likely the source of some of the problem.

The "object too deep for desired array" errors are probably coming from a separate issue where the returns from the functions controller and velocity are of the wrong shape; they need to be 1-d numpy arrays. You can fix this using the flatten() method:

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

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

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

But in order to fully debug, you need to fix the above error in the controller function.

Answer 3

This isn't really an answer, but a few comments on your code that I've noticed.

As @qmorgan points out, you've named a parameter in your controller function the same as another function: velocity Try to be consistent in your variable names to avoid such things. For example, all of your constants are short abbreviations. So, when you name parameters in functions, perhaps use words, this way you'll have a habit of knowing where you've used what.

You've made a few similar mistakes where you have a parameter for something in your function but you ignore it in the body of the function, and instead use a constant. For example, you've defined velocity to take (mass, desired, theta, t) but you pass it (m, v, theta, time) where v is your starting velocity. Be careful! You didn't notice this mistake because in fact, velocity ignores desired and instead just uses vr the global constant. Instead, this whole bit should have something like

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

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

---

To convert a list to a numpy array, you can just use np.array(vel, float) instead of np.array([x for x in vel], float) since [x for x in vel] is identical to vel itself.

---

In vderivs your loop through t can probably be completely eliminated, but at the very least I think it should be something more like:

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

Which doesn't ignore the input t.

---

As it stands now, uderivs takes a current velocity and a globally defined desired velocity and returns their difference:

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

But you always pass it vr (see definition of controller), so every time you integrate it, it will simply return a flat array of length t.size and value vr.

---

Instead of theta = 4 you probably want theta = np.radians(4)

---

There exists already in numpy the function np.sign, no need to implement your own.