Faster, Higher, and Greener: Vehicular Optimal Control

Vehicular optimal control problems have been studied extensively since the early part of the 20th century. Progress in solving these problems has been driven primarily by applications in space and atmospheric flight, including launch vehicles, Earth-based and interplanetary space orbital transfer, and high-performance supersonic aircraft. In all of these applications, the ability to solve increasingly complex optimal control problems has been made possible by advances in high-speed computing. The mathematical and computing techniques being developed are now so diverse, and the range of applications of optimal control so broad, that a comprehensive review of the entire scope of vehicular optimal control is an impossible task. Instead, we will attempt only to provide a flavor of the scope and variety of these problems.


1066-033X/15©2015Ieee
Digital Object Identifier 10.110910. /MCS.2014 Date of publication: 17 March 2015 ehicular optimal control problems have been studied extensively since the early part of the 20th century. Progress in solving these problems has been driven primarily by applications in space and atmospheric flight, including launch vehicles, Earth-based and interplanetary space orbital transfer, and high-performance supersonic aircraft. In all of these applications, the ability to solve increasingly com-plex optimal control problems has been made possible by advances in high-speed computing. The mathematical and computing techniques being developed are now so diverse, and the range of applications of optimal control so broad, that a comprehensive review of the entire scope of vehicular optimal control is an impossible task. Instead, we will attempt only to provide a flavor of the scope and variety of these problems.
The albatross uses a variant of dynamic soaring, which is a flying strategy that can be optimized by exploiting wind gradients in the atmospheric boundary layer [1]. Dynamic gives them an energy-free loitering and traveling capability that prolongs their endurance [2]. Energy-optimal trajectories for UAVs equipped with solar cells are considered in [3]. These vehicles may be used to travel between given positions within an allowed time or for loitering indefinitely in a target neighborhood. Another problem relates to the optimal level turn of a solar-powered UAV flying in the atmosphere [4]. Readers who do not have optimal control as one of their main interests may wish to study one of the iconic fuel-optimal problems in vehicular optimal control. The famous Goddard Problem, posed by Robert H. Goddard (1882Goddard ( -1945 [5], is reviewed in "The Goddard Rocket." The minimal fuel thrust programming for the terminal phase of a lunar soft-landing mission is studied in [6]. The optimal thrust program consists of a period of zero thrust (free-fall) followed by full thrust until touchdown, which is a bang-bang strategy. Controlled thrusters and moment gyroscopes are used to maintain and correct the orbit and attitude of the International Space Station. The moment gyroscopes, which are zero-propellant actuators, can be used to re-align the structure. Optimal guidance trajectories can be designed using optimal control [7], which also plays a key role in the design of low-energy, low-thrust transfers to the moon [8]. The result is trajectories that use less propellant than required for standard low-energy transfers. Optimal low-thrust transfers to geosynchronous and Molniya orbits are studied in [9]. Another iconic optimal control problem, which moved to center stage with the introduction of the jet engine, is the minimum-time to climb problem. As with many of these problems, it comes in a variety of flavors. The optimal usage of a particular propulsion system to arrive at a predetermined altitude and speed in minimum time is studied in "Minimum Time to Climb." The solution is not obvious! Hybrid vehicles are equipped with at least two energy sources and powertrains that combine at least two modes of propulsion. While hybrid powertrains increase the complexity of the drivetrain control problem, they also increase its flexibility and in some sense can outperform anything achievable with a single-source drivetrain [10], [11]. Optimal control is used in [12] to optimize the powertrain energy flow and racing line of a hybrid electric vehicle (HEV) around a closed race track. The utility of optimal control in the solution of minimum-fuel problems in HEVs is demonstrated in [13]. Fuel-optimal control strategies for HEVs are also studied in [14]. The optimal control of hybrid vehicles using direct tran-scription is studied in [15], where the problem is discretized and the resulting finite-dimensional optimization problem is solved using a nonlinear programming code.
New technologies, as well as political, environmental, and marketing pressures, are driving the rapid development of hybrid vehicles. As a result, pollution thresholds are becoming increasingly more stringent, and hydrocarbonbased fuels are becoming more expensive and less socially acceptable. At the same time, the road-vehicle customer base is demanding the maintenance of customary levels of vehicle performance. More recently, turbo-compounded internal combustion (IC) engines and thermal energy recovery systems (ERSs) have received attention. These systems facilitate the capture of thermal energy from the notoriously inefficient IC engine's exhaust gases. The introduction of these technologies inevitably brings new technological complications, with the deployment of kinetic and thermal ERSs increasing the complexity of the associated control problem. The optimal control of multiple energy-related resources in real driving conditions, when the future trajectory is unknown, is a difficult problem.
The computation of numerical solutions to optimal control problems in motor racing originated in the 1950s, when simple heuristic arguments were employed to estimate lap times and optimize setup parameters [16]. In a development of the motorcycle work in [17], an indirect method is used to examine minimum-time maneuvers for a singletrack car model through a U-turn [18]. An early example of a direct numerical method for time-optimal control being applied to road vehicles is given in [19], which considers a lane-change maneuver for a simple nonlinear single-mass car model, and the optimal control problem was solved using a gradient projection algorithm. Direct multiple shooting is used in a minimum lap-time setting in [20]. This approach was deemed a sensible compromise between the sensitivity problems associated with shooting methods and the much larger nonlinear programs (NLPs) associated with direct collocation techniques [21]. Full-lap optimal control calculations for the Barcelona track using a simple four-wheeled car model required between 28 and 60 h to converge on a Sun Sparc workstation [20]. Some related work is described in [22]. The method reported simultaneously finds the optimal racing line and the controls required to follow this line in minimum time.
Optimal control calculations for closed-circuit racing must ensure that the car remains within the track boundaries while its wheels remain in contact with the road.

IEEE CONTROL SYSTEMS MAGAZINE » APRIL 2015
The Goddard Rocket T he Goddard rocket problem and its variants determine the thrust program that maximizes the altitude achievable by a rocket that expends a fixed amount of propellant [5]. The atmospheric drag plays a significant role in the switching structure of the optimal thrust for this problem. It was shown in 1951 [S1] that, for the low-drag case, a simple bang-bang structure results, where the first subarc is flown at maximum thrust while the second subarc consists of a coast. When atmospheric drag is introduced, the optimal thrust program may, under certain conditions, contain a singular interval that consists of a variable thrust subarc. If a realistic transonic drag model is used (in which the coefficient of drag transitions quickly from a relatively low value in the subsonic region to a higher value in the supersonic region), the optimal thrust program may be bang-singular-bang-singular-bang [S2].
To convey some of the technical detail of this problem, we follow the treatment given in [34]. Specifically, suppose it is desired to determine the thrust program that maximizes the terminal altitude of a one-dimensional sounding rocket subject to the expenditure of a fixed amount of propellant. The equations of motion for this problem are where h is the altitude, v is the velocity, m is the mass, T is the thrust (the control), ( , ) D D v h = is the drag force, c is the exhaust velocity of the propellant, and g is the acceleration due to gravity. The thrust is assumed to be constrained to lie in the interval The objective is to maximize ( ) h tf such that the final mass is in which all the new quantities are constants. In this case, (S5) becomes the nonlinear In the case that (S5) is never satisfied, a bang-bang thrust program results and the vehicle expends all of its propellant at the maximum allowed rate and then coasts to its maximum altitude.
Otherwise, upon completion of the maximum thrust subarc, the singular arc characterized by (S5) is encountered and the feedback law (S6) is used to propel the vehicle until its propellant depleted.
On the third subarc the vehicle coasts to its maximum altitude.
The performance of a supersonic ballistic missile, based loosely on the German V-2 illustrated in Figure S1, is shown in Figure S2. The maximum thrust phase propels the vehicle to an altitude of almost 20 km at a supersonic speed. The altitude of the vehicle increased to almost 35 km along a reduced thrust singular subarc. The vehicle then coasts to its maximum altitude of over 80 km.

Figure S2
An optimized thrust profile for the V-2 ballistic missile. Part (a) shows the classical bang-singular-bang thrust program, (b) shows the missile rising to an apogee of approximately 80 km, while (c) shows the missile's speed reaching a peak during the singular subarc. The full-thrust subarc is in red, the singular subarc is in blue, and the coasting subarc is in green. The apparently unusable computation times associated with high-fidelity optimal control calculations motivated a search for faster optimization methods. One hypothesis was that the computation time could be reduced by supplying the racing line rather than computing it as part of the laptime optimization problem. Another was the use of quasisteady-state models, which are now widespread within the Formula One racing community. Both of these ideas were investigated in [23], and, not surprisingly, much shorter computation times were achieved. The reported difference between the predicted lap times using time-optimal and quasi-steady-state analysis for the Barcelona circuit was 2.19 s, which is a significant discrepancy.
Another approach to minimizing the burden associated with minimum-lap-time studies is the use of a linear-quadratic preview to follow a prescribed driving line at a fixed speed [24]. The proposed method makes use of multiple linear models and a gain scheduling scheme, and it is demonstrated on three track segments. Ideas based on linear approximations have been further developed in [25], where (linear time-varying) model predictive control rather than linear preview control is exploited. This approximation allows the problem of finding a racing line and speed profile to be defined as a sequence of convex optimization problems. Prior work by the authors on minimum-lap-time optimal control is reported in [26], where the optimization of the car setup is also considered. The optimal control of hybrid powertrains is considered in [27], with extensions to three-dimensional tracks given in [28] and [29].
To improve the "green" credentials of motor sport, and accelerate road-car-relevant powertrain developments, the Fédération Internationale de l'Automobile (FIA) set out, in 2006, a "green agenda" for Formula One racing, announcing its intention to turn a sport in which cars consume over 150 kg of fuel per race into a catalyst for the development of "green technologies" for road cars. Thermal-electric hybrid powertrains were introduced into Formula One in 2014 [30]. In the lead up to the 2008-2017 Formula One engine development freeze, it was estimated that an average of 4 ms per lap were gained for every million dollars spent on engine development [31]. The hope is that the introduction of hybrid thermal-kinetic ERSs will produce gains that are far more significant and cost-effective than further engine development alone. The 2014 Formula One technical regulations mandate a fuel mass flow limit on all competing engines as well as a total fuel consumption limit of 100 kg per race. The results of recent research show that contemporary lap times can be maintained using less powerful turbo-charged engines, which require only two-thirds of the fuel required by contemporary cars, when thermalelectric ERSs are simultaneously employed [27].
The focus of this article is on optimal control applied to high-performance road vehicles, together with recently developed orthogonal collocation methods for the numerical solution of continuous optimal control problems. These topics were chosen as focal points for two reasons. First, the solution of optimal control problems in road vehicles has lagged behind solutions to air and space vehicles, and only recently has this gap begun to close in exciting and innovative ways. Second, the recently developed orthogonal collocation methods, together with methods for solving large, sparse NLPs, are only now being widely accepted as a viable solution method for complex, constrained optimal control problems. High-performance road vehicle optimal control problems that were previously challenging are now becoming computationally tractable.

HYBRID POWERTRAIN
From the beginning of the 2014 Formula One season, the race cars will be equipped with hybrid thermal-kinetic ERSs, which will be governed by a complex set of rules. These include fuel consumption limits, fuel mass-flow rate restrictions, energy storage limits, and ERS power constraints. These changes represent the biggest technological upgrade to have occurred over the entire history of the sport.
The kinetic ERS (KERS), which is part of the ERS, comprises an electric machine that is coupled directly to the crank shaft of an IC engine. This machine can be used to recover kinetic energy under braking conditions or augment the engine's drive power when the car is accelerating or traveling at high speed. The KERS motor is called the motor-generator-unit-kinetic (MGU-K). Any energy recovered, either thermally or under braking, is stored in the energy store (ES), which is typically a Li-ion battery pack. The 2014 technical regulations also mandate the use of turbo-compounded engines [32]. The pressure-boosting system includes a turbine, a compressor, and a motor-generator on a common drive shaft. This motor-generator is referred to as the motor-generator-unit-heat (MGU-H); see Figure 1. The turbine is bypassed by a controllable waste gate. Closing the waste gate increases the power generated by the MGU-H but decreases the power generated by the engine. The key features of the thermal ERS are summarized in Figure 2, which shows a plot of power from the IC engine PIC (under maximum fueling) versus power supplied by the MGU-H .
Ph Figure 2 shows that when the waste gate is fully open, the MGU-H requires 60 kW to operate the compressor and the IC engine output power is boosted by 20 kW. When the waste gate is closed, 40 kW is recovered from the MGU-H and the IC engine output power drops to approximately 440 kW. It is clear that operating the engine with an open waste gate is very inefficient from an energy-usage point of view. In summary, the powertrain controls are the fuel mass-flow rate, the waste gate opening, and the MGU-K demand. It is assumed here that the engine full-speed rotational losses are a constant 40 kW.
The 2014 Formula One energy-recovery rules are more complex than anything seen previously, with the optimal control problem correspondingly intricate. The key control-related features of the technical regulations [30], which must be translated into control constraints, are as follows: 1) The car is restricted to a 100 kg of fuel per race.
2) The fuel mass flow rate must not exceed 100 kg/h.
3) The car's ES may be a maximum of 4 MJ. 4) The ES can accept up to 2 MJ per lap from the MGU-K. 5) The MGU-K can draw up to 4 MJ per lap from the ES. 6) Power flows to and from the MGU-K are restricted to !120 kW. 7) The power and energy flow to and from the MGU-H is unrestricted. The asymmetry in rules 4) and 5) encourages the use and development of thermal energy recovery, which is a relatively untried technology.
We will now explain how these various constraints can be set up in an optimal control setting. The power available at the rear wheels Pm is constrained so that The Pmax term in (1) represents the power generated by the IC engine under full fueling when the waste gate is closed; we will use kW.
The normalized fuel mass flow rate is [ , ].
The second term in (1) represents the power boost resulting from opening the waste gate; [ , ] W 0 1 g ! is the waste gate control with kW. P 20 Wg = The third term represents the engine's rotational losses and is set at a constant 40 kW, Pk is the power delivered to the MGU-K ( Pk is positive when the MGU-K is motoring), and Pm represents the mechanical power delivered to the rear wheels (positive values accelerate the car). The power generated/absorbed by the MGU-H is  Figure S3). Kaiser's key observation was that the slow dynamics of the aircraft could be captured in a single state variable, which was called the "gesa-  Finally, it is assumed that c is small and that ( , ) .
Implicit in the energy-state model is the assumption that kinetic and potential energy can be exchanged instantaneously, and is a function of the waste gate opening and the fuel mass flow rate; kW P 40 min h =-and kW. P 60 max h = To monitor the resource constraints associated with the engine fueling, and the ES and the MGU-K usage, four auxiliary state variables are introduced. The first of these states is used to monitor the fuel consumption and is found by integrating the fuel mass flow rate (as a function of elapsed distance) around the circuit It is often convenient to use the elapsed distance rather than the elapsed time as the independent variable. Elapsed distance s and time t are related by , ds vdt = where v is the vehicle's speed, which is assumed positive. The consumed fuel mass is constrained by / , F n 0 100 kg l # # in which nl is the number of laps in the race. As pointed out in regulation 2) above, the fuel mass flow rate , F o when normalized, is constrained to the interval [0, 1]. The second auxiliary state is described by This state is used to monitor the energy stored in the ES and is constrained by .
The third auxiliary state is described by This state is used to monitor the energy supplied to the MGU-K from the ES and is constrained by .
$ The first alternative corresponds to the case when both electrical machines are operating as motors. In this case, it is only the MGU-K energy usage that is monitored for compliance purposes. The second alternative corresponds to the with es representing the slow dynamics, while changes in h and c occur quickly (at constant es ). In reality, kinetic-potential energy trades can be made relatively quickly using zoom-climb and zoom-dive maneuvering. Figure S4 shows the Kaiser diagram for Airplane 1, as given in [S8], using a four-state, flight-mechanics model. The initial zero-altitude speed of the aircraft is 130 m/s, and it is required to climb to 20 km and Mach 1 in minimum time. In accordance with the principles described, the aircraft gathers speed (and specific energy) through an increase in kinetic energy at zero altitude until it reaches 285 m/s. The reason for this initial ma-neuver is explained by the excess specific thrust contours (that  case in which the MGU-K is motoring but the MGU-H is generating. In this case, the power generated by the MGU-H is offset against the MGU-K power requirements, with only the difference monitored for compliance purposes. The fourth auxiliary state is used to monitor energy flows from the MGU-K into the ES and is described by This state is constrained by The first alternative corresponds to the case when both electrical machines are operating as generators. In this case, it is only the MGU-K energy generation that is monitored. The second alternative corresponds to the case in which the MGU-K is generating but the MGU-H is motoring. In this case, the power required by the MGU-H is offset against the power generated by the MGU-K, with the difference monitored. Differentiable approximations to the Heaviside step function are used to enforce constraints (5) and (6). The inequality constraint (1) and the four auxiliary states (3)-(6) and their associated box constraints are used to model the power and energy constraints required by the 2014 Formula One rules [30].

VEHICLE MODEL
Optimal control calculations for closed-circuit racing must ensure that the car remains within the track boundaries while its wheels remain in contact with the road. To enforce the associated path constraints, it is necessary to model the circuit and monitor the position and orientation of the vehicle on this circuit. As explained in [26], track and car kinematics are modeled using elementary ideas from classical differential geometry. The orientation of the car is described in terms of the angle between its center line and the center line of the track. The lateral position of the vehicle on the track is described by the orthogonal distance between the geometric center of the vehicle and the spine (or center line) of the track. The dynamic model of the car follows from Lagrange's equations and will be based on the geometric quantities shown in Figure 3. The external forces acting on the vehicle come from the tires and the aerodynamic loads; both of these influences will be described briefly. Several more subtle influences must also be modeled relating to tire load transfer effects, the braking system, the differential system, and the powertrains. Each of these modeling components will be described in the remainder of this section. The car parameters are available in [26] and [27]. In the interest of clarity, the track will be assumed to be flat (the interested reader can find a detailed treatment of three-dimensional track modeling in [28], with the associated three-dimensional vehicle dynamics and optimal control dealt with in [29]).

Track Model
The track is modeled using a curvilinear coordinate system that follows the vehicle using the track center-line position as the curvilinear abscissa [17]. Referring to Figure 4, the location  of the vehicle's mass center is described in terms of the curvilinear abscissa ( ) s t and the vector ( ( )) . n s t The former quantity defines the distance traveled along the track spine, while the latter gives the position of the vehicle's center of mass in a direction perpendicular to the track spine tangent vector ( ( )) . t s t It is assumed that the traveled distance s t h is an increasing function of time, and that "time" and "distance" are alternate independent variables. The standard dot notation will be used to signify derivatives with respect to time. At any point , s the track's curvature is given by C (where the s-dependence is implied), and its radius of curvature is given by .
R The track spine tangent vector t will be described in terms of the trackorientation angle .
i with the track's left-and right-hand half widths given by Nl and , Nr respectively. The yaw angle of the vehicle is given by } and the angle between the vehicle and the track by ; p . } i p = + In this coordinate system, constraints on the track width are easily expressed in terms of constraints on the magnitude of .
n The right-and left-hand half widths need not be equal.
It follows by routine calculation [26] that , cos sin s n in which u and v are the longitudinal and lateral components of the velocity of the car. This equation describes the way in which the point s on the track spine is "dragged along" by the moving car. The rate of change of n is given by

Independent Variable
The distance traveled will be used as the independent variable. This (rather than time) has the advantage of maintaining an explicit connection with the track position, as well as reducing (by one) the number of problem state variables. Suppose where S f comes from (7) as follows: The quantity S f is the reciprocal of the component of the vehicle velocity in the track-tangent direction (on the spine at s). It follows from (8) that and from (9) that , ds where ~} = o is the vehicle yaw rate.

