"Finding Ultimate Limits of Performance for Hybrid Electric Vehicles"[1]

An exercise in recreating the results from Tate and Boyd's Paper

Problem Description
Proposed Solution
Data Sources
Solution for Hybrid Vehicle
Visualization for Hybrid Vehicle
Solution for Internal Combustion Vehicle(without battery)
Visualization for Internal Combustion Vehicle(without battery)
Analysis and Conclusions
References

By Sindhura Buggaveeti and Allison Waters

Problem Description

  In this project, the objective was to reproduce the results from Tate
and Boyd's paper, "Finding ultimatelimits of performance for hybrid
electric vehicles" [1]. The article proses problem of finding optimal
engine operation in a pure series hybrid vehicle over a fixed drive cycle
subject to a number of practical contraints such as nonlinear fuel/power
maps, minimum and maximum battery charge, battery efficiency, nonlinear
vehicle dynamics and losses, drive train efficiency, engine slew rate
limits.

Nonlinear fuel/power maps
Min and max battery charge
Battery efficiency
Nonlinear vehicle dynamics and losses
Drive train efficiency
Engine slew rate limits

The problem of optimizing fuel efficiency is initially formulated as a
nonlinear convex optimization problem. This convex problem is then
accurately approximated as a large linear program. This optimal solution
is the lower limit of fuel consumption that any control law can achieve
for the given drive cycle and vehicle.

Proposed Solution

Objective function for optimal fuel economy can be described as:

$min \sum_{j=k_0}^{k_f} {f(k)}$

Following are the equality and inequality constraints for the proposed optimization problem.

$$f(k)=f_e(P_e(k))$$

$$E(k+1)-E(k)=f_b(P_s(k))$$

$$E(k_0)=E(k_f)$$

$$P_e(k)=P_s(k)+P_m(k)+P_b(k)$$

$P_e(k+1)-P_e(k)\leq Max Engine Slew Up$

$P_e(k+1)-P_e(k)\geq Max Engine Slew Down$

$P_e(k)\geq 0$

$P_e(k)\leq Max Engine Power$

$P_b(k)\geq 0$

$P_b(k)\leq Max Braking Power$

$P_s(k)\geq Max Discharge Rate$

$P_s(k)\leq Max Charge Rate$

$E(k)\geq Min Battery Energy$

$E(k)\leq Max Battery Energy$

The parameters and variable described in the above algorithm are:

* Variables: $$f(k), E(k), P_e(k), P_b(k), P_s(k)$$
* Parameters: Max Engine Slew Up, Max Engine Slew Down, Max Engine
Power, Max Braking Power, Max Discharge and charge rate, Min and Max
BatteryEnergy

   The optimization problem above was converted into a convex linear
program and solved using Grant and Boyd's convex optimization software
for Matlab, CVX [2][3].

Data Sources

  To illustrate the results obtained using this method, a numerical
example was performed. The drive cycle is first 1371 seconds of FTP-75
Urban Dynamometer Driving Schedule'(UDDS) [4]. The vehicle mass is 1072kg
with a coefficient of drag(CD) of 0.3 and frontal area of 1.96m^2. The
coefficient of rolling resistance is 0.015. Air density is 1.22kg/m^3.
The vehicle is modeled with a 2kW battery load. The fuel consumption
curve for engine is approximated by

$f(t)\geq P_e(t)*0.059+0.075$

$f(t)\geq P_e(t)*0.096-1.041$

where f(t) is in unit of grams per second and P_e is in units of W.

The battery is lead acid(PbA) with 0.60 kW-hr maximum energy capcity and
a reserve of 0.12 KW-hrs. The maximum charge and discharge rates are
9.54kW. The battery has a charge efficiency of 80% and discharge
efficiency of 100% as given in the paper.

Solution for Hybrid Vehicle

Input Parameters

The mechanical power generated by the vehicle is calculated by the
equation below using the following inputs given in the paper.

$P_m=(1/2)*p*C_d*A*v^2+m*g*(C_r+a/g)*v/1000+2$

