# Power Management Problem for Civil Aircraft under More Electric Environment

Authors: Yuanchao Yang; Zichen Gao
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

The civil aviation industry is moving toward the more electric aircraft (MEA) which is to use electrical power to meet the load demands on multiple aircraft subsystems which are conventionally driven by other power resources. Thus, there will be introduced a large amount of new electrical power demands which are safety-critical for aircraft’s flight and this may lead the challenge for a reliable and efficient power management problem (PMP): the balance between the aircraft power supply and demands while minimizing the operation costs. To cope with the PMP for civil aircraft under more electric environment, in this paper, we explicitly give a detailed and complete modeling of all power supply resources (fuel and battery) and safety-critical electrical loads and cast the PMP as a mixed-integer nonlinear programming problem; we develop a practical solution methodology for the application on the real civil MEA. The proposed formulation and solution algorithm can give an efficient power schedule result with the minimal fuel and battery operation cost through a smart codispatch between the gas turbine generator, storage devices, and all electrical loads of MEA. Numerical testing results based on one real civil aircraft case demonstrate the economic and operational effectiveness.


---

## Body

## 1. Introduction

In a conventionally civil aircraft, except that the commercial and avionics loads are provided by the electric power [1], the main aircraft’s power extraction, such as fuel pumping into engine combustion chamber, hydraulic pump drive for primary flight control actuation, compressor spinning for cabin environmental control, and warming of the leading-edge wings/elevators for ice protection, is all resorted to the nonelectric resources, which include hydraulic, mechanical, and pneumatic power [2]. These resources usually introduce lower energy efficiency and cause more unfriendly environmental impact than electric power and have thus brought severe burden on operation cost and environment issue to the civil aviation industry. For that matter, over the last few years, there has been tremendous progress in the efforts to move toward more electric aircraft (MEA) environment [3–5]: many aircraft subsystems that previously used nonelectric power have been fully or partially replaced with electrical systems. This revolution of employing more electricity leads to MEA having improved energy extraction efficiency and reliability and reduced weight, fuel consumption, and cost. The MEA integrated with the associated more electric techniques become the tendency for the next generation of civil aircraft.However, the integration of the more electric concept into the civil aircraft will introduce several challenges for the operation of its power system [3–6], especially the power management problem (PMP) which is an important scheduling problem before an aircraft executing a flight mission, in order to ensure there generates enough power in any instant in time to be equal to the consumed power within the whole flight duration [6]. The PMP with its dispatched power generation level and load shedding results is one of the critical prerequisites for a conventional aircraft making a reliable and safe flight [7]. As shown in [8], considering the total amount of the load demand for a conventional aircraft power system is relatively small and the load shape is almost “flat” since the commercial and avionics electric loads are usually fixed within the whole flight duration, the PMP of the conventional aircraft is just to select the rated power of each possible load and calculate the sum value as the expected power generation level for the gas turbine generators. The conventional PMP is designed to protect electrical generators from overloads through prioritizing the loads as essential or nonessential one and will shed the loads based on the priority list in case of the abrupt generator outage. However, with the development of the MEA, the conventional PMP cannot make use of all generation resources and controllable loads in fulfilling a reliable power balance [7]. One of the significant reasons is that the MEA will introduce a large amount of electrical power demands and these electric loads are usually highly variable accompanied with the aircraft’s different flight phases represented by the altitude, velocity, and attitude [6–8]. The majority of increased electric loads for the MEA are also closely corresponding to the flight safety subsystems, such as flight control, engine control, cabin pressure hold, and leading-edge wing/elevator ice protection. This will lead the significant challenge for the conventional power management scheme in order to fulfill the balance between the aircraft power supply and variable demands in a reliable and efficient way [9].Moreover, another salient characteristic for MEA is that the power system will integrate more storage devices such as lithium/nickel-cadmium battery pack with large amount of capacity [10], in contrast to the conventional aircraft power system where the storage resource capacity is relatively small and only used as the emergency backup power source for the flight safety computer loads. The reasons of the large-scale integration of storage devices into MEA power system come from the following [11]. (1) Since the conventional electric power and flight thrust power for an aircraft are both collected from a common Brayton thermodynamics cycle implemented within the aircraft engine [12], the electric power generation from the gas turbine generator on aircraft is highly correlated with the expected thrust value predetermined by the given flight phase. The gas turbine generator thus has not enough capability to cope with the increased variability of the electric loads for a MEA. The storage device is thus an efficient power injection/sink source to provide the necessary load following capability for PMP. (2) The storage devices convert the chemical energy directly into electric energy without the combustion in a hydrocarbon fuel-consuming engine. They are emission-free and thus cleaner for the environment. However, the large-scale integration of storage resources into the aircraft power system brings key challenge to the PMP where there needs to make an efficiently integrated scheduling between the gas turbine generator, storage devices, and all controllable electrical loads so as to realize the minimal power generation and operation costs. This will require the PMP for MEA to have an explicit representation of the storage devices’ charging/discharging behavior with its corresponding operation cost and also reflecting the physical relationship between the power generation by the gas turbine generator and its associated fuel cost.The previous publications for enhancing the power management scheme of civil MEA power system focused on the following two ways. The first class such as [7, 13, 14] presents a platform for modeling and architectural exploration of PMP. The main focus of these publications is on the selection of the generator capacity and the design of storage devices rather than optimal design of schedule policy between different generation resources and load demands. While some optimization techniques have been reported for the storage capacity and shedding policy of the noncritical loads [7, 13], the optimizing of the overall PMP for MEA has received scant attention in these literatures. The second class including [15, 16] tries to fill the above gap and provides a mixed integer linear programming formulation for the optimal power management in which the decision variables are the generator power generation level, load connection/disconnection, and battery charging schedules. The global objectives are to minimize the number of load shedding behaviors which follows some preset priorities and ensure the generators to operate within their optimal operation ranges. However, the key challenge of this kind of approaches is that they are only aimed at cutting and reconnecting loads regarding their priority so as to prevent overloads. Thus, it is usually limited to switchable loads and cannot deal sufficiently with continuously controllable loads which are more common under the MEA environment. Furthermore, in most cases, the safety-critical load demand of the MEA is not constant during a flight. It can depend on the flight phase and using the fixed priorities as done in previous publications thus have not enough flexibility to cope with the variable characteristics for different loads. The scheme of PMP mentioned above cannot thus be capable of optimizing the overall efficiency and enhancing the reliability of the power system for MEA. A more practical way is coming back to analyze in deep the electrical load behaviors for the different safety-critical subsystems of MEA and connect the relationship between the electrical load variation and the basic power requirement of these subsystems with the flight phase change through an explicit mathematical representation. To the best of our knowledge, few researches have been made to deal with the above challenges. Our recent work [17] has been made to do a power optimization for the aircraft environmental control system considering the variable characteristics of its associated thermal power loads; however, the work in [17] only focused on the power management of the environmental control system but did not shed light on the larger and more important power management problem, i.e., the complete aircraft power optimization including the engine, environmental control system, flight control system, and anti-icing system.In this paper, we address the above challenges and propose a completely novel formulation and efficient method for PMP of the civil MEA. The main contributions of the paper are that (1) rather than just dividing the loads by shedding and nonshedding policy, we present a detailed and complete analysis of all safety-critical electrical load behaviors with the flight phase change through an explicit mathematical formulation and provide the model representing the relationship between the load variation and the basic power consumption change of all safety-critical subsystems. (2) We formulate the PMP of the MEA power system as a nonlinear constrained optimization problem which encapsulates a smart codispatch between the gas turbine generator, storage devices, and all controllable electrical loads with an objective for a minimal fuel cost and battery operation cost. (3) We provide a solution algorithm to solve the above problem for the real application on a civil MEA power system. We firstly divide the mathematical formulation of PMP into a two-level (outer and inner subproblems) optimization problem and then obtain the efficient real power schedule results through the two-level Generalized Benders Decomposition (GBD) techniques. (4) Finally, the numerical testing results based on the real civil aircraft case demonstrate the effectiveness of the proposed method in terms of its operation efficiency.The remainder of this paper is organized as follows. In Section2, we provide the safety-critical load modeling and present the formulation of PMP for a civil MEA power system with the gas turbine generator, storage devices, and all controllable electrical loads as a mixed-integer nonlinear constrained programming problem. In Sections 3 and 4, we separately discuss the proposed solution algorithm and demonstrate the application of the methodology on one real civil aircraft power system case. We present a summary of our conclusions in Section 5.


## 2. The PMP Formulation for Civil MEA Power System

All of the safety-critical electrical loads for a MEA are fundamentally used for providing the critical power to the multiple aircraft subsystems which are designed to keep a safe, comfortable, and economically efficient flight. These electrical loads for a MEA are thus closely corresponding to the MEA’s flight behavior, which basically causes the associated subsystems’ power demands. The flight behavior is generally stated by a series of flight phases representing a route on which the aircraft will make a flight, and the flight phases are usually predetermined and loaded into the aircraft’s flight management computer before the aircraft executing this flight [18]. If we take a single-aisle and two-engine civil aircraft as an example (such as Boeing 737/Airbus 320 aircraft), the definition of the flight phases of this kind of aircraft can be given in Table 1. If we connect these flight phases sequentially, the flight behavior of the aircraft can be demonstrated by Figure 1.Table 1
Definition of the flight phases.
Phase no.
Flight phase
Definition
1
Initial climb
Periods from aircraft lift-up to 450 m with constant speedV0; the pitch angle is 4°.
2
Stage 1 for climb
Periods from 450 m to 3000 m with constant speedV1; the pitch angle is 9.5°.
3
Stage 2 for climb
Periods for constant flight altitude 3000 m with speed increased fromV1 to V2; the pitch angle is 0°.
4
Stage 3 for climb
Periods from 3000 m to cruise altitude with constant speedV2; the pitch angle is 7°~10°.
5
Stage 1 for cruise
Periods for keeping cruise altitude with speed increased fromV2 to cruise speed; the pitch angle is 0°.
6
Stage 2 for cruise
Periods for keeping cruise altitude with cruise speed; the pitch angle is 0°.
7
Stage 3 for cruise
Periods for keeping cruise altitude with speed decreased from cruise speed toV2; the pitch angle is 0°.
8
Stage 1 for descent
Periods from cruise altitude to 3000 m with speedV2; the pitch angle is -10°~-7°.
9
Stage 2 for descent
Periods for constant flight altitude 3000 m with speed decreased fromV2 to V1; the pitch angle is 0°.
10
Stage 3 for descent
Periods from 3000 m to 1800 m with constant speedV1; the pitch angle is -9°.
11
Stage 4 for descent
Periods for constant flight altitude 1800 m with speed decreased fromV1 to V3; the pitch angle is 0°.
12
Stage 5 for descent
Periods from 1800 m to 450 m with constant speedV3; the pitch angle is -4°.
13
Approach and landing
Periods from 450 m to touchdown with constant speedV4; the pitch angle is -1.5°.Figure 1
Flight behavior represented by flight phases for a typical civil aircraft.We define the average flight duration from the lift-up to touchdown for a civil aircraft isS minutes and the smallest scheduling snapshot for the PMP is one minute, then the PMP will have S scheduling periods. Before addressing the PMP formulation, we make the following assumptions.Assumption 1.
Within a given time period, the aircraft is assumed to fly at a fixed pitch angle and zero roll and yaw angles as shown in Table1. In other words, there has only one kind of flight phases within a period and the horizontal turn maneuvers or the pitch attitude changes are assumed to be instantaneous and occur only at points which connect the consecutive two periods. The aircraft is thus modeled to fly along the phases represented by the piece-wise linear segments, which are constructed by sequences of points (Figure 1). The above assumption comes from the previous publications [19, 20] and is reasonable without loss of generality by the following facts: (1) a typical civil aircraft trajectory is usually to keep straight flight with constant pitch angle during the most of flight time and only make a horizontal/vertical maneuver at the predetermined waypoints and (2) the time for an aircraft’s maneuver is usually very short, such as 5-20 seconds. Considering the smallest scheduling snapshot for a PMP is one minute, it is reasonable to assume that the aircraft’s maneuver is instantaneous, and in this paper, we thus do not consider the impact by the aircraft’s maneuver on the electrical load modeling in the PMP formulation.Assumption 2.
In this paper, we only consider the two-engine type of civil aircraft because it is currently the main civil aircraft type [21]. However, the main modeling and results of this paper can also apply to the other types of civil aircraft. Considering that the two-engine aircraft always operates its engines with the same technical parameters under a same work mode except for the emergency condition, in the rest of this paper, we assume the two engines are “integrated” into a single engine with double-sized capacity so as to help make the statement of this paper more clear and concise.The altitude and speed values where the aircraft locates at the initial and end time in periods, s=1,2,⋯,S, are separately represented as hinits and Vinits and hends and Vends, while the pitch angle within this period is addressed by θs. We define that hs and Vs are separately the average flight altitude and speed for the aircraft within the period s and are equal to hinits+hends/2 and Vinits+Vends/2. Moreover, within the period s, the air’s average static pressure p0s and average density ρ0s at altitude hs can be obtained based on the transformation relationship between these values and hs [22]. The air’s average static temperature T0s at altitude hs can thus be determined by the ideal gas equation of state: p0s=R⋅ρ0s⋅T0s where R is the gas constant of the air. Thus, if the flight phases for all time periods 1,2,⋯,S are previously given, the parameters hs,Vs,θs, s=1,2,⋯,S, can completely describe the whole flight behavior of the aircraft. In this section, we describe the modeling of the safety-critical electrical loads explicitly represented by these parameters.


### 2.1. Safety-Critical Electrical Load Modeling of PMP for Civil MEA Power System