Car Model
Each tire produces longitudinal and lateral forces that are responsive to the tire's normal loads and slip; the tire forces together with the steer and yaw angle definitions are given in Figure 5. The tire's longitudinal slip is described by a longitudinal slip coefficient , l while the lateral slip is described by a slip angle a [33]. Following standard conventions we use where R is the wheel radius and w the wheel's spin velocity. The quantities uw and vw are the absolute speed components of the wheel center in a wheel-fixed coordinate system. The four tire-slip angles are given by    Figure 6(a) shows the load dependence of the longitudinal tire force under pure longitudinal slip conditions; larger longitudinal forces are available when the normal load is increased. The main reason for introducing aerodynamic down-force generating systems on Formula One cars is to exploit this effect. Figure 6(b) shows how the longitudinal force is compromised when the tire is also side slipping. As the side-slip angle increases, the longitudinal peak force reduces and moves toward higher slip values. The normal-load-dependent inverted witch's hat in Figure 7 gives a three-dimensional representation of the combined slip characteristics of the tire model at a fixed normal load of 2000 N. The tire's pure slip characteristics are obtained by taking vertical cuts through the lines 0 a = and . 0 l = A more detailed account of tire modeling can be found in [26].
Balancing forces in the longitudinal and lateral directions, while also balancing the yaw moments, gives  (18) in which Fax is the aerodynamic drag force. These equations can be expressed in terms of the independent variable s

