昨天我在这里发布了一个问题:ValueError and odepack.error using integrate.odeint(),我认为已经成功回答了.但是,此后我注意到了两件事.
>运行此程序时,它不会趋向于所需的速度vr
>当运行程序并随着时间变化而改变角度(代表道路的坡度)时,它并不总是返回所需的速度vr,甚至不会返回到先前的稳定状态.
该程序旨在为整体控制系统(特别是巡航控制)建模.它当前从速度v0开始,以该速度运行一段时间,然后启用巡航控制.在这一点上,我们应该看到速度的变化(确实如此),该变化最终稳定在所需的速度vr上.没有.由于未知的原因,它达到了其他一些值,并且该值根据骑乘者的坡度而有所不同.不管初始速度如何,它仍然无法达到所需的速度
我一直在尝试不同的参数和变量,但无济于事.我认为问题在于控制器未通过正确的当前速度,但是我不确定如何解决该问题.
如果您需要其他信息,请告诉我.如果我应该编辑上一个问题,请告诉我,我会做的,请提前对不起.
这是我的代码:
import matplotlib.pylab as plt
import numpy as np
import scipy.integrate as integrate
##Parameters
kp = .5 #proportional gain
ki = .1 #integral gain
vr = 25 #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(61, 500, 5000) #time
theta = np.radians(4) #Theta
def torque(v):
return Tm * (1 - B*(an*v/wm - 1)**2)
def vderivs(v, t):
v1 = an * controller(v, t) * 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*(t>=200)
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)
t0l = [i for i in range(61)]
vf=[v0 for i in range(61)]+[v for v in velocity(v0,theta,t)]
tf=t0l+[time for time in t]
plt.plot(tf, vf, 'k-', label=('V(0) = '+str(v0)))
v0=35.
vf=[v0 for i in range(61)]+[v for v in velocity(v0,theta,t)]
plt.plot(tf, vf, 'b-', label=('V(0) = '+str(v0)))
v0=vr
vf=[v0 for i in range(61)]+[v for v in velocity(v0,theta,t)]
plt.plot(tf, vf, 'g-', label=('V(0) = Vr'))
plt.axhline(y=vr, xmin=0, xmax=1000, color='r', label='Desired Velocity')
plt.legend(loc = "upper right")
plt.axis([0,500,18,36])
plt.show()
这是使用一定范围的初始速度按原样运行程序的图
速度的第一个急剧变化是启动巡航控制时的速度,第二个是骑乘者的坡度变化时的速度
解决方法:
你的猜测是正确的.目标系统和PI控制器已集成在一起,您无法将其分为两个节点.我修改了您的代码,该系统具有两个状态变量:一个用于系统的速度,一个用于控制误差的积分:
import matplotlib.pylab as plt
import numpy as np
import scipy.integrate as integrate
##Parameters
kp = .5 #proportional gain
ki = .1 #integral gain
vr = 25 #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(61, 500, 5000) #time
theta = np.radians(4) #Theta
def torque(v):
return Tm * (1 - B*(an*v/wm - 1)**2)
def vderivs(status, t):
v, int_err = status
err = vr - v
control = kp * err + ki * int_err
v1 = an * control * 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*(t>=200)
return vtot/m, err
def velocity(desired, theta, time):
return integrate.odeint(vderivs, [desired, 0], time)[:, 0]
t0l = [i for i in range(61)]
vf=[v0 for i in range(61)]+[v for v in velocity(v0,theta,t)]
tf=t0l+[time for time in t]
plt.plot(tf, vf, 'k-', label=('V(0) = '+str(v0)))
v0=35.
vf=[v0 for i in range(61)]+[v for v in velocity(v0,theta,t)]
plt.plot(tf, vf, 'b-', label=('V(0) = '+str(v0)))
v0=vr
vf=[v0 for i in range(61)]+[v for v in velocity(v0,theta,t)]
plt.plot(tf, vf, 'g-', label=('V(0) = Vr'))
plt.axhline(y=vr, xmin=0, xmax=1000, color='r', label='Desired Velocity')
plt.legend(loc = "upper right")
plt.axis([0,500,18,36])
plt.show()
输出: