Fuel Minimisation for a Vehicle Equipped with a Flywheel and Battery on a ThreeDimensional Track
Abstract
An optimal control based methodology is proposed for minimising the combustible fuel consumption of a hybrid vehicle equipped with an internal combustion engine, a highspeed flywheel and a battery. The threedimensionality of the road is recognised by the optimal control calculations. Fuel efficiency is achieved by optimally exploiting the primary and secondary energy sources and controlling the vehicle so that the fuel consumption is minimised for a given, but arbitrary threedimensional route. A timeofarrival constraint rather than a driving cycle is used. The benefits of using multiple auxiliary storage systems are demonstrated and a lowerbound estimate of the fuel consumption is presented.
I Introduction
The low efficiency and the airpollution side effects associated with internal combustion engine (ICE) usage are well known. As is now widely appreciated, several of the problems endemic to ICEs can be mitigated using a variety of secondary energy storage technologies such as fuel cells, lithiumbased battery systems, highspeed flywheels and supercapacitors. In broad terms, fuel cells and lithiumbased battery systems have good energy storage properties, while highspeed flywheels and supercapacitors can be utilised for their high power density characteristics [Doucette2011, Guzzella2007]. As a result of their poor power density properties, battery packs in commercial electric vehicles tend to be oversized. Other batteryrelated issues include long charging times and shortened life expectancy, especially when they are cycled at high charge and discharge rates. Supercapacitors and flywheelbased systems are an ideal complement to batteries, because of their high power density characteristics under both charging and discharging. The fundamentals and applications of fuel cells, including the main reactions, are reviewed in [Carrette2000]. A comprehensive review of lithiumion battery technologies is given in [Horiba2014]. The physical structure of some flywheelbased systems is reviewed in a vehicle context in [Dhand2013]. Hybrid energy storage systems (HESS) based on batteries and supercapacitors are reviewed in [Ju2014] and the references therein.
Aside from the technologies themselves, energy storage modelling as well as their optimal deployment are important issues. The modelling and control of hybrid electric vehicles (HEVs) are reviewed in [Shen2011], where various powertrain topologies and control strategies are discussed. It is pointed out that global optimization can be used as a ‘what’s possible’ benchmark for evaluating energy management strategies. Another good survey paper reviewing 180 papers on the optimal energy management of HEVs and plugin HEVs can be found in [Panday2014]. Power management control strategies can be divided into offline and online methodologies. Online management strategies include methods such as lookup tables, state machines, thermostat control, Equivalent Consumption Minimisation Strategies (ECMS) [Paganelli2002, Rousseau2008, Sciarretta2004], Neural networks [Baumann1998], particle swarm optimisation [Wang2006], model predictive control (MPC) [West2003] and fuzzy control [Lee1998]. Offline strategies that often focus on global optimisation include dynamic programming (DP) [Lin2003, Rousseau2008], linear programming [Tate2000], nonlinear programming [Perez2009], Stochastic DP [liu2008, Moura2011, Kim2007] and genetic algorithms [Piccolo2001]. An evolutionary algorithm was applied in [Qi2016] over a sliding window to allow realtime power management and a datadriven reinforcement learning algorithm was proposed in [Qi2016data]. ECMS does not guarantee charge sustainability and hence an Adaptive ECMS (AECMS) algorithm was introduced in [Musardo2005] to update the equivalence factor ‘onthefly’ on the basis of past and predicted driving conditions.
The majority of the work in the literature uses driving cycles as benchmarks for performance evaluation. As was recognised in [Roy2014], optimising vehicle parameters over one drive cycle does not necessarily mean that the vehicle will perform well on other drive cycles. In [Roy2014], a method was proposed for optimising an HEV over a range of drive cycles with different levels of driving aggressiveness and traffic conditions, in order to reduce the fuel economy variability with respect to drive cycle changes. While seeking to address the limitations associated with single driving cycle usage, this method is still restricted by the driving cycle combination used and there were cases where a single driving cycle resulted in a lower fuel variability compared with the proposed multicycle method.
The (combustible) fuel used by any vehicle will depend on the vehicle’s speed, with higher average speeds typically resulting in a higher fuel usage. A contribution of this work is proposing a method of minimising fuel consumption over different driving conditions without using drive cycles. The idea is to simultaneously optimise the powertrainâs energy deployment and driving strategy over a given (but arbitrary) route. Key in this procedure is the selection of a timeofarrival constraint, which acts as an ‘aggressiveness’ surrogate. A short travel time corresponds to aggressive driving and a generally higher fuel consumption. An optimal control algorithm then seeks to minimise the fuel consumption, while ensuring a ‘justintime’ arrival. The optimal control calculation makes use of a realistic vehicle model, with three degrees of freedom, and a nonlinear tire model. It is shown that flywheels are an excellent means of reducing fuel consumption in manoeuvres where high levels of braking are involved, such as extraurban driving on fast roads.
Another distinguishing feature of this work over the majority of the existing literature is that unlike rulebased methods, system dynamics and nonlinear constraints can be included explicitly as part of the optimal control problem and optimal rather than near optimal results are obtained. A pseudospectral method will be used to solve the optimal control problem, which makes use of first and secondorder gradient information for fast convergence. Compared to DP/SDP, models of much higher complexity can be solved using this approach.
In this paper, we will take the roadâs threedimensional geometric profile into consideration. In mechanical systems with comparable potential and kinetic energies, the optimal control strategy is strongly dependent on tradeoffs between the two energy sources, as was demonstrated by the famous ‘Minimum Time to Climb’ problem associated with jetpowered aircraft [Kaiser1944]. In Section II a simple motivating example is given that highlights the importance of changes in the road gradient. In Section III the vehicle and track models employed in this study are described. In Section IV the optimal control problem is cast in a standard form and the numerical method used to solve it is described. Numerical results are presented and discussed in Section V. The conclusions are given in Section VI and the vehicle parameters used in the study are provided in the Appendix.
Ii An Illustrative Example
Imagine a pathconstrained point mass (bead) that is required to reach a given destination within a prespecified time. Suppose the trajectory’s starting point is in absolute Cartesian coordinates and that the path is parabolic and defined by . If the destination point is given, there holds and so the path is specified by a single free parameter. We would like to minimise the energy supplied to the bead in order that it arrives at the destination ‘just in time’. If the effects of air resistance and friction forces are neglected, the bead’s height is indicative of its speed. If the bead has unity mass, conservation of energy dictates
(1) 
The bead’s speed is , its instantaneous height above the origin is and is the amount of external energy supplied. The speed of the bead is thus
(2) 
An infinitesimal path segment can be described as
(3) 
The thrust programme that minimises the total external energy supplied is the optimal control problem of minimising
(6) 
subject to state dynamics
(7) 
where is the externally applied force. The boundary conditions are
(8) 
The control Hamiltonian of this system is
(9) 
Since is free, . For the quadratic path the control Hamiltonian becomes
(10) 
with costate equations
(11) 
where
(12) 
Since the drive force is bounded, Pontryagin’s minimum principle determines that
(13) 
In a manner reminiscent of simple variants of the Goddard rocket problem, the optimal thrust programme turns out to be impulsive [Tsien1951], with all the external energy injected into the system at the beginning of the manoeuvre. To illustrate this the described optimal control problem was solved with GPOPS in the case that and ; the states, control and costates are shown in Fig. 1.
Nine cases are considered that highlight the significance of interactions between the kinetic and potential energy of the bead. These paths correspond to three values for the parameter and three values for the terminal point; the resulting paths are shown in Fig. 2.
For each of the nine cases, Table I shows the path lengths , the zerofuel times (in the cases that a zerofuel journey is possible), and the external energy required to complete the journey under a 1 s timeofarrival constraint.
Case  

