ElectricVehicleCruiseControl

Title: Modeling and PID Controller Example - Cruise Control for an Electric Vehicle

Author: Alex Dowling

Date Prepared: May 3, 2009

= Introduction =

Controls principles developed in this course can be applied to non-chemical engineering systems such as automobiles. Some companies, such as NAVTEQ, are developing adaptive cruise control products that use information about the upcoming terrain to shift gears in a more intelligent manner which improves speed regulation and fuel economy. This case study will examine the basics of develop a speed controller for an electric vehicle.

An electric vehicle was chosen for the following reasons:
 * Electric vehicles are interesting from an engineering perspective and may become a reality for consumers in the future
 * Torque produced by an electric motor is instantaneous (for all practical purposes). Thus actuator lag can be ignored, simplifying the development of said controller.
 * Some electric vehicles feature motors directly integrated into the hub of the drive wheel(s). This eliminates the need for a transmission and simplifies vehicle dynamics models.

= Forces =

As shown in the free body diagram below, there are six forces acting on the vehicle:
 * 1) Rolling Resistance
 * 2) Aerodynamic Drag
 * 3) Aerodynamic Lift
 * 4) Gravity
 * 5) Normal
 * 6) Motor



Figure 1: A free body diagram of the forces acting on the vehicle. $$-\theta\,\!$$ is used to denote the grade of the road such that a positive value of $$\theta\,\!$$ corresponds to the vehicle traveling uphill.

Rolling Resistance
Rolling resistance is due the tires deforming when contacting the surface of a road and varies depending on the surface being driven on. It can be model using the following equation:

$$ F_{RR} = C_{rr1} v + C_{rr2} F_{N}\,\!$$

The two rolling resistance constants can be determined experimentally and may be provided by the tire manufacture. $$F_{N}\,\!$$ is the normal force.

Aerodynamic Drag
Aerodynamic drag is caused by the momentum loss of air particles as they flow over the hood of the vehicle. The aerodynamic drag of a vehicle can be modeled using the following equation:

$$F_{drag} = (1/2) \rho CdA v^2\,\!$$
 * $$\rho\,\!$$ is the density of air. At 20°C and 101kPa, the density of air is 1.2041 kg/m3.
 * $$CdA\,\!$$ is the coefficient of drag for the vehicle times the reference area. Typical values for automobiles are listed here.
 * $$v\,\!$$ is the velocity of the vehicle.

Aerodynamic Lift
Aerodynamic lift is caused by pressure difference between the roof and underside of the vehicle. Lift can be modeled using the following equation:

$$F_{lift} = (1/2) \rho ClA v^2\,\!$$
 * $$\rho\,\!$$ is the density of air. At 20°C and 101kPa, the density of air is 1.2041 kg/m3.
 * $$ClA\,\!$$ is the coefficient of lift for the vehicle times the reference area.
 * $$v\,\!$$ is the velocity of the vehicle.

Gravity
In the diagram above, there is a component of gravity both in the dimension normal to the road and in the dimension the vehicle is traveling. Using simple trigonometry, the component in the dimension of travel can be calculated as follows:

$$F_{G, travel} = mg \sin(-\theta)\,\!$$
 * $$m\,\!$$ is the mass of the vehicle.
 * $$g\,\!$$ is the acceleration due to gravity.

Normal Force
The normal force is the force excerted by the road on the vehicle's tires. Because the vehicle is not moving up or down (relative to the road), the magnitude of the normal forces equals the magnitude of the force due to gravity in the direction normal to the road.

$$F_{N} = F_{G, norm.} - F_{lift} = mg\cos(-\theta) - (1/2) \rho ClA v^2\,\!$$


 * $$m\,\!$$ is the mass of the vehicle.
 * $$g\,\!$$ is the acceleration due to gravity.

Motor
The torque produced by an electric motor is roughly proportional to the current flowing through the stater of the motor. In this case study, the current applied to the motor will be controlled to regulate speed. Applying a negative current will cause the vehicle to regeneratively brake.

$$\tau = k_{motor} I\,\!$$

$$ F_{M} = \frac{\tau}{r} = \frac{k_{motor} I}{r}\,\!$$
 * $$\tau\,\!$$ is the torque produced by the motor.
 * $$I\,\!$$ is the current flowing through the motor.
 * $$r\,\!$$ is the radius of the tire.
 * $$k_{motor}\,\!$$ is a constant.

