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 packages
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 numpy as np
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.

Course Information

Assignments

Projects

Exams

Dynamic Modeling

Equipment Design

Control Design

Optimal Control

Related Courses

Admin

Streaming Chatbot
💬