1  1  0.1  0.7  5.29  1.57  6.02 
2  1  0  0.2  5.10  2.30  8.56 
3  1  0.1  0.3  5.29    13.51 
4  0  0.1  0.5  5.20  2.38  9.70 
5  0  0  0  5.00  12.50  
6  0  0.1  0.5  5.20    17.73 
7  +1  0.1  0.3  5.29    15.83 
8  +1  0  0.2  5.10    18.37 
9  +1  0.1  0.7  5.29    23.31 
The second column in Table I shows the elevation of the terminal point. To calculate the arrival time, (5) is solved numerically by running a bisection on until the journey time constraint is met.
Superior fuel consumption performance is achieved when the initial part of the journey is downhill, because the early conversion of potential energy into kinetic energy results in higher speeds throughout the journey. Shorter journey paths do not necessarily result in lower fuel consumption (e.g. compare Case 1 and 2). As one would expect, the journeys in Cases 7 to 9 are the most arduous interms of fuel consumption, due to elevated terminal points. This example demonstrates that the path elevation curvature can have a significant influence on fuel usage and sometimes in a counterintuitive manner.
Iii The Mathematical Model
A bicycle model of the car is used that has yaw, lateral and longitudinal freedoms, and nonlinear tires. The road is assumed threedimensional and is represented by a geometric construct called a ‘ribbon’, which is describable in terms of three curvature variables [Perantoni2015Track].
Iiia Track Model
A moving coordinate system (called a Darboux frame) is used to describe the track. As shown in Fig. 3 the origin of this moving system is ‘dragged along’ by the car. The independent variable in the track description is , which is the distance travelled by the car from some starting point, projected onto the track spine. The spine could be, but need not be, the track centre line.
The Darboux frame is described by the orthogonal moving triad , with normal to the road surface. The road is represented locally by the  plane, which ‘travels’ with the car down the road.
The orientation of this local ‘patch’ is described in terms of three Euler angles that are all functions of [Perantoni2015Track]. If the track orientation is described by the roll, yaw and pitch angles, respectively , and , then the track curvatures are given by
(14) 
The angular velocity of the Darboux frame is given by
(15) 
The next kinematic relationships we will require relates to the way in which the car progresses down the road. Suppose that the absolute velocity of the car in its bodyfixed coordinate system is ; the longitudinal velocity component is determined by the throttle/brakes, the lateral component is determined by the steering, while the vertical component is determined by the track characteristics. The car’s geometric centre (the car’s mass centre projected down on to the road) is given by in the Darboux frame. The absolute velocity of the car in the Darboux frame is:
(16) 
where
(17) 
represents the yawing of the car relative to the spine of the track. The first term in (16) derives from the angular velocity of the Darboux frame, while the second is the velocity of the car’s geometric centre expressed in Darboux coordinates. The first row of (16) gives the speed of the origin of the Darboux frame in its tangent direction and can be rewritten as
(18) 
since . The function
(19) 
transforms ‘time’ as the independent variable into the ‘elapsed distance’ as the independent variable; it is assumed that and its inverse are nonzero everywhere on the spine curve. The second row of (16) is
(20) 
which describes the way the vehicle moves normal to the spine. Transforming (20) into the distancetravelled domain gives
(21) 
in which is the derivative of with respect to . Expressing the absolute angular velocity of the car in its bodyfixed reference frame gives
(22) 
This will be used in the next section to derive the vehicle’s equations of motion. The car’s yaw angle in Darboux frame expressed in distance domain is deduced from the third row of (22) by integrating
(23) 
IiiB Dynamics
The equations describing the dynamics of the car are derived using standard vectorial methods. The absolute velocity of the car’s mass centre (expressed on the vehicle’s coordinate system) can be written as
(24) 
where denotes the cross product. The NewtonEuler equations for this system are given by
(25)  
(26) 
where the car’s inertia matrix is assumed to be diagonal and is given by , with and being the external force and moment. The last term in (25) is due to the gravitational acceleration of the car’s mass centre and can be expressed as
(27) 
The car’s equations of motion can now be assembled from (25), (26) and (27) as follows:
(28)  
(29)  
(30) 
in which , are the resultant longitudinal and lateral forces, and is the zaxis tire moment acting on the car. These quantities are given by
(31)  
(32)  
(33) 
The tire force system is illustrated in Fig. 4 and is discussed in Section IIID. The coefficient is the tire rolling resistance coefficient and the last two terms in (31) and (32) represent the rolling resistance forces. The aerodynamic drag force acts in negative axis direction and is given by
(34) 
where is the drag coefficient. The equations of motion (28), (29) and (30), expressed in terms of the elapsed arc length, are as follows
(35)  
(36)  
(37) 
The angular acceleration of the Darboux (road) frame is neglected in the dynamic equations for the examples given here as its influence is negligibly small.
IiiC Load Transfer
In order to compute the tire normal loads we balance the forces acting on the car normal to the road, and then balance moments around the bodyfixed axis (see Fig. 4). These calculations must recognise the gravitational, inertial, centripetal and aerodynamic forces acting on the car as well as the roadrelated gyroscopic moments. The vertical force balance gives
(38) 
in which the and are the front and rear tire normal loads, the second term derives from the axis component of the track’s angular acceleration, the third term is the acceleration due to gravity, while the last two terms are centripetal accelerations.
Balancing the pitching moments around the car’s mass centre gives
(39) 
The first two terms represent the pitching moments produced by the vertical tire forces, the third term is the pitching moment produced by the longitudinal force given in (31), the fourth is the inertial moment around the car’s pitch axis, whilst the fifth term is a gyroscopic moment acting in the car’s pitch direction.
IiiD Tire Forces
The tire forces have normal, longitudinal and lateral components that act on the vehicle chassis at the tire ground contact points and react on the inertial frame. The rearwheel tire force is expressed in the vehicle’s bodyfixed reference frame, while the front tire force is expressed in a steered reference frame; refer again to Fig. 4. We make use of the wellknown Magic Formula tire model [Pacejka2006], where these forces are a function of the normal load and the tire’s longitudinal slip coefficient and a lateral slip angle . The tire equations were also described in details in the Appendix of [Limebeer20153DOptimal]. The same tire parameters are used in this work except that the peak longitudinal and lateral friction coefficients have been scaled down by 30 %. Following standard conventions we use
(40)  
(41) 
where is the wheel radius and the wheel’s spin velocity. The quantities and are the absolute velocity components of the wheel centre in a wheelfixed coordinate system. The front and rear tire lateral slip angles are given by
(42)  
(43) 
IiiE Battery Model
We will use a simple ‘voltage behind output resistance’ model for the (lithiumion) battery, as shown in Fig. 5.
The terminal voltage is
(44) 
where is the battery opencircuit voltage, is the battery current and is the internal resistor. In this convention positive corresponds to discharging the battery, while negative corresponds to charging. The unloaded battery voltage has a dependency on the state of the charge as follows
(45) 
where the is defined as
(46) 
The maximum battery charge is given by
(47) 
where represents the maximum energy storage capacity of the battery. Using (44) and (45) the power delivered, or drawn from the battery is
(48) 
Again negative implies charging and positive implies discharging.
The battery charge can be modelled by the dynamic equation
(49) 
The power transmission between the battery and the rear wheel requires an electric motor, which is assumed to have an efficiency factor . The electric motor output power is thus
(50) 
IiiF Engine Map
The engine used in the work presented here is the 1.5 L Prius engine with maximum power output of 43 kW. The fuel consumption map (see Fig. 6) for this engine was obtained from the ADVISOR software [ADVISOR2]. One measure of fuel efficiency is brake specific fuel consumption (BSFC), which represents the fuel mass needed to release one unit of energy. A BSFC map can be calculated from the fuel consumption map using
(51) 
where is the fuelmass consumption rate and and are the engine torque and rotational speed respectively. We will use the BSFC to evaluate the efficiency performance of the engine. The internal combustion engine power is represented by and is given by
(52) 
In the optimal control calculations a quadratic multivariate polynomial was fitted to the engine map using a Linear Least Squares algorithm to speed up the fuel consumption calculation. The polynomial captures the shape of the map quite well and has an average absolute error of 2.6% over the entire engine operating range.
IiiG Flywheel
A highspeed flywheel is included in the drivetrain to provide a highpower energy redeployment capability, which complements the lowpower highenergy storage capability of the battery. While batteries can store energy for a relatively long time due to their low inherent losses, flywheel storage systems suffer from high losses especially when running at high speeds. A basic flywheel storage system can be represented by a spinning inertia with kinetic energy
(53) 
where is the moment of inertia and is the flywheel’s angular velocity. The dominant losses in the flywheel come from the friction in the bearings. These losses are modelled using the empirical relationship given on pg. 147 of [Doucette2013]
(54) 
where is given in . The flywheel dynamics are described by
(55) 
There are also losses in the continuous variable transmission (CVT) that can be lumped into an efficiency factor .
IiiH Power Transmission
The power transmitted to the rear wheel can be modelled by the constraint
(56) 
which ensures that the power delivered to the back wheels never exceeds the combined power delivery capability of the internal combustion engine, the battery and the flywheel. If the rearwheel tire force is positive, the vehicle is being driven, and the sum of powers delivered by the engine, the flywheel CVT and electric motor will match the mechanical power delivered to the rear wheel. Under braking, rearwheel tire force is negative, and the mechanical power at the back wheels is used to charge the battery and/or the flywheel, or else it is dissipated as heat.
At , is discontinuous and hence must be approximated using a smooth function. The approximation used in this study is
(57) 
in which is a constant. As is increased, the approximation (57) approaches .
Iv Optimal Control
The minimum fuel problem can be formulated as an optimal control problem that is now described.
The system is described by a set of equations of the form
(58) 
in which the statevector is given by
(59) 
The associated differential equations are given by (21), (23), (28), (29), (30), (49) and (55), respectively. The control vector is given by
(60) 
The problem is also subject to the path constraints (38), (39), (56) and bounds on the states and controls. There is also a constraint on maximum engine torque available as shown in Fig. 6.
The minimumfuel performance index to be minimised is given by
(61) 
The arrival time constraint can be written as an integral constraint as below
(62) 
where is the arrival time.
In practice, however, to avoid jerky controls and singular arcs, we actually control and minimise
(63) 
in which is an appropriate weighting matrix. We also impose slewrate limits on controls by placing hard constraints on .
Iva Numerical Optimal Control
An optimal control problem formulation general enough for our purposes is of Lagrange form. The aim is to determine states , controls and static parameters which minimise a cost functional
(64) 
subject to the state dynamics,
(65) 
path constraints
(66) 
and boundary conditions
(67) 
The normalised optimisation interval can be transformed into the general interval using the affine transformation .
The pseudospectral numerical optimal control solver GPOPSII [Patterson2013], which is based on the LegendreGaussRadau (LGR) collocation scheme, was used to solve the minimum lap time problem in this paper. In this scheme the state is approximated using a Lagrange polynomial of order
(68) 
where
(69) 
The statederivative approximation is thus given by
(70) 
Collocating the state dynamics at LGR points gives
(71) 
where are the elements of the LGR differentiation matrix. Note that is a noncollocated point. Using (71) we can discretise (65) and essentially transform the state dynamics given by ordinary differential equations into algebraic constraints.
The optimal control problem can then be approximated by an NLP problem. The cost function of this NLP is obtained by approximating the cost functional (64) using LGR quadrature. The NLP problem can be described by the task of finding ’s, ’s and which minimise
(72) 
in which the ’s are the quadrature weights [Abramowitz1965], subject to the following constraints
(73) 
(74)  
(75) 
(76) 
For clarity of exposition, the description provided here is for a singleinterval pseudospectral method (global collocation). GPOPSII uses a mesh and polynomial degree refinement scheme [Darby2011] so that error reduction can be achieved in the presence of nonsmooth problem features. The extension to multiple segments is straightforward, with the only requirement being the need to enforce continuity between each mesh interval.
The transcribed NLP problem is typically large but sparse. The IPOPT[Biegler2008] software library (based on interior point methods) was used to solve the NLP problem. Automatic Differentiation was used to provide IPOPT with accurate and computationally efficient first and second order derivatives [Weinstein2014].
V Results
In order to illustrate the concepts described, the motorracing Circuit de SpaFrancorchamps will be used as an exemplar track. This track is approximately 7 km long with an elevation change of approximately 110 m. The path is restricted to the neighbourhood of the track centre line in order to make comparisons easier; this was achieved by constraining the state , which is the perpendicular distance from the car mass centre to the track centre line, to be ‘small’. The threedimensional track, as well as a twodimensional projection on to a ground plane are shown in Fig. 7.
The route that the car is constrained to take together with the corner distances are shown in Fig. 8.
This figure will be useful in analysing the results to be presented later. In all the simulations described here, the vehicle will start at rest from the startfinish line (SF) (in practice the vehicle will start from a low speed in order to keep its reciprocal well defined), and will complete the circuit such that the combustible fuel usage is minimised. The car will come to a standstill at the end of its journey. The vehicle parameters are summarised in Table III given in the Appendix.
The minimum fuel problem was solved for a journey time of 240 s on the twodimensional and threedimensional track descriptions. These comparitive calculations are used to quantify the effects of three dimensionality. The optimal speed profile of the full hybrid vehicle for both the two and three dimensional tracks, along with the track elevation changes, are depicted in Fig. 9.
One can see that the track is on a slight incline when the car starts its journey. This results in the car’s speed at 200 m from the SF line being slightly higher on the 2D track. After the hairpin bend at 400 m the track falls away and the predicted speed on the 3D track exceeds that of the 2D track model. The 3D track speed advantage is then ‘given back’ as the car enters the uphill section between 1100 m and 2400 m. At the start of the incline at 1100 m, one also sees an increase in the fuel consumption rate as shown in Fig. 10.
The vehicle’s fuel consumption rate is high initially in order for it to accelerate from rest. As the vehicle approaches the hairpin bend at 400 m, the throttle is eased off at approximately 150 m, with the brake applied at approximately 200 m (see also Fig. 11). The fulllap fuel usage for the 2D case was 321.2 g, while that in the 3D case was 320.5 g. Although this difference is small, Fig. 10 shows that the fuel consumption strategy is different in the two cases; more fuel is used ascending hills, while fuel is saved coming down them.
Fig. 11 shows the power management strategy of the vehicle. An inspection of this plot brings a number of issues to light: (i) the battery power delivered and absorbed is relatively small compared to the flywheel; (ii) the controller tries to utilise the stored flywheel energy as fast as possible in order to avoid energy dissipation resulting from the rotational losses; (iii) the battery stored energy is only used when that in the flywheel has been depleted; (iv) often times the power delivery from the flywheel exceeds that of the engine; (v) the battery is discharged at relatively few isolated points on the circuit, and only when the flywheel energy has been depleted.
The power losses associated with various dissipation mechanisms are shown in Fig. 12. Predominant are the aerodynamic drag power losses that are proportional to the cube of the forward speed; see (34). The rolling resistance loss is given by the product of forward speed and the rolling resistance force, which is also relatively large in this case. The flywheel losses are the sum of internal rotational losses (54) and the CVT losses. The peaks which appear in flywheel loss plot are substantially due to transmission losses. The positive gradients appearing on some of the peaks are associated with the increasing losses as the flywheel is charged. The negative gradients, on the other hand, are associated with decreasing losses as the flywheel is discharged. The low loss regions, such as the one which appears between 3900 m and 4300 m, are due to selfdischarging. The battery and electric motor losses come from the battery’s output resistance and motor losses.
It can be seen that the optimal controller tries to operate the engine as close as possible to the optimum efficiency line; the controller chooses the engine speed and torque so that maximum mechanical energy is extracted per unit of fuel mass burnt.
The results presented thus far can be expanded to include realworld influences such as speed limits, enforced stops, changes in the roadsurface conditions and wind gusting. One approach to the solution of these problems is to set them up as multiphase optimal control problems [Rao2010]. In this formalism each phase can contain new models and/or new model parameters, new inputs, and new constraints. New parameters might include downgraded tyre parameters that particularise degraded road surface conditions, new inputs might include windrelated disturbances, and new constraints might include such things as speed limits and traffic controls. Fig. 14 demonstrates how the speed and fuel consumption rate vary when a 35 m/s speedlimit is imposed. At an elapsed distance of 1000 m it is evident that the fuel consumption rate reduces (below the unrestricted speed case) on entry to the speedrestricted section of road, and then increases above the unrestricted speed case in order to make up for the time lost. Similar variations in fuel consumption can be observed on the 5300 m to 6200 m road section.
It will be shown that the flywheel can provide significant fuel savings, especially when the journey times are low, and when the vehicle is required to complete the route aggressively. In these cases there will be many braking regions that will regenerate energy back into the flywheel that can then be quickly redeployed to accelerate the vehicle. For longer journey times the vehicle can complete the circuit at low speed with little or no braking. These ideas are illustrated in Figs. 15 and 16, which show the power and energy stored in the flywheel, respectively, for a vehicle with engine and flywheel only. For a journey time of 230 s, the flywheel frequently reaches high levels of stored energy. In contrast, when the journey time is increased to 265 s, the braking regions become less frequent, and lighter, resulting in fewer opportunities to scavenge energy to recharge the flywheel. Flywheel selfdischarging power losses are evident as negative gradients on the flywheel energy peaks in Fig. 16.
In order to analyse the battery management strategy, the vehicle with ICE and battery is considered in Fig. 17 for a journey time of 250 s.
The battery state of charge at the start and end of the lap is constrained to be 60 %. The battery is discharged in the acceleration regions for power and the recharged in braking phases. It is evident that unlike the flywheel, the discharge does not happen in an ‘onoff’ fashion. This is because the power losses in the battery are proportional to square of the current. This results in the optimal controller spreading the power delivery over longer time periods in order to maintain high efficiency. In the short charging phases all the available braking power is absorbed. In the regeneration phases, the electric motor losses result in less power than the power available at the wheels being delivered to the battery. These efficiency losses are also evident in the power assist mode when part of battery power is lost on its way to the wheels, and explains why the electric motor power is higher in braking and lower in acceleration.
In a final study, the minimum fuel problem was solved for a number of arrival times for four different energy storage combinations; the results are summarised in Fig. 18.
As one would expect, at lower journey times, when vehicle has to be driven more aggressively, the fuel consumption increases. Furthermore, at higher speeds the aerodynamic drag losses increase. The vehicle with both the flywheel and battery offers maximum advantage at low journey times. However, as the journey time increases, the vehicle with engine and flywheel performs the best as it does not have to carry the additional 100 kg of battery load. With increasing journey times the fuel usage drops significantly and the benefit of using an axillary storage system also diminishes, as energy regeneration opportunities decreases. The minimum journey times possible for the vehicle with engine only, engine and battery only, engine and flywheel only, and engine, battery and flywheel were 230.09 s, 225.76 s, 219.67 s and 217.35 s respectively. This highlights the power boost capability of the flywheel in aggressive manoeuvrings.
Some computational details
In this study the optimal control solver was initialised with 100 mesh segments. The number of collocation points was allowed to vary between 4 to 10 per mesh segment. The error tolerance across all meshes was and IPOPT tolerance was set to . For the vehicle with both battery and flywheel, when arrival time was set to 240 s, the error tolerance was reached after 18 mesh adaptation iterations. The final mesh had a total of 367 mesh segments and 2861 collocation points. The whole problem took under 2 hours to solve on an 8 core 3.5 GHz computer. For the engineonly case the error tolerance was reached with a similar number of collocations and meshes. However, the total solution time was only 20 minutes as only 10 mesh adaptation iterations were required and each iteration was faster to complete. The solution time is sensitive to the problem setup. For example, increasing the size of in (57) slows down convergence as the problem becomes ‘less continuous’. Avoiding the use of look up tables in the cost function speeds up the algorithm significantly. In general, even though the direct pseudospectral method shows great robustness, careful attention is required in problem setup when fast solution speeds are desired.
Global Optimality
As with any gradient based optimisation method, the optimal solution obtained from the algorithm might be a local rather than the global minimum. It is therefore important to quantify the sensitivity of the optimal solution to the initial guess. In this study, the problem was initialised using a set of sensible constant values across the entire solution space. A series of cases were then considered to examine whether the solution obtained can be trusted to be a global optimum. In all the cases, the vehicle was equipped with a battery and flywheel and the arrival time was set to 240 s.
Firstly, the initial state and control guesses were chosen constant, and given in terms of the bounds by and , with a weighting between zero and one. The optimal fuel usage for these cases are shown in Table II as ‘Linear x%’. For the 10%, 20% and 80% cases, the solution did not converge (marked as DNC). However, for all the other cases the solutions were essentially identical with the differences being less than the mesh tolerance error of .
Initial Guess  Fuel (g)  Initial Guess  Fuel (g) 

