Optimal energy management for hybrid electric aircraft

A convex formulation is proposed for optimal energy management in aircraft with hybrid propulsion systems consisting of gas turbine and electric motor components. By combining a point-mass aircraft dynamical model with models of electrical and mechanical powertrain losses, the fuel consumed over a planned future flight path is minimised subject to constraints on the battery, electric motor and gas turbine. The resulting optimisation problem is used to define a predictive energy management control law that takes into account the variation in aircraft mass during flight. A simulation study based on a representative 100-seat aircraft with a prototype parallel hybrid electric propulsion system is used to investigate the properties of the controller. We show that an optimisation-based control strategy can provide significant fuel savings over heuristic energy management strategies in this context.


INTRODUCTION
Aviation currently contributes around 2% of worldwide human-made CO2 emissions but the demand for air travel and transport is growing at a significant rate. The aviation industry is committed to realising this growth sustainably with a drastic reduction of CO2 emissions by 2050. One avenue identified to contribute to the required CO2 reduction is through hybridisation of aircraft propulsion systems. This refers to enabling technologies for boundary layer ingesting aircraft (Hall et al., 2017) as well as rotary/tilt wing aircraft configurations in the Urban Air Mobility markets (Thipphavong et al., 2018). Hybrid electric architectures require real-time dynamic power management in order to minimise CO2 output. This paper addresses an optimal energy management problem for a hybrid electric aircraft with a propulsion system consisting of a gas turbine and a battery-powered electric motor in a parallel configuration. Although we consider here a battery as a secondary energy source, the approach is equally applicable to other primary and secondary energy sources such as hydrogen powered reciprocating engines, fuel-cells and super-capacitors.
Any optimisation methodology for primary power management must satisfy the basic requirements of determinism, convergence in finite time and verifiability. We propose a solution based on model predictive control (MPC) employing convex optimisation. Predicted performance, expressed in terms of the fuel consumption over a given future flight path, is optimised subject to constraints on power flow and stored energy, and subject to the nonlinear aircraft dynamics, which include nonlinear losses in powertrain components. The proposed convex formulation of the optimisation problem is made suitable for a real time nonlinear 1 On part-time secondment from Rolls-Royce PLC.
MPC implementation by introducing several key simplifying assumptions on the characteristics of powertrain components. Specifically, the gas turbine and the electric motor are modelled via sets of convex quasi-static power maps, battery losses are modelled using a time-invariant equivalent circuit, and the available data on the future flight path is assumed sufficient to determine powertrain shaft speeds across the prediction horizon.
Supervisory control methodologies for energy management have been proposed in the context of hybrid electric ground vehicles (e.g. Sciarretta and Guzzella, 2007). Several approaches have been proposed for this problem, including methods based on indirect optimal control (Kim et al., 2011;Onori and Tribioli, 2015), Dynamic Programming (Lin et al., 2003) and MPC (Koot et al., 2005;East and Cannon, 2019). Optimal control of hybrid propulsion systems in aircraft is a new application area that poses a number of distinct challenges, perhaps the most significant of which are complex nonlinear flight dynamics and the effects of the time-varying aircraft weight due to the burning of fuel during flight. On the other hand, the future power demand is likely to be more reliably predictable in aircraft than in cars since a pre-planned flight path is generally available for aircraft, whereas the driving cycle is subject to greater uncertainty in route and traffic conditions (Di Cairano et al., 2014;Josevski et al., 2017). The contribution of this paper is to demonstrate that the optimal energy management problem for hybrid aircraft can be posed as a convex optimisation problem. To the authors' knowledge, this is the first attempt to address an important new application area of energy management.
The paper is organised as follows. Section 2 derives a continuous time hybrid electric aircraft model. This model is the basis of the discrete-time model and the MPC strategy that are proposed in Section 3. Section 4 shows that the minimisation of fuel consumption can be expressed as a convex problem. Section 5 describes simulation results and conclusions are drawn in Section 6.

MODELLING
We assume a parallel hybrid electric aircraft propulsion system in which a gas turbine producing power P gt (t) is combined with an electric motor with power output P em (t) (Fig. 1). The net power output of the propulsion system, P drv (t), is produced by combining these two power sources via the relation P drv (t) = P gt (t)+P em (t) (assuming 100% efficiency in drivetrain components). When the drive power is negative, which may occur for example while the aircraft is descending, it is assumed that the same powertrain could be used to generate electrical energy (i.e. it is capable of operating in a "windmilling" mode) in order to recharge the battery. In practice, a variable-pitch fan would be required, which increases complexity. The gas turbine and electric motor shaft rotation speeds are ω gt (t) and ω em (t) respectively.