= Newton's Second Law =

Using Newton's Second Law, a differential for the vehicle's speed can be obtained.

$$ ma = \sum F = F_{M} - F_{drag} - F_{RR} + F_{G, travel}$$

Substituting in the expressions for various forces detailed above yields the following:

$$ ma = \frac{\tau}{r} - (1/2) \rho CdA v^2 - C_{rr1}v - C_{rr2} F_{N} + mg\sin(-\theta)$$

Further substituting the expression for normal forces yields the following:

$$ ma = \frac{\tau}{r} - (1/2) \rho CdA v^2 - C_{rr1}v - C_{rr2}[mg\cos(-\theta) - (1/2) \rho ClA v^2 ] + mg\sin(-\theta)$$

Substituting $$\frac{dv}{dt}$$ for $$a\,\!$$ results in the following:

$$ m \frac{dv}{dt} = \frac{\tau}{r} - (1/2) \rho CdA v^2 - C_{rr1}v - C_{rr2}[mg\cos(-\theta) - (1/2) \rho ClA v^2 ] + mg\sin(-\theta)$$

Grouping like turns results in the following:

$$ m \frac{dv}{dt} = \frac{\tau}{r} - [(1/2) \rho CdA - C_{rr2}(1/2) \rho ClA] v^2 - C_{rr1}v - C_{rr2} mg\cos(-\theta) + mg\sin(-\theta)$$

In order to simply the remaining analysis, several constants are defined as follows:

$$\alpha = (1/2) \rho r (CdA - C_{rr2} ClA)\,\!$$

$$\beta = r C_{rr1}\,\!$$

$$\gamma = r C_{rr2} mg\cos(-\theta) - rmg\sin(-\theta)\,\!$$

Substituting these into the differential equation results in the following expression:

$$ \frac{dv}{dt} = \frac{\tau - \alpha v^2 - \beta v - \gamma(\theta)}{m r}$$

It is important that $$\theta\,\!$$ (and thus also $$\gamma\,\!$$) is a function of vehicle position.

= PID Controller = One way to regulate vehicle speed is to control the torque generated by the electric motor. If the motor torque is greater than the resistive torque acting on the vehicle (a summation of aerodynamic drag, rolling resistance, etc.) the vehicle will accelerate. If the motor torque is less than the resistive torque, the vehicle will slow down.

Expression for Phase Current
For an electric motor, the phase current flowing through the motor is proportional to the torque produced. Thus one strategy for controlling the vehicles speed is to controller the motor phase current.

Using a PID controller architecture, the expression for motor current is the following:

$$I = K_{c}(v_{set}-v) + \frac{1}{\tau_{I}} \int {(v-vset)dt} + \tau_{D} \frac{d(v_{set}-v)}{dt} + C_{offset}$$

Differential Equation for Velocity
Substituting this expression into the differential equation for vehicle position results in the following:

$$ \frac{dv}{dt} = \frac{k_{motor}[K_{c}(v_{set}-v) + \frac{1}{\tau_{I}} \int {(v_{set}-v)dt} + \tau_{D} \frac{d(v_{set}-v)}{dt} + C_{offset}]- \alpha v^2 - \beta v - \gamma(\theta)}{m r}$$

Defining another variable $$x_{1}\,\!$$ allows for the removal of the integral from the expression.

$$ \frac{dx_{1}}{dt} = v_{set}-v$$

$$ \frac{dv}{dt} = \frac{k_{motor}[K_{c}(v_{set}-v) + \frac{x_{1}}{\tau_{I}} + \tau_{D} \frac{d(v_{set}-v)}{dt} + C_{offset}]- \alpha v^2 - \beta v - \gamma(\theta)}{m r}$$

If all changes in $$v_{set}\,\!$$ are gradual then $$\frac{dv_{set}}{dt} \approx 0$$. Applying this simplification results in the following expression:

$$ \frac{dv}{dt} = \frac{k_{motor}[K_{c}(v_{set}-v) + \frac{x_{1}}{\tau_{I}} - \tau_{D} \frac{dv}{dt} + C_{offset}]- \alpha v^2 - \beta v - \gamma(\theta)}{m r}$$

