Energy Optimization Analysis of Turbofan Engine for More-Electric Civil Aircraft

The civil aviation industry is moving into the more-electric environment where the civil aircraft uses electricity to meet the multiple energy demands of the associated aircraft subsystems. The civil aircraft’s turbofan engine, which is the largest energy supply system of civil aircraft, will thus utilize more fuel resource to provide increased electric energy besides of its conventional responsibility of maintaining the fundamental thrust requirements by aircraft. This will introduce new challenge for the energy optimization analysis of aircraft turbofan engine: it was nearly a simple optimal setting of engine’s component parameters for keeping required thrust, but under the more-electric environment it will become an optimization problem in order to minimize the fuel consumption while obeying the multiple constraints by thermodynamic limits of turbofan engine and by varying electric power demands associated with the flight profile of aircraft. We present a complete modeling in this paper for the energy optimization analysis of turbofan engine for more-electric civil aircraft and formulate it as a nonlinear programming form. We propose an algorithm based on Benders decomposition method to solve this problem; the numerical results demonstrate the economic effectiveness of the proposed modeling and algorithm.


I. INTRODUCTION
The main tasks of the turbofan engine for a conventional civil aircraft are basically to provide sufficient thrust power for keeping aircraft's safe flight and also to supply thermal energy through its ram air for the operation of aircraft's other subsystems, such as the environmental control system [1]. To realize the above tasks, the turbofan engine for a conventional civil aircraft is designed to be ''bleeding'' small parts of ram air from the fan/compressor(s) into the aircraft's other subsystems so as to achieve the required amounts of thermal energy of the subsystems, and using the rest of the ram air to produce the required thrust levels through standard Brayton thermodynamic cycle of the engine while maintaining the engine operated in an economic and stable manner [2]. However, since the way of driving the ram air from turbofan engine into aircraft's other subsystems usually introduces a large amount of energy waste, the energy efficiency of the conventional turbofan engine with air bleeding function is rel-The associate editor coordinating the review of this manuscript and approving it for publication was Shuaihu Li . atively low [3]. Thus, it is necessary to implement an energy optimization analysis for turbofan engine, which is basically a performance evaluation of engine's operation parameters based on prescribed evaluation indexes such as fuel burn, emission and exergy under different levels of thrust and ram air flow rate [4]- [6]. With considering the fact that, in most of the case in real engineering environment, the amount of energy driven by the bleed air is usually fixed (in order to keep the turbofan engine operated in a reliable mode, the temperature and pressure of the bleed air are usually set to be at high and constant values which are oversized than the level of those required by aircraft's other subsystems; i.e., the air bleeding function within the conventional turbofan engine is designed only under the ''worst case'' which is to have the largest energy utilization by aircraft's other subsystems under some particular situation(s) during flight), the conventional way of energy optimization analysis for turbofan engine is nearly a simple optimal setting of component parameters for turbofan engine with only considering keeping required thrust through some heuristic and quantitative analysis [4]- [6].
In recent years, driven by the increased environmental and economic awareness [7], the techniques on more-electric civil aircraft (MECA) has introduced large interest of the civil aviation industry [8]. The MECA will resort to more electrical power to supply thermal energy demands for the aircraft subsystems and can basically improve the energy efficiency for aircraft and reduce the associated pollution and fuel usage [9]. However, the implementation of MECA will introduce significant technique changes for civil aircraft and one of them will be for the turbofan engine. Under the more-electric environment, the turbofan engine eliminates the conventional pneumatic bleed manifold and converts the power extraction source from the ''bleed air'' to the electric power which is generated by the electrical motors driven by the axial power of turbofan engine [10]. This ''no-bleed'' engine architecture can thus provide the benefits which include the improvements the turbofan engine's efficiency of the energy usage [11]. However, the introduction of the ''no-bleed'' engine into the civil aircraft also raises new challenge for the energy optimization analysis of turbofan engine [12]: the turbofan engine will utilize additional fuel resource to produce large amount of axial power to make a balance with the increased electric energy demands besides of its conventional responsibility of maintaining the fundamental thrust requirements by aircraft. However, since the increased electric energy demands basically own a time-varying characteristic which is highly correlated with the flight profile of aircraft [13], then an efficient energy optimization analysis of turbofan engine shall be able to realize an integrated scheduling for a minimal amount of fuel consumption while maintaining the engine operated under a reliable way and the energy demand balance by the thrust and electricity. Obviously, the conventional technique of energy optimization analysis for turbofan engine, which is based on the design under the ''worst case'' mentioned above under some typical flight phases, does not consider the new characteristic of ''no-bleed'' turbofan engine and cannot be considered as an efficient way for energy optimization analysis under the more-electric environment. It is thus necessary to redefine the technique of energy optimization analysis for turbofan engine.
Previous research works have been proposed to cope with the energy optimization analysis of turbofan engine for moreelectric environment. The papers in [14], [15] separately present an evaluation method based on exergy analysis and a minimal fuel consumption to the modeling of an optimal sizing of more-electric turbofan engine (including combination of compressor/fan, combustion chamber, turbine and electric generator and its associated parameter settings), in order to realize a minimal energy consumption under the ''worst case''. These works finally verifies its performance with different levels of thrust and electric power demands. In contrast, the papers in [16]- [18] provide a power management/control model to co-optimize the energy output of the more-electric turbofan engine and electric loads. The above approaches for energy optimization analysis are novel but have the following challenges: 1) fixed electric loads provided by turbofan engine are assumed, i.e. the quantification of relationship between the electric loads by aircraft's other subsystems and flight's profile is not considered and thus the previous papers cannot fully quantify the energy optimization analysis of turbofan engine under MECA; 2) they do not consider the basic thermodynamic scheme within the turbofan engine, especially the thermodynamics constraints driven by Brayton cycle of turbofan engine [19] so as to reach a stable and economic work mode while making an efficient energy balance between the thrust requirement and electric demands. In fact, rather than just to resolve the ''optimal'' sizing of the components for turbofan engine, or rather than just a simple energy balance scheduling between the generations and loads, an efficient energy optimization analysis for turbofan engine under more-electric environment should address an integrated scheduling to minimize the fuel consumption with considering electric energy supply and demand balance and thermodynamic schematics of turbofan engine under every flight phase of its associated flight profile. Few research efforts have been made to deal with the above problem for energy optimization analysis of turbofan engine for MECA.
In order to cope with the above challenges, we present a new modeling and algorithm for energy optimization analysis of turbofan engine for MECA. The main contributions of this paper are that: 1) to give an explicit modeling of the energy optimization analysis of turbofan engine under the more-electric environment as a nonlinear programming problem with an objective for a minimal fuel consumption and with the constraints stipulated by engine thermodynamic limits and electric energy demands associated with aircraft's varying flight phases; 2) to present an efficient algorithm based on Benders decomposition scheme to solve the above problem: at first separate the formulated model as a two-layer optimization problem (upper and lower sub-problems) and then obtain the schedule results through a two-layer decomposition method; 3) the numerical analysis results based on a test case of turbofan engine of MECA validate the operation efficiency of the proposed method. Section 2 will give the complete modeling of the energy optimization analysis of turbofan engine for MECA. Section 3 will propose the proposed algorithm. Section 5 will present the numerical analysis results under a real turbofan engine environment. A conclusion will be given in Section 5.