Battery Electric Bus
Motor/ Generator Fuel Gas Turbine The electric motor is powered by a battery with state of charge (SOC) E(t) and rate of change of energy P b , and the state of charge dynamics are given bẏ The battery is modelled as an equivalent circuit with internal resistance R and open-circuit voltage U , so that where U and R are assumed constant (East and Cannon, 2018). The map relating the mechanical power output of the electric motor P em to electrical input power P c is P c = h(P em ). We assume that, for fixed ω em , h(P em , ω em ) is non-decreasing and differentiable with respect to P em and h(·) is determined empirically from electric motor loss map data as h(P em , ω em ) = κ 2 (ω em )P 2 em + κ 1 (ω em )P em + κ 0 (ω em ) for some functions κ 2 (·), κ 1 (·), κ 0 (·). with κ 2 (ω em ) ≥ 0 and κ 1 (ω em ) > 0 for all ω em in the operating range.
The aircraft motion is constrained by its dynamic equations. Assuming a point-mass model (Stevens et al., 2016) and referring to Figure 2, the equilibrium of forces yields Using the polar coordinates parametrisation (v,γ), the drive power is given as follows where v is the magnitude of the velocity vector, S is the wing area, ρ the density of air, g the acceleration due to gravity, γ the flight path angle, C D = C D (α) the drag coefficient, α the angle of attack. Similarly, projecting equation (1) along the lift vector The contribution of the thrust in the vertical direction being very small, the term T sin (α) can be neglected (it can be checked a posteriori that α is small).
The rate of change of the aircraft mass is given bẏ , ω gt (t)) whereṁ fuel is the rate of fuel consumption and f (P gt , ω gt ) is assumed to be convex, differentiable and non-decreasing with respect to P gt for fixed ω gt . We assume that f (·) can be determined empirically from fuel map data in the form f (P gt , ω gt ) = β 2 (ω gt )P 2 gt + β 1 (ω gt )P gt + β 0 (ω gt ), with β 2 (ω gt ) ≥ 0, β 1 (ω gt ) > 0 in the operating range of ω gt .
The problem at hand is to find the real-time optimal power split between the gas turbine and electric motor that minimises while satisfying constraints on the battery SOC, limits on power flows throughout the powertrain, and producing sufficient power to follow a prescribed flight path.

DISCRETE-TIME OPTIMAL CONTROL
This section describes a discrete-time model that enables the optimisation of the power split between the electric motor and the gas turbine over a given future flight path to be formulated as a finite-dimensional optimisation problem. For a fixed sampling interval δ, we consider a predictive control strategy that minimises, online at each sampling instant, the predicted fuel consumption over the remaining flight path. The minimisation is performed subject to the discrete-time dynamics of the aircraft mass and the battery SOC. The problem is also subject to bounds on the stored energy in the battery (to prevent deep discharging or overcharging), as well as limits on power flows that correspond to physical and safety constraints.
The optimal solution to the fuel minimisation problem at the kth sampling instant is computed using estimates of the battery SOC E(kδ) and the aircraft mass m(kδ). The control law at time kδ is defined by the first time step of this optimal solution. The notation {x 0 , x 1 , . . . x N −1 } is used for the sequence of future values of a variable x predicted at the kth discrete-time step, so that x i denotes the predicted value of x (k +i)δ . The horizon N is chosen so that N = T /δ − k, and hence N shrinks as k increases and kδ approaches T .
The discrete-time approximation of the objective (2) is (5) where the forward Euler approximation has been used for derivatives. Using the same approach, the discrete-time battery model is for The problem to be solved at each time step k is therefore: Here (E, E) are the bounds on SOC that are required for normal battery operation, (P gt , P gt ) and (P em , P em ) are the bounds on gas turbine power and electric motor power respectively, and (ω em , ω em ) and (ω gt , ω gt ) are the bounds on the gas turbine and electric motor shaft rotation speeds.