Vel=xlsread('FTP.xlsx');    %FTP drive cycle data [4]
m=1072;     %Mass of vehicle
Cd=0.3;     %Coefficient of drag
A=1.96;     %Frontal Area
Cr=0.015;   %Rolling resistance coefficient
p=1.22;     %Density of air
g=9.81;     %Gravity
k=1371;     %Number of discrete time steps (each time step is 1 second)
V=0;        %Initial velocity
g1=0.8;     %Battery charge efficiency
g2=1;       %Battery discharge efficiency
ve=0.44704*Vel;  %Converting FTP velocity from mph to m/s
for i=1:k
    %Vehicle acceleration using derivative approximation
    a(i,1)=ve(i)-V;
    %Updating V to be the previous velocity for the next acceleration
    %calculation
    V=ve(i);
    %Calculating mechanical power required to drive the vehicle
    Pm(i,1)=(0.5*p*Cd*A*(ve(i,1)^2)+m*g*(Cr+a(i,1)/g))*ve(i,1)/1000+2;
    time(i)=i;          %Time in seconds
end

Matrices and arrays

one=ones(1,k);              %lxk matrix of ones
none=one(1,1:k-1);         %1xk-1 matrix of ones
%D is an "identity" matrix with -1 along diagonal and 1 along diagonal one
%above the center diagonal
D=-diag(one)+diag(none,1);
Z=zeros(k);                 %kxk matrix of zeros
T=eye(k);                   %kxk identity matrix
Onit=ones(k,1);             %kx1 matrix of ones
ED=D;               %ED is a kxk identity matrix with the element (k,k)=0
ED(k,k)=0;
c=zeros(1,8*k);
c(1,1:k)=one;

Generating equality constraint matrices

The matrix equality constraint Ax==b is made from the equality equations
where A is a (4*k+1)x11*k matrix and b is a (4*k+1)x1 matrix. Slack and surplus
variables were added to 2 equations to convert the inequality
constraints that could cause the set of feasible solutions to degenerate
the rank of the constraints into equality constraints. The storage power
of the battery (Ps) was split into positive (Psp) and negative values
(Psn) so that the variables, Psp and Psn, could be positive. It is
required that all variables in the linear program be greater than zero.
The battery power was determined from the equation below.

$P_s=P_s^+-P_s^-$

By substituting the equation for battery power into the constraints,
Ax==b consisted of the following 5 equations:

$f-0.059P_e-s_1=0.075$

$$-f+0.096P_e+s_2=1.041667$$

$P_e+P_s^--P_s^+-P_b=P_m$

$E(k+1)-E(k)-0.8P_s^+(k)+P_s^-(k)=0$

$E(1)-E(1371)=0$

A1=[T,-0.059.*T,-T,Z,Z,Z,Z,Z;-T,0.096.*T,Z,T,Z,Z,Z,Z;Z,T,Z,Z,T,-T,-T,Z;Z,Z,Z,Z,T,-0.8.*T,Z,ED];
A=zeros(4*k+1,8*k);
A(4*k+1,7*k+1)=1;
A(4*k+1,8*k)=-1;
A(1:4*k,1:8*k)=A1+A(1:4*k,1:8*k);

b=zeros(4*k+1,1);               %Initializing b
b(1:k,1)=0.075.*Onit;           %Equation 1 right side
b(k+1:2*k,1)=1.041667.*Onit;    %Equation 2 right side
b(2*k+1:3*k,1)=Pm;              %Equation 3 right side

Generating inequality constraint matrices

The matrix inequality constraint Gx<=H was made from the inequality
equations where G is a 7*kx11*k matrix and H is a 7*kx1 matrix. The
inequality constraints consist of the upper bound limits on the four
parameters Pe, Psp, Psn, Pb, E, and the upper bound limits on the rate
of change in engine supply power (Pe). In total, 7 equations were
combined to make the matrix inequality constraint.

$P_e\leq50$

$P_s^+\leq9.54$

$P_s^-\leq9.54$

$P_b\leq13.3625$

$E\leq2160$

$P_e(k+1)-P_e(k)\leq10$

$P_e(k)-P_e(k+1)\leq20$