II. THE NUMERICAL MODELING OF THE PROBLEM
The turbofan engine for civil aircraft are mainly to realize the stable thrust supply and sufficient thermal energy for aircraft's other subsystems by ram air (the subsystem which requires the largest amount of thermal energy from turbofan engine is the environmental control system) [1]. Since all of the energy used for the above functions are basically absorbed from the combustion of the fuel source, the civil aircraft usually relies a complicated hydraulic-mechanical fuel pump of fuel supply system to provide fuel into combustion chamber for the conventional turbofan engine [1]. As for the ''no-bleed'' turbofan engine for MECA [20], however, a significant change in contrast to the conventional turbofan engine of civil aircraft is that an electricity-controlled fuel pump has replaced the hydraulic-mechanical facilities to continuously drive the fuel into the combustion chamber of turbofan engine so as to keep its safe and stable operation; an electricity-controlled fuel pump will be driven by a direct current (DC) motor and the provided power load by DC motor is equal to the value of the fuel volume flow rate into the combustion chamber multiply with the chamber's inner pressure. Besides of that, since the turbofan engine for MECA eliminates the conventional pneumatic bleed manifold and converts the thermal energy supply for aircraft's other subsystems directly into the electric energy, the turbofan engine will utilize more mechanical axial power to drive its co-axis electrical motor so as to make an energy balance with the increased electric energy demand (the increased electrical demand also includes the energy provided for electricitycontrolled fuel pump of fuel supply system).
In Fig. 1, a typical system architecture of the turbofan engine for MECA is presented (here to note Fig. 1 only gives half of demonstration of the cross section for a turbofan engine which is actually symmetrical along the direction of free stream) [10], [21]. The system architecture of this ''nobleed'' turbofan engine can be divided as two modules: 1) the first module is the thermodynamic assemblies which in basic provide all of energy from fuel sources, while eliminating the conventional hydraulic and pneumatic systems. The free stream is driven into the inlet of turbofan engine and compressed by the fan (a ducted propeller); after that, the compressed air is distributed into two branches (ducts): one branch of air (within the inner duct) is used for the inner Brayton thermodynamic cycle [19] (equivalent to a two-axis turbojet), where the air is sequentially compressed by lowpressure booster and high-pressure compressor, combusted with injected fuel at the chamber and expanded by highpressure and low-pressure turbines (HPT and LPT) and core nozzle, so as to provide the basic thrust; meantime, the other branch of air (within the outer duct) distributed from the fan will be completely expanded at the bypass nozzle in order to provide an additional thrust; i.e., the total thrust is produced both by the inner and outer ducts; 2) the second module includes the electrical units of a turbofan engine and its associated interconnection facilities with the aircraft's other subsystems. The electric energy generation of the turbofan engine is provided by a generator which is mechanically and co-axially connected with the fan, the booster and LPT. The generator transmits the electricity through the aircraft's DC distribution system to three kinds of mainly electrical energy demands: 1) the DC motor of fuel supply system which drives the fuel pump continuously to inject fuel into combustion chamber; 2) environmental control system which utilizes the motor to realize the stable and comfortable air conditioning and cabin pressurization control, and is also the system with largest non-propulsion thermal energy consumption; 3) the avionics and commercial service systems which require the electricity to maintain the operation of their multiple computers.
Since the electrical demands driven by the fuel pump and environmental control system are usually time-varying and (indirectly or directly) highly correlated with the change of the different flight phases, the electrical energy provided by the generator of turbofan engine shall thus be correlated with the flight profile [2], [22]. A typical flight phase usually belongs to one of the three typical flight stages [23]: climb, cruise and descent which are differentiated by the aircraft's altitude, velocity and flight-path attitude associated with the flight stages. We define the duration between the take-off at the climb stage to landing at the descent stage for a MECA is T minutes. The smallest scheduling snapshot for the energy optimization analysis of turbofan engine is 1min. Then, the problem is a T -periods scheduling problem. We define the average altitude, velocity and flight-path attitude at time period t = 1, 2, . . . , T as (h (t), v (t), ϑ (t)) and if the flight phases for t = 1, 2, . . . , T are given, (h (t), v (t), ϑ (t)) for t = 1, 2, . . . , T will completely describe the whole flight profile and the values of the air static pressure P 0 (t) and static temperature T 0 (t) [23]. Besides of that, we have given the station numbering of the turbofan engine in Fig. 1 and will address the total pressure/temperature of different components of turbofan engine associated with the constructed numbering (i.e., P x (t) and T x (t) for x = 1, 2, . . . , 8 are separately defined as the total pressure and total temperature at the position of associated numbering in Fig. 1).
At first, the modeling of the energy optimization analysis for MECA is needed to explicitly present the objective function with a minimal energy utilization. Since all of the energy demands of MECA are driven from the fuel source, the objective function is given in (1) to realize a minimal fuel consumption: (1) q in (t), is a decision variable which represents the mass flow rate of the air in the inner duct of turbofan engine at period t, and is equal to all of the associated parameters in this formulation are defined in Nomenclatures in which Ma (t) represents the Mach number before the turbofan engine inlet at period t and is equal to v (t) √ γ i ·g c ·R·T 0 (t) , and P 1 (t) and T 1 (t) separately depict the total pressure and temperature before the inlet of turbofan engine at period t and are separately equal to P 0 (t) 2 ).A in (t), which defines the area of the inner duct of turbofan engine, is the decision control variable associated with q in (t)) [2]. c f represents the transformation coefficient between the fuel mass flow rate and the fuel cost and f is the parameter which represents the optimal ratio of fuel to air mass flow amount within the combustion chamber.
Eq. (2) stipulates the isentropic change of the total temperatures and total pressures associated with the fan.P 2 (t) and T 2 (t), which are separately the total pressure and temperature before the fan at period t, have the following relationships P 2 (t) = π i · P 1 (t) and T 2 (t) = T 1 (t) where π i is the fixed pressure ratio in the inlet [24].
Eq. (3) stipulates the isentropic change of the total temperatures and total pressures associated with the low-pressure booster.
Eq. (4) stipulates the isentropic change of the total temperatures and total pressures associated with the high-pressure compressor.
Eq. (5.a) and (5.b) describe: there has an isobaric change in the combustion chamber where the pressure keeps 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 within the inner duct plus to the combustion energy released by the fuel.
Eq. (6) stipulates the isentropic change of the total temperatures and total pressures associated with the high-pressure turbine.
Eq. (7) stipulates the isentropic change of the air state variables associated with the low-pressure turbine.
Eq. (8.a) and (8.b) separately depict the relationships between the rotation speed and the total temperature variation within the fan/booster at the low-pressure axis; (8.c) stipulates the relationship between the rotation speed and the total temperature variation within the compressor at the high-pressure axis (in Fig. 1).
q in (t) · η m h · (1 + f ) · c pt · (T 6 (t) − T 7 (t)) = q in (t)·c pc ·(T 5 (t)−T 4 (t)), t = 1, 2, . . . , T (9.a) Since the compressor and HPT are co-axial and the fan, the booster and LPT are mechanically co-axial, (9.a) and (9.b) separately describe the power balance constraints for the two mechanical axes: the reduced amount of the total enthalpy for the HPT at period t is equal to the increased amount of the total enthalpy for the high-pressure compressor at period t; the reduced amount of the total enthalpy for the LPT at period t is equal to the increased amount of the total enthalpies for the booster and the fan at period t plus to the electrical power output driven by the co-axial generator within this period (p g (t)). η m h and η m l (both < 1) are separately the mechanical efficiency for the high-pressure and low-pressure axis. B in (9.b) is defined to be the bypass ratio of the turbofan engine.
We define the required thrust at period t to be F (t) and its value can be determined by (h (t), v (t), ϑ (t)) through the standard flight mechanics analysis [23]. We assume the two branches of ram air into the turbofan engine are completely expanded at core and bypass nozzles separately (i.e., the static pressures after the core and bypass nozzles are both equal to P 0 (t)) [24]. Thus, Eq. (10) will represent the flight thrust balance relationship obtained from aircraft engine theory [2], [24] and the first and second terms at the right side of the equation are separately the thrust values provided by the inner and outer ducts.
Eq. (11) states electrical power balance at period t between the generator and the different electric demands by three main kinds of aircraft's subsystems, which are: p 1 (t), the electric power at period t by fuel pump of fuel supply system; p 2 (t), the electric power at period t by environmental control system which requires the largest thermal power in aircraft; p 3 (t), the electric power at period t by avionics and commercial service systems. p 1 (t) is equal to q in (t)·f ρ f · P 5 (t)(ρ f is the fuel density) and is thus indirectly correlated with (h (t), v (t), ϑ (t)) through the decision variables q in (t) and P 5 (t);p 2 (t) can be determined through the power optimization on environmental control system and can be finally quantified as a representation of parameters (h (t), v (t), ϑ (t)) [22], [25], [26]. To conclude, the above two electric demands can be both modeled as a formulation (indirectly or directly) correlated with (h (t), v (t), ϑ (t)). In contrast, since the computers in avionics and commercial service systems are usually operated under stable power level, p 3 (t) is usually taken as a fixed value for all t = 1, 2, . . . , T .
Eq. (12) represents that all of the associated decision variables are non-negative and their correspondingly upper bounds are stipulated by the associated physical and stable operation limits. The above modeling in (1)-(12) is a complex nonconvex mathematical optimization problem. We will propose an efficient algorithm in the next section to solve this problem.