The conventional electrical load modeling method for PMP, which is mostly based on the load shedding idea that detects overloading and disconnects/connects relevant loads according to their priority, can only cope with the “flat” commercial loads (i.e., the electric power requirements without impact on the flight safety; they usually contain cabin illumination system, entertainment system, and galley refrigerators) and avionics loads (i.e., the electric power for the operation of the avionics-related computers located on the equipment cabinet in the cockpit) but cannot efficiently reflect the increased safety-critical electrical loads which are highly variable associated with the aircraft’s flight phases represented by the change of the average altitude, velocity, and attitudehs,Vs,θs under the MEA environment. The increased electrical load behaviors for the MEA cannot be simply modeled through the rule-of-thumb shedding policy but should be explicitly represented with the basic power requirement of their associated safety-critical subsystems. The main increased safety-critical electrical loads are usually the following four kinds:
(1)
The first kind of the increased electrical loads for the MEA is used for continuously pumping the fuel into the engine combustion chamber so as to keep the safe and stable operation for the engine [12, 23]. The pump is electrically driven by a direct current (DC) motor of which the load demand is equal to the value of the fuel volume flow rate injecting into the chamber multiply with the chamber inner pressure. The fuel volume (or mass) flow rate and the inner pressure of chamber can be determined through the stable operation equations of the engine [23] based on the required generator power output and the expected engine thrust given by its associated hs,Vs,θs. Although the current civil aircraft have multiple engine types including turbofan and turboprop, the general operation characteristics of the engines always obey the standard Brayton thermodynamic cycle and can all be simplified by a single compressor-turbine axis work mode [12]. The general operation flowchart of an aircraft engine and its required electrical load for fuel pumping is addressed in Figure 2. Its associated Brayton cycle is depicted in Figure 3Figure 2
The general operation flowchart of an aircraft engine.Figure 3
The Brayton cycle for an aircraft engineFrom [23], if the values hs,Vs,θs are given, the Mach number Mas is Vs/γi⋅gc⋅R⋅T0s; at the altitude hs, the air’s average total pressure is p1s=p0s⋅1+Mas2γi−1/2γi/γi−1 and the average total temperature is T1s=T0s⋅1+Mas2γi−1/2. The air mass flow rate into the engine inlet in period s ṁ0s is
(1)γigcR⋅Mas⋅1+Mas2γi−12−γi+1/2γi−1⋅p1sT1s⋅A1s,where Ais is the cross-sectional area of the engine inlet in period s which is taken as a decision variable controlled by the Full Authority Digital Engine Control (FADEC) computer [23]. Based on the aircraft inlet theory [12], we have p2s=πi⋅p1s and T2s=T1s. From Figures 2 and 3, the stable operation equations of the aircraft engine can be represented by equations (2), (3), (4), (5), (6), (7), and (8):
(2)p3sp2s=1+ηcT3sT2s−1γc/γc−1,s=1,2,⋯,S.Equation (2) represents the isentropic change characteristic of the air state variables within the engine compressor. ηc<1 is the fixed operation efficiency parameter of the engine compressor.
(3)p4s=p3s,s=1,2,⋯,S,(4)fs=cpt⋅T4s−cpc⋅T3sηb⋅hpr−cpt⋅T4s,s=1,2,⋯,S.Equations (3) and (4) describe that there has an isobaric characteristic of the air thermodynamic process in the combustion chamber where the pressure stays constant; furthermore, the total energy taken by the exported mixed gas from the combustion chamber is equal to the total energy of the imported air plus the combustion energy released by the fuel.
(5)T5sT4s=1−ηt1−p5sp4sγt−1/γt,s=1,2,⋯,S.Equation (5) represents the isentropic change characteristic of the air state variables within the engine turbine. ηt<1 is the fixed operation efficiency parameter of the engine turbine.
(6)T3s−T2s=Kc⋅ns2,s=1,2,⋯,S.Equation (6) states the relationships between the rotation speed and the change of total temperature for the engine compressor. Kc is the associated transfer coefficient for the engine compressor.
(7)ṁ0sηm1+fscptT4s−T5s=ṁ0scpcT3s−T2s+Pgs,s=1,2,⋯,S.Since the compressor and turbine are mechanically coaxial, (7) describes the reduced amount of the total enthalpy for the turbine in period s that is equal to the energy output of the generator within this period plus the increased amount of the total enthalpy for the compressor in period s. ηm<1 is the mechanical efficiency.Moreover, if we define the required thrust in periods to be Fs, its value can be determined through the flight mechanics analysis as shown in Appendix A; thus, Fs is represented by the parameters hs,Vs,θs. Then, based on aircraft engine theory [23], we can have
(8)Fs=ṁ0sgc1+fs2cptgcT5s1−p0sp5sγt−1/γt−Vs,s=1,2,⋯,S.Finally, since all decision variables in (2), (3), (4), (5), (6), (7), and (8) are positive and have their associated upper bounds (especially, Ais is also required to be no less than a minimum cross-sectional area A¯i, so as to keep the enough air mass flow rate into the engine for its stable operation), we also have
(9)0≤p3s≤p¯3,0≤T3s≤T¯3,0≤p4s≤p¯4,0≤T4s≤T¯4,A¯i≤Ais≤A¯i,0≤p5s≤p¯5,0≤T5s≤T¯5,0≤ns≤n¯,s=1,2,⋯,S.From the above equations, we can find the load demand of fuel pumping for the engine operation in periods P1s is equal to ṁ0s/ρf⋅fs⋅p3s where ṁ0s and p3s are constrained by (2), (3), (4), (5), (6), (7), (8), and (9).
(2)
The second kind of the increased electrical loads for MEA is the electric drive of the compressor spinning for aircraft environmental control system (ECS) in order to control the temperature and pressure of the aircraft pressure vessel which includes the cockpit, passenger cabin, and interior cargo compartments [24–26]. The ECS contains two main functions: air conditioning and vessel pressurization control. The air conditioning function conditions the fresh air from the outside of the aircraft and supplies it to the pressure vessel at the requested temperature. Here, to note that the air conditioning function in this paper only refers to the regulation of temperature, we do not consider the removal of humidity by the air conditioning since the function of controlling the humidity in the pressure vessel is usually implemented only when the aircraft is on ground [26]. Meanwhile, the pressurization control function regulates the pressure within the vessel by controlling the outflow of air supplied by the air conditioning function and also via an electrically controlled valve which controls the pressure level through the rate of air change in the vessel so as to provide satisfactory pressure values for comfort and safety for all the passengers and crew on the civil aircraft. A typical operation flowchart of an aircraft electrical ECS [9] including its associated electrical load for spinning the compressor is addressed in Figure 4Figure 4
The operation flowchart of an aircraft electrical ECS.In Figure4, at first, since the initial total temperature of the imported air is relatively low, it is usually preheated where the energy is provided by a DC motor. After that, ECS compressor 1 is electrically driven by the same motor so as to realize a higher pressure and temperature for the imported air. Through releasing part of the energy while keeping the pressure unchanged by heat exchanger 1, the imported air is then conditioned by a typical Brayton inverse thermodynamic cycle so as to reach an appropriate pressure and low temperature which can be injected into the vessel [27]. Finally, an electrical control valve at position A will control the air mass flow rate exhausted outside the vessel so as to manage the pressure of vessel to the requested level; the imported air with low temperature will balance the temperature of vessel at a comfortable level.From the previous deductions,Mas, p1s, and T1s can be determined if the values hs,Vs,θs are given for period s. Moreover, since the cross-sectional area of the inlet for ECS compressor 1 A0e is usually taken as a fixed value, the mass flow rate of the imported air in period s ṁ0es is equal to
(10)γigcR⋅Mas⋅1+Mas2γi−12−γi+1/2γi−1⋅p1sT1s⋅A0e.In this paper, we assume the pressure ratio driven by the inlet for ECS compressor 1 is 1; then, the stable operation of the ECS can be represented by equations (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), and (24) [26]:
(11)p2esp1s=1+ηc1eT2esT1es−1γce/γce−1,s=1,2,⋯,S,(12)T2es−T1s=Kc1e⋅nc1es2,s=1,2,⋯,S.Equation (11) represents the isentropic change characteristic of the air state variables within ECS compressor 1. ηc1e<1 is the operation efficiency parameter of ECS compressor 1. Equation (12) states the relationships between the rotation speed and the change of total temperature for ECS compressor 1. Kc1e is defined to be the associated transfer coefficient for ECS compressor 1.Based on the heat transfer theory [27], there has an isobaric characteristic of the air thermodynamic process in heat exchanger 1 where the pressure stays constant; the total energy taken by the cold side air (the air imported through the inlet of the heat exchanger so as to provide the low-temperature air resource) is equal to the total energy released by the hot side air (the high-temperature air from ECS compressor 1); furthermore, the cold side air and hot side air reach a common temperature after the heat transfer within heat exchanger 1. The above characteristics are represented by
(13)p2es=p3es,s=1,2,⋯,S,(14)ṁ1es⋅T3es−T1s=ṁ0es⋅T2es−T3es,s=1,2,⋯,S.In (13) and (14), the mass flow rate of the cold side air into heat exchanger 1 in period s, ṁ1es, can be stated as
(15)γigcR⋅Mas⋅1+Mas2γi−12−γi+1/2γi−1⋅p1sT1s⋅A1e.where A1e is the cross-sectional area of the inlet at the cold side of heat transfer 1.
(16)p4esp3es=1+ηc2eT4esT3es−1γce/γce−1,s=1,2,⋯,S.Eq. (16) represents the isentropic change characteristic of the air state variables within the ECS compressor 2. ηc2e<1 is the operation efficiency parameter of ECS compressor 2.
(17)p5es=p4es,s=1,2,⋯,S,(18)ṁ1es⋅T5es−T3es=ṁ0es⋅T4es−T5es,s=1,2,⋯,S.Similar with equations (13) and (14), the isobaric characteristic of the air thermodynamic process in heat exchanger 2 is represented by equations (17) and (18). The cold side air flow into heat exchanger 2 in period s is the exported air flow at the cold side of heat exchanger 1.
(19)T6esT5es=1−ηte1−p6esp5esγte−1/γte,s=1,2,⋯,S,Equation (19) represents the isentropic change characteristic of the air state variables within the ECS turbine. ηte<1 is the operation efficiency parameter of the ECS turbine.
(20)T4es−T3es=Kc2e⋅nc2es2,s=1,2,⋯,S.Equation (20) states the relationships between the rotation speed and the change of total temperature for ECS compressor 2. Kc2e is defined to be the associated transfer coefficient for ECS compressor 2.
(21)T4es−T3es=ηme⋅T5es−T6es,s=1,2,⋯,S.Since compressor 2 and the turbine are mechanically coaxial, (21) describes that the reduced amount of the total enthalpy in period s for the ECS turbine is equal to the increased amount of the total enthalpy for ECS compressor 2 in period s. ηme<1 is the associated mechanical efficiency.
(22)mvs60sec+ṁ0es−ṁouts⋅Tvs+1=Qv1s+Qv2s+Qv3s+Qv4scpe⋅60sec+ṁ0es⋅T6es−ṁouts⋅Tvs+mvs60sec⋅Tvs,s=1,2,⋯,S.Equation (22) describes the temperature control characteristics in the vessel where the enthalpy change amount of the air in vessel at period s is equal to the sum of the imported enthalpy by the upstream air, the exported enthalpy by the exhausted air from valve A, the enthalpy introduced by the passenger/pilot activity, the heat transfer between the aircraft’s skin and the air, the radiation by the sun, and avionics device operation.
(23)p6es=mvs⋅R⋅TvsDv,s=1,2,⋯,S.Equation (23) represents the pressure balance constraints in vessel where Dv is the fixed volume of the vessel.Finally, all decision variables in (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), and (23) are all positive values and have their associated upper bounds; besides that, during the whole flight, the pressure and temperature within the vessel are separately required to be larger than the value maxp~,p1s (p~ is the pressure value over the sea level by 1800 m) and the minimum comfortable temperature T¯v. Thus, we have
(24)T1s≤T1es≤T¯1e,0≤T2es≤T¯2e,0≤p2es≤p¯2e,0≤T3es≤T¯3e,0≤p3es≤p¯3e,0≤T4es≤T¯4e,0≤p4es≤p¯4e,0≤T5es≤T¯5e,0≤p5es≤p¯5e,0≤T6es≤T¯6e,maxp~,p1s≤p6es≤p¯6e,0≤nc1es≤n¯c1e,0≤nc2es≤n¯c2e,T¯v≤Tvs≤T¯v,0≤mvs≤m¯v,0≤ṁouts≤ṁ¯out,s=1,2,⋯,S.Based on the above equations, we can find that the load demand of the electrical DC motor for the ECS in periodsP2s is equal to ṁ0es⋅cpe⋅T2es−T1s. T2es for s=1,2,⋯,S are constrained by equations (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), and (24).
(3)
The third kind of the increased loads for the MEA is the electrical drive of the hydraulic pump for the primary flight control actuators, such as the electro-hydrostatic actuators for the surfaces (elevators, ailerons, and rudders) [28, 29]. In the electrically driven flight control actuation system, a DC motor energizes the pump operation of an actuator in order to maintain the required hydraulic oil mass flow rate and pressure to the actuation system. In this paper, since we have assumed the aircraft are to fly at fixed pitch angle and zero roll and yaw angles within a given time period, we will not consider the electrical loads introduced by actuating the elevators for a vertical maneuver and actuating the ailerons and rudders for a horizontal maneuver. From the flight control actuation theory [29], the only required power within a period so as to maintain the expected flight behaviors is to keep the elevators actuating with fixed deflection angle and it can thus be calculated according to the specified flight phase which is represented by hs,Vs,θs. We define the required air-dynamical force on elevators in period s to be Les, and its value can thus be determined through the flight mechanics analysis shown in Appendix B; in other words, Les can be represented by the parameters hs,Vs,θs. Finally, we have known that the electrical load required by the flight control system in period s P3s is Kleak⋅Les/Ae2/ρhy where Kleak is the leak coefficient addressing the relationship between the air-dynamical force on elevators and required hydraulic oil mass flow rate [29]. ρhy is the density of the elevators’ hydraulic oil and Ae is the area of the elevators
(4)
The fourth kind of the increased loads for the MEA is the electrical warming of the leading-edge wings and elevators for ice protection [30, 31]. The icing, which is usually caused by the supercooled water droplets in clouds where the aircraft is passing by, has been recognized as a potentially dangerous phenomenon for aircraft’s flight, since the buildup of ice on the leading edges of wings/elevators can cause an irregular change of the aircraft’s maneuverability characteristics and may finally result in the aircraft’s loss of control [30]. Anti-icing refers to the prevention of any buildup of ice on a wing/elevator during flight and can be achieved by the use of thermal energy. Under the MEA environment, the electrical energy is used for anti-icing via electrothermal heaters where the electrothermal pads are bonded to the outer surface in need of anti-icing. The heat generated by these electrothermal pads is conducted to the wing’s/elevator’s interface where a thin film of water is formed. This breaks the adhesive bond between the ice and the outer surface of the pad so that the ice can be swept away by aerodynamic forces. The anti-icing electrical power is usually implemented by cycling the load between the different electrothermal heaters such as the electrothermal heaters on the wings and elevators. This will reduce the electrical power need while realizing the minimal anti-ice requirements. In this paper, we define electrical loads for anti-icing the leading-edge wings and elevators in period s P4s to be separately equal to Pws and Pes which are activated iteratively in every 3 minutes after the aircraft’s altitude becomes higher than the minimum altitude with icing condition. Pws and Pes are separately the rated power values of the electrothermal heaters on the wings and elevators, respectively