CONVEX FORMULATION
The optimisation problem in (12) is nonconvex, which makes a real-time implementation of an MPC algorithm that relies on its solution computationally intractable. In this section a convex formulation is proposed that is suitable for an online solution. We assume that the aircraft speed v i and flight path angle γ i are chosen externally by a suitable guidance algorithm for i = 0, . . . , N − 1.
A convex formulation of the drive power is derived by expressing the drag and lift coefficients, C D and C L , as functions of the angle of attack α and combining the equations that constrain the aircraft motion in the forward and vertical directions. Over a restricted domain and for given Reynolds and Mach numbers, the drag and lift coefficients can be expressed respectively as a quadratic nondecreasing function and a linear non-decreasing function (Abbott et al., 1945): Combining (9), (10), (13) and (14), the angle of attack can be eliminated from the expression for drive power, so that P drv,i can be expressed as a quadratic function of the aircraft mass, m i , as follows where Here the flight path angles γ i and speeds v i are assumed to be fixed and are not optimisation variables. Since η 2,i > 0 for all i, the drive power is a convex function of m i . Note that there is no guarantee that satisfying equation (15) enforces equations (9) and (10) individually. In practice, assuming that we have full control over the eliminated variable α (via the elevator and fans), both individual dynamical equations can be satisfied a posteriori. The inequality constraint on α also has to be checked a posteriori.
For the given parallel hybrid configuration, we assume for simplicity that the electric motor and gas turbine share a common shaft rotation speed which is equal to the speed of rotation of the fan, i.e. ω gt,i = ω em,i for all i. If the shaft speed is known at each discrete-time step of the prediction horizon, then the coefficients in (4) and (8) can be estimated from a set of polynomial approximations of h(·) and f (·) at a pre-determined set of speeds. This allows h(P em,i , ω em,i ) and f (P gt,i , ω gt,i ) in (7) and (12) to be replaced by time-varying functions of the gas turbine power and electric motor power alone: h(P em,i , ω em,i ) = h i (P em,i ) = κ 2,i P 2 em,i + κ 1,i P em,i + κ 0,i , f (P gt,i , ω gt,i ) = f i (P gt,i ) = β 2,i P 2 gt,i + β 1,i P gt,i + β 0,i (17) with κ 2,i ≥ 0, κ 1,i > 0 and β 2,i ≥ 0, β 1,i > 0 for all i.
In order to estimate the shaft speed ω gt,i = ω em,i , and hence determine the coefficients κ 2,i , κ 1,i , κ 0,i , β 2,i , β 1,i , β 0,i in (16) and (17), we use a pre-computed look-up table relating the drive power to rotational speed of the fan, for a given altitude, Mach number and air conditions (temperature and specific heat at constant pressure). This enables the shaft speed to be determined as a function of the fan power output at each discrete-time step along the flight path. Although P drv,i depends on the aircraft mass m i , which is itself an optimisation variable, a prior estimate of the required power output can be obtained from (15) assuming a constant mass m i = m 0 for all i. The simulation results described in Section 5 show that this assumption has a negligible effect on solution accuracy.
We define g i (·) in terms of h i (·) as is necessarily a convex, non-decreasing, one-toone function if the lower bound on P em,i is redefined as since this bound ensures that h i (·) is a one-to-one nondecreasing convex function of P em,i .
The dynamic constraints (5), (6) and the power balance (11) can be expressed using (15), (16) and (17) as These constraints are nonconvex due to their quadratic dependence on the optimisation variables P gt,i , P em,i and m i . To convexify these constraints, we first eliminate P em,i from (19) and (20) using P b,i = g i (P em,i ) and P em,i = g −1 i (P b,i ). Then (19) becomes linear, E i+1 = E i − P b,i δ. Moreover, under the assumptions on g i (·) (convex, nondecreasing and one-to-one), the inverse mapping g −1 i (·) is a concave, increasing function (e.g. East and Cannon, 2018). Note that g −1 (·) is given explicitly as if κ 2,i > 0, and by at any time steps i such that κ 2,i = 0. Therefore, by relaxing the equality constraints in (18) and (20) to inequalities, a pair of convex constraints: is obtained since g −1 i (·) is concave and f i (·) is convex. With these modifications, and noting that the objective in (12) is equivalent to minimising m 0 −m N , the optimisation to be solved to determine the optimal power profile at the kth time step can be expressed as the convex problem: min where P b,i = g i (P em,i ), P b,i = g i (P em,i ), and the constraints are imposed for i = 0, . . . , N − 1. The form of the objective in (23) ensures that any feasible solution that does not satisfy the constraints in (21) and (22) with equality is suboptimal. Thus the solutions of (23) and (12) are necessarily equal if (12) is feasible.

NUMERICAL RESULTS
This section uses the optimisation problem (23) to construct an energy management case study involving a representative hybrid-electric passenger aircraft. Solutions of (23) were computed using the general purpose convex optimisation solver CVX (Grant and Boyd, 2008). Since the minimisation in (23) is convex, convergence of the solver to a global optimum is ensured.

Simulation scenario
The parameters of the model used in simulations are collected in Table 1. These are based on publicly available data for the BAe 146 aircraft. The propulsion system is assumed to consist of four gas turbines and electric motors, each with the hybrid-parallel configuration shown in Fig. 1 For the purposes of this study it is assumed that velocity and height profiles are known a priori as a result of the fixed flight plan entered prior to take-off. We consider The electric loss map coefficients κ i,2 , κ i,1 , κ i,0 can be estimated in two steps from these profiles. First, the shaft rotation speed, ω i (= ω gt,i = ω em,i ), is interpolated from a precomputed look-up table relating measured shaft rotation speed, altitude and drive power at a given Mach number (Fig. 5). Then, the coefficients are interpolated from a precomputed record of losses in the electric motor as a function of rotation speed. This procedure requires drive power P drv to be approximated a priori, e.g. by assuming constant aircraft mass for the duration of the flight. This assumption is supported by Figure 4, which shows that the electric map coefficients are almost identical for the estimated drive power profile and for the actual drive power profile computed retrospectively. We also find that the κ 2,i coefficients are negligible for all i.
The gas turbine fuel map used in this study is approximately linear (β 2,i ≈ 0 ∀i) and furthermore the fuel consumption does not depend significantly on shaft rotation speed. Therefore the fuel map coefficients are given in Table 1 as constants (i.e. β 1,i = β 1 , β 0,i = β 0 ∀i). √ T in where T in = T 0 (h) + v 2 /2c p is the temperature at inlet of the fan, c p = 1000 JK −1 kg −1 is the specific heat of air at constant pressure and T 0 (h) is the temperature of air at altitude h.