Aerodynamic Loads
The external forces acting on the car come from the tires and from aerodynamic influences. The aerodynamic force is applied at the center of pressure, which is located in the plane of symmetry of the vehicle. The drag and lift forces are given by and . , respectively. The negative sign in (22) indicates that the drag acts in the negative body-fixed x-direction. The speeddependent drag and down-force coefficients and the speeddependent location of the aerodynamic center of pressure are shown in Figure 8. The interpolating functions for , CD , CL and the aerodynamic center of pressure are all quadratic functions of the speed of the vehicle. In more elaborate aerodynamic models, the drag and down force coefficients are functions of the car's ride height, pitch, and yaw and steering angles.

Load Transfer
To compute the time-varying tire loads normal to the ground plane, the forces acting on the car are balanced in the nz direction and balance moments around the body-fixed xb-and yb-axes (see Figure 3). Next, balancing the vertical forces gives in which the s F z $$ l are the vertical tire forces for each of its four wheels (the first subscript, r or , f refers to the rear or front of the vehicle, and the second subscript, r or , l refers to the right or left of the vehicle), g is the acceleration due to gravity, and Faz is the aerodynamic down force acting on the car. Balancing moments around the car's body-fixed xb-axis gives in which Fy is the lateral inertial force acting on the car's mass center; see (18). Balancing moments around the car's body-fixed yb -axis gives where Fx is the longitudinal inertial force acting on the car's mass center [see (17)], while Faz is the aerodynamic down force. Equations (24)-(26) are a set of linear equations in four unknowns. A unique solution for the tire loads can be obtained by adding a suspension-related roll balance relationship, in which the lateral load difference across the front axle is some fraction of the whole where [ , ].
The first column gives the steady-state tire loading due to the vehicle's mass and aerodynamic down force; the second column derives from the lateral load transfer; while the third column is due to the vehicle's acceleration, braking, and drag force and produces the longitudinal load transfer effect.