### 2.2. Smart PMP Formulation

The PMP formulation firstly needs to explicitly account for the objective function with a minimal fuel cost and battery operation cost. The decision-making of a PMP is to implement a codispatch between the gas turbine generator, storage devices, and all controllable electrical loads in an optimal fashion. An effective quantification of the system operation costs in the PMP formulation is stated in the following:(25)min∑s=1SCgs+Cbs.Equation (25) describes that the total operation cost for the whole scheduling horizon is the sum of the fuel cost by generator output, battery operation costs by battery charging and discharging numbers. In (25), Cgs is equal to ag⋅ṁ0s⋅fs where ag is the linear coefficient between the fuel mass flow rate and fuel cost; Cbs is equal to hb⋅zbcs+zbds where hb is the operation cost coefficient.We next state each of the constraints we include in the problem formulation.


#### 2.2.1. Battery Charging and Discharging Behavior Constraints

(26)zbcs+zbds≤1,zbds⋅P¯bo−zbcs⋅P¯bi≤Pbs,Pbs≤zbds⋅P¯bo−zbcs⋅P¯bi,s=1,2,⋯,S,(27)E¯b≤Ebs≤E¯b,Ebs+1=Ebs−Pbs⋅60sec,s=1,2,⋯,S.Equation (26) stipulates that in every time period, the state of the battery will be one of the three cases: charging (the discrete variable zbcs∈0,1 is 1), discharging (the discrete variable zbds∈0,1 is 1), and shut-off (zbcs and zbds are both 0); the charge/discharge power in period s Pbs is within its upper and lower bound. Equation (27) addresses the upper and lower bound constraints for state-of-charge (SOC) energy in every time period and its dynamic characteristic.


#### 2.2.2. Load Balance Constraints

(28)Pgs+Pbs=P1s+P2s+P3s+P4s+P5s,s=1,2,⋯,S.In (28), P5s is defined to be the fixed power load for commercial and avionics services in period s. Equation (28) stipulates that in every time period, the total power generated by generator and battery is equal to the total load demands by engine fuel pumping, ECS, flight control actuation, anti-ice for leading edge of wings/elevators, and commercial and avionics services.


#### 2.2.3. Safety-Critical Controllable Load Constraints

Based on Section2.1, all associated constraints with the safety-critical controllable loads are listed in equations (2), (3), (4), (5), (6), (7), (8), (9), (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), and (24).


#### 2.2.4. Generator Power Output Constraints

(29)0≤Pgs≤P¯g,s=1,2,⋯,S.Equation (29) addresses that the power output by the generator in period s should be positive and less than the rated power capacity P¯g.


## 2.1. Safety-Critical Electrical Load Modeling of PMP for Civil MEA Power System

The conventional electrical load modeling method for PMP, which is mostly based on the load shedding idea that detects overloading and disconnects/connects relevant loads according to their priority, can only cope with the “flat” commercial loads (i.e., the electric power requirements without impact on the flight safety; they usually contain cabin illumination system, entertainment system, and galley refrigerators) and avionics loads (i.e., the electric power for the operation of the avionics-related computers located on the equipment cabinet in the cockpit) but cannot efficiently reflect the increased safety-critical electrical loads which are highly variable associated with the aircraft’s flight phases represented by the change of the average altitude, velocity, and attitudehs,Vs,θs under the MEA environment. The increased electrical load behaviors for the MEA cannot be simply modeled through the rule-of-thumb shedding policy but should be explicitly represented with the basic power requirement of their associated safety-critical subsystems. The main increased safety-critical electrical loads are usually the following four kinds:
(1)
The first kind of the increased electrical loads for the MEA is used for continuously pumping the fuel into the engine combustion chamber so as to keep the safe and stable operation for the engine [12, 23]. The pump is electrically driven by a direct current (DC) motor of which the load demand is equal to the value of the fuel volume flow rate injecting into the chamber multiply with the chamber inner pressure. The fuel volume (or mass) flow rate and the inner pressure of chamber can be determined through the stable operation equations of the engine [23] based on the required generator power output and the expected engine thrust given by its associated hs,Vs,θs. Although the current civil aircraft have multiple engine types including turbofan and turboprop, the general operation characteristics of the engines always obey the standard Brayton thermodynamic cycle and can all be simplified by a single compressor-turbine axis work mode [12]. The general operation flowchart of an aircraft engine and its required electrical load for fuel pumping is addressed in Figure 2. Its associated Brayton cycle is depicted in Figure 3Figure 2
The general operation flowchart of an aircraft engine.Figure 3
The Brayton cycle for an aircraft engineFrom [23], if the values hs,Vs,θs are given, the Mach number Mas is Vs/γi⋅gc⋅R⋅T0s; at the altitude hs, the air’s average total pressure is p1s=p0s⋅1+Mas2γi−1/2γi/γi−1 and the average total temperature is T1s=T0s⋅1+Mas2γi−1/2. The air mass flow rate into the engine inlet in period s ṁ0s is
(1)γigcR⋅Mas⋅1+Mas2γi−12−γi+1/2γi−1⋅p1sT1s⋅A1s,where Ais is the cross-sectional area of the engine inlet in period s which is taken as a decision variable controlled by the Full Authority Digital Engine Control (FADEC) computer [23]. Based on the aircraft inlet theory [12], we have p2s=πi⋅p1s and T2s=T1s. From Figures 2 and 3, the stable operation equations of the aircraft engine can be represented by equations (2), (3), (4), (5), (6), (7), and (8):
(2)p3sp2s=1+ηcT3sT2s−1γc/γc−1,s=1,2,⋯,S.Equation (2) represents the isentropic change characteristic of the air state variables within the engine compressor. ηc<1 is the fixed operation efficiency parameter of the engine compressor.
(3)p4s=p3s,s=1,2,⋯,S,(4)fs=cpt⋅T4s−cpc⋅T3sηb⋅hpr−cpt⋅T4s,s=1,2,⋯,S.Equations (3) and (4) describe that there has an isobaric characteristic of the air thermodynamic process in the combustion chamber where the pressure stays constant; furthermore, the total energy taken by the exported mixed gas from the combustion chamber is equal to the total energy of the imported air plus the combustion energy released by the fuel.
(5)T5sT4s=1−ηt1−p5sp4sγt−1/γt,s=1,2,⋯,S.Equation (5) represents the isentropic change characteristic of the air state variables within the engine turbine. ηt<1 is the fixed operation efficiency parameter of the engine turbine.
(6)T3s−T2s=Kc⋅ns2,s=1,2,⋯,S.Equation (6) states the relationships between the rotation speed and the change of total temperature for the engine compressor. Kc is the associated transfer coefficient for the engine compressor.
(7)ṁ0sηm1+fscptT4s−T5s=ṁ0scpcT3s−T2s+Pgs,s=1,2,⋯,S.Since the compressor and turbine are mechanically coaxial, (7) describes the reduced amount of the total enthalpy for the turbine in period s that is equal to the energy output of the generator within this period plus the increased amount of the total enthalpy for the compressor in period s. ηm<1 is the mechanical efficiency.Moreover, if we define the required thrust in periods to be Fs, its value can be determined through the flight mechanics analysis as shown in Appendix A; thus, Fs is represented by the parameters hs,Vs,θs. Then, based on aircraft engine theory [23], we can have
(8)Fs=ṁ0sgc1+fs2cptgcT5s1−p0sp5sγt−1/γt−Vs,s=1,2,⋯,S.Finally, since all decision variables in (2), (3), (4), (5), (6), (7), and (8) are positive and have their associated upper bounds (especially, Ais is also required to be no less than a minimum cross-sectional area A¯i, so as to keep the enough air mass flow rate into the engine for its stable operation), we also have
(9)0≤p3s≤p¯3,0≤T3s≤T¯3,0≤p4s≤p¯4,0≤T4s≤T¯4,A¯i≤Ais≤A¯i,0≤p5s≤p¯5,0≤T5s≤T¯5,0≤ns≤n¯,s=1,2,⋯,S.From the above equations, we can find the load demand of fuel pumping for the engine operation in periods P1s is equal to ṁ0s/ρf⋅fs⋅p3s where ṁ0s and p3s are constrained by (2), (3), (4), (5), (6), (7), (8), and (9).
(2)
The second kind of the increased electrical loads for MEA is the electric drive of the compressor spinning for aircraft environmental control system (ECS) in order to control the temperature and pressure of the aircraft pressure vessel which includes the cockpit, passenger cabin, and interior cargo compartments [24–26]. The ECS contains two main functions: air conditioning and vessel pressurization control. The air conditioning function conditions the fresh air from the outside of the aircraft and supplies it to the pressure vessel at the requested temperature. Here, to note that the air conditioning function in this paper only refers to the regulation of temperature, we do not consider the removal of humidity by the air conditioning since the function of controlling the humidity in the pressure vessel is usually implemented only when the aircraft is on ground [26]. Meanwhile, the pressurization control function regulates the pressure within the vessel by controlling the outflow of air supplied by the air conditioning function and also via an electrically controlled valve which controls the pressure level through the rate of air change in the vessel so as to provide satisfactory pressure values for comfort and safety for all the passengers and crew on the civil aircraft. A typical operation flowchart of an aircraft electrical ECS [9] including its associated electrical load for spinning the compressor is addressed in Figure 4Figure 4
The operation flowchart of an aircraft electrical ECS.In Figure4, at first, since the initial total temperature of the imported air is relatively low, it is usually preheated where the energy is provided by a DC motor. After that, ECS compressor 1 is electrically driven by the same motor so as to realize a higher pressure and temperature for the imported air. Through releasing part of the energy while keeping the pressure unchanged by heat exchanger 1, the imported air is then conditioned by a typical Brayton inverse thermodynamic cycle so as to reach an appropriate pressure and low temperature which can be injected into the vessel [27]. Finally, an electrical control valve at position A will control the air mass flow rate exhausted outside the vessel so as to manage the pressure of vessel to the requested level; the imported air with low temperature will balance the temperature of vessel at a comfortable level.From the previous deductions,Mas, p1s, and T1s can be determined if the values hs,Vs,θs are given for period s. Moreover, since the cross-sectional area of the inlet for ECS compressor 1 A0e is usually taken as a fixed value, the mass flow rate of the imported air in period s ṁ0es is equal to
(10)γigcR⋅Mas⋅1+Mas2γi−12−γi+1/2γi−1⋅p1sT1s⋅A0e.In this paper, we assume the pressure ratio driven by the inlet for ECS compressor 1 is 1; then, the stable operation of the ECS can be represented by equations (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), and (24) [26]:
(11)p2esp1s=1+ηc1eT2esT1es−1γce/γce−1,s=1,2,⋯,S,(12)T2es−T1s=Kc1e⋅nc1es2,s=1,2,⋯,S.Equation (11) represents the isentropic change characteristic of the air state variables within ECS compressor 1. ηc1e<1 is the operation efficiency parameter of ECS compressor 1. Equation (12) states the relationships between the rotation speed and the change of total temperature for ECS compressor 1. Kc1e is defined to be the associated transfer coefficient for ECS compressor 1.Based on the heat transfer theory [27], there has an isobaric characteristic of the air thermodynamic process in heat exchanger 1 where the pressure stays constant; the total energy taken by the cold side air (the air imported through the inlet of the heat exchanger so as to provide the low-temperature air resource) is equal to the total energy released by the hot side air (the high-temperature air from ECS compressor 1); furthermore, the cold side air and hot side air reach a common temperature after the heat transfer within heat exchanger 1. The above characteristics are represented by
(13)p2es=p3es,s=1,2,⋯,S,(14)ṁ1es⋅T3es−T1s=ṁ0es⋅T2es−T3es,s=1,2,⋯,S.In (13) and (14), the mass flow rate of the cold side air into heat exchanger 1 in period s, ṁ1es, can be stated as
(15)γigcR⋅Mas⋅1+Mas2γi−12−γi+1/2γi−1⋅p1sT1s⋅A1e.where A1e is the cross-sectional area of the inlet at the cold side of heat transfer 1.
(16)p4esp3es=1+ηc2eT4esT3es−1γce/γce−1,s=1,2,⋯,S.Eq. (16) represents the isentropic change characteristic of the air state variables within the ECS compressor 2. ηc2e<1 is the operation efficiency parameter of ECS compressor 2.
(17)p5es=p4es,s=1,2,⋯,S,(18)ṁ1es⋅T5es−T3es=ṁ0es⋅T4es−T5es,s=1,2,⋯,S.Similar with equations (13) and (14), the isobaric characteristic of the air thermodynamic process in heat exchanger 2 is represented by equations (17) and (18). The cold side air flow into heat exchanger 2 in period s is the exported air flow at the cold side of heat exchanger 1.
(19)T6esT5es=1−ηte1−p6esp5esγte−1/γte,s=1,2,⋯,S,Equation (19) represents the isentropic change characteristic of the air state variables within the ECS turbine. ηte<1 is the operation efficiency parameter of the ECS turbine.
(20)T4es−T3es=Kc2e⋅nc2es2,s=1,2,⋯,S.Equation (20) states the relationships between the rotation speed and the change of total temperature for ECS compressor 2. Kc2e is defined to be the associated transfer coefficient for ECS compressor 2.
(21)T4es−T3es=ηme⋅T5es−T6es,s=1,2,⋯,S.Since compressor 2 and the turbine are mechanically coaxial, (21) describes that the reduced amount of the total enthalpy in period s for the ECS turbine is equal to the increased amount of the total enthalpy for ECS compressor 2 in period s. ηme<1 is the associated mechanical efficiency.
(22)mvs60sec+ṁ0es−ṁouts⋅Tvs+1=Qv1s+Qv2s+Qv3s+Qv4scpe⋅60sec+ṁ0es⋅T6es−ṁouts⋅Tvs+mvs60sec⋅Tvs,s=1,2,⋯,S.Equation (22) describes the temperature control characteristics in the vessel where the enthalpy change amount of the air in vessel at period s is equal to the sum of the imported enthalpy by the upstream air, the exported enthalpy by the exhausted air from valve A, the enthalpy introduced by the passenger/pilot activity, the heat transfer between the aircraft’s skin and the air, the radiation by the sun, and avionics device operation.
(23)p6es=mvs⋅R⋅TvsDv,s=1,2,⋯,S.Equation (23) represents the pressure balance constraints in vessel where Dv is the fixed volume of the vessel.Finally, all decision variables in (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), and (23) are all positive values and have their associated upper bounds; besides that, during the whole flight, the pressure and temperature within the vessel are separately required to be larger than the value maxp~,p1s (p~ is the pressure value over the sea level by 1800 m) and the minimum comfortable temperature T¯v. Thus, we have
(24)T1s≤T1es≤T¯1e,0≤T2es≤T¯2e,0≤p2es≤p¯2e,0≤T3es≤T¯3e,0≤p3es≤p¯3e,0≤T4es≤T¯4e,0≤p4es≤p¯4e,0≤T5es≤T¯5e,0≤p5es≤p¯5e,0≤T6es≤T¯6e,maxp~,p1s≤p6es≤p¯6e,0≤nc1es≤n¯c1e,0≤nc2es≤n¯c2e,T¯v≤Tvs≤T¯v,0≤mvs≤m¯v,0≤ṁouts≤ṁ¯out,s=1,2,⋯,S.Based on the above equations, we can find that the load demand of the electrical DC motor for the ECS in periodsP2s is equal to ṁ0es⋅cpe⋅T2es−T1s. T2es for s=1,2,⋯,S are constrained by equations (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), and (24).
(3)
The third kind of the increased loads for the MEA is the electrical drive of the hydraulic pump for the primary flight control actuators, such as the electro-hydrostatic actuators for the surfaces (elevators, ailerons, and rudders) [28, 29]. In the electrically driven flight control actuation system, a DC motor energizes the pump operation of an actuator in order to maintain the required hydraulic oil mass flow rate and pressure to the actuation system. In this paper, since we have assumed the aircraft are to fly at fixed pitch angle and zero roll and yaw angles within a given time period, we will not consider the electrical loads introduced by actuating the elevators for a vertical maneuver and actuating the ailerons and rudders for a horizontal maneuver. From the flight control actuation theory [29], the only required power within a period so as to maintain the expected flight behaviors is to keep the elevators actuating with fixed deflection angle and it can thus be calculated according to the specified flight phase which is represented by hs,Vs,θs. We define the required air-dynamical force on elevators in period s to be Les, and its value can thus be determined through the flight mechanics analysis shown in Appendix B; in other words, Les can be represented by the parameters hs,Vs,θs. Finally, we have known that the electrical load required by the flight control system in period s P3s is Kleak⋅Les/Ae2/ρhy where Kleak is the leak coefficient addressing the relationship between the air-dynamical force on elevators and required hydraulic oil mass flow rate [29]. ρhy is the density of the elevators’ hydraulic oil and Ae is the area of the elevators
(4)
The fourth kind of the increased loads for the MEA is the electrical warming of the leading-edge wings and elevators for ice protection [30, 31]. The icing, which is usually caused by the supercooled water droplets in clouds where the aircraft is passing by, has been recognized as a potentially dangerous phenomenon for aircraft’s flight, since the buildup of ice on the leading edges of wings/elevators can cause an irregular change of the aircraft’s maneuverability characteristics and may finally result in the aircraft’s loss of control [30]. Anti-icing refers to the prevention of any buildup of ice on a wing/elevator during flight and can be achieved by the use of thermal energy. Under the MEA environment, the electrical energy is used for anti-icing via electrothermal heaters where the electrothermal pads are bonded to the outer surface in need of anti-icing. The heat generated by these electrothermal pads is conducted to the wing’s/elevator’s interface where a thin film of water is formed. This breaks the adhesive bond between the ice and the outer surface of the pad so that the ice can be swept away by aerodynamic forces. The anti-icing electrical power is usually implemented by cycling the load between the different electrothermal heaters such as the electrothermal heaters on the wings and elevators. This will reduce the electrical power need while realizing the minimal anti-ice requirements. In this paper, we define electrical loads for anti-icing the leading-edge wings and elevators in period s P4s to be separately equal to Pws and Pes which are activated iteratively in every 3 minutes after the aircraft’s altitude becomes higher than the minimum altitude with icing condition. Pws and Pes are separately the rated power values of the electrothermal heaters on the wings and elevators, respectively