Using a little algebra, an expression for $$\frac{dv}{dt}$$ can be obtained as follows:

$$\frac{dv}{dt} = \frac{k_{motor} K_c (v_{set}-v)+\frac{k_{motor}}{\tau_I} x_1+k_{motor} C_{offset}-\alpha v^2 -\beta v - \gamma(\theta)}{k_{motor} t_{D}+m r}$$

Find Fixed Point
The fixed point can be obtained by setting the derivatives to zero and solving the system of equation.

System of Equations:

$$\frac{dv}{dt} = \frac{k_{motor} K_c (v_{set}-v)+\frac{k_{motor}}{\tau_I} x_1+k_{motor} C_{offset}-\alpha v^2 -\beta v - \gamma(\theta)}{k_{motor} t_{D}+m r} = 0$$

$$\frac{x_1}{dt} = v_{set}-v = 0$$

Solution:

$$x_{1} = \frac{-\tau_{I} (k_{motor} C_{offset}- \alpha v_{set}^2-\beta v_{set}-\gamma(\theta)])}{k_{motor}}$$

$$v = v_{set}\,\!$$

As expected, a fixed point exists when the set velocity equals the actual velocity of the vehicle.

Linearize System of ODEs
Before the stability of the system can be (easily) examined, the system must be linearized around a fixed point.

Note: To simply the matrix expressions in section the following notation will be used:

$$y'_{i} = \frac{dy_{i}}{dt}$$

Overall a linearized system of ODEs has the following form:

$$\begin{bmatrix} y'_1 \\ y'_2 \\ \vdots \\ y'_n \\ \end{bmatrix} = \mathbf{J} \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \\ \end{bmatrix} + \begin{bmatrix} k_1 \\ k_2 \\ \vdots \\ k_n \\ \end{bmatrix}$$

The first step of linearizing any system of ODEs to calculate the Jacobian. For this particular system, the Jacobian can be calculated as follows:

$$\mathbf{J} = \begin{bmatrix} \frac{\partial v'}{\partial v} & \frac{\partial v'}{\partial x_1} \\ \frac{\partial x'_1}{\partial v} & \frac{\partial x'_1}{\partial x_1} \\ \end{bmatrix} = \begin{bmatrix} \frac{-k_{motor}K_c-mr[2\alpha v + \beta]}{mr+k_{motor} \tau_D} & \frac{k_{motor}}{t_I (mr+k_{motor} \tau_D)} \\ -1 & 0 \\ \end{bmatrix}$$

The Jacobian is then evaluated at the fixed point:

$$\mathbf{J} = \begin{bmatrix} \frac{-k_{motor}K_c-mr[2\alpha v_{set} + \beta]}{mr+k_{motor} \tau_D} & \frac{k_{motor}}{t_I (mr+k_{motor} \tau_D)} \\ -1 & 0 \\ \end{bmatrix}$$

The next step if to calculate the vector of constants. For this particular system, said vector can be calculated as follows:

$$\begin{bmatrix} k_v \\ k_{x_1} \\ \end{bmatrix} = - \mathbf{J} \begin{bmatrix} v \\ x_1 \\ \end{bmatrix} \Bigg|_{fixed point}

= - \begin{bmatrix} \frac{-k_{motor}K_c-2\alpha v_{set} - \beta]}{mr+k_{motor} \tau_D} & \frac{k_{motor}}{t_I (mr+k_{motor} \tau_D)} \\ -1 & 0 \\ \end{bmatrix} \begin{bmatrix} v_{set} \\ \frac{-\tau_{I} (k_{motor} C_{offset}-\alpha v_{set}^2-\beta v_{set}-\gamma(\theta)])}{k_{motor}} \\ \end{bmatrix} $$

$$\begin{bmatrix} k_v \\ k_{x_1} \\ \end{bmatrix} = \begin{bmatrix} \frac{-(-k_{motor} K_c -2 \alpha v_{set} - \beta] ) v_{set}}{(mr+k_{motor} \tau_{D})}+\frac{(k_{motor} C_{offset}-\alpha v_{set}^2-\beta v_{set} - \gamma(\theta)])}{(mr+k_{motor} \tau_D)} \\ v_{set} \\ \end{bmatrix} $$

Combining the Jacobian and vector of constants results in the following linearized system:

