Model Predictive Control Exercise
Objective: Implement a model predictive controller that automatically regulates vehicle velocity. Implement the controller in Python and tune the controller for acceptable performance. Discuss factors that may be important for evaluating controller performance.
The dynamic relationship between a vehicle gas pedal position (MV) and velocity (CV) is given by the following set of conditions and a single dynamic equation.
Constants m = 500 ! Mass (kg) Parameters b = 50 ! Resistive coefficient (N-s/m) K = 0.8 ! Gain (m/s-%pedal) p = 0 >= 0 <= 100 ! Gas pedal position (%) Variables v = 0 ! initial condition Equations m * $v = -v * b + K * b * p
Implement a model predictive controller that adjusts gas pedal position to regulate velocity. Start at an initial vehicle velocity of 0 m/s and accelerate to a velocity of 40 m/s.
Discuss the controller performance and how it could be tuned to meet multiple objectives including:
- minimize travel time
- remain within speed limits
- improve vehicle fuel efficiency
- discourage excessive gas pedal adjustments
- do not accelerate excessively
There is no need to implement these advanced objectives in simulation for this second part of the exercise, only discuss the possible competing objectives.
Solution
The Model Predictive Control both solves the differential equations that describe the velocity of a vehicle as well as minimizes the control objective function.
The multiple objectives can be implemented with variable constraints or alternative objective functions. For example, if the only objective is to minimize travel time, the solution would be to use full gas pedal position (100%) for the duration of the simulation. If the objective is to minimize travel time but stay within the speed limit, the solution would be to reach 40 m/s as fast as possible without overshoot. Each additional objective has the potential to adapt the solution to achieve an optimal tradeoff. It is a tradeoff because some of the desirable outcomes have conflicting objectives.
See the Dynamic Optimization Course for Excel, MATLAB, and Simulink source files to this problem. Python source in Python Gekko (newer) and APM Python is shown below.
import numpy as np
from random import random
from gekko import GEKKO
import matplotlib.pyplot as plt
#%% Build model
#initialize GEKKO model
m = GEKKO()
#time
m.time = np.linspace(0,20,41)
#constants
mass = 500
#Parameters
b = m.Param(value=50)
K = m.Param(value=0.8)
#Manipulated variable
p = m.MV(value=0, lb=0, ub=100)
#Controlled Variable
v = m.CV(value=0)
#Equations
m.Equation(mass*v.dt() == -v*b + K*b*p)
#%% Tuning
#global
m.options.IMODE = 6 #control
#MV tuning
p.STATUS = 1 #allow optimizer to change
p.DCOST = 0.1 #smooth out gas pedal movement
p.DMAX = 20 #slow down change of gas pedal
#CV tuning
#setpoint
v.STATUS = 1 #add the SP to the objective
m.options.CV_TYPE = 2 #L2 norm
v.SP = 40 #set point
v.TR_INIT = 1 #setpoint trajectory
v.TAU = 5 #time constant of setpoint trajectory
#%% Solve
m.solve()
#%% Plot solution
plt.figure()
plt.subplot(2,1,1)
plt.plot(m.time,p.value,'b-',LineWidth=2)
plt.ylabel('gas')
plt.subplot(2,1,2)
plt.plot(m.time,v.value,'r--',LineWidth=2)
plt.ylabel('velocity')
plt.xlabel('time')
plt.show()
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import time
# alpha = exp(-1/10) for tau=10 and time step of 1
alpha = 0.904837
# process gain
Kp = 0.8
def process(y0, u0, Kp=Kp, alpha=alpha):
ynext = alpha * y0 + Kp*(1-alpha) * u0
return ynext
# horizons
time_horizon = 200
control_horizon = 100
# set point
sp = np.zeros(time_horizon)
sp[10:50] = 50
sp[50:100] = 100
sp[100:150] = 25
sp[150:] = 75
# initial condition
y0 = 0
# predictive objective function
## "hat" is used for all predicted values
## "current" is used to track current values
## u(i) is used to calculate y(i+1), u[0] is predicted and used to calculate y[1]
### len(u) = control_horizon or time_horizon
### len(y) = control_horizon+1 or time_horizon+1
def objective(u_hat, y_current, sp_current):
y_hat_current = y_current
y_hat = np.zeros(control_horizon+1)
y_hat[0] = y_hat_current
for i in range(control_horizon):
y_hat_next = process(y_hat_current, u_hat[i])
y_hat[i+1] = y_hat_next
y_hat_current = y_hat_next
sp_hat = sp_current*np.ones(len(y_hat)) # optimize to current setpoint, not future setpoints
sqerror = np.sum((y_hat[1:] - sp_hat[1:])**2) # only optimize future states, not current state
return sqerror
# input bounds
u_lo = 0
u_hi = 100
# set current state values; current state = initial state
# i = 0
y_current = y0
# storage of real process values
u = []
y = []
trialsteps = time_horizon
for i in range(trialsteps):
print(f"\nTimestep {i+1} of {trialsteps}")
# find optimal u_hat values
u_hat_guess = np.zeros(control_horizon)
sp_current = sp[i]
print(f'Initial Objective Squared Error: {objective(u_hat_guess, y_current, sp_current)}')
start = time.time()
sol = minimize(objective, u_hat_guess, args=(y_current, sp_current), bounds=[(u_lo, u_hi)]*len(u_hat_guess))
u_hat = sol.x
end = time.time()
print(f"Optimization Time: {end - start} s")
print(f'Final Objective Squared Error: {objective(u_hat, y_current, sp_current)}')
u_pred = u_hat[0]
print(f"Next u: {u_pred}")
u.append(u_pred)
y.append(y_current)
# PLOT CURRENT STATE
# Clear previous plots
plt.clf()
# Subplot 1 of 2 - Output
plt.subplot(2,1,1)
# Subplot 1: Plot y-values
plt.plot(y, 'b-', label='y(t)')
# Subplot 1: Plot predicted y values based on predicted u_hat inputs
y_pred = [y_current] # Start from current y value
y_temp = y_current # temporary y value is used to calculate y_pred
for u_val in u_hat:
y_temp = process(y_temp, u_val)
y_pred.append(y_temp)
t_pred = np.arange(i, i+len(y_pred))
plt.plot(t_pred, y_pred, 'g--', label='Predicted y')
# Subplot 1: Plot setpoint
plt.plot(sp[:len(y)], 'r-', label='Setpoint')
plt.plot(t_pred, sp_current*np.ones(len(y_pred)), 'r--')
plt.axvline(x=i, color='k', linestyle='-', linewidth=0.75) # Add vertical line at current timestep
plt.legend()
plt.grid(True)
plt.ylabel('Output')
# Subplot 2 of 2 - Input
plt.subplot(2,1,2)
# Subplot 2: Plot u-values
plt.step(u, 'b-', label='u(t)')
# Subplot 2: Plot predicted u_hat starting from current timestep
t_pred = np.arange(i, i+len(u_hat))
plt.step(t_pred, u_hat, 'g--', label='Predicted u')
plt.axvline(x=i, color='k', linestyle='-', linewidth=0.75) # Add vertical line at current timestep
plt.legend()
plt.grid(True)
plt.ylabel('Input')
plt.xlabel('Timesteps')
plt.draw()
plt.pause(0.1)
# update current state
y_next = process(y_current, u_pred)
y_current = y_next
plt.ioff() # Turn off interactive mode
plt.show() # Keep final plot visible
Thanks to Jonathan Pershing for generating the Scipy code.