## 2.2. Smart PMP Formulation

The PMP formulation firstly needs to explicitly account for the objective function with a minimal fuel cost and battery operation cost. The decision-making of a PMP is to implement a codispatch between the gas turbine generator, storage devices, and all controllable electrical loads in an optimal fashion. An effective quantification of the system operation costs in the PMP formulation is stated in the following:(25)min∑s=1SCgs+Cbs.Equation (25) describes that the total operation cost for the whole scheduling horizon is the sum of the fuel cost by generator output, battery operation costs by battery charging and discharging numbers. In (25), Cgs is equal to ag⋅ṁ0s⋅fs where ag is the linear coefficient between the fuel mass flow rate and fuel cost; Cbs is equal to hb⋅zbcs+zbds where hb is the operation cost coefficient.We next state each of the constraints we include in the problem formulation.


### 2.2.1. Battery Charging and Discharging Behavior Constraints

(26)zbcs+zbds≤1,zbds⋅P¯bo−zbcs⋅P¯bi≤Pbs,Pbs≤zbds⋅P¯bo−zbcs⋅P¯bi,s=1,2,⋯,S,(27)E¯b≤Ebs≤E¯b,Ebs+1=Ebs−Pbs⋅60sec,s=1,2,⋯,S.Equation (26) stipulates that in every time period, the state of the battery will be one of the three cases: charging (the discrete variable zbcs∈0,1 is 1), discharging (the discrete variable zbds∈0,1 is 1), and shut-off (zbcs and zbds are both 0); the charge/discharge power in period s Pbs is within its upper and lower bound. Equation (27) addresses the upper and lower bound constraints for state-of-charge (SOC) energy in every time period and its dynamic characteristic.


### 2.2.2. Load Balance Constraints

(28)Pgs+Pbs=P1s+P2s+P3s+P4s+P5s,s=1,2,⋯,S.In (28), P5s is defined to be the fixed power load for commercial and avionics services in period s. Equation (28) stipulates that in every time period, the total power generated by generator and battery is equal to the total load demands by engine fuel pumping, ECS, flight control actuation, anti-ice for leading edge of wings/elevators, and commercial and avionics services.


### 2.2.3. Safety-Critical Controllable Load Constraints

Based on Section2.1, all associated constraints with the safety-critical controllable loads are listed in equations (2), (3), (4), (5), (6), (7), (8), (9), (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), and (24).


### 2.2.4. Generator Power Output Constraints

(29)0≤Pgs≤P¯g,s=1,2,⋯,S.Equation (29) addresses that the power output by the generator in period s should be positive and less than the rated power capacity P¯g.


## 2.2.1. Battery Charging and Discharging Behavior Constraints

(26)zbcs+zbds≤1,zbds⋅P¯bo−zbcs⋅P¯bi≤Pbs,Pbs≤zbds⋅P¯bo−zbcs⋅P¯bi,s=1,2,⋯,S,(27)E¯b≤Ebs≤E¯b,Ebs+1=Ebs−Pbs⋅60sec,s=1,2,⋯,S.Equation (26) stipulates that in every time period, the state of the battery will be one of the three cases: charging (the discrete variable zbcs∈0,1 is 1), discharging (the discrete variable zbds∈0,1 is 1), and shut-off (zbcs and zbds are both 0); the charge/discharge power in period s Pbs is within its upper and lower bound. Equation (27) addresses the upper and lower bound constraints for state-of-charge (SOC) energy in every time period and its dynamic characteristic.


## 2.2.2. Load Balance Constraints

(28)Pgs+Pbs=P1s+P2s+P3s+P4s+P5s,s=1,2,⋯,S.In (28), P5s is defined to be the fixed power load for commercial and avionics services in period s. Equation (28) stipulates that in every time period, the total power generated by generator and battery is equal to the total load demands by engine fuel pumping, ECS, flight control actuation, anti-ice for leading edge of wings/elevators, and commercial and avionics services.


## 2.2.3. Safety-Critical Controllable Load Constraints

Based on Section2.1, all associated constraints with the safety-critical controllable loads are listed in equations (2), (3), (4), (5), (6), (7), (8), (9), (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), and (24).


## 2.2.4. Generator Power Output Constraints

(29)0≤Pgs≤P¯g,s=1,2,⋯,S.Equation (29) addresses that the power output by the generator in period s should be positive and less than the rated power capacity P¯g.


## 3. The Proposed Solution Algorithm

### 3.1. Mathematical Structure Analysis for the PMP Formulation and the Architecture of the Proposed Two-Level Decomposition Algorithm