III. THE ALGORITHM
The nonlinearity of the above problem in (1)-(12) can be generalized as follows (x 1 , x 2 represent the variables): 1) (2)-(4), (6), (7) and (10) have a non-convex form x 1 / x 2 ; 2) (8.a)-(8.c) have a convex form (x 1 ) 2 ; 3) (9.b), (10) and (11) have a bilinear form x 1 · x 2 . An efficient way to tackle the problem in (1)-(12) can use (Benders) decomposition method and select the complicating variables which separate the original problem as the outer-layer and inner-layer problems; when the complicating variables are fixed, it can make the inner-layer problem more tractable [27]. As for the problem in (1)-(12), the outer-layer problem will cope with the complicating variables which are taken as fixed parameters in the inner-layer problem. The inner-layer problem usually has a linear/nonlinear convex form which can finally help the outerlayer problem formulate the Benders cuts [28]. (The Benders cuts usually have a linearly approximated formulation about the solution point found in the inner-layer problem and can reduce the search region of the outer-layer problem.) The outer-layer problem and the inner-layer problem are solved successively until they both reach the required convergence criterion [27]. We at first select P 3 (t), T 4 (t), P 5 (t), P 6 (t), T 7 (t), P 8 (t), q in (t) and p g (t) as the complicating variables. Thus, with given values of the complicating variables, all of the nonlinear constraints in (1)-(12) will thus have a convex form. The outer-layer problem will thus be formulated and easily linear-approximately restated by a mixed integer and linear programming (MILP) problem which can be solved efficiently by the advanced commercialized software. The above two-layer Benders decomposition method can help to obtain an optimal fuel consumption result with only considering the complicating variables at the outer layer, and to find the feasible schedule results for all of the convex constraints at the inner layer. If any violations of the bounds on the variables exist at the outer layer, the Benders cuts will be formulated and sent back into the next round of the outer-layer problem [28].
The proposed algorithm is given as follows: 1) x * 0 is defined to be the optimal solution for the following problem: Then, solve the inner-layer optimization problem with given x * 0 (the solution method will be given later in section 3.2); y * 0 is defined to be the optimal solution of the inner-layer problem with given x * 0 . Set the lower bound at outer-layer L O−layer = −∞, upper bound at outer-layer U O−layer = +∞, the iteration number at outer-layer k = 1 and the convergence tolerance level at outer-layer ζ = 0.001.
2) k ≥ 1: Step 1) Solve the outer-layer problem. The outer-layer problem is: α k and β k are the introduced positive slack variables and M is a column vector in which the elements are very large positive real value.y * k−1 is obtained at the (k-1)-th round of inner-layer problem and is the fixed parameter for (19). From (14) and (16), we can find parts of the formulations in G (x k , y * k−1 ) are nonlinear corresponding to x k with given y * k−1 (the constraints (3), (4), (6), (7), (9.b) and (10)) and parts of the formulation in H 1 (x k ) are nonlinear corresponding to x k (the constraints (11)); for that matter, in order to solve the outer-layer problem, we at first recast the problem in (19) as a linear-approximated MILP formulation and then resort to the commercial optimization software to finally solve it. In this paper we use the linearization method in [29] to represent the strictly decreasing/increasing convex formulations in (3), (4), (6), (7) and (10) and use the linearization method in [30] to linearly approximate the bilinear function in (9.b) and (11). Let (x * k , α * k , β * k ) be the optimal solution for the k-th round of outer-layer problem and Step 2) Check the outer-layer convergence. If U O−layer − L O−layer < ζ , stop and return (x * k , y * k ); Otherwise, set U O−layer = c T x * k + M T α * k + M T β * k and go to the step 3.
Step 3) Solve inner-layer problem. We will address the inner-layer problem in section 3.2. Here we only define y * k to be the optimal solution of the inner-layer problem given x * k and set y * k as parameters for the (k+1)-th round of the outer-layer problem. Set k = k + 1 and go to the step 1.

