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.
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 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.
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.
We define the average flight duration from the lift-up to touchdown for a civil aircraft is minutes and the smallest scheduling snapshot for the PMP is one minute, then the PMP will have 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 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 , , are separately represented as and and and , while the pitch angle within this period is addressed by . We define that and are separately the average flight altitude and speed for the aircraft within the period and are equal to and . Moreover, within the period , the air’s average static pressure and average density at altitude can be obtained based on the transformation relationship between these values and [22]. The air’s average static temperature at altitude can thus be determined by the ideal gas equation of state: where is the gas constant of the air. Thus, if the flight phases for all time periods are previously given, the parameters , , 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 attitude 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 . 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 are given, the Mach number is ; at the altitude , the air’s average total pressure is and the average total temperature is . The air mass flow rate into the engine inlet in period is where is the cross-sectional area of the engine inlet in period 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 and . 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. 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. 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 compressor. is the associated transfer coefficient for the engine compressor.
Since the compressor and turbine are mechanically coaxial, (7) describes the reduced amount of the total enthalpy for the turbine in period 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 . is the mechanical efficiency.
Moreover, if we define the required thrust in period to be , its value can be determined through the flight mechanics analysis as shown in Appendix A; thus, is represented by the parameters . 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, is also required to be no less than a minimum cross-sectional area , 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 is equal to where and 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 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 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, , , and can be determined if the values are given for period . Moreover, since the cross-sectional area of the inlet for ECS compressor 1 is usually taken as a fixed value, the mass flow rate of the imported air in period 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. 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. 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
In (13) and (14), the mass flow rate of the cold side air into heat exchanger 1 in period , , can be stated as where is the cross-sectional area of the inlet at the cold side of heat transfer 1.
Eq. (16) represents the isentropic change characteristic of the air state variables within the ECS compressor 2. is the operation efficiency parameter of ECS compressor 2.
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 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. 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. 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 for the ECS turbine is equal to the increased amount of the total enthalpy for ECS compressor 2 in period . 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 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 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 ( is the pressure value over the sea level by 1800 m) and the minimum comfortable temperature . 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 is equal to . for 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 . We define the required air-dynamical force on elevators in period to be , and its value can thus be determined through the flight mechanics analysis shown in Appendix B; in other words, can be represented by the parameters . Finally, we have known that the electrical load required by the flight control system in period is where is the leak coefficient addressing the relationship between the air-dynamical force on elevators and required hydraulic oil mass flow rate [29]. is the density of the elevators’ hydraulic oil and 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 to be separately equal to and which are activated iteratively in every 3 minutes after the aircraft’s altitude becomes higher than the minimum altitude with icing condition. and 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:
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), is equal to where is the linear coefficient between the fuel mass flow rate and fuel cost; is equal to where 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
Equation (26) stipulates that in every time period, the state of the battery will be one of the three cases: charging (the discrete variable is 1), discharging (the discrete variable is 1), and shut-off ( and are both 0); the charge/discharge power in period 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
In (28), is defined to be the fixed power load for commercial and avionics services in period . 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 Section 2.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
Equation (29) addresses that the power output by the generator in period should be positive and less than the rated power capacity .
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:
The variable is a vector of the decisions including , , , , , , , , , , , , , , , , , and for each time period () of the whole PMP scheduling horizon. The variables and 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 variable based on (25), the objective function can be stated as the linear expression in (30). The linear constraints in (31), involving only the variable , 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 , include (12) and (28). The nonlinear constraints (33) and (34), which couple the variables with and with , separately include (2), (4), (5), (6), (7), (8), (11), (16), (17), (18), (19), (20), (21), (22), and (23). The compact sets , , and 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 as the complicating variables under the GBD architecture, (33) and (34) will separately have a convex form for and with given . 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 . 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.
3.2. Algorithm Description
The description of the proposed algorithm is given as follows.
3.2.1. Outer Level: GBD Algorithm
(1) Initialization. Let be the optimal solution for the optimization problem: subject to , , . Solve the two inner problems based on (the solution methodology is given after in this subsection) and let and be separately the optimal solutions of the two inner problems based on . Set the outer level lower bound , upper bound , and the number of iteration . Choose an outer level convergence tolerance level (>0).
(2) Iteration .
Step 1. Solve the GBD master problem. The master problem is the following MILP program: subject to , , , , and . and are the introduced slack variables which have positive values, and is a column vector where all of its elements are ; and are determined in the (-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 and is nonlinear corresponding to ; moreover, has also a nonlinear form for . 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 be the optimal solution for the th round of master problem. We set .
Step 2. Check the outer level convergence. If , stop and return . Otherwise, set 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 and be separately the optimal solutions of the two inner problems based on and input them as parameters into ()th round of the master problem. Let , 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 () passed from the th iteration of the outer level master problem. Set the two inner level lower bounds , upper bounds , and number of iteration . Choose the inner level convergence tolerance level (>0).
(2) Iteration .
Step 1. We at first present the two slave problems as follows. SP I: find a feasible solution of within , ; SP II: find a feasible solution of within , . We can equivalently restate the problems SP I and SP II separately as the following optimization problems: subject to , , and subject to , . , , , and are the introduced slack continuous variables.
Step 2. The th round of iteration for solving the slave problems. From the above, we can find that, if and 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 and , can be determined analytically as and . Thus, we define to be the linear approximation point at the th round of SLP for the above two optimization problems. Let be equal to the initial solution of : and we stipulate and where and are separately determined by solving the following two convex quadratic programming problems ( and are separately strictly decreasing direction for the value of and as increases; the components of vector are all very small positive real values):
Step 3. Check the inner level convergence. At first, let and . If and or and , then stop and output the current solutions , , , and as , , , and into the master problem; and are used to construct the associated Benders cuts for the master problem. Otherwise, let and and further set 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 Tables 2 and 3. Based on Tables 2 and 3, we can obtain the values for period and thus obtain , , and . The aircraft aerodynamic data are presented in Table 4. From Tables 2–4, we can obtain and for .
4.1.2. The Engine Data
These data which include all technical and performance parameters of the aircraft engine are demonstrated in 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 , , , , and for .
4.1.3. The ECS Data
These data which include all technical and performance parameters of the aircraft ECS are demonstrated in Table 6. From Table 6, we can obtain and for . Meanwhile, based on the values , we can also obtain the values of , , , and for by the heat payload analysis of the aircraft environment control theory [26]. In this paper, we assume and are all fixed parameters under the different time periods and and are linearly and time period varying separately with and . 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. , 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). and , 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, , the fixed power load for commercial and avionics services in period , 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 and the convergence tolerance for the inner level SLP algorithm to be . 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 position (10%), and mass of the aircraft (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 and the battery power charging/discharging power for , are given in Figure 6. The power load for engine fuel pumping , the power load for ECS compressor spinning , the power load for flight control actuation , the power load for electrical anti-icing , and the power load for commercial and avionics services for , are depicted in Figure 7.
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 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 by the primary flight control actuators. Since is one of main part of the total controllable electric load and is relatively fixed within the whole flight profile, a small value of 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.
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.
Appendix
A. The Aircraft Force Balance Analysis
Here, we provide the approach to determine the required thrust value in period based on the given values . From the force balance analysis in Figure 11(a), we have where is the aircraft total lift force and is the total drag force. is the mass of the aircraft and is the acceleration of gravity. Based on the flight mechanics theory, we also have where is the area of the wings and is the total lift coefficient. and are the transfer coefficients from to .
(a) Force balance
(b) Pitch moment balance
From equations (A.1) and (A.2), we can thus have
B. The Aircraft Pitch Moment Balance Analysis
The aircraft total lift force in Figure 11(a) is the sum of the lift forces created by the wings and the elevators which are demonstrated in Figure 11(b). We thus have where and are separately the fixed lift coefficient with the wings and the elevators. is the elevators’ area, and and are separately the aircraft’s angle of attack and elevators’ deflection angle.
From equations (A.1) and (A.2), we thus have
Moreover, from the pitch moment balance analysis in Figure 11(b), is equal to , and from equation (B.1), we can have
In equation (B.3), is the distance between the aerodynamic center on the wing and the aircraft’s center of gravity ; is the distance between the aerodynamic center on the elevator’s and the aircraft’s center of gravity .
Through solving the above two linear equations with two unknowns (, ) represented by equations (B.2) and (B.3), we can obtain the values of and . Finally, from equation (B.1), we have
Parameters
| : | Time period () |
| sec: | Second |
| rpm: | Revolutions per minute |
| : | Very large positive number |
| : | Altitude value where the aircraft locates at the initial time in period (m) |
| : | Altitude value where the aircraft locates at the end time in period (m) |
| : | Speed value where the aircraft locates at the initial time in period (m/sec) |
| : | Speed value where the aircraft locates at the end time in period (m/sec) |
| : | Average speed value in period (m/sec) |
| : | Constant pitch angle in period (°) |
| : | Average altitude value in period (m) |
| : | Air’s average density before the engine and ECS inlets in period (kg/m3) |
| : | Air’s average static pressure before the engine and ECS inlets in period (Pa) |
| : | Air’s average static temperature before the engine and ECS inlets in period (K) |
| : | Air’s average total pressure before the engine and ECS inlets in period (Pa) |
| : | Air’s average total temperature before the engine and ECS inlets in period (K) |
| : | Gravitational conversion factor |
| : | Mach number in period |
| : | Ratio of specific heats in engine and ECS inlets |
| : | Ratio of specific heats in engine compressor |
| : | Ratio of specific heats in engine turbine |
| : | Ratio of fuel to air mass flow amount within the combustion chamber in period |
| : | Air’s average total pressure before the engine compressor in period (Pa) |
| : | Air’s average total temperature before the engine compressor in period (K) |
| : | Pressure ratio driven by the engine inlet |
| : | Operation efficiency of the engine compressor |
| : | Combustion efficiency of the chamber |
| : | Low heating value of the fuel (J/kg) |
| : | Fuel density before the fuel is pumped into the chamber (kg/m3) |
| : | Operation efficiency of the engine turbine |
| : | Mechanical efficiency between the engine compressor and turbine |
| : | Specific heat at constant pressure in the engine compressor (J/(kg·K)) |
| : | Specific heat at constant pressure in the engine turbine (J/(kg·K)) |
| : | Transfer coefficient between the rotation speed and the change of total temperature for the engine compressor (K/(rpm)2) |
| : | Required thrust value in period (N) |
| : | Air mass flow rate into the inlet of ECS compressor 1 in period (kg/sec) |
| : | Air mass flow rate into the cold side of ECS heat exchangers 1 and 2 in period (kg/sec) |
| : | Cross-sectional area of the inlet for ECS compressor 1 (m2) |
| : | Cross-sectional area of the inlet at the cold side of ECS heat exchangers 1 and 2 (m2) |
| : | Ratio of specific heats in ECS compressors 1 and 2 |
| : | Ratio of specific heats in ECS turbine |
| : | Operation efficiency of ECS compressor 1 |
| : | Operation efficiency of ECS compressor 2 |
| : | Operation efficiency of ECS turbine |
| : | Mechanical efficiency between ECS compressor 2 and turbine |
| : | Specific heat at constant pressure in ECS compressors 1 and 2 (J/(kg·K)) |
| : | Transfer coefficient between the rotation speed and the change of total temperature for ECS compressor 1 (K/(rpm)2) |
| : | Transfer coefficient between the rotation speed and the change of total temperature for ECS compressor 2 (K/(rpm)2) |
| : | Heat energy imported into the vessel by the heat transfer between the aircraft’s skin and the air outside in period (J) |
| : | Heat energy imported into the vessel by the sun radiation in period (J) |
| : | Heat energy imported into the vessel by the passenger/pilot activity in period (J) |
| : | Heat energy imported into the vessel by the commercial and avionics loads in period (J) |
| : | Volume of the vessel (m3) |
| : | Area of the elevators (m2) |
| : | Elevators’ aerodynamic payload in period (N) |
| : | Leak coefficient addressing the relationship between the air-dynamical force on elevators and required hydraulic oil mass flow rate () |
| : | Density of the elevators’ hydraulic oil (kg/m3) |
| : | Rated power value of the electrothermal heater on the wings in period (W) |
| : | Rated power value of the electrothermal heater on the elevators in period (W) |
| : | Coefficient of the linear relationship between the fuel mass flow rate and the fuel cost ($/(kg/s)) |
| : | Operation cost coefficient of the battery ($) |
| : | Minimum permissible charging power for the battery in every period (W) |
| : | Maximum permissible charging power for the battery in every period (W) |
| : | Minimum permissible discharging power for the battery in every period (W) |
| : | Maximum permissible discharging power for the battery in every period (W) |
| : | Lower bound of the state-of-charge energy for the battery (J) |
| : | Upper bound of the state-of-charge energy for the battery (J) |
| : | Generator’s rated power capacity (W) |
| : | column vector of the decision variables in the master problem |
| : | column vector of the decision variables in the slave problem I for the engine |
| : | column vector of the decision variables in the slave problem II for the ECS. |
| : | Air’s average total pressure before the combustion chamber in period (Pa) |
| : | Air’s average total temperature before the combustion chamber in period (K) |
| : | Air’s average total pressure before the engine turbine in period (Pa) |
| : | Air’s average total temperature before the engine turbine in period (K) |
| : | Air’s average total pressure before the nozzle in period (Pa) |
| : | Air’s average total temperature before the nozzle in period (K) |
| : | Air mass flow rate into the engine inlet in period (kg/sec) |
| : | Cross-sectional area of the engine inlet in period (m2) |
| : | Stable rotation speed of the engine compressor and turbine in period (rpm) |
| : | Generator power output in period (W) |
| : | Air’s average total temperature before ECS compressor 1 in period (K) |
| : | Air’s average total pressure before the hot side of ECS heat exchanger 1 in period (Pa) |
| : | Air’s average total temperature before the hot side of ECS heat exchanger 1 in period (K) |
| : | Air’s average total pressure before ECS compressor 2 in period (Pa) |
| : | Air’s average total temperature before ECS compressor 2 in period (K) |
| : | Air’s average total pressure before the hot side of ECS heat exchanger 2 in period (Pa) |
| : | Air’s average total temperature before the hot side of ECS heat exchanger 2 in period (K) |
| : | Air’s average total pressure before ECS turbine in period (Pa) |
| : | Air’s average total temperature before ECS turbine in period (K) |
| : | Air’s average total pressure after ECS turbine in period (Pa) |
| : | Air’s average total temperature after ECS turbine in period (K) |
| : | Temperature value within the vessel in period (K) |
| : | Air mass within the vessel in period (kg) |
| : | Air mass flow rate exhausted outside the vessel from valve A in period (kg/sec) |
| : | Stable rotation speed of ECS compressor 1 in period (rpm) |
| : | Stable rotation speed of ECS compressor 2 and turbine in period (rpm) |
| : | The charge/discharge power of the battery in period (W) |
| : | State-of-charge energy of the battery in period (J) |
| : | ; discrete variable describing the state of battery 1 if the battery is charging, 0 if the battery is discharging or shut-off in period |
| : | ; discrete variable describing the state of battery 1 if the battery is discharging, 0 if the battery is charging or shut-off in period |
| : | column vector of the parameters for the objective function |
| : | left side parameter matrix for the linear constraints with |
| : | right side column vector of the parameters for the linear constraints with |
| : | nonlinear vector-valued function expression for the constraints with and |
| : | nonlinear vector-valued function expression for the constraints with and |
| : | right side column vector of the parameters for the constraints with and |
| : | right side column vector of the parameters for the constraints with and |
| : | nonlinear vector-valued function expression for nonlinear constraints with |
| : | right side column vector of the parameters for the nonlinear constraints with . |
Data Availability
All of the test data associated with this paper has been given in Tables 2–8 of Section 4.
Conflicts of Interest
The authors declare that they have no conflicts of interest.