The mathematical presentation of the PMP has been stated by (2), (3), (4), (5), (6), (7), (8), (9), (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), (24), (25), (26), (27), (28), and (29), and this proposed model is a large-scale, mixed-integer and nonlinear programming problem (MINLP) which can be constructed by two parts: (25), (26), (27), (28), and (29) address a standard real power dispatch problem which has a mixed-integer linear program (MILP); (2), (3), (4), (5), (6), (7), (8), and (9) and (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), and (24) separately state the associated two kinds of constraints which stipulate the safety-critical controllable loads within the aircraft engine and ECS stable operation regions. Because the subclass of the real power dispatch problem in (25), (26), (27), (28), and (29) as a MILP is NP-hard and the subclass of the constraints in (2), (3), (4), (5), (6), (7), (8), and (9) and (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), and (24) as the nonconvex and nonlinear equations are also NP-hard, the generalized MINLP problem in (2), (3), (4), (5), (6), (7), (8), (9), (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), (24), (25), (26), (27), (28), and (29) combines the difficulties of both subclasses and is thus very challenging to solve. In this paper, we tackle the above MINLP problem by Generalized Benders method [32–34] which is developed with exploiting the special structure of large-scale mathematical programming problems by selecting the complicating variables which divide the original problem into the master and slave problems and, when fixed, render the slave problem considerably more tractable [32–34]. As for a mixed-integer program such as the problem in (2), (3), (4), (5), (6), (7), (8), (9), (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), (24), (25), (26), (27), (28), and (29), the master problem—a mixed-integer optimization problem—deals with only the discrete variables and part of the continuous variables. Then, the solutions of the variables in the master problem are taken as the complicating variables and fixed as the parameters in the slave problem which thus usually has a linear or nonlinear convex optimization form where the solutions of the slave problem can provide the information so as to help the master problem formulate the associated Benders cuts. The Benders cuts reduce the search region for the master problem and usually have a linear or linearly approximated formulation about the solution point found in the nonlinear programming (NLP) slave problem. Therefore, the master problem is solved successively, and the slave problem as well, until they both reach the required convergence criterion.As analyzed above, we at first present a compact matrix formulation for the PMP in (2), (3), (4), (5), (6), (7), (8), (9), (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), (24), (25), (26), (27), (28), and (29) as follows:
(30)minx,y,zcTx(31)s.t.Fx≤f,(32)Hx=h,(33)B1x,y=b1,(34)B2x,z=b2,(35)x∈X,(36)y∈Y,(37)z∈Z.The variablex is a vector of the decisions including p3s, p4s, T4s, Ais, Pgs, p2es, T2es, p3es, T3es, T5es, P6es, nc1es, mvs, ṁouts, Pbs, Ebs, Zbcs, and Zbds for each time period s (s=1,2,⋯,S) of the whole PMP scheduling horizon. The variables y and z are separately the vector of the rest of decision variables related to maintain a stable operation for the engine and ECS.The objective function is to minimize the sum of the fuel cost by the generator output and the battery operation costs caused by its charging and discharging numbers. Since the above two kinds of costs are only corresponding to the variablex based on (25), the objective function can be stated as the linear expression in (30). The linear constraints in (31), involving only the variable x, contain the constraints which are stated by (3), (13), (14), (26), and (27). The nonlinear constraints in (32), which represent part of the constraints stipulating the stable operation for the engine and ECS while corresponding only to the variables x, include (12) and (28). The nonlinear constraints (33) and (34), which couple the variables x with y and x with z, separately include (2), (4), (5), (6), (7), (8), (11), (16), (17), (18), (19), (20), (21), (22), and (23). The compact sets X, Y, and Z in (35), (36), and (37) represent the upper and lower bound constraints stated by (9) and (24).The structure for the PMP formulation stated by (30), (31), (32), (33), (34), (35), (36), and (37) has the following special characteristics: (1) the variables have a clear partition structure which has a merit on using the decomposition method; (2) although the nonlinear expression in (30), (31), (32), (33), (34), (35), (36), and (37) has a nonconvex basics, if we take x as the complicating variables under the GBD architecture, (33) and (34) will separately have a convex form for y and z with given x. Meanwhile, the master problem can be represented by a simplified MINLP formulation and be easily linear-approximately restated by a MILP problem which is finally solved by the mature commercialized software package. The above key characteristics will permit us to devise a GBD-based two-level decomposition scheme for the PMP: the outer level is to find an optimal real power generation, load management dispatch, and battery charge/discharge decision result with only considering the complicating variables x. The inner level is then devised to find the feasible schedule results for the constraints in (33) and (36) and (34) and (37). In other words, the constraints in (33) and (36) and (34) and (37) are not taken into account in the initial master problem while the adjustments of the associated continuous variables for the engine and ECS are separately solved in the two slave problems at the inner level. The above two-level algorithm is then depicted in Figure 5: the outer level employs the GBD-type cutting plane algorithm to obtain the optimal real power outputs, load management, and battery charge/discharge decisions using the cuts information generated from the inner level, which in parallel and iteratively solves two nonlinear convex optimization problems through the successive linear programming (SLP) techniques [35]. If the two inner slave problems cannot converge or any violations of the bounds on the variables exist, the corresponding Benders cuts will be active and fed back into the next round of calculation for the master problem. These cuts represent the coupled information on the power output, load management, and battery charge/discharge decisions with the associated operation constraints for the engine and ECS, and this interaction is in essence represented by the Benders coefficients.Figure 5
Two-level architecture of the proposed decomposition algorithm for the PMP formulation.


### 3.2. Algorithm Description

The description of the proposed algorithm is given as follows.


#### 3.2.1. Outer Level: GBD Algorithm

(1) Initialization. Let x0∗ be the optimal solution for the optimization problem: minxcTx subject to Fx≤f, Hx=h, x∈X. Solve the two inner problems based on x0∗ (the solution methodology is given after in this subsection) and let y0∗ and z0∗ be separately the optimal solutions of the two inner problems based on x0∗. Set the outer level lower bound LO−level=−∞, upper bound UO−level=+∞, and the number of iteration k=1. Choose an outer level convergence tolerance level ζ (>0).(2) Iterationk≥1.Step 1.
Solve the GBD master problem. The master problem is the following MILP program:minxk,αk,βkcTxk+MTαk+MTβk subject to αk≥B1xk,yk‐1∗−b1, βk≥B2xk,zk‐1∗−b2, Fxk≤f, Hxk=h, and xk∈X. αk and βk are the introduced slack variables which have positive values, and M is a column vector where all of its elements are M; yk‐1∗ and zk‐1∗ are determined in the (k-1)th round of inner level problems and are thus taken as given parameters. From (33) and (34), we can find that part of the functions in B1xk,yk‐1∗ and B2xk,zk‐1∗ is nonlinear corresponding to xk; moreover, Hxk has also a nonlinear form for xk. Thus, to solve the master problem, it can be at first recast by its linear-approximated MILP formulation [36, 37] and then resort to be solved by the commercial optimization software package. Let xk∗,αk∗,βk∗ be the optimal solution for the kth round of master problem. We set LO‐level=cTxk∗+MTαk∗+MTβk∗.Step 2.
Check the outer level convergence. IfUO−level−LO−level<ζ, stop and return xk∗,yk∗,zk∗. Otherwise, set UO−level=cTxk∗+MTαk∗+MTβk∗ and go to Step 3.Step 3.
Solve two slave problems. We will discuss the methodology to solve the inner level slave problems in the following. Letyk∗ and zk∗ be separately the optimal solutions of the two inner problems based on xk∗ and input them as parameters into (k+1)th round of the master problem. Let k=k+1, go to Step 1.


#### 3.2.2. Inner Level: SLP Technique

A SLP algorithm is used to solve the two inner optimization problems where we at first present their formulations. All of the nonlinear terms in the inner problems are then successively linearized around intermediate solution points, and the Benders cuts generated by solving the inner problems are finally added into the master problem. Since the two inner level problems are both convex,δ optimal solution of the inner problems is guaranteed by the SLP algorithm where δ is the inner level convergence tolerance level and this solution can thus generate a valid cut for the outer level master problem.(1) Initialization. Fix the decision of xk∗ (k≥0) passed from the kth iteration of the outer level master problem. Set the two inner level lower bounds L1I−level=L2I−level=−∞, upper bounds U1I−level=U2I−level=+∞, and number of iteration j=1. Choose the inner level convergence tolerance level δ (>0).(2) Iterationj≥1.Step 1.
We at first present the two slave problems as follows. SP I: find a feasible solution ofykj within B1xk∗,ykj=b1, ykj∈Y; SP II: find a feasible solution of zkj within B2xk∗,zkj=b2, zkj∈Z. We can equivalently restate the problems SP I and SP II separately as the following optimization problems: R1jxk∗=minykj∈Y,s1,kj,s2,kj1Ts1,kj+1Ts2,kj subject to B1xk∗,ykj−s1,kj+s2,kj=b1, s1,kj≥0,s2,kj≥0., and R2jxk∗=minzkj∈Z,s3,kj,s4,kj1Ts3,kj+1Ts4,kj subject to B2xk∗,zkj−s3,kj+s4,kj=b2,, s3,kj≥0,s4,kj≥0. s1,kj, s2,kj, s3,kj, and s4,kj are the introduced slack continuous variables.Step 2.
Thejth round of iteration for solving the slave problems. From the above, we can find that, if ykj and zkj are given, the above two optimization problems can be both reformulated as the linear programming problem and their optimal values, which are separately defined as G1jxk∗,ykj and Gj2xk∗,zkj, can be determined analytically as G1jxk∗,ykj=1TB1xk∗,ykj−b1 and Gj2xk∗,zkj=1TB2xk∗,zkj−b2. Thus, we define y~kj,z~kj to be the linear approximation point at the jth round of SLP for the above two optimization problems. Let yk−1∗,zk−1∗ be equal to the initial solution of y~kj,z~kj : y~k1,z~k1 and we stipulate y~kj+1=y~kj+d^kj and z~kj+1=z~kj+d¯kj where d^kj and d¯kj are separately determined by solving the following two convex quadratic programming problems (d^kj and d¯kj are separately strictly decreasing direction for the value of G1jxk∗,y~kj and Gj2xk∗,z~kj as j increases; the components of vector Δ are all very small positive real values):
(38)maxd^kjd^kjd^kjT(39)subjecttoB1xk∗,y~kjμ−b1μ⋅∇B1xk∗,y~kjd^kjμ<0,∀μ=1,2,⋯,6⋅S∇B1xk∗,y~kjd^kj<B1xk∗,y~kj−b1,−Δ≤d^kj≤Δ,y~kj+d^kj∈Y,(40)maxd¯kjd¯kjd¯kjT(41)subjecttoB2xk∗,z~kjμ−b2μ⋅∇B2xk∗,z~kjd¯kjμ<0,∀μ=1,2,⋯,9⋅S∇B2xk∗,z~kjd¯kj<B2xk∗,z~kj−b2,−Δ≤d¯kj≤Δ,z~kj+d¯kj∈Z.Step 3.
Check the inner level convergence. At first, letL1I−level=G1jxk∗,y~kj and L2I−level=G2jxk∗,z~kj. If U1I−level=0 and U2I−level=0 or U1I−level−L1I−level<δ and U2I−level−L2I−level<δ, then stop and output the current solutions y~kj, z~kj, G1jxk∗,y~kj, and G2jxk∗,z~kj as yk∗, zk∗, G1xk∗,yk∗, and G2xk∗,zk∗ into the master problem; B1xk+1,yk∗−b1 and B2xk+1,zk∗−b2 are used to construct the associated Benders cuts for the master problem. Otherwise, let U1I−level=G1jxk∗,y~kj and U2I−level=G2jxk∗,z~kj and further set j=j+1 and go to Step 2 of the SLP algorithm.


## 3.1. Mathematical Structure Analysis for the PMP Formulation and the Architecture of the Proposed Two-Level Decomposition Algorithm

The mathematical presentation of the PMP has been stated by (2), (3), (4), (5), (6), (7), (8), (9), (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), (24), (25), (26), (27), (28), and (29), and this proposed model is a large-scale, mixed-integer and nonlinear programming problem (MINLP) which can be constructed by two parts: (25), (26), (27), (28), and (29) address a standard real power dispatch problem which has a mixed-integer linear program (MILP); (2), (3), (4), (5), (6), (7), (8), and (9) and (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), and (24) separately state the associated two kinds of constraints which stipulate the safety-critical controllable loads within the aircraft engine and ECS stable operation regions. Because the subclass of the real power dispatch problem in (25), (26), (27), (28), and (29) as a MILP is NP-hard and the subclass of the constraints in (2), (3), (4), (5), (6), (7), (8), and (9) and (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), and (24) as the nonconvex and nonlinear equations are also NP-hard, the generalized MINLP problem in (2), (3), (4), (5), (6), (7), (8), (9), (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), (24), (25), (26), (27), (28), and (29) combines the difficulties of both subclasses and is thus very challenging to solve. In this paper, we tackle the above MINLP problem by Generalized Benders method [32–34] which is developed with exploiting the special structure of large-scale mathematical programming problems by selecting the complicating variables which divide the original problem into the master and slave problems and, when fixed, render the slave problem considerably more tractable [32–34]. As for a mixed-integer program such as the problem in (2), (3), (4), (5), (6), (7), (8), (9), (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), (24), (25), (26), (27), (28), and (29), the master problem—a mixed-integer optimization problem—deals with only the discrete variables and part of the continuous variables. Then, the solutions of the variables in the master problem are taken as the complicating variables and fixed as the parameters in the slave problem which thus usually has a linear or nonlinear convex optimization form where the solutions of the slave problem can provide the information so as to help the master problem formulate the associated Benders cuts. The Benders cuts reduce the search region for the master problem and usually have a linear or linearly approximated formulation about the solution point found in the nonlinear programming (NLP) slave problem. Therefore, the master problem is solved successively, and the slave problem as well, until they both reach the required convergence criterion.As analyzed above, we at first present a compact matrix formulation for the PMP in (2), (3), (4), (5), (6), (7), (8), (9), (11), (12), (13), (14), (16), (17), (18), (19), (20), (21), (22), (23), (24), (25), (26), (27), (28), and (29) as follows:
(30)minx,y,zcTx(31)s.t.Fx≤f,(32)Hx=h,(33)B1x,y=b1,(34)B2x,z=b2,(35)x∈X,(36)y∈Y,(37)z∈Z.The variablex is a vector of the decisions including p3s, p4s, T4s, Ais, Pgs, p2es, T2es, p3es, T3es, T5es, P6es, nc1es, mvs, ṁouts, Pbs, Ebs, Zbcs, and Zbds for each time period s (s=1,2,⋯,S) of the whole PMP scheduling horizon. The variables y and z are separately the vector of the rest of decision variables related to maintain a stable operation for the engine and ECS.The objective function is to minimize the sum of the fuel cost by the generator output and the battery operation costs caused by its charging and discharging numbers. Since the above two kinds of costs are only corresponding to the variablex based on (25), the objective function can be stated as the linear expression in (30). The linear constraints in (31), involving only the variable x, contain the constraints which are stated by (3), (13), (14), (26), and (27). The nonlinear constraints in (32), which represent part of the constraints stipulating the stable operation for the engine and ECS while corresponding only to the variables x, include (12) and (28). The nonlinear constraints (33) and (34), which couple the variables x with y and x with z, separately include (2), (4), (5), (6), (7), (8), (11), (16), (17), (18), (19), (20), (21), (22), and (23). The compact sets X, Y, and Z in (35), (36), and (37) represent the upper and lower bound constraints stated by (9) and (24).The structure for the PMP formulation stated by (30), (31), (32), (33), (34), (35), (36), and (37) has the following special characteristics: (1) the variables have a clear partition structure which has a merit on using the decomposition method; (2) although the nonlinear expression in (30), (31), (32), (33), (34), (35), (36), and (37) has a nonconvex basics, if we take x as the complicating variables under the GBD architecture, (33) and (34) will separately have a convex form for y and z with given x. Meanwhile, the master problem can be represented by a simplified MINLP formulation and be easily linear-approximately restated by a MILP problem which is finally solved by the mature commercialized software package. The above key characteristics will permit us to devise a GBD-based two-level decomposition scheme for the PMP: the outer level is to find an optimal real power generation, load management dispatch, and battery charge/discharge decision result with only considering the complicating variables x. The inner level is then devised to find the feasible schedule results for the constraints in (33) and (36) and (34) and (37). In other words, the constraints in (33) and (36) and (34) and (37) are not taken into account in the initial master problem while the adjustments of the associated continuous variables for the engine and ECS are separately solved in the two slave problems at the inner level. The above two-level algorithm is then depicted in Figure 5: the outer level employs the GBD-type cutting plane algorithm to obtain the optimal real power outputs, load management, and battery charge/discharge decisions using the cuts information generated from the inner level, which in parallel and iteratively solves two nonlinear convex optimization problems through the successive linear programming (SLP) techniques [35]. If the two inner slave problems cannot converge or any violations of the bounds on the variables exist, the corresponding Benders cuts will be active and fed back into the next round of calculation for the master problem. These cuts represent the coupled information on the power output, load management, and battery charge/discharge decisions with the associated operation constraints for the engine and ECS, and this interaction is in essence represented by the Benders coefficients.Figure 5
Two-level architecture of the proposed decomposition algorithm for the PMP formulation.


