Power Management Problem for Civil Aircraft under More Electric Environment

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.


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][4][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][4][5][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][7][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 safetycritical 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 2 International Journal of Aerospace Engineering 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 Section 2, 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.

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. We define the average flight duration from the lift-up to touchdown for a civil aircraft is S 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. Stage 2 for climb Periods for constant flight altitude 3000 m with speed increased from V 1 to V 2 ; the pitch angle is 0°. 4 Stage 3 for climb Periods from 3000 m to cruise altitude with constant speed V 2 ; the pitch angle is 7°~10°. 5 Stage 1 for cruise Periods for keeping cruise altitude with speed increased from V 2 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 to V 2 ; the pitch angle is 0°. 8 Stage 1 for descent Periods from cruise altitude to 3000 m with speed V 2 ; the pitch angle is -10°~-7°. 9 Stage 2 for descent Periods for constant flight altitude 3000 m with speed decreased from V 2 to V 1 ; the pitch angle is 0°. 10 Stage 3 for descent Periods from 3000 m to 1800 m with constant speed V 1 ; the pitch angle is -9°. 11 Stage 4 for descent Periods for constant flight altitude 1800 m with speed decreased from V 1 to V 3 ; the pitch angle is 0°. 12 Stage 5 for descent Periods from 1800 m to 450 m with constant speed V 3 ; the pitch angle is -4°. 13 Approach and landing Periods from 450 m to touchdown with constant speed V 4 ; the pitch angle is -1.5°.  3 International Journal of Aerospace Engineering 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 Table 1. 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 period s, s = 1, 2, ⋯, S, are separately represented as h init ðsÞ and V init ðsÞ and h end ðsÞ and V end ðsÞ, while the pitch angle within this period is addressed by θðsÞ. We define that hðsÞ and VðsÞ are separately the average flight altitude and speed for the aircraft within the period s and are equal to ðh init ðsÞ + h end ðsÞÞ/2 and ðV init ðsÞ + V end ðsÞÞ/2. Moreover, within the period s, the air's average static pressure p 0 ðsÞ and average density ρ 0 ðsÞ at altitude hðsÞ can be obtained based on the transformation relationship between these values and hðsÞ [22]. The air's average static temperature T 0 ðsÞ at altitude hðsÞ can thus be determined by the ideal gas equation of state: p 0 ðsÞ = R ⋅ ρ 0 ðsÞ ⋅ T 0 ðsÞ 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 ðhðsÞ, VðsÞ, θð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.

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 shed-ding 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 attitude ðhðsÞ, VðsÞ, θð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 safetycritical 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 ðhðsÞ, VðsÞ, θð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 3 From [23], if the values ðhðsÞ, VðsÞ, θðsÞÞ are given, the Mach number MaðsÞ is VðsÞ/ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where A i ðsÞ is the cross-sectional area of the engine inlet in period s which is taken as a decision variable controlled by 4 International Journal of Aerospace Engineering the Full Authority Digital Engine Control (FADEC) computer [23]. Based on the aircraft inlet theory [12], we have p 2 ðsÞ = π i ⋅ p 1 ðsÞ and T 2 ðsÞ = T 1 ðsÞ. 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): 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.
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.
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.
Equation (6) states the relationships between the rotation speed and the change of total temperature for the engine 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 period s to be FðsÞ, its value can be determined through the flight mechanics analysis as shown in Appendix A; thus, FðsÞ is represented by the parameters ðh ðsÞ, V ðsÞ, θ ðsÞÞ. Then, based on aircraft engine theory [23], we can have Finally, since all decision variables in (2), (3), (4), (5), (6), (7), and (8) are positive and have their associated upper bounds (especially, A i ðsÞ 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 From the above equations, we can find the load demand of fuel pumping for the engine operation in period s P 1 ðsÞ is equal to ð _ m 0 ðsÞ/ρ f Þ ⋅ f ðsÞ ⋅ p 3 ðsÞ where _ m 0 ðsÞ and p 3 ðsÞ 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][25][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 4 In Figure 4, 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  International Journal of Aerospace Engineering 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, MaðsÞ, p 1 ðsÞ, and T 1 ðsÞ can be determined if the values ðhðsÞ, VðsÞ, θðsÞÞ are given for period s. Moreover, since the cross-sectional area of the inlet for ECS compressor 1 A e 0 is usually taken as a fixed value, the mass flow rate of the imported air in period s _ m e 0 ðsÞ is equal to 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]: Equation (11) represents the isentropic change characteristic of the air state variables within ECS compressor 1. η e c1 ð<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. K e c1 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 lowtemperature 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 In (13) and (14), the mass flow rate of the cold side air into heat exchanger 1 in period s, _ m e 1 ðsÞ, can be stated as where A e 1 is the cross-sectional area of the inlet at the cold side of heat transfer 1.
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.
Equation (19) represents the isentropic change characteristic of the air state variables within the ECS turbine. η e t ð<1Þ is the operation efficiency parameter of the ECS turbine.
Equation (20) states the relationships between the rotation speed and the change of total temperature for ECS compressor 2. K e c2 is defined to be the associated transfer coefficient for ECS compressor 2.
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. η e m ð<1Þ is the associated mechanical efficiency.
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.
Equation (23) represents the pressure balance constraints in vessel where D v 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 max ðp, p 1 ðsÞÞ (p is the pressure value over the sea level by 1800 m) and the minimum comfortable temperature T v . Thus, we have Based on the above equations, we can find that the load demand of the electrical DC motor for the ECS in period s P 2 ðsÞ is equal to _ (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 electrohydrostatic 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 main-tain 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 ðhðsÞ, VðsÞ, θðsÞÞ. We define the required airdynamical force on elevators in period s to be L e ðsÞ, and its value can thus be determined through the flight mechanics analysis shown in Appendix B; in other words, L e ðsÞ can be represented by the parameters ðhðsÞ, VðsÞ, θðsÞÞ. Finally, we have known that the electrical load required by the flight control system in period s P 3 ðsÞ is K leak ⋅ ðL e ðsÞ/A e Þ 2 /ρ hy where K leak 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 A e 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 antiicing the leading-edge wings and elevators in period s P 4 ðsÞ to be separately equal to P w ðsÞ and P e ðsÞ which are activated iteratively in every 3 minutes after the aircraft's altitude becomes higher than the 8 International Journal of Aerospace Engineering minimum altitude with icing condition. P w ðsÞ and P e ðsÞ 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 decisionmaking 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: 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), C g ðsÞ is equal to a g ⋅ _ m 0 ðsÞ ⋅ f ðsÞ where a g is the linear coefficient between the fuel mass flow rate and fuel cost; where h b is the operation cost coefficient. We next state each of the constraints we include in the problem formulation.

Battery Charging and Discharging Behavior Constraints.
Equation (26) stipulates that in every time period, the state of the battery will be one of the three cases: charging (the discrete variable z c b ðsÞ ∈ f0, 1g is 1), discharging (the discrete variable z d b ðsÞ ∈ f0, 1g is 1), and shut-off (z c b ðsÞ and z d b ðsÞ are both 0); the charge/discharge power in period s P b ðsÞ 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.

Load Balance Constraints.
In (28), P 5 ðsÞ 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 actu-ation, anti-ice for leading edge of wings/elevators, and commercial and avionics services.
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 GBDbased 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.

Algorithm
Description. The description of the proposed algorithm is given as follows.  x ∈ X. Solve the two inner problems based on x * 0 (the solution methodology is given after in this subsection) and let y * 0 and z * 0 be separately the optimal solutions of the two inner problems based on x * 0 . Set the outer level lower bound L O−level = −∞, upper bound U O−level = +∞, and the number of iteration k = 1. Choose an outer level convergence tolerance level ζ (>0).
Step 1. Solve the GBD master problem. The master problem is the following MILP program: min Fx k ≤ f, Hðx k Þ = h, and x k ∈ 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; y * k-1 and z * k-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 B 1 ðx k , y * k-1 Þ and B 2 ðx k , z * k-1 Þ is nonlinear corresponding to x k ; moreover, Hðx k Þ has also a nonlinear form for x k . 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 ðx * k , α * k , β * k Þ be the optimal solution for the kth round of master problem. We set Step 2. Check the outer level convergence. If U O−level − L O−level < ζ, stop and return ðx * k , y * k , z * k Þ. Otherwise, set U O−level = c T x * k + M T α * k + M T β * 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. Let y * k and z * k be separately the optimal solutions of the two inner problems based on x * k and input them as parameters into (k + 1)th round of the master problem. Let k = k + 1, go to Step 1.

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. (2) Iteration j ≥ 1.
Step 1. We at first present the two slave problems as follows.

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.  Tables 2 and 3. Based on Tables 2 and 3, we can obtain the values ðhðsÞ, VðsÞ, θðsÞÞ for period s = 1, 2, ⋯, 120 and thus obtain p 0 ðsÞ, ρ 0 ðsÞ, and T 0 ðsÞ. The aircraft aerodynamic data are presented in Table 4. From Tables 2-4, we can obtain F ðsÞ and L e ðsÞ for s = 1, 2, ⋯, 120.  Table 5. 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 MaðsÞ, p 1 ðsÞ, p 2 ðsÞ, T 1 ðsÞ, and T 2 ðsÞ for s = 1, 2, ⋯, 120.

The ECS Data.
These data which include all technical and performance parameters of the aircraft ECS are demonstrated in Table 6. From Table 6 [26]. In this paper, we assume Q v 3 ðsÞ and Q v 4 ðsÞ are all fixed parameters under the different time periods and Q v 1 ðsÞ and Q v 2 ðsÞ are linearly and time period varying separately with VðsÞ and hðsÞ. The associated parameters are given in Table 7.
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 Table 8.
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. K leak , 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). P w ðsÞ and P e ðsÞ, 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, P 5 ðsÞ, 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.

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 position X wg /ðX wg + X eg Þ (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 P g ðsÞ and the battery power charging/discharging power P b ðsÞ for s = 1, 2, ⋯, 120, are given in Figure 6. The power load for engine fuel pumping P 1 ðsÞ, the power load for ECS compressor spinning P 2 ðsÞ, the power load for flight control actuation P 3 ðsÞ, the power load for electrical anti-icing P 4 ðsÞ, and the power load for commercial and avionics services P 5 ðsÞ for s = 1, 2, ⋯, 120, are depicted in Figure 7. International Journal of Aerospace Engineering 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 Figure 8 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 Table 3: Data for the flight phase.
1500 700 1500 5 2.5 500 1500 15000      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 P 3 ðsÞ by the primary flight control actuators. Since P 3 ðsÞ is one of main part of the total controllable electric load and is relatively fixed within the whole flight profile, a small value of P 3 ðsÞ 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.

Conclusions
In this paper, to cope with the PMP for the civil MEA efficiently, we at first present a clear modeling of all safetycritical 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     14 International Journal of Aerospace Engineering 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.

A. The Aircraft Force Balance Analysis
Here, we provide the approach to determine the required thrust value FðsÞ in period s based on the given values ðhðsÞ, VðsÞ, θðsÞÞ. From the force balance analysis in Figure 11(a), we have where LðsÞ is the aircraft total lift force and DðsÞ is the total drag force. M is the mass of the aircraft and g is the acceleration of gravity. Based on the flight mechanics theory, we also have

B. The Aircraft Pitch Moment Balance Analysis
The aircraft total lift force LðsÞ in Figure 11(a) is the sum of the lift forces created by the wings L W ðsÞ and the elevators L e ðsÞ which are demonstrated in Figure 11(b). We thus have