Brakes
We approximate equal brake caliper pressures with equal braking torques when neither wheel on an axle is locked. If a wheel locks up, the braking torque applied to the locked wheel may be lower than that applied to the rolling wheel. For the front wheels, this constraint is modeled as in which fr and fl are the angular velocities of the front-right and front-left wheel, respectively. If either road wheel "locks up," the corresponding angular velocity will be nonpositive and the braking torque constraint (29) becomes inactive.

Differential
The drive torque is delivered to the rear wheels through a limited-slip differential, which is modeled by  Speed (m/s) Figure 8 Car aerodynamic maps. The drag coefficient CD is the solid blue curve, which is a function of speed. The down-force coefficient , CL also a function of speed, is given by the dot-dash red curve. The speed-dependent aerodynamic center of pressure location is given by the dashed magenta curve in meters from the front axle. The + symbols represent measured data points. , in which lr and rr are the rear-wheel angular velocities, R is the wheel radius, and kd is a torsional damping coefficient. The special cases of an open-and a locked-differential correspond to k 0 d = and kd arbitrarily large, respectively. Limited slipping occurs between these extremes.

OPTIMAL CONTROL
Optimal control has been described extensively in a wide array of excellent books including [34]- [39], while IEEE Control Systems Magazine articles on optimal control are found in [40] and [41]. In addition, excellent survey articles on numerical methods for optimal control are found in [42] and [43].
Most often, optimal control problems are posed in a general Bolza form [44]. The objective is to determine the state ( ) , subject to the dynamic constraint the path constraint , The Bolza optimal control problem given in (31)   is the control Hamiltonian, o is the Lagrange multiplier associated with the boundary conditions, and U is the admissible control set. Equations (35) and (36) are called a Hamiltonian system because they come from differentiating the (control) Hamiltonian [35], [38], [47]. Equation (37) is known as Pontryagin's minimum principle [51] and is a classical result to determine the optimal control. The conditions on the initial and final costate given in (38) are called transversality conditions [34]- [38], [48]- [50], [52], while the conditions on the Lagrange multipliers of the path constraints given in (40) are called complementary slackness conditions [53]- [55]. The Hamiltonian system, together with its boundary conditions, transversality conditions, and complementary slackness conditions, is called a Hamiltonian boundary-value problem (HBVP) [35], [38], [56]. Any solution ( ( ), ( ), ( ), ( ), ) x u t t t t m n o is called an extremal solution and consists of the state, the costate, the controls, any Lagrange multipliers that satisfy the boundary conditions, and any interior-point constraints on the state and costate.