G=[Z,T,Z,Z,Z,Z,Z,Z;Z,Z,Z,Z,T,Z,Z,Z;Z,Z,Z,Z,Z,T,Z,Z;Z,Z,Z,Z,Z,Z,T,Z;Z,Z,Z,Z,Z,Z,Z,T;Z,D,Z,Z,Z,Z,Z,Z;Z,-D,Z,Z,Z,Z,Z,Z];
H=[50.*Onit;9.54.*Onit;6.875.*Onit;13.585.*Onit;2160.*Onit;10.*Onit;20.*Onit];

CVX Code

With the additions of slack and surplus variables, the total number of
variables in the convex optimization problem is 10968 with 8 vectors
each with 1371 elements (the vectors are fuel consumption (f), engine
power (Pe), surplus variable (s1), slack variable (s2), negative storage
power (Psn), positive storage power (Psp), brake power (Pb), battery
energy (E)). The problem is a linear program optimization problem,
therefore, all of the variables must be positive. The battery energy has
a minimum value constraint that it must be greater than 432
kiloWatt-seconds.

cvx_begin
    variables x(8*k,1)
    minimize c*x
    subject to
        A*x==b
        G*x<=H
        x>=0
        x(7*k+1:8*k,1)>=432
cvx_end

 
Calling SDPT3 4.0: 23305 variables, 6853 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints = 6853
 dim. of linear var  = 21935
 dim. of free   var  = 1370
 *** convert ublk to linear blk
********************************************************************************************
   SDPT3: homogeneous self-dual path-following algorithms
********************************************************************************************
 version  predcorr  gam  expon
    NT      1      0.000   1
it pstep dstep pinfeas dinfeas  gap     mean(obj)    cputime    kap   tau    theta
--------------------------------------------------------------------------------------------
 0|0.000|0.000|8.7e+03|5.6e+03|2.7e+11| 1.868927e+08| 0:0:00|2.7e+11|1.0e+00|1.0e+00| spchol 1   1
 1|0.000|0.000|8.7e+03|5.6e+03|2.7e+11| 1.868928e+08| 0:0:00|2.7e+11|1.0e+00|1.0e+00| spchol 1   1
 2|0.000|0.000|8.7e+03|5.6e+03|2.7e+11| 1.868930e+08| 0:0:00|2.7e+11|1.0e+00|1.0e+00| spchol 1   1
 3|0.000|0.000|8.7e+03|5.6e+03|2.7e+11| 1.868942e+08| 0:0:00|2.7e+11|1.0e+00|1.0e+00| spchol 1   1
 4|0.000|0.000|8.7e+03|5.6e+03|2.7e+11| 1.868992e+08| 0:0:00|2.7e+11|1.0e+00|1.0e+00| spchol 1   1
 5|0.000|0.000|8.7e+03|5.6e+03|2.7e+11| 1.869138e+08| 0:0:00|2.7e+11|1.0e+00|1.0e+00| spchol 1   1
 6|0.000|0.000|8.7e+03|5.6e+03|2.7e+11| 1.869579e+08| 0:0:00|2.7e+11|1.0e+00|1.0e+00| spchol 1   1
 7|0.001|0.001|8.6e+03|5.6e+03|2.7e+11| 1.870899e+08| 0:0:00|2.6e+11|1.0e+00|1.0e+00| spchol 1   1
 8|0.003|0.003|8.6e+03|5.6e+03|2.6e+11| 1.874855e+08| 0:0:00|2.6e+11|1.0e+00|1.0e+00| spchol 1   1
 9|0.010|0.010|8.6e+03|5.5e+03|2.6e+11| 1.886679e+08| 0:0:00|2.6e+11|1.0e+00|9.9e-01| spchol 1   1
