ValueError and odepack.error using integrate.odeint()

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:

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!

Solution

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:

short time

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

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