ValueError и odepack.error с использованием integrate.odeint()
Я пытаюсь написать уравнение для моделирования и затем построить интегральную систему управления (особенно в отношении круиз-контроля). Тем не менее, я получаю две ошибки всякий раз, когда я запускаю его:
ValueError: объект слишком глубокий для нужного массива odepack.error: Результат вызова функции не является правильным массивом с плавающей точкой.
Я прочитал эти вопросы:
- scipy curve_fit error: результат от вызова функции не является правильным массивом с плавающей точкой
- Как решить это дифференциальное уравнение используя scipy odeint?
- Слишком глубокий объект для требуемого массива - 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)
Пожалуйста, дайте мне знать, если вам что-то нужно от меня!
3 ответа
Публикация этого как отдельного, потому что я получил ваш код на работу. Ну, чтобы запустить и произвести вывод: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()
Это на самом деле не ответ, а несколько комментариев к вашему коду, которые я заметил.
Как указывает @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
Но ты всегда проходишь это vr (см. определение controller), поэтому каждый раз, когда вы его интегрируете, он просто возвращает плоский массив длины t.size и значение vr,
Вместо theta = 4 ты наверное хочешь theta = np.radians(4)
Там уже в NumPy функция np.signНет необходимости реализовывать свои собственные.
Я вижу одну ошибку в функции controller; вы пытаетесь вычесть функцию (velocity) из целого числа (vr) что, вероятно, является источником некоторых проблем.
Ошибки "объект слишком глубокий для нужного массива", вероятно, происходят из-за отдельной проблемы, когда возвращаются функции controller а также velocity имеют неправильную форму; они должны быть 1-ю массивами. Вы можете исправить это используя flatten() метод:
In [82]: z.shape
Out[82]: (50, 1)
In [83]: z_flat=z.flatten()
In [84]: z_flat.shape
Out[84]: (50,)
Но для полной отладки необходимо исправить вышеуказанную ошибку в функции контроллера.