NUMERICAL METHODS USED IN OPTIMAL CONTROL
Most practical optimal control problems must be solved numerically, and these solution techniques are separated along two different lines. The first distinction relates to the actual problem being solved, while the second difference pertains to the type of numerical method used to solve the problem. The problem being solved is either the one that arises from the first-order necessary conditions or the optimal control problem itself. A method that attempts to find a solution to the first-order necessary conditions is called an indirect method, while a method that attempts to solve the optimal control problem itself is called a direct method. In an indirect approach, a multiple-point BVP is solved with the boundary conditions at the endpoints and/or interior points only partially known. In a direct method, an optimization problem is solved numerically by approximating the original optimal control problem with a finite-dimensional optimization problem. Once the solution approach has been chosen, the problem is solved using either explicit or implicit simulation. In explicit simulation, the dynamics are integrated using a time-marching approach, where the solution at a given time step is obtained from the solution of the dynamics at one or more previous time steps. This process of time marching is then repeated by making progressively better approximations to the unknown boundary conditions until the known boundary conditions are satisfied. In implicit simulation, the solution of the dynamics at each time step along with the boundary conditions is obtained simultaneously without any notion of time marching. Both methods were developed in the context of a general BVP ( ( ), ), ( ( ), , ( ), ) , where neither the initial conditions ( ) x t0 nor the terminal conditions ( ) x t f are, in general, completely known. As a result, (41) must be solved, while seeking simultaneously to determine any missing components of ( ) x t0 and/or ( ).
x t f

Explicit Simulation (Time-Marching)
In explicit simulation, typically performed by time-marching, the solution of the differential equation at a future time step is obtained using the solution at current and/or previous time steps. The most basic explicit simulation method is the shooting method [57]. In shooting, an initial guess is made of the unknown boundary conditions at one end of the interval. Using this guess, together with the known conditions at that endpoint, the differential equation in (41) (41) are evaluated and an assessment is made as to how to update the unknown boundary conditions. The process of integrating the dynamics and updating the unknown conditions is repeated until the boundary conditions on (41) are satisfied to within a specified accuracy tolerance. The standard shooting method can present numerical difficulties due to instabilities in the dynamics. To overcome this difficulty, the multiple-shooting method [58] was developed. In a multiple-shooting method, the time interval [ , ] t tf 0 is divided into K subintervals Despite being of higher dimension, the multiple-shooting method represents an improvement over the standard shooting method because the shorter integration intervals reduce the sensitivity to errors of the unknown terminal conditions.

Implicit Simulation (Collocation)
In implicit simulation, the BVP (41) This parameterization is used to integrate (41) from tk 1 -to one or more stage times ; , , , The integrals given in (43) are then replaced by quadrature approximations of the form where ; , , are the values of the function approximation ( ) X t at the stage points ; , , t j I 1 jk f = and ; , , , is the integration matrix associated with the particular integration (quadrature) rule that is used to integrate the dynamics in mesh interval .
Sk Equation (44) can then be rearranged in the form which are called defect conditions that must be satisfied at all stage times in each mesh interval. In addition, assuming that the solution is continuous, continuity conditions