B. INNER-LAYER
The inner-layer problem has a convex formulation, an optimal solution of the inner-layer problem can thus be efficiently obtained by the successively linear programming technique [27] and this solution can generate a valid cut for the outerlayer problem. All of the nonlinear terms in the inner-layer problem are successively linearized [27], [31] around intermediate solution points and the Benders cuts generated by solving the inner-layer problem are finally added into the outer-layer problem.
1) Fix x * k (k ≥ 0), set the lower bound at inner-layer L I −layer = −∞, upper bound at inner-layer U I −layer = +∞, iteration number at inner-layer j = 1 and choose the convergence tolerance level at inner-layer δ = 0.0001.
2) j ≥ 1: Step 1) Formulate the inner-layer problem. The inner-layer problem is basically to find a solution of y j k within G (x * k , y j k ) = g, H 2 (y j k ) = h 2 , y j k ∈ Y. We thus restate the inner-layer problem as follows: Step 2) j-th round of iteration at inner-layer. We can find, the problem in (20) can be equivalently recast as two independent linear programming (LP) problems if y j k VOLUME 8, 2020 is given, and the optimal values of the two problems, which are defined as D j 1 (x * k , y j k ) and D j 2 (y j k ), can be determined analytically through separately solving the following two problems: where ∀µ 1 ∈ [1, 2, . . . , 10·T ] and the constraint set 1 (j, k) is given as: where ∀µ ∈ [1, 2, . . . , 2 · T ] and the constraint set 2 (j, k) is given as follows: Moreover, we can find, the problem in (21) can be separated as 10 · T of independent ''two-variables'' LP problems and the optimal value D j 1 (x * k , y j k ) can thus be easily analytically determined as 1 T G (x * k , y j k ) − g ; similarly, the problem in (22) can be separated as 2 · T of independent ''twovariables'' LP problems and the optimal value D In (23), the are defined to be a vector with very small positive values. From (23), we can findd j k is a strictly Step