Results
The mission is simulated with sampling interval δ = 10 s over a one-hour shrinking horizon by solving the optimisation problem (23) at each time step and implementing the first element of the optimum power split sequence as an MPC law. The closed loop energy management control strategy is shown in Figure 6, which gives the power split for a single coupled gas turbine and electric motor. Clearly the constraints on the gas turbine and electric motor power are respected. The evolution of the battery SOC and fuel consumption are shown in Figure 7. The upper plot illustrates that the constraints on SOC are respected and that the SOC reaches a minimum when the drive power becomes negative, as expected. The lower plot in Fig. 7 shows that, as expected, the rate of fuel consumption is greater during the initial climb phase when the gas turbine power output is high. The fuel consumption recorded for this simulation is F * = 1799 kg. In comparison, a fully gas turbine-powered flight with the same initial total aircraft weight has a fuel consumption of F gt = 2034 kg. We note however that this reduction is achieved at the expense of reduced available payload as a result of the weight of the electric components of the powertrain (battery storage and electric motors).
In order to evaluate the optimality of the power split solution, we compare it with the strategy of supplementing the gas turbine with the maximum electric motor power (P em ) until the battery is fully depleted, then switching to a sustaining mode in which only the gas turbine operates. In hybrid vehicles this is known as a Charge-Depleting-Charge-Sustaining (CDCS) strategy (Onori and Tribioli, 2015). Using this strategy the power split is as shown in Figure 8 and the fuel consumption is F CDCS = 1858kg.
To investigate the potential for windmilling (energy recovery when the net power demand is negative), the lower bound on electric motor power is set to P em = −2 MW, to allow transmission of power from the fan to the battery  results shown in Figures 9 and 10. The windmilling effect can be seen at the end of the flight and is characterised by negative electrical power and battery recharge.

Discussion
Intuitively, the optimal power split strategy might be expected to consume as much fuel as possible at the beginning of the flight so as to reduce the aircraft mass, and thus reduce the drive power needed during level flight and descent. However, the MPC strategy maintains an almost constant electric power over the whole flight (Fig. 6). This is explained by the relatively short flight duration and the characteristics of the aircraft model, as a result of which the change in total mass of the aircraft is relatively small (less than 5%). Despite this, the MPC strategy achieves a non-negligible reduction in fuel consumption (3.2%) over the CDCS strategy.
More radical optimal power split solutions are obtained if the change in aircraft mass during flight is more significant.
In particular, the MPC strategy allocates more electrical power at the end of the flight if the gas turbine fuel consumption is increased. For example, Figure 11 shows the power split solution for a situation in which the rate of fuel consumption is increased so that the change in aircraft total mass during flight is 15% (with all other simulation parameters unchanged). The fuel consumption for the CDCS strategy in this case is about 4% higher than that of the MPC strategy.
We next consider the case that the upper limit on the gas turbine power output is reduced to P gt = 3 MW. The MPC energy management strategy for this case is shown in Figure 12. Here the power demand is such that the gas turbine is at maximum power while the aircraft climbs. As a result, the electric motor is needed to meet the total power output requirement while the gas turbine power output is saturated. The fuel consumption for this scenario is increased slightly (by 0.1%) since the electrical power is mostly used at the beginning of the flight to compensate for the limit on the gas turbine power output.

CONCLUSIONS
This paper proposes a model predictive control law for energy management in hybrid-electric aircraft. The main contribution of the work is a convex formulation of the problem of minimising fuel consumption for a given future flight path. We provide a simulation study to illustrate the approach, and demonstrate that significant fuel savings can be achieved relative to heuristic strategies. The convexity of the formulation is crucial for computational tractability and is expected to be a basic requirement for verification by the aviation industry. Future work will consider the design of bespoke solvers. In particular, first order solution methods are expected to provide computational savings by exploiting the high degree of separability in the problem, while also being suitable for real-time implementation. The modelling approach described in this paper provides a framework for optimising system design, and future work will explore flight path optimisation and evaluate alternative hybrid propulsion configurations.