10|0.030|0.030|8.4e+03|5.4e+03|2.6e+11| 1.921746e+08| 0:0:00|2.6e+11|9.9e-01|9.6e-01| spchol 1   1
11|0.088|0.088|7.9e+03|5.1e+03|2.6e+11| 2.023113e+08| 0:0:00|2.4e+11|9.8e-01|8.9e-01| spchol 1   1
12|0.253|0.253|6.5e+03|4.2e+03|2.2e+11| 2.288910e+08| 0:0:00|1.9e+11|9.4e-01|7.0e-01| spchol 1   1
13|0.630|0.630|2.9e+03|1.8e+03|1.1e+11| 2.731481e+08| 0:0:00|7.7e+10|9.0e-01|3.0e-01| spchol 1   1
14|0.966|0.966|1.0e+02|6.8e+01|4.6e+09| 2.682438e+08| 0:0:00|1.4e+09|9.1e-01|1.1e-02| spchol 1   1
15|0.764|0.764|2.3e+01|1.5e+01|8.5e+08| 6.016079e+07| 0:0:00|3.7e+07|1.1e+00|2.9e-03| spchol 1   1
16|0.790|0.790|5.0e+00|3.2e+00|1.6e+08| 1.204795e+07| 0:0:00|1.0e+06|1.3e+00|7.4e-04| spchol 1   1
17|0.808|0.808|1.1e+00|6.2e-01|2.6e+07| 4.318146e+06| 0:0:01|2.5e+04|1.5e+00|1.7e-04| spchol 1   1
18|0.853|0.853|1.0e+00|9.8e-02|3.7e+06| 9.697342e+05| 0:0:01|6.4e+02|1.8e+00|3.1e-05| spchol 1   1
19|0.966|0.966|3.3e-01|5.3e-03|2.1e+05| 5.125682e+04| 0:0:01|2.2e+02|2.0e+00|1.7e-06| spchol 1   1
20|0.960|0.960|1.3e-02|6.6e-04|8.6e+03| 4.064272e+03| 0:0:01|2.3e+01|2.0e+00|7.1e-08| spchol 1   1
21|0.905|0.905|4.1e-03|5.1e-04|4.5e+03| 3.207877e+03| 0:0:01|2.8e+00|2.0e+00|3.2e-08| spchol 1   1
22|1.000|1.000|1.3e-03|4.0e-04|1.5e+03| 2.473243e+03| 0:0:01|3.6e-01|2.0e+00|9.4e-09| spchol 1   1
23|0.966|0.966|5.4e-04|3.5e-04|8.2e+02| 2.331539e+03| 0:0:01|1.3e-01|2.0e+00|5.2e-09| spchol 1   1
24|1.000|1.000|1.3e-04|3.1e-04|1.9e+02| 2.197099e+03| 0:0:01|6.6e-02|2.0e+00|1.2e-09| spchol 1   1
25|0.938|0.938|9.5e-05|2.8e-04|1.1e+02| 2.185700e+03| 0:0:01|1.9e-02|2.0e+00|6.8e-10| spchol 1   1
26|1.000|1.000|2.7e-05|2.5e-04|2.7e+01| 2.169412e+03| 0:0:01|8.6e-03|2.0e+00|1.8e-10| spchol 1   1
27|1.000|1.000|6.2e-06|2.2e-04|8.0e+00| 2.165330e+03| 0:0:01|2.2e-03|2.0e+00|5.1e-11| spchol 1   1
28|0.900|0.900|1.5e-06|3.0e-05|1.8e+00| 2.162092e+03| 0:0:01|8.0e-04|2.0e+00|1.1e-11| spchol 1   1
29|0.971|0.971|3.5e-07|1.2e-06|3.1e-01| 2.161492e+03| 0:0:01|1.6e-04|2.0e+00|1.9e-12| spchol 1   1
30|0.983|0.983|3.1e-08|5.5e-08|5.9e-03| 2.161411e+03| 0:0:01|2.7e-05|2.0e+00|3.3e-14| spchol 1   1
31|0.994|0.994|1.3e-08|3.8e-09|1.2e-04| 2.161409e+03| 0:0:01|6.3e-07|2.0e+00|7.6e-16| spchol 1   1
32|0.997|0.997|1.3e-08|3.6e-10|2.2e-06| 2.161409e+03| 0:0:01|1.2e-08|2.0e+00|2.7e-17|
  Stop: max(relative gap,infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 32
 primal objective value =  2.16140926e+03
 dual   objective value =  2.16140927e+03
 gap := trace(XZ)       = 2.17e-06
 relative gap           = 1.00e-09
 actual relative gap    = -1.01e-09
 rel. primal infeas     = 1.27e-08
 rel. dual   infeas     = 3.61e-10
 norm(X), norm(y), norm(Z) = 3.7e+01, 4.6e+04, 6.5e+04
 norm(A), norm(b), norm(C) = 1.2e+07, 3.7e+01, 8.2e+04
 Total CPU time (secs)  = 0.88  
 CPU time per iteration = 0.03  
 termination code       =  0
 DIMACS: 1.3e-08  0.0e+00  3.6e-10  0.0e+00  -1.0e-09  5.0e-10
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +463.909

Renaming variables

f=x(1:k,1)*3.6;         %Fuel consumption [kg/hr]
Pen=x(k+1:2*k,1);       %Engine power [kW]
s1=x(2*k+1:3*k,1);      %Surplus variable
s2=x(3*k+1:4*k,1);      %Slack variable
Psn=x(4*k+1:5*k,1);     %Negative component of battery power [kW]
Psp=x(5*k+1:6*k,1);     %Positive component of battery power [kW]
Ps=Psp-Psn;             %Battery power [kW]
Pb=x(6*k+1:7*k,1);      %Brake power [kW]
E=x(7*k+1:8*k,1)/3600;  %Battery energy [kW-hr]

Visualization for Hybrid Vehicle

Generating plots

for i=1:k
    if Pm(i,1)<-6.875       %Inverter has a lower limit of -6.875 kW
        Pm(i,1)=-6.875;
    end
end
figure;
plot(Pm)
title('Drive Inverter Power for Hybrid Vehicle')
ylabel('Power [kW]')
xlabel('Time [s]')
axis([0 1400 -50 50])
figure;
plot(Pen)
title('Spark Ignition Power for Hybrid Vehicle')
ylabel('Power [kW]')
xlabel('Time [s]')
axis([0 1400 0 25])
figure;
plot(f)
title('Spark Ignition Fuel Use for Hybrid Vehicle')
ylabel('Fuel Usage [kg/hr]')
xlabel('Time [s]')
axis([0 1400 0 6])
figure;
plot(Pb)
title('Brake Power for Hybrid Vehicle')
ylabel('Power [kW]')
xlabel('Time [s]')
axis([0 1400 -50 50])
figure;
plot(Ps)
title('PbA Battery Power for Hybrid Vehicle')
ylabel('Power [kW]')
xlabel('Time [s]')
axis([0 1400 -10 10])
figure;
plot(E)
title('PbA Battery Energy for Hybrid Vehicle')
ylabel('Power [kW-hr]')
xlabel('Time [s]')
axis([0 1400 0.1 0.4])

Calculating optimal global fuel economy

%Integrating fuel consumption [g/s] to total fuel [g]
fuel=trapz(time,x(1:k,1));
%Integrating velocity [m/s] to total distance travelled [m]
distance=trapz(time,ve);
%Calculating distance travelled per total fuel consumed and converting to
%miles per gallon.
MPG=(distance*1000/fuel)*0.734/(0.2642*1609.34);
display(MPG)

MPG =

   44.6356

Solution for Internal Combustion Vehicle(without battery)

Generating equality constraint matrices

The solution for the internal combustion vehicle uses the same input
parameters and matrices and arrays as initialized above. The matrix
equality constraint Ax==b is made from the equality equations where A is
a 3*kx11*k matrix and b is a 3*kx1 matrix. Slack and surplus variables
were added to 2 equations to convert the inequality constraints into
equality constraints. Ax==b consists of 3 equations:

$f-0.059*P_e-s_1=0.075$

$$-f+0.096*P_e+s_2=1.041667$$

$P_e-P_b=P_m$

A=[T,-0.059.*T,-T,Z,Z;-T,0.096.*T,Z,T,Z;Z,T,Z,Z,-T];

b=zeros(3*k,1);                 %Initializing b
b(1:k,1)=0.075.*Onit;           %Equation 1 right side
b(k+1:2*k,1)=1.041667.*Onit;    %Equation 2 right side
b(2*k+1:3*k,1)=Pm;              %Equation 3 right side

c=zeros(1,5*k);
c(1,1:k)=one;