IV. NUMERICAL ANALYSIS A. DATASET
The test data are obtained from a conventional turbofan engine model which is refitted to be a ''no-bleed'' one and used for a single-aisle MECA. Since a single-aisle MECA usually owns two turbofan engines operated as same mode, to make the statement more concisely in the numerical analysis, we revised the associated data of one turbofan engine through doubling the values of the data and can thus refit it into a ''larger'' turbofan engine where we assume this single turbofan engine can provide sufficient thrust and electric energy for the MECA. In the numerical tests of this paper, the flight duration is set to be 1.5 hours and thus the scheduling horizon T will be 90. The definition of the three kinds of flight phases with the associated time period allocation are presented in Table 1 (here we want to clarify that we assume the phases of the climb and descent for an aircraft will have a fixed flight-path angle, with considering the fact that the future tendency of next generation air transportation will adopt the continuous climb/descent procedures [32], [33]). We then obtain (h (t), v (t), ϑ (t)) for t = 1, 2, . . . , T based on Table 1 and they will be utilized to determine P 0 (t) and T 0 (t) (P 0 (t) can be determined through the relationship between P 0 (t) and h (t); T 0 (t) is obtained through ideal gas equation of state [23]) and further Ma (t), P 1 (t) and T 1 (t). Finally, based on flight mechanics analysis [23], we can also obtain the thrust value F (t) based on (h (t), v (t), ϑ (t)).
The revised data of the integrated single turbofan engine are given in Table 2. Meanwhile, with considering the thermal power loads required by environmental control system can usually be divided by four main kinds [22], [25] (conductive, solar, internal electric and occupants' heat loads associated with the pressurization cabin of aircraft), the total thermal power load by environmental control system p 2 (t) will thus  be equal to the sum of the above four kinds of loads. A typical combination of these four kinds of thermal loads is given in Table 3 which indirectly reflects the relationship between p 2 (t) and (h (t), v (t), ϑ (t)) for t = 1, 2, . . . , T (T s (t) is defined as the total temperature of the aircraft skin and is equal to T 1 (t) · (1 + 0.18 · (Ma (t)) 2 ) [26]; T cabin (t), which is the required temperature within the aircraft cabin so as to maintain the crews/passengers comfortable and electric facilities operated, is set to be 25 • C). The electric power at period t by avionics and commercial service systems p 3 (t) will be set to be a fixed value (20kw) for all time periods.
The proposed two-layer Benders decomposition algorithm is operated under Visio Studio 2016 environment [34] and is programmed by Gurobi 8.1.1 [35] on a PC laptop with a 2.50 GHz CPU and 4 GB RAM. The MILP gap for the outerlayer problem is 0.0001.

V. NUMERICAL RESULTS
We firstly obtain the computational results of the total fuel cost under the different combinations of the cruise altitudes (8000m-12000m) and cruise velocities (700km/h, 750km/h, 800km/h, 850km/h, 900km/h) (shown in Fig. 2). From Fig. 2, we generalize the two findings: 1) under a fixed cruise altitude, a larger cruise velocity introduces more fuel consumption. It is because that a larger cruise velocity will introduce more drag which is required to be balanced by the thrust and also generate more heat through the friction between the air and aircraft skin such that more cooling power by environmental control system will be required to reach balanced temperatures in the pressurization cabin of the aircraft; 2) in contrast, under a fixed cruise velocity, a higher cruise altitude generally needs less fuel consumption because a higher cruise altitude requires less thrust. However, with the cruise altitude increasing, the thermal power required by environmental control system will also be increased; this will introduce the fact that, under some particular cases such as the cruise velocities 700km/h and 750km/h, the optimal cruise altitudes with minimal fuel cost are not their largest cruise altitudes but separately 11000m and 11500m.
Secondly, besides of the above results, the curves in Fig. 3 make a comparison analysis between the thrust values, electric power loads and fuel costs under different time periods under the typical cruise altitude (10km) and cruise velocity (800km/h). Typically, the variation of the fuel cost curve is generally in accordance with those of thrust and electric power load curves (at climb stage, the three curves increase; at cruise stage, the three curves keep constant; at descent stage, the three curves decrease). However, in contrast with the curve of electric power load, the fuel cost curve is more correlated with the curve of thrust value since the supply of thrust usually requires a larger amount of energy than that of the electric load in aircraft.
Finally, the economic efficiency is also analyzed under different time periods. In Fig. 4, we give the numerical results of the fuel cost under typical combinations of cruise altitude and cruise velocity. In Fig. 4, ''HA'', ''LA'', ''LV'' and ''SV'' separately represent high cruise altitude (12000m), low cruise altitude (8000m), large cruise velocity (900km/h) and small cruise velocity (700km/h). From the numerical results we can find: the minimal and maximal amount of fuel will be separately consumed at HA & SV and LA & LV case; if the cruise altitude is fixed, the total fuel cost is positively and linearly correlated with the cruise velocity and if the cruise velocity is fixed, the total fuel cost will be nearly negatively and linearly correlated with the cruise altitude.