## 3.2. Algorithm Description

The description of the proposed algorithm is given as follows.


### 3.2.1. Outer Level: GBD Algorithm

(1) Initialization. Let x0∗ be the optimal solution for the optimization problem: minxcTx subject to Fx≤f, Hx=h, x∈X. Solve the two inner problems based on x0∗ (the solution methodology is given after in this subsection) and let y0∗ and z0∗ be separately the optimal solutions of the two inner problems based on x0∗. Set the outer level lower bound LO−level=−∞, upper bound UO−level=+∞, and the number of iteration k=1. Choose an outer level convergence tolerance level ζ (>0).(2) Iterationk≥1.Step 1.
Solve the GBD master problem. The master problem is the following MILP program:minxk,αk,βkcTxk+MTαk+MTβk subject to αk≥B1xk,yk‐1∗−b1, βk≥B2xk,zk‐1∗−b2, Fxk≤f, Hxk=h, and xk∈X. αk and βk are the introduced slack variables which have positive values, and M is a column vector where all of its elements are M; yk‐1∗ and zk‐1∗ are determined in the (k-1)th round of inner level problems and are thus taken as given parameters. From (33) and (34), we can find that part of the functions in B1xk,yk‐1∗ and B2xk,zk‐1∗ is nonlinear corresponding to xk; moreover, Hxk has also a nonlinear form for xk. Thus, to solve the master problem, it can be at first recast by its linear-approximated MILP formulation [36, 37] and then resort to be solved by the commercial optimization software package. Let xk∗,αk∗,βk∗ be the optimal solution for the kth round of master problem. We set LO‐level=cTxk∗+MTαk∗+MTβk∗.Step 2.
Check the outer level convergence. IfUO−level−LO−level<ζ, stop and return xk∗,yk∗,zk∗. Otherwise, set UO−level=cTxk∗+MTαk∗+MTβk∗ and go to Step 3.Step 3.
Solve two slave problems. We will discuss the methodology to solve the inner level slave problems in the following. Letyk∗ and zk∗ be separately the optimal solutions of the two inner problems based on xk∗ and input them as parameters into (k+1)th round of the master problem. Let k=k+1, go to Step 1.


### 3.2.2. Inner Level: SLP Technique

A SLP algorithm is used to solve the two inner optimization problems where we at first present their formulations. All of the nonlinear terms in the inner problems are then successively linearized around intermediate solution points, and the Benders cuts generated by solving the inner problems are finally added into the master problem. Since the two inner level problems are both convex,δ optimal solution of the inner problems is guaranteed by the SLP algorithm where δ is the inner level convergence tolerance level and this solution can thus generate a valid cut for the outer level master problem.(1) Initialization. Fix the decision of xk∗ (k≥0) passed from the kth iteration of the outer level master problem. Set the two inner level lower bounds L1I−level=L2I−level=−∞, upper bounds U1I−level=U2I−level=+∞, and number of iteration j=1. Choose the inner level convergence tolerance level δ (>0).(2) Iterationj≥1.Step 1.
We at first present the two slave problems as follows. SP I: find a feasible solution ofykj within B1xk∗,ykj=b1, ykj∈Y; SP II: find a feasible solution of zkj within B2xk∗,zkj=b2, zkj∈Z. We can equivalently restate the problems SP I and SP II separately as the following optimization problems: R1jxk∗=minykj∈Y,s1,kj,s2,kj1Ts1,kj+1Ts2,kj subject to B1xk∗,ykj−s1,kj+s2,kj=b1, s1,kj≥0,s2,kj≥0., and R2jxk∗=minzkj∈Z,s3,kj,s4,kj1Ts3,kj+1Ts4,kj subject to B2xk∗,zkj−s3,kj+s4,kj=b2,, s3,kj≥0,s4,kj≥0. s1,kj, s2,kj, s3,kj, and s4,kj are the introduced slack continuous variables.Step 2.
Thejth round of iteration for solving the slave problems. From the above, we can find that, if ykj and zkj are given, the above two optimization problems can be both reformulated as the linear programming problem and their optimal values, which are separately defined as G1jxk∗,ykj and Gj2xk∗,zkj, can be determined analytically as G1jxk∗,ykj=1TB1xk∗,ykj−b1 and Gj2xk∗,zkj=1TB2xk∗,zkj−b2. Thus, we define y~kj,z~kj to be the linear approximation point at the jth round of SLP for the above two optimization problems. Let yk−1∗,zk−1∗ be equal to the initial solution of y~kj,z~kj : y~k1,z~k1 and we stipulate y~kj+1=y~kj+d^kj and z~kj+1=z~kj+d¯kj where d^kj and d¯kj are separately determined by solving the following two convex quadratic programming problems (d^kj and d¯kj are separately strictly decreasing direction for the value of G1jxk∗,y~kj and Gj2xk∗,z~kj as j increases; the components of vector Δ are all very small positive real values):
(38)maxd^kjd^kjd^kjT(39)subjecttoB1xk∗,y~kjμ−b1μ⋅∇B1xk∗,y~kjd^kjμ<0,∀μ=1,2,⋯,6⋅S∇B1xk∗,y~kjd^kj<B1xk∗,y~kj−b1,−Δ≤d^kj≤Δ,y~kj+d^kj∈Y,(40)maxd¯kjd¯kjd¯kjT(41)subjecttoB2xk∗,z~kjμ−b2μ⋅∇B2xk∗,z~kjd¯kjμ<0,∀μ=1,2,⋯,9⋅S∇B2xk∗,z~kjd¯kj<B2xk∗,z~kj−b2,−Δ≤d¯kj≤Δ,z~kj+d¯kj∈Z.Step 3.
Check the inner level convergence. At first, letL1I−level=G1jxk∗,y~kj and L2I−level=G2jxk∗,z~kj. If U1I−level=0 and U2I−level=0 or U1I−level−L1I−level<δ and U2I−level−L2I−level<δ, then stop and output the current solutions y~kj, z~kj, G1jxk∗,y~kj, and G2jxk∗,z~kj as yk∗, zk∗, G1xk∗,yk∗, and G2xk∗,zk∗ into the master problem; B1xk+1,yk∗−b1 and B2xk+1,zk∗−b2 are used to construct the associated Benders cuts for the master problem. Otherwise, let U1I−level=G1jxk∗,y~kj and U2I−level=G2jxk∗,z~kj and further set j=j+1 and go to Step 2 of the SLP algorithm.


## 3.2.1. Outer Level: GBD Algorithm

(1) Initialization. Let x0∗ be the optimal solution for the optimization problem: minxcTx subject to Fx≤f, Hx=h, x∈X. Solve the two inner problems based on x0∗ (the solution methodology is given after in this subsection) and let y0∗ and z0∗ be separately the optimal solutions of the two inner problems based on x0∗. Set the outer level lower bound LO−level=−∞, upper bound UO−level=+∞, and the number of iteration k=1. Choose an outer level convergence tolerance level ζ (>0).(2) Iterationk≥1.Step 1.
Solve the GBD master problem. The master problem is the following MILP program:minxk,αk,βkcTxk+MTαk+MTβk subject to αk≥B1xk,yk‐1∗−b1, βk≥B2xk,zk‐1∗−b2, Fxk≤f, Hxk=h, and xk∈X. αk and βk are the introduced slack variables which have positive values, and M is a column vector where all of its elements are M; yk‐1∗ and zk‐1∗ are determined in the (k-1)th round of inner level problems and are thus taken as given parameters. From (33) and (34), we can find that part of the functions in B1xk,yk‐1∗ and B2xk,zk‐1∗ is nonlinear corresponding to xk; moreover, Hxk has also a nonlinear form for xk. Thus, to solve the master problem, it can be at first recast by its linear-approximated MILP formulation [36, 37] and then resort to be solved by the commercial optimization software package. Let xk∗,αk∗,βk∗ be the optimal solution for the kth round of master problem. We set LO‐level=cTxk∗+MTαk∗+MTβk∗.Step 2.
Check the outer level convergence. IfUO−level−LO−level<ζ, stop and return xk∗,yk∗,zk∗. Otherwise, set UO−level=cTxk∗+MTαk∗+MTβk∗ and go to Step 3.Step 3.
Solve two slave problems. We will discuss the methodology to solve the inner level slave problems in the following. Letyk∗ and zk∗ be separately the optimal solutions of the two inner problems based on xk∗ and input them as parameters into (k+1)th round of the master problem. Let k=k+1, go to Step 1.


## 3.2.2. Inner Level: SLP Technique

A SLP algorithm is used to solve the two inner optimization problems where we at first present their formulations. All of the nonlinear terms in the inner problems are then successively linearized around intermediate solution points, and the Benders cuts generated by solving the inner problems are finally added into the master problem. Since the two inner level problems are both convex,δ optimal solution of the inner problems is guaranteed by the SLP algorithm where δ is the inner level convergence tolerance level and this solution can thus generate a valid cut for the outer level master problem.(1) Initialization. Fix the decision of xk∗ (k≥0) passed from the kth iteration of the outer level master problem. Set the two inner level lower bounds L1I−level=L2I−level=−∞, upper bounds U1I−level=U2I−level=+∞, and number of iteration j=1. Choose the inner level convergence tolerance level δ (>0).(2) Iterationj≥1.Step 1.
We at first present the two slave problems as follows. SP I: find a feasible solution ofykj within B1xk∗,ykj=b1, ykj∈Y; SP II: find a feasible solution of zkj within B2xk∗,zkj=b2, zkj∈Z. We can equivalently restate the problems SP I and SP II separately as the following optimization problems: R1jxk∗=minykj∈Y,s1,kj,s2,kj1Ts1,kj+1Ts2,kj subject to B1xk∗,ykj−s1,kj+s2,kj=b1, s1,kj≥0,s2,kj≥0., and R2jxk∗=minzkj∈Z,s3,kj,s4,kj1Ts3,kj+1Ts4,kj subject to B2xk∗,zkj−s3,kj+s4,kj=b2,, s3,kj≥0,s4,kj≥0. s1,kj, s2,kj, s3,kj, and s4,kj are the introduced slack continuous variables.Step 2.
Thejth round of iteration for solving the slave problems. From the above, we can find that, if ykj and zkj are given, the above two optimization problems can be both reformulated as the linear programming problem and their optimal values, which are separately defined as G1jxk∗,ykj and Gj2xk∗,zkj, can be determined analytically as G1jxk∗,ykj=1TB1xk∗,ykj−b1 and Gj2xk∗,zkj=1TB2xk∗,zkj−b2. Thus, we define y~kj,z~kj to be the linear approximation point at the jth round of SLP for the above two optimization problems. Let yk−1∗,zk−1∗ be equal to the initial solution of y~kj,z~kj : y~k1,z~k1 and we stipulate y~kj+1=y~kj+d^kj and z~kj+1=z~kj+d¯kj where d^kj and d¯kj are separately determined by solving the following two convex quadratic programming problems (d^kj and d¯kj are separately strictly decreasing direction for the value of G1jxk∗,y~kj and Gj2xk∗,z~kj as j increases; the components of vector Δ are all very small positive real values):
(38)maxd^kjd^kjd^kjT(39)subjecttoB1xk∗,y~kjμ−b1μ⋅∇B1xk∗,y~kjd^kjμ<0,∀μ=1,2,⋯,6⋅S∇B1xk∗,y~kjd^kj<B1xk∗,y~kj−b1,−Δ≤d^kj≤Δ,y~kj+d^kj∈Y,(40)maxd¯kjd¯kjd¯kjT(41)subjecttoB2xk∗,z~kjμ−b2μ⋅∇B2xk∗,z~kjd¯kjμ<0,∀μ=1,2,⋯,9⋅S∇B2xk∗,z~kjd¯kj<B2xk∗,z~kj−b2,−Δ≤d¯kj≤Δ,z~kj+d¯kj∈Z.Step 3.
Check the inner level convergence. At first, letL1I−level=G1jxk∗,y~kj and L2I−level=G2jxk∗,z~kj. If U1I−level=0 and U2I−level=0 or U1I−level−L1I−level<δ and U2I−level−L2I−level<δ, then stop and output the current solutions y~kj, z~kj, G1jxk∗,y~kj, and G2jxk∗,z~kj as yk∗, zk∗, G1xk∗,yk∗, and G2xk∗,zk∗ into the master problem; B1xk+1,yk∗−b1 and B2xk+1,zk∗−b2 are used to construct the associated Benders cuts for the master problem. Otherwise, let U1I−level=G1jxk∗,y~kj and U2I−level=G2jxk∗,z~kj and further set j=j+1 and go to Step 2 of the SLP algorithm.