$$\begin{bmatrix} v' \\ x'_1 \\ \end{bmatrix} = \mathbf{J} \begin{bmatrix} v \\ x_1 \\ \end{bmatrix} + \begin{bmatrix} k_v \\ k_{x_1} \\ \end{bmatrix}$$

$$\begin{bmatrix} v' \\ x'_1 \\ \end{bmatrix} = \begin{bmatrix} \frac{-k_{motor}K_c-2\alpha v_{set} - \beta}{m+k_{motor}] \tau_D} & \frac{k_{motor}}{t_I (mr+k_{motor} \tau_D)} \\ -1 & 0 \\ \end{bmatrix} \begin{bmatrix} v \\ x_1 \\ \end{bmatrix} + \begin{bmatrix} \frac{-(-k_{motor} K_c -2 \alpha v_{set} - \beta] ) v_{set}}{(mr+k_{motor} \tau_{D})}+\frac{(k_{motor} C_{offset}-\alpha v_{set}^2-\beta v_{set}-\gamma(\theta)])}{(mr+k_{motor} \tau_D)} \\ v_{set} \\ \end{bmatrix}$$

Stability Analysis
To assess the stability of the controller, the eigenvalues of the Jacobian in the linearized systems of ODEs can be examined. In general, an eigenvalue ($$lambda$$) is the solution to the following equation:

$$\big| \mathbf{J} - \lambda \mathbf{I} \big| = 0$$

Using a computer to solve said equation, the eigenvalues of this particular system can be found to be the following:

$$\lambda = \frac{k_{motor} K_c \tau_I+\beta \tau_I+2 \alpha v_{set} \tau_I \pm \sqrt{ \tau_I (k_{motor}^2 K_c^2 \tau_I+2 k_{motor} K_c \tau_I \beta+4 k_{motor} K_c \tau_I \alpha v_{set}+\beta^2 \tau_I+4 \beta \tau_I \alpha v_{set}+4 \alpha^2 v_{set}^2 \tau_I-4 k_{motor} mr-4 k_{motor}^2 t_D)}}{2 \tau_I (mr+k_{motor} \tau_D)}$$

For the system to be stable, the real component of all eigenvalues must be non-positive. The following inequality must be true for a stable controller:

$$k_{motor} K_c \tau_I+\beta \tau_I+2 \alpha v_{set} \tau_I < - Real(\tau_I (k_{motor}^2 K_c^2 \tau_I+2 k_{motor} K_c \tau_I \beta+4 k_{motor} K_c \tau_I \alpha v_{set}+\beta^2 \tau_I+4 \beta \tau_I \alpha v_{set}+4 \alpha^2 v_{set}^2 \tau_I-4 k_{motor} mr-4 k_{motor}^2 t_D))$$

For the system to not oscillate, the imaginary component of all eigenvalues must be zero. The following inequality must be true for a non-oscillating controller: $$ 0 > \tau_I (k_{motor}^2 K_c^2 \tau_I+2 k_{motor} K_c \tau_I \beta+4 k_{motor} K_c \tau_I \alpha v_{set}+\beta^2 \tau_I+4 \beta \tau_I \alpha v_{set}+4 \alpha^2 v_{set}^2 \tau_I-4 k_{motor} mr-4 k_{motor}^2 t_D) $$

Interestingly, neither of these criteria depend on the grade of the road ($$\theta$$). However, during the analysis, it was assumed that $$\theta$$ is constant. For most roads, this is not the case; $$\theta$$ is actually a function of vehicle position. In order to add this additional level of detail, the original system of ODEs needs to be revised:

$$\frac{dv}{dt} = \frac{k_{motor} K_c (v_{set}-v)+k_{motor} x_1+k_{motor} C_{offset}-\alpha v^2-\beta v -\gamma(\theta)}{k_{motor} t_{D}+mr}$$

$$\frac{dx_1}{dt} = v_{set} - v$$

$$\frac{ds}{dt} = v $$

$$\theta = f(s)\,\!$$

Unfortunately for any normal road, the grade is not a simple (or even explicate) function of position (s). This prevents an in depth analytical analysis of stability. However for a very smooth road with very gradual changes in grade, the stability of the controller should be unaffected.

Example Electric Vehicle
Simulating the system also for other properties of the controller to be examined. This section presents the results from simulating a specified fictitious electric vehicle.