INDIRECT METHODS
The first-order necessary conditions given in (35)-(40) constitute a two-point BVP. This BVP can be solved using either explicit or implicit simulation as described above. A typical explicit simulation method would be indirect shooting or indirect multiple shooting, while a typical simultaneous indirect method would be indirect collocation. In an indirect-shooting method [57], an initial guess is made for the unknown boundary conditions at one boundary. Using this guess, together with the known initial conditions, the Hamiltonian system in (35) and (36) is integrated to the other boundary (that is, either forward from t0 to t f or backward from t f to t0 ) using a time-marching method. Upon reaching , t f the terminal conditions obtained from the numerical integration are compared to the known terminal conditions given in (34) and (38). If the integrated terminal conditions differ from the known terminal conditions by more than a specified tolerance, the unknown initial conditions are adjusted, and the process is repeated until the difference between the integrated terminal conditions and the required terminal conditions is less than some specified threshold. While the simplicity of simple shooting is appealing, this technique suffers from significant numerical difficulties due to ill-conditioning, which is examined in more detail in "Uncertainty Principle." The shooting method poses particularly poor characteristics when the optimal control problem is hypersensitive [59]- [63], that is, when time interval of interest is long in comparison with the time scales of the Hamiltonian system in a neighborhood of the optimal solution. In an indirect collocation method, the state and costate can be parameterized using piecewise polynomials, as shown in (42). The collocation procedure leads to a large-scale root-finding problem ( ) F z 0, = where the vector of unknown coefficients z consists of the coefficients of the piecewise polynomial.

Uncertainty Principle
2n-dimensional phase space can be associated with the HBVP given in (35) and (36). If these equations are well behaved, each point in the phase space evolves along a unique path, and we can think of these paths as flows in phase space.
The new elemental volume in phase space is where J is the n n 2 2 # Jacobian matrix . As shown in Figure S5, the Hamiltonian phase flow ( , ) x m U is volume preserving but not shape preserving. Suppose a region of the phase space has its boundaries determined by particles representing initial conditions for the HBVP given in (35) and (36). As the positions of these particles evolve in time, so will the boundaries of the associated control volume. In accordance with Liouville's theorem, any evolution of the control volume in the phase space, with its moving boundaries as specified, will maintain a constant volume. Its shape, however, will, in general, distort, becoming dendritic so as to disperse over different parts of phase space. It is the formation of these dendrites that causes numerical stability problems in the solution of the HBVP. Unlike the quantum mechanical setting, the set of starting points for the HBVP is a collection of discrete predetermined points rather than a probability distribution. It is this critical difference between quantum mechanical systems and classical shooting-type schemes that allow us to avoid the consequences of Liouville's theorem if we are careful about our choice of solution techniques. Recognizing that volume must be preserved, in a classical system, we can "push" the starting points (classical particles) closer together, while forcing the spaces between the points to the periphery of the control volume. The obvious way to do this is to use an implicit integration scheme, rather than time matching, to solve the HBVP because this essentially converts the HBVP problem into a large rootsolving problem on an integration mesh that can be solved to high accuracy.

Figure S5
The conservation of volume in phase space.
This system of nonlinear equations is then solved using an appropriate root-finding technique.

DIRECT METHODS
The most basic direct explicit simulation method for solving optimal control problems is the direct-shooting method. The control is parameterized as are the parameters to be determined by optimization. The dynamics are then satisfied by integrating the differential equations using a time-marching integration method. Similarly, the cost function (31) is determined using a quadrature approximation that is consistent with the numerical integrator used to solve the differential equations. The NLP that arises from direct shooting then minimizes the cost subject to any path and interior-point constraints.
Direct shooting suffers from issues similar to those associated with indirect shooting in that instabilities arise due to the long time interval over which the time-marching method is applied. Thus, in a manner similar to that for indirect methods, an alternate direct explicit simulation method is the direct multiple-shooting method. In this case, the time interval [ , ] t tf 0 is divided into K subintervals. The aforementioned direct-shooting method is then used over each subinterval [ , ] t t k k 1 -with the values of the state at the beginning of each subinterval and the unknown coefficients in the control parameterization being unknowns in the optimization. To enforce continuity, the following conditions are enforced at the interface of each subinterval The continuity conditions (47) result in a vector root-finding problem, where it is desired to drive the values of the difference between ( ) ( ) x x t t k k --+ to zero. It is seen that the direct multiple-shooting method increases the size of the optimization problem because the values of the state at the beginning of each subinterval are parameters in the optimization. Despite the increased size of the problem due to these extra variables, the direct multiple-shooting method is an improvement over the standard direct-shooting method because the sensitivity to errors in the unknown initial conditions is reduced due to integration performed over significantly smaller time intervals.

DIRECT COLLOCATION METHODS
Over the past two decades, implicit simulation direct collocation methods have become the de facto standard for the numerical solution of optimal control problems. A direct collocation method is an implicit simulation approach where both the state and control are parameterized in terms of basis functions (generally piecewise polynomials) for implicit simulation. In the context of optimal control, fixedorder integration methods such as Runge-Kutta and Hermite-Simpson methods [21], [64]- [73] have been employed.
In direct collocation methods, the following NLP arises from the approximation of the optimal control problem: As a rule, this NLP is large, containing tens or even hundreds of thousands of variables and constraints. The key feature of these problems, which makes them tractable, is the fact that the NLP is sparse. A variety of NLP solvers, such as SNOPT [74], IPOPT, and [75], have been developed specifically for solving large, sparse NLPs. In a typical fixed-order collocation method (such as a Runge-Kutta method), accuracy is improved by increasing the number of mesh intervals. A limitation of a fixed-order method is that, for complex problems, the only way to achieve convergence to the true optimal solution is to increase significantly the size of the mesh.