## 4. Numerical Test Results

### 4.1. Numerical Test Design and Test Data

We present a numerical study to evaluate the performance of the proposed PMP formulation and two-level GBD-based algorithm. We test on a revised power system operated on a real civil aircraft, and the revision is mainly aimed at making the power system compatible to the more electric environment for the civil aircraft. We select a two-engine and single-aisle civil aircraft which is the currently main civil aircraft type as an example. We define that the average flight duration from the lift up to touchdown for a civil aircraft is two hours (in the real environment, a typical and also an economical flight duration for this kind of aircraft is around two hours) and the smallest scheduling snapshot for the PMP is one minute, then the PMP will have 120 scheduling periods. The tests require the following data of the aircraft.


#### 4.1.1. The Flight Plan and Aerodynamic Data

The flight plan information which is constructed by the flight phases is listed in Tables2 and 3. Based on Tables 2 and 3, we can obtain the values hs,Vs,θs for period s=1,2,⋯,120 and thus obtain p0s, ρ0s, and T0s. The aircraft aerodynamic data are presented in Table 4. From Tables 2–4, we can obtain Fs and Les for s=1,2,⋯,120.Table 2
Time period allocation for the flight phase.
1
2
3
4
5
6
7
8
9
10
11
12
13
Flight phase
Initial climb
Stage 1 for climb
Stage 2 for climb
Stage 3 for climb
Stage 1 for cruise
Stage 2 for cruise
Stage 3 for cruise
Stage 1 for descent
Stage 2 for descent
Stage 3 for descent
Stage 4 for descent
Stage 5 for descent
Approach and landing
Time period no.
1
2-3
4-13
14-18
19-36
37-72
73-90
91-95
96-105
106
107-112
113-116
117-120Table 3
Data for the flight phase.
V0 (km/h)
V1 (km/h)
V2 (km/h)
V3 (km/h)
V4 (km/h)
Cruise altitude (km)
Cruise speed (km/h)
375
465
600
300
250
10
850Table 4
Aircraft aerodynamic data.
Ae (m2)
AW (m2)
g (m/s2)
M (kg)
cD0
cD1
Xwg/Xwg+Xeg
25
125
9.8
60000
0.026
0.043
10%


#### 4.1.2. The Engine Data

These data which include all technical and performance parameters of the aircraft engine are demonstrated in Table5. Here, to note the engine mentioned here is the integrated “larger” one by the real two engines of the aircraft. Based on Table 5, we can obtain Mas, p1s, p2s, T1s, and T2s for s=1,2,⋯,120.Table 5
Aircraft engine data.
(a)
R (kJ/(kg·K))
gc
γi
γc
γt
fs
πi
hpr (kJ/kg)
ρf (g/cm3)
0.286
1
1.4
1.36
1.3
0.026
1.2
43000
0.7
(b)
ηc
ηb
ηt
ηm
cpc (kJ/(kg·K))
cpt (kJ/(kg·K))
Kc (K/(rpm)2)
P¯g (kw)
0.87
0.95
0.9
0.9
1.079
1.234
11.52
500
(c)
p¯3 (kPa)
T¯3 (K)
p¯4 (kPa)
T¯4 (K)
A¯i (m2)
A¯i (m2)
p¯5 (kPa)
T¯5 (K)
n¯ (rpm)
700
1500
700
1500
5
2.5
500
1500
15000


#### 4.1.3. The ECS Data

These data which include all technical and performance parameters of the aircraft ECS are demonstrated in Table6. From Table 6, we can obtain ṁ0es and ṁ1es for s=1,2,⋯,120. Meanwhile, based on the values hs,Vs,θs, we can also obtain the values of Q1vs, Q2vs, Q3vs, and Q4vs for s=1,2,⋯,120 by the heat payload analysis of the aircraft environment control theory [26]. In this paper, we assume Q3vs and Q4vs are all fixed parameters under the different time periods and Q1vs and Q2vs are linearly and time period varying separately with Vs and hs. The associated parameters are given in Table 7.Table 6
Aircraft ECS data.
(a)
Dv (m3)
γce
γte
A0e (m2)
A1e (m2)
cpe (kJ/(kg·K))
p~ (kPa)
200
1.4
1.3
0.002
0.0001
1.004
80.8
(b)
ηc1e
ηc2e
ηte
ηme
Kc1e (K/(rpm)2)
Kc2e (K/(rpm)2)
0.9
0.9
0.9
0.8
51.84
51.84
(c)
T¯1e (K)
T¯2e (K)
p¯2e (kPa)
T¯3e (K)
p¯3e (kPa)
T¯4e (K)
p¯4e (kPa)
T¯5e (K)
p¯5e (kPa)
500
800
500
800
500
800
500
800
500
(d)
T¯6e (K)
p¯6e (kPa)
n¯c1e (rpm)
n¯c2e (rpm)
T¯v (K)
T¯v (K)
m¯v (kg)
ṁ¯out (kg/min)
500
200
6000
6000
300
293
500
100Table 7
Heat energy generated within the vessel.
Q1vs (kJ)
Q2vs (kJ)
Q3vs (kJ)
Q4vs (kJ)
‐2000⋅Vs/950
1920⋅hs/12
1284
1200


#### 4.1.4. The Battery and Generator Data

These data which include the fuel and battery operation cost information, the power output upper and lower bounds for the battery, and the state-of-charge upper and lower bounds of the battery are all presented in Table8.Table 8
Aircraft battery and generator data.
ag ($/(kg/s))
hb ($)
P¯bi (kw)
P¯bi (kw)
P¯bo (kw)
P¯bo (kw)
E¯b (kJ)
E¯b (kJ)
60
1.5
200
0
200
0
15000
0


#### 4.1.5. The Other Electrical Load Data

These data include the information corresponding to the electrical loads driven by the flight control system, electrical anti-ice protection system, and the commercial and avionics systems.Kleak, the coefficient addressing the relationship between the air-dynamical pressure on elevators and the required hydraulic oil mass flow rate, is set to be 0.01 (m·s). Pws and Pes, separately the rate power values of the electrothermal heater in the wings and elevators, are set to be 52.5 kW and 12.5 kW. The electrothermal heaters in the wings and elevators are iteratively activated for every 3 minutes after the aircraft’s altitude is higher than the minimum altitude with icing condition, 3000 m [30]. Finally, P5s, the fixed power load for commercial and avionics services in period s=1,2,⋯,120, are all set to be 50 kW.The proposed two-level GBD-based decomposition algorithm for the PMP formulation is implemented in Gurobi 7.0 [38] where the mixed integer and linear program in the master problem is solved by the branch and bound method on a personal computer laptop with a 2.50 GHz CPU and 4 GB RAM. We set the convergence tolerance for the outer level GBD algorithm to be ζ=0.001 and the convergence tolerance for the inner level SLP algorithm to be δ=0.0001. The MILP gap for the GBD master problem is set to 0.0001. We want to analyze the performance characteristics of the proposed PMP formulation and decomposition method on their computational performance and economic efficiency.


### 4.2. Computational Performance and Economic Efficiency

At first, we will collect the cost, computation time, power generation, and load information under the data of the standard cruise altitude (10 km), cruise speed (850 km/h), gravity center positionXwg/Xwg+Xeg (10%), and mass of the aircraft M (60000 kg) which are demonstrated in Tables 3 and 4. The computed fuel cost and battery operation cost are separately 55491$ and 136.5$; the computation time is 1685 sec. The generator power output Pgs and the battery power charging/discharging power Pbs for s=1,2,⋯,120, are given in Figure 6. The power load for engine fuel pumping P1s, the power load for ECS compressor spinning P2s, the power load for flight control actuation P3s, the power load for electrical anti-icing P4s, and the power load for commercial and avionics services P5s for s=1,2,⋯,120, are depicted in Figure 7.Figure 6
Generator and battery power generation results under standard cruise altitude speed, mass, and gravity center position.Figure 7
All electrical load results under standard cruise altitude, speed, mass, and gravity center position.Besides the above results, the performance with respect to the economic efficiency and computational time is also compared under different values of cruise altitude and cruise speed, gravity center position, and mass of the aircraft. (1) With standard aircraft mass (60000 kg) and gravity center position (10%), we present the total cost (fuel cost + battery operation cost) analysis under different combinations of cruise altitude and speed; the associated results are given in Figure8 where we can find that the total cost is nearly positive linear-correlated with the cruise speed when the cruise altitude is fixed; furthermore, the total cost will decrease as the cruise altitude increases under the same cruise speed. (2) With standard cruise altitude (10 km), we address the “cruise speed and total cost” analysis under different combinations of aircraft mass and gravity center position; the associated results are demonstrated in Figure 9 where “LM,” “SM,” “FC,” and “BC” separately represent large mass (75000 kg), small mass (45000 kg), forward gravity center position (5%), and backward gravity center position (15%), respectively. Although the curves in Figure 9 do not demonstrate a unified curve shape, they all reveal a rising tendency of the total cost as the cruise speed increases. (3) With standard cruise speed (850 km/h), we give the “cruise altitude and total cost” analysis under different combinations of aircraft mass and gravity center position; the associated results are shown in Figure 10 where we can find that the total cost is nearly negatively linear-correlated with the cruise altitude when the cruise speed is fixed under the different combinations of the aircraft mass and gravity center position. Moreover, we want to point out that in Figures 9 and 10, under the scenarios of forward gravity center position, the curve shapes are obviously different with the others. This is because a forward gravity center position usually associates with a relatively small power demand P3s by the primary flight control actuators. Since P3s is one of main part of the total controllable electric load and is relatively fixed within the whole flight profile, a small value of P3s will make the total electric load more variable (in Figure 7) and thus generate the more “volatile” total cost results under different cruise speeds (in Figure 9) or cruise altitudes (in Figure 10). Finally, the average computation time of the above three sets of test samples is 2336 sec.Figure 8
Total cost analysis under different combinations of cruise altitude and speed with standard aircraft mass and gravity center position.Figure 9
“Cruise speed and total cost” analysis under different combinations of aircraft mass and gravity center position with standard cruise altitude.Figure 10
“Cruise altitude and total cost” analysis under different combinations of aircraft mass and gravity center position with standard cruise speed.


## 4.1. Numerical Test Design and Test Data

We present a numerical study to evaluate the performance of the proposed PMP formulation and two-level GBD-based algorithm. We test on a revised power system operated on a real civil aircraft, and the revision is mainly aimed at making the power system compatible to the more electric environment for the civil aircraft. We select a two-engine and single-aisle civil aircraft which is the currently main civil aircraft type as an example. We define that the average flight duration from the lift up to touchdown for a civil aircraft is two hours (in the real environment, a typical and also an economical flight duration for this kind of aircraft is around two hours) and the smallest scheduling snapshot for the PMP is one minute, then the PMP will have 120 scheduling periods. The tests require the following data of the aircraft.


### 4.1.1. The Flight Plan and Aerodynamic Data

The flight plan information which is constructed by the flight phases is listed in Tables2 and 3. Based on Tables 2 and 3, we can obtain the values hs,Vs,θs for period s=1,2,⋯,120 and thus obtain p0s, ρ0s, and T0s. The aircraft aerodynamic data are presented in Table 4. From Tables 2–4, we can obtain Fs and Les for s=1,2,⋯,120.Table 2
Time period allocation for the flight phase.
1
2
3
4
5
6
7
8
9
10
11
12
13
Flight phase
Initial climb
Stage 1 for climb
Stage 2 for climb
Stage 3 for climb
Stage 1 for cruise
Stage 2 for cruise
Stage 3 for cruise
Stage 1 for descent
Stage 2 for descent
Stage 3 for descent
Stage 4 for descent
Stage 5 for descent
Approach and landing
Time period no.
1
2-3
4-13
14-18
19-36
37-72
73-90
91-95
96-105
106
107-112
113-116
117-120Table 3
Data for the flight phase.
V0 (km/h)
V1 (km/h)
V2 (km/h)
V3 (km/h)
V4 (km/h)
Cruise altitude (km)
Cruise speed (km/h)
375
465
600
300
250
10
850Table 4
Aircraft aerodynamic data.
Ae (m2)
AW (m2)
g (m/s2)
M (kg)
cD0
cD1
Xwg/Xwg+Xeg
25
125
9.8
60000
0.026
0.043
10%


### 4.1.2. The Engine Data

These data which include all technical and performance parameters of the aircraft engine are demonstrated in Table5. Here, to note the engine mentioned here is the integrated “larger” one by the real two engines of the aircraft. Based on Table 5, we can obtain Mas, p1s, p2s, T1s, and T2s for s=1,2,⋯,120.Table 5
Aircraft engine data.
(a)
R (kJ/(kg·K))
gc
γi
γc
γt
fs
πi
hpr (kJ/kg)
ρf (g/cm3)
0.286
1
1.4
1.36
1.3
0.026
1.2
43000
0.7
(b)
ηc
ηb
ηt
ηm
cpc (kJ/(kg·K))
cpt (kJ/(kg·K))
Kc (K/(rpm)2)
P¯g (kw)
0.87
0.95
0.9
0.9
1.079
1.234
11.52
500
(c)
p¯3 (kPa)
T¯3 (K)
p¯4 (kPa)
T¯4 (K)
A¯i (m2)
A¯i (m2)
p¯5 (kPa)
T¯5 (K)
n¯ (rpm)
700
1500
700
1500
5
2.5
500
1500
15000


### 4.1.3. The ECS Data