Parameters
For the fictitious electric vehicle simulated in this analysis, the following parameters were used. These parameters are roughly based on the parameters one would expect to see in a typical electric automobile.

$$m = 1400 ~ kg\,\!$$

$$CdA = 0.58 ~ m^2\,\!$$

$$ClA = 0 ~ m^2\,\!$$

$$C_{rr1} = 2.75 ~ \frac{N s}{m}0.018\,\!$$

$$C_{rr2} = 0.018\,\!$$

$$k_{motor} = 1.2 ~ \frac{Nm}{A_{rms}}$$

Root Locus Plots
Using the parameters above, Root Locus Plots of the system were constructed to numerically explore the stability. The following controller constants were used when constructing the plots:

$$v_{set} = 25 ~ m/s\,\!$$

$$K_c = 55\,\!$$

$$\tau_I = 0.5\,\!$$

$$\tau_D = 1\,\!$$

$$C_{offset} = 0\,\!$$







These Root Locus plots show that the controller with said constants is both stable and does not oscillate.

Phase Portrait
A phase portrait of the system with the following parameters was also constructed:

$$v_{set} = 25 ~ m/s\,\!$$

$$K_c = 55\,\!$$

$$\tau_I = 0.5\,\!$$

$$\tau_D = 20\,\!$$

$$C_{offset} = 0\,\!$$



The phase portrait shows that the system is both stable and does not oscillate, as predicted by the Root Locus plots.

Driving Simulation on Level Terrain
The vehicle was simulated starting at 10 m/s and accelerating to 25 m/s via cruise control on level terrain ($$\theta=0\,\!$$). For this simulation, the following constants were used:

$$v_{set} = 25 ~ m/s\,\!$$

$$K_c = 55\,\!$$

$$\tau_I = 0.5\,\!$$

$$\tau_D = 20\,\!$$

$$C_{offset} = 146 \approx \frac{\alpha v_{set}^2 + \beta v_{set} + \gamma(\theta = 0)}{k_{motor}}$$



Interestingly this graph shows oscillation, despite the Root Locus plots and phase diagrams. It is important to remember, however, that the Root Locus plot and stability methods involve linearizing the system. It is possible the linearized system is not a good approximation of the system.

It is also important to remember that this particular example involves a large set point change, which can induce oscillations in certain systems.

Driving Simulation on Unlevel Terrain
In order to explore how the controller behaves on a road with a non-zero grade, a route with hills was constructed from the following equation, where h is elevation and s is position.:

$$h = 100 \sin(\frac{s}{250 \pi}) + \frac{s}{2000}$$



The vehicle was simulated driving said road starting with the following initial conditions and controller constants:

Initial Conditions

$$s(t=0) = 0\,\!$$

$$x_1(t=0) = 0\,\!$$

$$v(t=0) = v_{set} = 25 ~ m/s\,\!$$

Controller Constants

$$k_c = 55\,\!$$

$$\tau_I = 0.5\,\!$$

$$\tau_D = 20\,\!$$

$$C_{offset} = 146\,\!$$



The current controller tuning is inadequate for this road. There are vary large variations in velocity. In order to reduce these variation, the propartional gain Kc was increased by a factor of 5. Below are the results:



Using the Optimization Toolbox in MATLAB, the controller was optimized by minimizing the sum of the errors (difference between vehicle velocity and set velocity) on this particular segment of road. Below are the optimized controller constants and a plot of the velocity profile:

$$k_c = 232.58\,\!$$

$$\tau_I = 0.001\,\!$$

$$\tau_D = 0\,\!$$

$$C_{offset} = 114.62\,\!$$



This example shows the power of optimization and model predictive control.

= Summary =

This example demonstrates the following:


 * 1) Modeling a simple, non-chemical engineering system.
 * 2) Develop a PID controller for said physical dynamical system.
 * 3) Manipulate said system to develop an system of differential equations.
 * 4) Find fixed point(s) for said system.
 * 5) Linearize said system of differential equations.
 * 6) Find the eigenvalues of said linearized system.
 * 7) Construct Root Locus plots for said linearized system.
 * 8) Construct a phase portrait for said dynamical system.
 * 9) Simulate said dynamical system under various conditions.
 * 10) Demonstrate the idea of model predictive control by optimize said controller under a specific scenario.