Linear 10%  DNC  Linear 20%  DNC 
Linear 30%  320.53  Linear 40%  320.52 
Linear 50%  320.53  Linear 60%  320.59 
Linear 70%  320.52  Linear 80%  DNC 
Random 1  320.51  Random 2  320.59 
Random 3  320.59  Random 4  320.60 
Random 5  320.59  Random 5  320.59 
Random 6  320.59  Random 7  320.59 
Random 8  320.58  Random 9  320.57 
Battery min  320.56  Battery max  320.53 
In an alternative test, a random value of between 20% and 80% of the total bound range, for each state and control, was chosen as the initial guess. This experiment was repeated 10 times and the results are shown as ‘Random n’ in Table II. All the runs converged successfully to the same solution. Finally, the problem was initialised with the lowest charging current and , and then with highest discharge current and for the battery. The results are shown as ‘Battery min’ and ‘Battery max’. Identical solutions were obtained once more, demonstrating that the optimal control problem has a good ‘radius of convergence’ when solved using a direct pseudospectral method. To ensure that the same minimal cost was not obtained from different trajectories the cases with the highest cost difference were selected (Case Random 1 and Random 2). The state and controls for the two cases where then compared. The worst case difference was found to be in the Engine Torque, , as plotted in Fig. 19. As can be seen, the trajectories are nearly identical, lending confidence to the idea that the solution obtained is indeed globally optimal.
Vi Conclusions
We have presented a novel method of evaluating the fuelconsumption performance of hybrid vehicles with multiple secondary energy sources. Rather than utilising standardised driving cycles, we use a specific route and then drive the vehicle so as to we reach the destination within a given journey time, while controlling the vehicle in order to minimise the combustible fuel consumption. In this framework driving aggressiveness can be systematically controlled by changing the arrival time. Another thrust of this work was to quantify the effectiveness of flywheels and batteries in reducing fuel consumption. Finally, the minimum fuel usage strategy for the hybrid vehicle over a three dimensional route was evaluated and the effects of different combinations of auxiliary energy storage systems were studied.
In future work it might be interesting to consider the effect of route optimisation. Adding traffic information to the model will also make the simulations more realistic. Component sizing of the combustion engine and auxiliary storage system can be formulated easily in the optimal control problem framework by including static parameters.
Symbol  Description  Value 
Vehicle total mass  1400  
Vehicle only mass  1280  
Battery and motor mass  100  
Flywheel mass  20  
x moment of inertia  500  
y moment of inertia  1000  
z moment of inertia  1000  
Mass centre from front axle  1.35  
Mass centre from rear axle  1.35  
Centre of mass height  0.5  
Drag coefficient  0.3  
Vehicle frontal area  1.8  
Air density  1.2  
Rolling resistance coefficient  0.009  
Max battery energy capacity  5  
Max battery power  25  
Min battery voltage  240  
Max battery voltage  210  
Battery internal resistance  0.5  
Min state of charge  40 %  
Max state of charge  80 %  
Max flywheel power  60  
Electric motor efficiency  85%  
Max flywheel energy  400  
Flywheel spinning inertia  0.02  
CVT efficiency  85%  

Mehdi Imani Masouleh received the M.Eng. degree in electrical and electronic engineering with firstclass honours from the Imperial College London, London, U.K., in 2011. He is currently pursuing the D.Phil. degree in engineering science at the University of Oxford. His current research interests include optimal control, vehicle dynamics, nonlinear stability and hybrid vehicles. 
David J N Limebeer received a B.Sc.(Eng) degree from the University of the Witwatersrand in 1974, MSc(Eng) and PhD degrees from the University of Natal in 1977 and 1980, respectively, and the DSc (Eng) from the University of London in 1992. He was a postdoc researcher at the University of Cambridge between 1980 and 1984. He then joined the Electrical and Electronic Engineering Department at Imperial College as a lecturer. He was promoted to Reader in 1989, Professor in 1993, Head of the Control Group in 1996, and Head of Department 19992009. In 2009 he moved to Oxford as Professor of Control Engineering and Professorial Fellow at New College Oxford. His research interests include applied and theoretical problems in control systems and engineering dynamics. He is a Fellow of the IEEE (1992), a Fellow of the IET (1994), and a Fellow of the Royal Academy of Engineering (1997). 