These data which include all technical and performance parameters of the aircraft ECS are demonstrated in Table6. From Table 6, we can obtain ṁ0es and ṁ1es for s=1,2,⋯,120. Meanwhile, based on the values hs,Vs,θs, we can also obtain the values of Q1vs, Q2vs, Q3vs, and Q4vs for s=1,2,⋯,120 by the heat payload analysis of the aircraft environment control theory [26]. In this paper, we assume Q3vs and Q4vs are all fixed parameters under the different time periods and Q1vs and Q2vs are linearly and time period varying separately with Vs and hs. The associated parameters are given in Table 7.Table 6
Aircraft ECS data.
(a)
Dv (m3)
γce
γte
A0e (m2)
A1e (m2)
cpe (kJ/(kg·K))
p~ (kPa)
200
1.4
1.3
0.002
0.0001
1.004
80.8
(b)
ηc1e
ηc2e
ηte
ηme
Kc1e (K/(rpm)2)
Kc2e (K/(rpm)2)
0.9
0.9
0.9
0.8
51.84
51.84
(c)
T¯1e (K)
T¯2e (K)
p¯2e (kPa)
T¯3e (K)
p¯3e (kPa)
T¯4e (K)
p¯4e (kPa)
T¯5e (K)
p¯5e (kPa)
500
800
500
800
500
800
500
800
500
(d)
T¯6e (K)
p¯6e (kPa)
n¯c1e (rpm)
n¯c2e (rpm)
T¯v (K)
T¯v (K)
m¯v (kg)
ṁ¯out (kg/min)
500
200
6000
6000
300
293
500
100Table 7
Heat energy generated within the vessel.
Q1vs (kJ)
Q2vs (kJ)
Q3vs (kJ)
Q4vs (kJ)
‐2000⋅Vs/950
1920⋅hs/12
1284
1200


### 4.1.4. The Battery and Generator Data

These data which include the fuel and battery operation cost information, the power output upper and lower bounds for the battery, and the state-of-charge upper and lower bounds of the battery are all presented in Table8.Table 8
Aircraft battery and generator data.
ag ($/(kg/s))
hb ($)
P¯bi (kw)
P¯bi (kw)
P¯bo (kw)
P¯bo (kw)
E¯b (kJ)
E¯b (kJ)
60
1.5
200
0
200
0
15000
0


### 4.1.5. The Other Electrical Load Data

These data include the information corresponding to the electrical loads driven by the flight control system, electrical anti-ice protection system, and the commercial and avionics systems.Kleak, the coefficient addressing the relationship between the air-dynamical pressure on elevators and the required hydraulic oil mass flow rate, is set to be 0.01 (m·s). Pws and Pes, separately the rate power values of the electrothermal heater in the wings and elevators, are set to be 52.5 kW and 12.5 kW. The electrothermal heaters in the wings and elevators are iteratively activated for every 3 minutes after the aircraft’s altitude is higher than the minimum altitude with icing condition, 3000 m [30]. Finally, P5s, the fixed power load for commercial and avionics services in period s=1,2,⋯,120, are all set to be 50 kW.The proposed two-level GBD-based decomposition algorithm for the PMP formulation is implemented in Gurobi 7.0 [38] where the mixed integer and linear program in the master problem is solved by the branch and bound method on a personal computer laptop with a 2.50 GHz CPU and 4 GB RAM. We set the convergence tolerance for the outer level GBD algorithm to be ζ=0.001 and the convergence tolerance for the inner level SLP algorithm to be δ=0.0001. The MILP gap for the GBD master problem is set to 0.0001. We want to analyze the performance characteristics of the proposed PMP formulation and decomposition method on their computational performance and economic efficiency.


## 4.1.1. The Flight Plan and Aerodynamic Data

The flight plan information which is constructed by the flight phases is listed in Tables2 and 3. Based on Tables 2 and 3, we can obtain the values hs,Vs,θs for period s=1,2,⋯,120 and thus obtain p0s, ρ0s, and T0s. The aircraft aerodynamic data are presented in Table 4. From Tables 2–4, we can obtain Fs and Les for s=1,2,⋯,120.Table 2
Time period allocation for the flight phase.
1
2
3
4
5
6
7
8
9
10
11
12
13
Flight phase
Initial climb
Stage 1 for climb
Stage 2 for climb
Stage 3 for climb
Stage 1 for cruise
Stage 2 for cruise
Stage 3 for cruise
Stage 1 for descent
Stage 2 for descent
Stage 3 for descent
Stage 4 for descent
Stage 5 for descent
Approach and landing
Time period no.
1
2-3
4-13
14-18
19-36
37-72
73-90
91-95
96-105
106
107-112
113-116
117-120Table 3
Data for the flight phase.
V0 (km/h)
V1 (km/h)
V2 (km/h)
V3 (km/h)
V4 (km/h)
Cruise altitude (km)
Cruise speed (km/h)
375
465
600
300
250
10
850Table 4
Aircraft aerodynamic data.
Ae (m2)
AW (m2)
g (m/s2)
M (kg)
cD0
cD1
Xwg/Xwg+Xeg
25
125
9.8
60000
0.026
0.043
10%


## 4.1.2. The Engine Data

These data which include all technical and performance parameters of the aircraft engine are demonstrated in Table5. Here, to note the engine mentioned here is the integrated “larger” one by the real two engines of the aircraft. Based on Table 5, we can obtain Mas, p1s, p2s, T1s, and T2s for s=1,2,⋯,120.Table 5
Aircraft engine data.
(a)
R (kJ/(kg·K))
gc
γi
γc
γt
fs
πi
hpr (kJ/kg)
ρf (g/cm3)
0.286
1
1.4
1.36
1.3
0.026
1.2
43000
0.7
(b)
ηc
ηb
ηt
ηm
cpc (kJ/(kg·K))
cpt (kJ/(kg·K))
Kc (K/(rpm)2)
P¯g (kw)
0.87
0.95
0.9
0.9
1.079
1.234
11.52
500
(c)
p¯3 (kPa)
T¯3 (K)
p¯4 (kPa)
T¯4 (K)
A¯i (m2)
A¯i (m2)
p¯5 (kPa)
T¯5 (K)
n¯ (rpm)
700
1500
700
1500
5
2.5
500
1500
15000


## 4.1.3. The ECS Data

These data which include all technical and performance parameters of the aircraft ECS are demonstrated in Table6. From Table 6, we can obtain ṁ0es and ṁ1es for s=1,2,⋯,120. Meanwhile, based on the values hs,Vs,θs, we can also obtain the values of Q1vs, Q2vs, Q3vs, and Q4vs for s=1,2,⋯,120 by the heat payload analysis of the aircraft environment control theory [26]. In this paper, we assume Q3vs and Q4vs are all fixed parameters under the different time periods and Q1vs and Q2vs are linearly and time period varying separately with Vs and hs. The associated parameters are given in Table 7.Table 6
Aircraft ECS data.
(a)
Dv (m3)
γce
γte
A0e (m2)
A1e (m2)
cpe (kJ/(kg·K))
p~ (kPa)
200
1.4
1.3
0.002
0.0001
1.004
80.8
(b)
ηc1e
ηc2e
ηte
ηme
Kc1e (K/(rpm)2)
Kc2e (K/(rpm)2)
0.9
0.9
0.9
0.8
51.84
51.84
(c)
T¯1e (K)
T¯2e (K)
p¯2e (kPa)
T¯3e (K)
p¯3e (kPa)
T¯4e (K)
p¯4e (kPa)
T¯5e (K)
p¯5e (kPa)
500
800
500
800
500
800
500
800
500
(d)
T¯6e (K)
p¯6e (kPa)
n¯c1e (rpm)
n¯c2e (rpm)
T¯v (K)
T¯v (K)
m¯v (kg)
ṁ¯out (kg/min)
500
200
6000
6000
300
293
500
100Table 7
Heat energy generated within the vessel.
Q1vs (kJ)
Q2vs (kJ)
Q3vs (kJ)
Q4vs (kJ)
‐2000⋅Vs/950
1920⋅hs/12
1284
1200


## 4.1.4. The Battery and Generator Data

These data which include the fuel and battery operation cost information, the power output upper and lower bounds for the battery, and the state-of-charge upper and lower bounds of the battery are all presented in Table8.Table 8
Aircraft battery and generator data.
ag ($/(kg/s))
hb ($)
P¯bi (kw)
P¯bi (kw)
P¯bo (kw)
P¯bo (kw)
E¯b (kJ)
E¯b (kJ)
60
1.5
200
0
200
0
15000
0


## 4.1.5. The Other Electrical Load Data

These data include the information corresponding to the electrical loads driven by the flight control system, electrical anti-ice protection system, and the commercial and avionics systems.Kleak, the coefficient addressing the relationship between the air-dynamical pressure on elevators and the required hydraulic oil mass flow rate, is set to be 0.01 (m·s). Pws and Pes, separately the rate power values of the electrothermal heater in the wings and elevators, are set to be 52.5 kW and 12.5 kW. The electrothermal heaters in the wings and elevators are iteratively activated for every 3 minutes after the aircraft’s altitude is higher than the minimum altitude with icing condition, 3000 m [30]. Finally, P5s, the fixed power load for commercial and avionics services in period s=1,2,⋯,120, are all set to be 50 kW.The proposed two-level GBD-based decomposition algorithm for the PMP formulation is implemented in Gurobi 7.0 [38] where the mixed integer and linear program in the master problem is solved by the branch and bound method on a personal computer laptop with a 2.50 GHz CPU and 4 GB RAM. We set the convergence tolerance for the outer level GBD algorithm to be ζ=0.001 and the convergence tolerance for the inner level SLP algorithm to be δ=0.0001. The MILP gap for the GBD master problem is set to 0.0001. We want to analyze the performance characteristics of the proposed PMP formulation and decomposition method on their computational performance and economic efficiency.


## 4.2. Computational Performance and Economic Efficiency

At first, we will collect the cost, computation time, power generation, and load information under the data of the standard cruise altitude (10 km), cruise speed (850 km/h), gravity center positionXwg/Xwg+Xeg (10%), and mass of the aircraft M (60000 kg) which are demonstrated in Tables 3 and 4. The computed fuel cost and battery operation cost are separately 55491$ and 136.5$; the computation time is 1685 sec. The generator power output Pgs and the battery power charging/discharging power Pbs for s=1,2,⋯,120, are given in Figure 6. The power load for engine fuel pumping P1s, the power load for ECS compressor spinning P2s, the power load for flight control actuation P3s, the power load for electrical anti-icing P4s, and the power load for commercial and avionics services P5s for s=1,2,⋯,120, are depicted in Figure 7.Figure 6
Generator and battery power generation results under standard cruise altitude speed, mass, and gravity center position.Figure 7
All electrical load results under standard cruise altitude, speed, mass, and gravity center position.Besides the above results, the performance with respect to the economic efficiency and computational time is also compared under different values of cruise altitude and cruise speed, gravity center position, and mass of the aircraft. (1) With standard aircraft mass (60000 kg) and gravity center position (10%), we present the total cost (fuel cost + battery operation cost) analysis under different combinations of cruise altitude and speed; the associated results are given in Figure8 where we can find that the total cost is nearly positive linear-correlated with the cruise speed when the cruise altitude is fixed; furthermore, the total cost will decrease as the cruise altitude increases under the same cruise speed. (2) With standard cruise altitude (10 km), we address the “cruise speed and total cost” analysis under different combinations of aircraft mass and gravity center position; the associated results are demonstrated in Figure 9 where “LM,” “SM,” “FC,” and “BC” separately represent large mass (75000 kg), small mass (45000 kg), forward gravity center position (5%), and backward gravity center position (15%), respectively. Although the curves in Figure 9 do not demonstrate a unified curve shape, they all reveal a rising tendency of the total cost as the cruise speed increases. (3) With standard cruise speed (850 km/h), we give the “cruise altitude and total cost” analysis under different combinations of aircraft mass and gravity center position; the associated results are shown in Figure 10 where we can find that the total cost is nearly negatively linear-correlated with the cruise altitude when the cruise speed is fixed under the different combinations of the aircraft mass and gravity center position. Moreover, we want to point out that in Figures 9 and 10, under the scenarios of forward gravity center position, the curve shapes are obviously different with the others. This is because a forward gravity center position usually associates with a relatively small power demand P3s by the primary flight control actuators. Since P3s is one of main part of the total controllable electric load and is relatively fixed within the whole flight profile, a small value of P3s will make the total electric load more variable (in Figure 7) and thus generate the more “volatile” total cost results under different cruise speeds (in Figure 9) or cruise altitudes (in Figure 10). Finally, the average computation time of the above three sets of test samples is 2336 sec.Figure 8
Total cost analysis under different combinations of cruise altitude and speed with standard aircraft mass and gravity center position.Figure 9
“Cruise speed and total cost” analysis under different combinations of aircraft mass and gravity center position with standard cruise altitude.Figure 10
“Cruise altitude and total cost” analysis under different combinations of aircraft mass and gravity center position with standard cruise speed.


## 5. Conclusions

In this paper, to cope with the PMP for the civil MEA efficiently, we at first present a clear modeling of all safety-critical electrical loads and formulate the PMP of the MEA power system as a nonlinear constrained optimization problem which encapsulates a smart codispatch between the gas turbine generator, storage devices, and all controllable electrical loads with an objective for a minimal fuel cost and battery operation cost. Then, we provide a solution algorithm to solve the above problem for the real application on a civil MEA power system through a two-level Generalized Benders Decomposition method. Finally, we present the numerical testing results based on the real civil aircraft case so as to demonstrate the effectiveness of the proposed method in terms of its operation efficiency.



---

*Source: 8474375-2020-01-24.xml*