ORTHOGONAL COLLOCATION METHODS
In recent years, a great deal of research has been done on orthogonal collocation methods for solving optimal control problems. In these methods, the optimization horizon is subdivided into several mesh intervals, with the state in each interval approximated by a linear combination of Lagrange polynomials. The collocation is performed at Gaussian quadrature collocation points. A Gaussian quadrature method converges exponentially when used to approximate the integral of a smooth function. In fixedorder algorithms, such as Runge-Kutta and Hermite-Simpson methods, accuracy can only be improved by adjusting the number and width of the mesh intervals. In an orthogonal collocation method, accuracy can also be improved by adjusting the degree of the polynomial approximations. In an orthogonal collocation method, by increasing the approximating polynomial degrees in segments where the solution is smooth, it is possible to achieve convergence to the optimal solution by exploiting the exponential convergence of the Gaussian quadrature. Mesh Orthogonal collocation methods have been used to solve optimal control problems since the early 1980s.
boundaries can also be added at strategic locations where smoothness in the solution is lost.
To see how a Gaussian quadrature method works in practice, we consider the most-studied orthogonal collocation method that employs collocation at Legendre-Gauss-Radau (LGR) points. In the LGR discretization, the time interval subject to the dynamic constraints ( ) ( ( ), ( ), ( , , )), ( , , ), the path constraints and the boundary conditions As the state must be continuous at each interior mesh point, it is required that the condition ( ) ( ), ) .
The multiple-interval form of the LGR collocation method [76]- [85] is then stated as follows. First, the state of the continuous-time Bolza optimal control problem is approximated in , Then, using the integral form of the dynamics given in (52) and discretizing the result with an LGR quadrature on the interval [ , ] .
The following NLP then arises from the LGR discretization. , where , ( , and the discretized boundary conditions, b( , , , ) , must also be recognized. The continuity in the state at the interior mesh points [ , , The results presented show that the racing performance achieved in 2013 with 150 kg of fuel can now be achieved with 100 kg of fuel when an optimized hybrid powertrain is used.
Computationally, the constraint (60) is eliminated from the problem by using the same variable for both X The NLP that arises from the LGR collocation method is then to minimize the cost function (56) subject to the algebraic constraints (57)- (59).

MESH REfINEMENT
A key aspect of efficiently generating accurate solutions to optimal control problems using collocation methods is the proper placement of the mesh segment boundaries and the collocation points within each mesh segment. The algorithmic determination of the mesh boundaries and collocation points is called mesh refinement. In a traditional fixedorder collocation method (trapezoid, Hermite-Simpson, or Runge-Kutta method), the degree of the polynomial in each interval is fixed and the mesh refinement determines the location of the mesh boundaries. Unlike a traditional fixed-order method, an orthogonal collocation method has the flexibility to change the degree of the approximating polynomial (that is, the order of the method). As seen from the prior discussion of the LGR method, the order of the method is changed by modifying the number of LGR points in any mesh interval. The freedom to change the degree of the approximating polynomial in a mesh interval, combined with the ability to modify the width and location of a mesh interval, are referred to, respectively, as the p and h parts of the mesh refinement.
Fixed-order ĥ h collocation methods have been used extensively because they are typically implicit simulation forms of fixed-order, time-marching methods, which were computationally efficient in the early days of computing. In the last decade, methods that employ only p collocation (that is, a method that employs a single mesh interval and only varies the degree of the approximating polynomial) showed promise for problems whose solutions are smooth and can be accurately approximated using relatively lowdegree (that is, lower than a degree 10 to 15) polynomials. The use of a pure , h or a pure p method has limitations. Achieving a desired accuracy tolerance may require an extremely fine mesh (in the case of an h method), or it may require the use of an unreasonably large degree polynomial approximation (in the case of a p method).
To reduce significantly the size of the finite-dimensional approximation, and thus improve computational efficiency of solving the NLP arising from orthogonal collocation methods, a new class of hp orthogonal collocation mesh refinement methods has been recently developed. As described above, in an hp mesh refinement scheme, the number of mesh intervals, the width of each mesh interval, and the degree of the approximating polynomial within each mesh interval are allowed to vary simultaneously. The hp approach has led to several new algorithms for mesh refinement [80]- [83]. This recent research has shown that an accurate solution to the continuous-time optimal control problem can be obtained in a much more efficient manner than by using a fixed-order method, while concurrently taking advantage of the key features of a p method in segments where the solution is smooth. Various hp mesh refinement methods have been implemented in the Matlab orthogonal collocation software GPOPS-II [87] that is used here.

Circuit de Spa-Francorchamps
An optimal control problem can be formulated to identify the track curvature (as a function of the elapsed distance) from noisy GPS measurement data [28]. In effect, this optimal control problem is an optimal parameter identifier for function-valued parameters (the track curvatures). Boundary constraints, which are part of the problem formulation, ensure the closure of the track. Figure 9 illustrates the corner positions on the Circuit de Spa-Francorchamps (Spa) in Belgium with the distances to each corner given in Table 1.   Figure 10 shows the geodesic curvature and elevation changes for Spa that were estimated using the algorithm described in [28]. The normal curvature and relative torsion are not shown but are "small" in comparison with the geodesic curvature.

Mesh Refinement
The elapsed distance along the track spine (rather than time) is used as the independent variable. This choice of independent variable has the advantage of maintaining a direct link with the track boundary constraints and the position of the car. Figure 11 illustrates the mesh refinement process used when solving the optimal control problem. Figure 11 shown as black dots. Figure 11(b) shows the refined mesh after 16 mesh-refinement steps. After refinement, the mesh intervals are of unequal length as are the orders of the interpolating polynomials. Figure 12(a) shows the speed profile and racing line of the optimally controlled car. The computed lap time was 1 min, 46.96 s, which compares favorably with the record time of 1 min, 47.26 s, which was set in 2009 in a car with a significantly more powerful engine. This preliminary result shows that the hybrid powertrain with its associated fuel usage restrictions is not compromising the car's performance. Also marked on this figure are the approximate locations of the track's 19 corners. Figure 12(b) shows the vehicle's racing line. To maintain an optimal racing line, a race car driver will usually take up a position toward the outside edge of a corner before entering it and then turn in following a path through the apex. After "hitting the apex" and upon exiting the corner, the driver will allow the vehicle to return to the outer edge of the track, thus allowing the car to follow a wide radius of curvature path through the bend. The racing line affects the speed into and out of every corner and onto the following straight. A late apex is used to maximize the acceleration onto the following straight. In contrast, an early apex maximizes the use of speed from the incoming straight. This technique is used when approaching a corner, for which the following straight is significantly shorter than that prior to the corner. Hitting the apex mid-corner is a good way to take a corner, while also maximizing the mid-corner minimum speed. Although this is a perennial discussion point for competition drivers, it is "just physics," and an optimal control algorithm knows how to best exploit the track's layout.

Optimal Racing
Solution times are influenced by many things, including initial conditions, scaling strategies, approximation methods for nonsmooth functions, mesh refinement strategies, gradient calculation techniques, and error tolerances. In our case, these calculations were performed in approximately 20 min each on a standard tower computer based on a 3.5-GHz Intel Core i7 processor using algorithmic derivatives and the NLP solver IPOPT. Figure 13 shows the optimal utilization of the vehicle's energy resources. The normalized speed (speed/100) is shown (in cyan) as a convenient reference. The first point to note (shown in green) is that the car utilizes fully its fuel quota of 2.273 kg (Spa is a 44-lap race). The MGU-K (shown in blue) generates its maximum quota of 2 MJ per lap, while the maximum quota of 4 MJ that can be transferred from the ES to the MGU-K (shown in magenta) is not fully exploited. This asymmetry in the MGU-K's motoring and generating quotas was introduced to force the car to generate energy with the MGU-H. The state of charge of the energy storage (shown in red) is constrained to be the same at the start and end of the lap so that the optimal lap can, in principle, be repeated if tire degradation and fuel usage are ignored. Interestingly, the energy storage capability of the ES is only partially utilized. If the plots are examined in detail, it is possible to see the MGU-K driving the car on the fast sections of the track, while it generates under braking on entry into most of the corners. It is not possible to identify easily the contributions being made by the MGU-H, but they will be reflected in discrepancies between the ES's state of charge and MGU-K usage. Figure 14 shows the power transfers occurring within the powertrain. The green trace shows the fuel usage that operates in a bang-bang manner between its maximum and minimum values. Bang-bang type control behavior is to be expected because the optimal control problem's performance index seeks to minimize the elapsed time and is not a function of the controls. In addition, the controls enter the system dynamics in a linear fashion. If the performance index and the system dynamics are linear in the controls, bang-bang control strategies are likely to result. The bang-bang principle follows immediately from the PMP [35, p. 246-247]. The red trace shows the bang-bang-type operation of the MGU-K, which operates as a generator on the entry to corners, and as a motor on the exit from corners and on fast sections of the track. The car behavior at individual corners can be identified by referring to the distance markers for the Spa circuit given in Table 1. The blue trace, which is the power delivered to the rear wheels, is also bang-bang in character, with a total power delivery of 560 kW possible. Due to the energy efficiencies associated with its use, the waste gate is closed throughout an optimal lap and so is not shown. Under racing conditions, the waste gate is predominantly a safety device that can be used to prevent compressor over speeding under fault conditions. During heavy braking, the mechanical power delivered to the rear wheels is negative, the MGU-K will be generating, and the fuel flow rate will be cut to zero. As shown in Figure 13, this strategy allows the ES to be recharged. During high-speed driving, such as the section between 400 and 2000 m, both the engine and the MGU-K are used at full capacity, which drains the ES. On the section between 2000 and 2200 m, the MGU-K is unused and the ES is recharged from the MGU-H.

CONCLUSION
Orthogonal collocation methods have been used to solve optimal control problems since the early 1980s. These methods can now solve complex, high-dimension problems by combining the economy of Gaussian quadratures with state-of-theart, sparse, large-scale, NLP algorithms such as IPOPT, SNOPT and KNITRO. While the majority of these applications have thus far related to aerospace problems, we have demonstrated that they are also well suited to ground vehicular optimal control problems that may include hybrid powertrains, ERSs, and complex operating constraints.
The results in this article focus on the control of the thermal-kinetic ERSs that have recently been introduced into Formula One racing. In this application, the emphasis is on maintaining established levels of racing performance while simultaneously significantly reducing fuel consumption. The results presented show that the racing performance achieved in 2013 with 150 kg of fuel can now be achieved with 100 kg of fuel when an optimized hybrid powertrain is used. If a modest performance degradation can be accommodated, the fuel consumption can be reduced even further. In racing applications, the powertrain control is essentially a classical bang-bang strategy applied to both the fueling and MGU-K power delivery. A strategy reminiscent of bang-bang control is also employed during braking, although regenerative braking will be limited by the vehicular kinetic energy available, the capacity of the tires, and the restrictions placed on MGU-K usage. The potential for redeploying these technologies in everyday commercial road cars is self-evident.
control, and the control of aeroelastic phenomena in large structures. He is also interested in a variety of problems in multibody mechanics and the dynamics of two-and fourwheeled road vehicles. He is a Fellow of the IEEE (1992), the IET (1994), the Royal Academy of Engineering (1997), and the City and Guilds of London Institute (2002). His consultancy activities have been in a variety of vehicle-related matters; product liability litigation; and patent disputes in optical recording systems, drilling equipment, and highspeed packing machines. He can be contacted at the Department of Engineering Science, University of Oxford, Parks Road, Oxford, OX1 3PJ, United Kingdom.
Anil V. Rao