VI. CONCLUSION
Under more-electric environment, the civil aircraft's turbofan engine utilizes more fuel resource to provide increased electric energy demands besides of its conventional responsibility of maintaining the fundamental thrust. This recast the energy optimization analysis of aircraft turbofan engine, which was a simple optimal setting of fuel usage for keeping required thrust, into a complex scheduling problem which is to minimize the fuel consumption while obeying the multiple thermodynamic constraints on operation limits by turbofan engine and electric power supply and demand balance associated with the varying flight profiles. In order to solve this problem under more-electric environment, we proposed a new mathematical modeling for the energy optimization analysis of turbofan engine and formulate it as a nonlinear programming form. We provided an algorithm based on Benders decomposition method to tackle this problem and illustrated the numerical analysis to demonstrate the operational effectiveness of the proposed algorithm on a turbofan engine test case. ρ f : fuel density before the fuel is pumped into the chamber (g/cm 3 ) P 3 : upper bound of total pressure after fan (Pa) T 3 : upper bound of total temperature after fan (K) P 4 : upper bound of total pressure after low-pressure booster (Pa) T 4 : upper bound of total temperature after lowpressure booster (K) P 5 : upper bound of total pressure after high-pressure compressor (Pa) T 5 : upper bound of total temperature after highpressure compressor (K) P 6 : upper bound of total pressure after the combustion chamber (Pa) T 6 : upper bound of total temperature after the combustion chamber (K) P 7 : upper bound of total pressure after the highpressure turbine (Pa) T 7 : upper bound of total temperature after the highpressure turbine (K) P 8 : upper bound of total pressure after the lowpressure turbine (Pa) T 8 : upper bound of total temperature after the lowpressure turbine (K) n l : upper bound of rotation speed of the lowpressure axis (rpm) n h : upper bound of rotation speed of the highpressure axis (rpm) A in : upper bound of cross section area of the inner duct of inlet (m 2 ) p : upper bound of generator power output (w) B. DECISION VARIABLES P 3 (t) : total pressure after fan at period t (Pa) T 3 (t) : total temperature after fan at period t (K) P 4 (t) : total pressure after low-pressure booster at period t (Pa) T 4 (t) : total temperature after low-pressure booster at period t (K) P 5 (t) : total pressure after high-pressure compressor at period t (Pa) T 5 (t) : total temperature after high-pressure compressor at period t (K) P 6 (t) : total pressure after the combustion chamber at period t (Pa) T 6 (t) : total temperature after the combustion chamber at period t (K) P 7 (t) : total pressure after the high-pressure turbine at period t (Pa) T 7 (t) : total temperature after the high-pressure turbine at period t (K) P 8 (t) : total pressure after the low-pressure turbine at period t (Pa) T 8 (t) : total temperature after the low-pressure turbine at period t (K) n l (t) : rotation speed of the low-pressure axis at period t (rpm) n h (t) : rotation speed of the high-pressure axis at period t (rpm) q in (t) : mass flow rate of the air imported through the inner duct at period t (kg/sec) A in (t) : cross section area of the inner duct of inlet at period t (m 2 ) p g (t) : generator power output at period t (w)