Multi-Level Modeling Methodology for Optimal Design of Electric Machines Based on Multi-Disciplinary Design Optimization

The transportation sector is undergoing electrification to gain advantages such as lighter weight, improved reliability, and enhanced efficiency. As contributors to the safety of embedded critical functions in electrified systems, better sizing of electric machines in vehicles is required to reduce the cost, volume, and weight. Although the designs of machines are widely investigated, existing studies are mostly complicated and application-specific. To satisfy the multi-level design requirements of power systems, this study aims to develop an efficient modeling method of electric machines with a background of aircraft applications. A variable-speed variable-frequency (VSVF) electrically excited synchronous generator is selected as a case study to illustrate the modular multi-physics modeling process, in which weight and power loss are the major optimization goals. In addition, multi-disciplinary design optimization (MDO) methods are introduced to facilitate the optimal variable selection and simplified model establishment, which can be used for the system-level overall design. Several cases with industrial data are analyzed to demonstrate the effectiveness and superior performance of the modeling method. The results show that the proposed practices provide designers with accurate, fast, and systematic means to develop models for the efficient design of aircraft power systems.


Introduction
Electrified transportation reveals significant benefits in terms of increased reliability, improved energy efficiency, reduced emissions, and enhanced passenger comfort. With the breakthrough of power electronics, the adoption of the More Electric Aircraft (MEA) and electric propulsion aircraft concepts have become a major trend in the air industry [1,2]. As a result, many hydraulic and pneumatic power-driven systems are being replaced by their electrical counterparts in existing MEAs, e.g., Boeing 787 [3] and Airbus 380 [4]. Moreover, the hybrid/all-electric propulsion electric grid is now being investigated [5][6][7]. In electrified aircraft, the electric machine is the key equipment to realize the conversion and transformation of multi-physical energies; therefore, the sizing of the machine has become one of the leading topics in the study of transport electrification.
In addition, with the development of artificial intelligence and information technologies, model-based system engineering (MBSE) [8][9][10] and optimization-assisted design [11] have become the focus of researchers, becoming increasingly mature, and will be particularly powerful for the design of aircraft power systems. In this context, models of the electrical environment control system (ECS) [12,13],

Modeling Needs Analysis
Traditional machine simulation models calculate variables of power and energy used for the component selections [14]. Both steady-state and dynamic simulations are mandatory due to the complex couplings of multi-physical characteristics. These simulations are usually based on lumped circuit models or field models, and run in environments such as MATLAB/Simulink, Dymola, and ANSOFT. To run this kind of model, it is necessary that all component parameters should be known. Meanwhile, it is highly time-consuming to extract features from models. Therefore, it is suitable for detailed design, analysis, and verification, whereas time-independent symbolic analytical models are preferred for assisting with system integration.
From the perspective of design optimization, design requirements can be translated into design, constraint, and target variables of different levels. Then, the design is conducted hierarchically by algebraic solvers or intelligent algorithms to achieve the best goals with the satisfaction of constraints by well-determined design variables; see Figure 1. As the impacts of variations of variables can be propagated both vertically (from bottom to top, e.g., components, machine, and system) and horizontally (e.g., among components), every design is a combinational optimization problem based on both the underlying design results and the synthesis forms of underlying-level designs, e.g., the selected component technologies and the selected machine structure in the machine design process. Therefore, iterative calculations may occur because of the improper design or integration of the underlying units; in Figure 1, see the red and green-colored directed lines. Therein, the green-colored lines can be avoided by optimization algorithms, whereas the red-colored lines cannot be avoided and may lead to blind designs, due to the invisible and complex cross-links among variables.  Accordingly, the hierarchical designs require a well-structured modeling approach:  The primary-level model handles the design of the machine for the integration of More-electric systems, linking dimensions and characteristics of selected materials and technologies to multiphysical characteristics of components with full consideration of their energy conversion features [29]. Thus, couplings of variables can be illustrated clearly.
 The secondary-level model allows proper simplification of the primary-level multi-physics model, which focuses on the performance evaluation of the system, and should be consistent with the primary-level model.
In this way, the developed modeling approach covers enough information needed for different levels of design and evaluation, and appears to be the most appropriate means of accelerating the design process, which is the major task and objective of this paper. Accordingly, the hierarchical designs require a well-structured modeling approach:

Modular Multi-Physics Model
• The primary-level model handles the design of the machine for the integration of More-electric systems, linking dimensions and characteristics of selected materials and technologies to multi-physical characteristics of components with full consideration of their energy conversion features [29]. Thus, couplings of variables can be illustrated clearly.

•
The secondary-level model allows proper simplification of the primary-level multi-physics model, which focuses on the performance evaluation of the system, and should be consistent with the primary-level model. In this way, the developed modeling approach covers enough information needed for different levels of design and evaluation, and appears to be the most appropriate means of accelerating the design process, which is the major task and objective of this paper.

Modular Multi-Physics Model
To demonstrate the propagations of the impacts of variables, both multi-physical characteristics and combinatorial explosions of components are required in the multi-physics model. On the other hand, the modular approach in electrical machines denotes that their parts (the stator, the rotor, or both) are made up of distinctive segments, which can increase flexibility of design. Therefore, the structure of the improved multi-physics model is illustrated in Figure 2. To demonstrate the propagations of the impacts of variables, both multi-physical characteristics and combinatorial explosions of components are required in the multi-physics model. On the other hand, the modular approach in electrical machines denotes that their parts (the stator, the rotor, or both) are made up of distinctive segments, which can increase flexibility of design. Therefore, the structure of the improved multi-physics model is illustrated in Figure 2. Every component model comprises three parts, namely, technical parameters, multi-physics coupled constraints, and the design goal, in which the multi-physics coupled constraints represent the core linking the other two parts. Technical parameters comprise the design codes from design experiences, standards, and protocols (i.e., densities and strengths of materials, technology readiness level (TRL), etc.), which are used as pre-specified design requirements. According to the selected configurations in the technical parameters, the multi-physics coupled constraints describe the relationship between design variables and performance according to the principle of energy conversion and transmission; structure parameters describing the geometry and material information can be then obtained. Based on these two parts, the design goal calculates evaluation indexes depending on design requirements. The component integration model has a similar structure to component models. The difference is in the components coupled constraints describing the cross-links of components. Differently from the traditional machine design method, the machine here is designed based on the direct approach. The detailed process is as follows: according to the required performance (i.e., the rated power, operating points of components, etc.), main radial and axial dimensions are initially calculated by the basic equations of the machine; then, detailed parameters of components (e.g., dimensions of the rotor, stator, yoke, or windings) are determined by optimal technology selections; next, the whole Every component model comprises three parts, namely, technical parameters, multi-physics coupled constraints, and the design goal, in which the multi-physics coupled constraints represent the core linking the other two parts. Technical parameters comprise the design codes from design experiences, standards, and protocols (i.e., densities and strengths of materials, technology readiness level (TRL), etc.), which are used as pre-specified design requirements. According to the selected configurations in the technical parameters, the multi-physics coupled constraints describe the relationship between design variables and performance according to the principle of energy conversion and transmission; structure parameters describing the geometry and material information can be then obtained. Based on these two parts, the design goal calculates evaluation indexes depending on design requirements. The component integration model has a similar structure to component models. The difference is in the components coupled constraints describing the cross-links of components. Differently from the traditional machine design method, the machine here is designed based on the direct approach. The detailed process is as follows: according to the required performance (i.e., the rated power, operating points of components, etc.), main radial and axial dimensions are initially calculated by the basic equations of the machine; then, detailed parameters of components (e.g., dimensions of the rotor, stator, yoke, or windings) are determined by optimal technology selections; next, the whole machine is integrated in aspects of structural, electromagnetic, and thermal features; finally, constraints and initial assumptions are checked to verify the design. Calculations in [35] show that the differences between initially assumed parameters and calculated results based on the direct approach are small enough to form the closed loop of design.

Component Models
In this way, functional and performance characteristics are included and couplings of parameters are demonstrated explicitly in the multi-physics model. With changes in component structures and their integration forms, the trade-off of different machines can be achieved easily. Moreover, design codes are integrated into the model in this paper. Therefore, significant expertise for each component is not required for integrators to implement and shrink the component range to a reduced number of candidate technologies in system-level studies [24].

MDO-Based Model Simplification
To satisfy the system-level design demands, the MDO method is first used in the simplification of the multi-physics model to generate the simplified model, so that secondary-level models with characteristics of accuracy, simplicity, and consistency can be obtained. This process is illustrated in Figure 3. the differences between initially assumed parameters and calculated results based on the direct approach are small enough to form the closed loop of design. In this way, functional and performance characteristics are included and couplings of parameters are demonstrated explicitly in the multi-physics model. With changes in component structures and their integration forms, the trade-off of different machines can be achieved easily. Moreover, design codes are integrated into the model in this paper. Therefore, significant expertise for each component is not required for integrators to implement and shrink the component range to a reduced number of candidate technologies in system-level studies [24].

MDO-Based Model Simplification
To satisfy the system-level design demands, the MDO method is first used in the simplification of the multi-physics model to generate the simplified model, so that secondary-level models with characteristics of accuracy, simplicity, and consistency can be obtained. This process is illustrated in Figure 3.
By manipulating model inputs according to pre-specified design space, design of experiments (DOE) technology is applied to generate the data sample for the analysis. To describe the whole performance of the multi-physics model effectively, sufficient design points evenly spread in the design space are required. Then the model-based sensitivity analysis method is introduced to explore the design space and find cause-and-effect relationships between the inputs and outputs of the multiphysics model. From the sensitivity analysis, an optimal parameter combination for system-level machine models can be concluded. Based on the data and the optimal selected variable combination, surrogate models, which have been widely used in engineering design due to the low computational cost and high fidelity [36,37], are built for the system-level study. The approximation model can be improved by using better learners and adding additional sample points to achieve better precision.

Sensitivity Analysis
Classification Model

Model Approximation
Model Initialization

Model Validation
Simplified Models Figure 3. Multi-disciplinary design optimization (MDO)-based simplification of the multi-physics model.

Modular Multi-Physics Model
The variable-speed variable-frequency (VSVF) electrically excited synchronous generator is chosen as the case study to illustrate the improved multi-physics modeling process. The 2D structure of the machine is shown in Figure 4. In this study, the weight and power loss are naturally selected as the objectives. According to the proposed method and [29,35], the improved multi-physics model of machines includes the sub-models as discussed in the following. By manipulating model inputs according to pre-specified design space, design of experiments (DOE) technology is applied to generate the data sample for the analysis. To describe the whole performance of the multi-physics model effectively, sufficient design points evenly spread in the design space are required. Then the model-based sensitivity analysis method is introduced to explore the design space and find cause-and-effect relationships between the inputs and outputs of the multi-physics model. From the sensitivity analysis, an optimal parameter combination for system-level machine models can be concluded. Based on the data and the optimal selected variable combination, surrogate models, which have been widely used in engineering design due to the low computational cost and high fidelity [36,37], are built for the system-level study. The approximation model can be improved by using better learners and adding additional sample points to achieve better precision.

Modular Multi-Physics Model
The variable-speed variable-frequency (VSVF) electrically excited synchronous generator is chosen as the case study to illustrate the improved multi-physics modeling process. The 2D structure of the machine is shown in Figure 4. In this study, the weight and power loss are naturally selected as the objectives. According to the proposed method and [29,35], the improved multi-physics model of machines includes the sub-models as discussed in the following.

Multi-Physics Coupled Constraints
The shaft model involves mechanical, structural, and cooling features. The shaft is assumed to be a hollow cylinder, the diameter of which is determined by the electromagnetic power Pem, rotational speed n, and the mechanical properties of the material: where Ds is the outer diameter; Ash is the coefficient of shear stress; Bsh is the stiffness coefficient of material; [σ-1]b is the admissible bending stress; αs is the ratio of the inner to outer diameter, and F is the equivalent moment. Pem can be calculated according to Equations (2) and ( where P is the rated power; p represents the pole pairs; η is the assumed efficiency; cosφ is the power factor, and KE is the per-unit value of full-load potential, namely, the ratio of induced electromotive force to the output voltage at full loads. The gravity center of the rotor is assumed to be at the center of the shaft, and to be at the same position as the circle center of the stator. Thus, F can be calculated according to Equation ( where lef is the efficient length of machine; Fw is the bending stress caused by the rotor weight Wr, and Ft is the rated torque.

Multi-Physics Coupled Constraints
The shaft model involves mechanical, structural, and cooling features. The shaft is assumed to be a hollow cylinder, the diameter of which is determined by the electromagnetic power P em , rotational speed n, and the mechanical properties of the material: where D s is the outer diameter; A sh is the coefficient of shear stress; B sh is the stiffness coefficient of material; [σ −1 ] b is the admissible bending stress; α s is the ratio of the inner to outer diameter, and F is the equivalent moment. P em can be calculated according to Equations (2) and (3) [38]: P em = K E P/ cosϕ, synchronous generator K E P/ cosϕ/η, synchronous/induction motor , where P is the rated power; p represents the pole pairs; η is the assumed efficiency; cosϕ is the power factor, and K E is the per-unit value of full-load potential, namely, the ratio of induced electromotive force to the output voltage at full loads. The gravity center of the rotor is assumed to be at the center of the shaft, and to be at the same position as the circle center of the stator. Thus, F can be calculated according to Equation (4): where l ef is the efficient length of machine; F w is the bending stress caused by the rotor weight W r , and F t is the rated torque.

Design Goals
The shaft weight W sh is shown in Equation (5): where V sh is the shaft volume, and ρ sh is the material density. For the VSVF system, the influences of the variable rotational speed should be taken into consideration to ensure security, e.g., the stresses of shaft diameters.

Multi-Physics Coupled Constraints
The yoke model involves magnetical, structural, and thermal features. As one of the main dimensions of a hollow cylinder, the yoke's thickness h y can be calculated according to the magnetic flux Φ y and the selected induction B y . The magnetic potential F y can be obtained according to the B-H curve of the selected material: where k c is the lamination factor [39]; D y is the average diameter of the yoke; ξ is the length coefficient of the yoke magnetic circuit, and H y is the yoke magnetic intensity corresponding to B y . It is assumed that B y0 is the yoke induction when the machine is in an idle run, whereas B yN is the yoke induction when the machine runs at full load. The relationship between B y0 and B yN is shown in Equation (7). By Equation (6) and (7), magnetic potentials of the yoke at no load and full load can be calculated. Combined with the magnetic potentials of other components calculated in a similar way, the magnetic circuit can be integrated, and the design requirements of the excitation windings can be obtained.
where |E i | is the amplitude of the per-unit value of rated induced electromotive force, and the calculation is introduced in the winding component. The thermal circuit of the yoke satisfies Equation (8). Note that other components have similar thermal models.
where T y is the temperature rise of the yoke; P iron_y is the produced heat by the yoke iron loss; R y is the thermal resistance of the yoke, which is dependent on the material and proportional to the surface area [20]; P h _ y is the stored thermal power of the yoke; C y is the yoke thermal capacity, which is the product of the specific heat SpC y and mass m y ; ∆T y is the fluctuation of T y , and P y_o is the thermal power transferring from the yoke to other components. At the steady heat balance state, as the yoke temperature is fixed, namely ∆T y = 0, P y_o equals P iron_y , which indicates that the branch of C y is ignored in the thermal circuit in the design phase.

Design goals
The yoke weight W y and the iron loss P y can be calculated according to Equation (9) where ρ y is the density of the yoke material; p y is the yoke iron loss per kilogram; f is the operating frequency; p 10/50_y is the yoke iron loss when B = 1 T and f = 50 Hz, and k a_y is the empirical coefficient, usually determined by Equation (10) [35] k a_y = 3.6, DC machine or P< 100 kVA AC machine 1.3, P> 100 kVA AC machine .

Multi-Physics Coupled Constraints
The slot model involves structural, magnetical, and thermal features. The shape of slots is selected based on the rated power P and armature winding voltage E by experience; see Equation (11) [38]: where it is assumed that the slot induction B s equals the induction at a one-third height of the slot from the dedendum; see Figure 4. According to Equation (12), the magnetic potential of slots F s can be calculated: where Z is the slot number; b s and H s are, respectively, the tooth width and magnetic intensity corresponding to B s , and h s is the slot height.
The slot inductions at no load and full load B s0 and B sN satisfy

Design Goals
Sections of every slot are assumed to be trapezoidal; thus, the slot weight W s and the iron loss P s can be calculated according to Equation (14): where ρ s is density of the yoke material; b z and b d are the widths of the addendum and dedendum respectively; N Z is the slot number per pole per phase; m is the phase number; D a is the bore diameter; p s is the slot iron loss per kilogram; p 10/50_s is the yoke iron loss when B = 1 T and f = 50 Hz, and k a_s is the empirical coefficient, usually determined by Equation (15) [35]:

Multi-Physics Coupled Constraints
The air-gap model involves structural, magnetical, and thermal features. The air-gap thickness δ, coefficient k δ , and magnetic potential F δ are illustrated in Equation (16) [38]: where A is the electric loading; τ is the pole pitch; x * d is the relative direct-axis inductance, ranging from 1.1 to 1.4; Φ δ is the air-gap magnetic flux; µ 0 is the permeability of vacuum, and B δ is the air-gap induction: where B δ0 and B δN are the air-gap inductions at no load and full load.
As the oil spray cooling method is chosen for this machine, the equivalent thermal circuit of cooling oil is included in the thermal model of the air gap. To simplify the calculation, it is assumed that the oil is coated evenly on the surfaces of the stator and rotor.

Multi-Physics Coupled Constraints
The pole model involves magnetical, thermal, and structural features. The shape of the poles is determined by the machine type.
In this article, as the study case is a salient pole synchronous machine, the pole shape is illustrated in Figure 4. It is assumed that the pole shoes are of the same circle center as the rotor yoke, and its width is assumed to be equal to τ. Then, the pole-core width b p can be calculated by exciting flux Φ p and selected pole induction B p . The magnetic potential of pole F p can be calculated by the corresponding magnetic intensity H p and pole height h p . Therein, B p has the same meaning as the afore-mentioned B y ; see Equation (19): where B p0 and B pN are the pole inductions at no load and full load, respectively. For a PMSG, the volume of the permanent magnet V m can be estimated by the Larionoff method; see Equation (20) [38]. For the non-salient pole synchronous generator, the volume is calculated via the unsaturated vectogram, which is omitted here.
where K bh is a coefficient related to the material and selected operating point of the poles, and pole-arc coefficient, etc.; P| cosϕ = 0 is the rated power when the power factor is zero; K σ0 is the leakage flux coefficient at no load, and is dependent on the rotor structure, material, and pole pairs; K adF is the equivalent coefficient of armature potential, which is dependent on the pole material; B r and H C are the remanent induction and coercive force of the pole material, respectively; K K is the short-circuit current ratio. Detailed dimensions can be calculated according to the pole structure (shape and magnetization direction); the process is omitted here but can be found in [35].
The selection of the operating point of poles is based on the material, including its magnetic and thermal characteristics. It is assumed that the demagnetization curve of the selected material is linear; thus, the operating point of poles (B m ,H m ) can be calculated by Equations (21) and (22): where B rT and H CT are the remanent induction and coercive forces at temperature T, respectively; B r0 and H C0 are the remanent induction and coercive forces at standard temperature T 0 , respectively; α B is the temperature coefficient of remanence; T is the actual operating temperature; T h and T l are the upper and lower limits of the allowed temperature, respectively, and B rh and B rl are the inductions at T h and T l , respectively. The influence of temperature on the irreversible magnetic degradation is ignored here.
where K σ is the leakage flux coefficient.

Design Goals
The pole weight W p and the iron loss P p can be calculated according to Equation (23): where ρ p is the density of the pole material; p p is the pole iron loss per kilogram, and p 10/50_p is the pole iron loss when B = 1 T and f = 50 Hz.

Multi-Physics Coupled Constraints
The winding model involves electromagnetic and thermal features. The configurations of the armature winding are determined by optimally selecting the magnet wire based on the carrying capacity, which is calculated by the rated power P, the armature winding voltage E, the pitch ratio β, and the current density J a , whereas those of the excitation winding is determined by the excitation magnetic potential F f and current density J f .
In detail, the number of winding turns in series N is determined by E, and the number of paralleled conductors a can be calculated according to the standard MIL-W-5088L and the wire gauge standard [40]; see Equation (24): where f is the frequency of the output voltage; k W is the winding coefficient; Φ is the main magnetic flux; β is the ratio of the winding pitch to the polar pitch τ; q is the number of slots per pole per phase; α is the electrical angle; m is the phase number; J is the current density of the selected wire gauge; d is the diameter of the selected wire gauge; k sh is the bundling coefficient of conductors; k g is the altitude correction coefficient of conductors, and k T is the temperature coefficient of conductors.
According to the requirements of the wire insulation [41], the diameter of the selected electromagnetic wire d ew satisfies Equation (25): where d max is set as the maximum wire diameter allowed for an inlay technology reason, and ε ew is the dielectric strength of the insulation.
According to the first-order model of the machine [42], the resistance and inductance can be calculated; see Equations (26) and (27): where R a is the resistance of the armature winding; k r is the coefficient considering the skin effect, which is dependent on the shape, configuration, and frequency of conductor; ρ c is the electrical resistivity of the conductor, which is dependent on the material (electrical resistivity at 20 • C ρ 0 and temperature coefficient α c ) and temperature rise ∆T, and l c and s c are the average conductor length and area per turn, respectively.
where X m , X M , X σ , and X a are the main, mutual, leakage, and total inductances of the armature winding, respectively; the specific permeance coefficients λ m , λ M , λ s , λ δ , λ t , and λ E are dependent on the slot shape and configurations of windings; detailed expressions are omitted here due to the existence of multiple kinds of combined constructs [35].
To integrate the magnetic circuit, the per-unit value of the rated induced electromotive force E i , which is a gain coefficient to obtain the rated operating points of every component, is calculated according to the impedance of armature windings [29]; see Equation (28): where R * a and X * σ are the per-unit value of R a and X σ , respectively; W and Q are the real and imaginary parts of E i , respectively, and ε is the phase angle of E i .
In addition, the winding thermal model includes the thermal circuit of both wire and insulation.

Design Goals
According to the determined strands of wires, the weight W aw and power loss P aw can be calculated; see Equation (29): where W aw_s is the weight of single-strand wire, and ρ ew and ρ ins are the densities of conductor and insulation materials, respectively.

Components Coupled Constraints
The whole machine model is integrated into aspects of structural, magnetical, and thermal features. The structural integration is mainly focused on the radial dimensions of components, including the main geometrical parameter shown in Equation (30), constraints of lengths shown in Equation (31), and the space-filling of windings shown in Equation (32): a l e f /(P em /n) = 6.1/(α k W AB δ ) where CA is the machine constant; D a is the diameter of the armature; λ is the dimension ratio; the pole-arc coefficient α' is assumed to be 2/π, which results in a sinusoidal magnetic field; D sy is the stator outer diameter, and K p is an empirical coefficient related to p.
As component dimensions are calculated by their respective inductions, when components are integrated into the structure, implicit constraints exist among inductions of components to make sure that they can be put in a limited space. Therefore, component inductions should be properly selected to satisfy Equations (31) and (32): where h sy and h ry are the thicknesses of the stator and rotor yoke, respectively; a w and b w are the tangential and radial lengths of excitation windings, respectively; r r is the radius of the rotor; θ is one-quarter of the mechanical angle of machine, and K is the slot fill factor. The magnetical integration calculates the magnetic circuit, and the excitation requirement is obtained by adding the no-load and full-load magnetic potentials of components; see Equation (33): where F ry and F sy are the magnetic potentials of rotor yoke and stator yoke, respectively; F a is the magnetic motive force of armature winding, and F ad is the direct-axis magnetic motive force of armature winding; detailed expressions are omitted here as this is a complex equation. The thermal integration illustrates the thermal behavior of every component (slot, yoke, etc.) with respect to their resistances, capacities, and ambient conditions, by a thermal circuit, which is coupled with electromagnetic characteristics through losses; see Figure 5. Requirements of the cooling system, e.g., the flow rate of the cooling oil [43], can be deduced from Equation (6) by calculating the thermal resistance and capacity of cooling oil R co and C co , which offers a resting interface to the fuel pump system. In addition, the thermal load AJ a is also given as a constraint here. The thermal integration illustrates the thermal behavior of every component (slot, yoke, etc.) with respect to their resistances, capacities, and ambient conditions, by a thermal circuit, which is coupled with electromagnetic characteristics through losses; see Figure 5. Requirements of the cooling system, e.g., the flow rate of the cooling oil [43], can be deduced from Equation (6) by calculating the thermal resistance and capacity of cooling oil Rco and Cco, which offers a resting interface to the fuel pump system. In addition, the thermal load AJa is also given as a constraint here.  Figure 5. The simplified lumped thermal model of the electrically excited synchronous generator.

Design Goals
The total weight WG and power loss PG are calculated by Equation (34): where Wry, Wew, Waw, and Wsy are the weights of the rotor yoke, excitation windings, armature windings, and stator yoke, respectively; Pew and Paw are the copper loss of excitation and armature windings, respectively; Pry and Psy are the iron loss of the rotor and stator yoke, respectively; PΔ is the stray loss, and is assumed to be one percent of the rated power, and Pf is the total friction loss; see Equation (35) [44]. Therein, the rotor is assumed to be a cylinder when the windage friction loss is calculated: where Psb and Pair are the friction losses of bearing and windage, respectively; ρair is the air density at temperature T (℃) and pressure press (mmHg); ω is the angular velocity of the rotor, and Kf1 and Kf2 are the friction coefficients of the cylinder and disk, and are dependent on the Couette Reynolds Number Re.
It is assumed that Kf1 and Kf2 are equal to Kf, and Kf can be calculated by Equation (36) where Cr is the relative roughness coefficient of the rotor, and is assumed to be 1, and Re is obtained according to Equation (37):

Design Goals
The total weight W G and power loss P G are calculated by Equation (34): where W ry , W ew , W aw , and W sy are the weights of the rotor yoke, excitation windings, armature windings, and stator yoke, respectively; P ew and P aw are the copper loss of excitation and armature windings, respectively; P ry and P sy are the iron loss of the rotor and stator yoke, respectively; P ∆ is the stray loss, and is assumed to be one percent of the rated power, and P f is the total friction loss; see Equation (35) [44]. Therein, the rotor is assumed to be a cylinder when the windage friction loss is calculated: where P sb and P air are the friction losses of bearing and windage, respectively; ρ air is the air density at temperature T ( • C) and pressure p ress (mmHg); ω is the angular velocity of the rotor, and K f1 and K f2 are the friction coefficients of the cylinder and disk, and are dependent on the Couette Reynolds Number Re. It is assumed that K f1 and K f2 are equal to K f , and K f can be calculated by Equation (36) [45]: where C r is the relative roughness coefficient of the rotor, and is assumed to be 1, and Re is obtained according to Equation (37): where ν is the kinematic viscosity of the airflow.

Problem Definition
The machine design problem can be sorted into the forms illustrated in Equation (38) by the equations above.
Min [W G , P G ] s.t. Explicit : multi-physics couplings of components Implicit : structural couplings among components initial assumptions verification (38) where the weight, power loss, physical dimensions, and cost of the generator are usually the major design goals of a machine in aircraft applications. In this study, the weight and power loss are naturally selected as the objectives, and the physical dimension of the machine is given in the process of modeling. As the cost can be calculated by multiplying the material price by the weight [44], it is omitted here. Detailed calculations are included in the design goal of every modeling module. As illustrated in the constraints of every modeling module, the machine has electromagnetic, structural, and thermal features, and the constraints of the problem include multi-physics couplings inside and among components. They can be divided into two kinds, namely, explicit and implicit. The explicit constraints are those expressions in which the outputs can be represented by inputs without iterative calculations, e.g., Equation (1), whereas implicit constraints are those expressions that cannot, e.g., Equation (31). Explicit constraints can be substituted into the goal functions directly, whereas implicit ones cannot. Therefore, considering that the search algorithm has the ability to achieve optimal solutions, if the implicit constraints are not satisfied, the model outputs high penalty factors to break the re-design loop to simplify the calculation.

Model Improvements
Compared with other models, the developed multi-physics model is improved as follows: • The model is restructured and divided into several replaceable modules to illustrate the couplings of components. Thus, other machine types can be evaluated easily just by changes in shape, material, and radial positions of component modules, which also follows the concurrent design concept. • Design codes, e.g., the wire gauge, shape selection experiences, and empirical formulas are embedded to realize easier designs for system integrators.

•
Inputs are well selected according to the insulation requirements, which offer comprehensive links with practical design requirements, whereas penalty coefficients are used to break the closed loops of the design, which simplifies the calculation.
However, the dimensionality of the built model is greater than that of the model in [29]. This is because some design variables are not independent, e.g., component inductions. For the purpose of power system integration, this high-dimensional design space clearly leads to difficulties in the trade-off. Therefore, simplification is needed to improve the design efficiency.

Model Simplification
To further reduce the dimensionality of the built multi-physics model, this section explores its design space by sensitivity analysis methods and then establishes the simplified models used for system integration based on the model approximation technologies. To illustrate this process clearly, materials and empirical coefficients are pre-specified (the 45 carbon steel made shaft, the D41 cold-rolled silicon steel sheet, the fluoropolymer coating insulation, the ratio of the inner and outer diameter, and relative direct-axis inductance, etc.). Variation ranges of other undetermined design variables are listed in Table 1, according to the requirements of Class F insulation.
To implement the simplification, a data set, which can efficiently describe the whole performance of the multi-physics model, is required. To obtain the data set, the optimal Latin hypercube design (LHD) algorithm [46][47][48] is applied and the discrepancy of point sets is selected as the optimized criterion objective to generate a design matrix evenly spread in the design space presented in Table 1.
To facilitate an accurate analysis, a matrix of 10,000 points is used here. The process is demonstrated in Figure 6. With the increase of design points, the discrepancy of point sets declines evidently first and then remains steady after 800 iterations. Then, by substituting the design points into Equation (38), corresponding outputs of the multi-physics model are obtained. As constraints are not always satisfied, 990 unsuitable points are eliminated before the analysis. As every iteration generates nine points on average, the remaining 9010 design points can be regarded as the results at the 1001st iteration in Figure 6, which verifies the distribution performance of this filtered data set.

Sensitivity Analysis
To obtain the optimal parameter combination, this paper analyzes trends between inputs and outputs to identify key factors of the multi-physics model by the multiple quadratic regression method.
Equation (39) illustrates the principle of this method with a binary quadratic function. The differential of y can be divided into main and interaction effects of x1 and x2. For example, coefficients c1 and 2c3 can reflect the main effects of x1 on y; c2 and 2c4 can reflect the main effects of x2 on y, and c5 can reflect the interaction effect of x1 and x2 on y. Therefore, based on the multiple quadratic regression model, the influences of design variables on the weight and power loss can be evaluated: The determination coefficient (R 2 ) is used to evaluate the imitative effect of the multiple quadratic regression model. The R 2 of the weight and power loss models are 0.9705 and 0.97848, respectively, which indicates the accuracy of the analysis. Figure 7 illustrates the main factor plots of weight and power loss. The main factors of the weight WG are the rotational speed n, rated power P, pole pairs p, and electric loading A, whereas those of the power loss PG are the rated power P, electric loading A, pole pairs p, rotational speed n, and air-

Sensitivity Analysis
To obtain the optimal parameter combination, this paper analyzes trends between inputs and outputs to identify key factors of the multi-physics model by the multiple quadratic regression method.
Equation (39) illustrates the principle of this method with a binary quadratic function. The differential of y can be divided into main and interaction effects of x 1 and x 2 . For example, coefficients c 1 and 2c 3 can reflect the main effects of x 1 on y; c 2 and 2c 4 can reflect the main effects of x 2 on y, and c 5 can reflect the interaction effect of x 1 and x 2 on y. Therefore, based on the multiple quadratic regression model, the influences of design variables on the weight and power loss can be evaluated: The determination coefficient (R 2 ) is used to evaluate the imitative effect of the multiple quadratic regression model. The R 2 of the weight and power loss models are 0.9705 and 0.97848, respectively, which indicates the accuracy of the analysis. Figure 7 illustrates the main factor plots of weight and power loss. The main factors of the weight W G are the rotational speed n, rated power P, pole pairs p, and electric loading A, whereas those of the power loss P G are the rated power P, electric loading A, pole pairs p, rotational speed n, and air-gap induction B δ . Both weight and power loss rise linearly with the increase of the rated power. Moreover, with the rise of pole pairs, the weight declines whereas the power loss increases. In addition, it is notable that the voltage has few effects on weight and power loss. Due to different dimensions and ranges of variables, dimensionless processing is executed to achieve fair identification of couplings; the top 20 contributors of design variables and couplings are shown in Figure 8. The weight attributes 83.05% of itself to the couplings of the rated power P, pole pairs p, rotational speed n, electric loading A, air-gap flux density Bδ, and pitch ratio β, whereas the power loss attributes 79.35% of itself to the couplings of the rated power P, pole pairs p, rotational speed n, electric loading A, air-gap flux density Bδ, pitch ratio β, and power factor cosφ.  β-n Due to different dimensions and ranges of variables, dimensionless processing is executed to achieve fair identification of couplings; the top 20 contributors of design variables and couplings are shown in Figure 8. The weight attributes 83.05% of itself to the couplings of the rated power P, pole pairs p, rotational speed n, electric loading A, air-gap flux density B δ , and pitch ratio β, whereas the power loss attributes 79.35% of itself to the couplings of the rated power P, pole pairs p, rotational speed n, electric loading A, air-gap flux density B δ , pitch ratio β, and power factor cosϕ. Due to different dimensions and ranges of variables, dimensionless processing is executed to achieve fair identification of couplings; the top 20 contributors of design variables and couplings are shown in Figure 8. The weight attributes 83.05% of itself to the couplings of the rated power P, pole pairs p, rotational speed n, electric loading A, air-gap flux density Bδ, and pitch ratio β, whereas the power loss attributes 79.35% of itself to the couplings of the rated power P, pole pairs p, rotational speed n, electric loading A, air-gap flux density Bδ, pitch ratio β, and power factor cosφ. β-n According to the analysis above, some design variables can be ignored, and model inputs before and after the dimensionality reduction are demonstrated in Table 2. The dimensionalities of the weight and power loss models can be reduced from 15 to 6 and 7, respectively. As the input variables of the simplified model are the subset of those of the multi-physics model, they have a consistent interface. For a more comprehensive simplification, the pre-specified parameters can be added into the analysis.

Model Approximation
According to the reduced input variables in Table 2, the methods of radial basis function (RBF) and quartic response surface model (RSM) are used for the approximations of the data sample generated from DOE. Ten-fold cross-validations are executed by randomly partitioning the sample into nine training subsets to train the model and one test subset to evaluate it.
Comparisons between the built models and the test subsets are illustrated in Table 3. The RSMbased weight model has a relatively low accuracy and high degree of uncertainty, whereas other models achieve satisfactory results. As only a small amount of data is required by the approximation of the reduced-order quartic RSM-based model, the sample containing 9010 points cannot be covered. Accordingly, although ten-fold cross-validation helps keep models accurate by calculating multiple surrogate models of different subsets, the performance of reduced-order RSM-based models is worthless for large sample data. In other words, RBF-based surrogate models are preferred. Table 3. Cross-validation error results of surrogate models.

Models
Evaluation Indexes RBF RSM cosϕ

Effects on Power Loss P G (%)
A-n B δ p P-n According to the analysis above, some design variables can be ignored, and model inputs before and after the dimensionality reduction are demonstrated in Table 2. The dimensionalities of the weight and power loss models can be reduced from 15 to 6 and 7, respectively. As the input variables of the simplified model are the subset of those of the multi-physics model, they have a consistent interface. For a more comprehensive simplification, the pre-specified parameters can be added into the analysis. Table 2. Candidate variables for simplified models.

Models Alternative Inputs
Before W G P, p, A, n, B δ , β, cosϕ, B p , B ry , B st , B sy , J a , J f , λ, and E P G After W G P, p, A, n, B δ , and β P G P, p, A, B δ , n, β, and cosϕ

Model Approximation
According to the reduced input variables in Table 2, the methods of radial basis function (RBF) and quartic response surface model (RSM) are used for the approximations of the data sample generated from DOE. Ten-fold cross-validations are executed by randomly partitioning the sample into nine training subsets to train the model and one test subset to evaluate it.
Comparisons between the built models and the test subsets are illustrated in Table 3. The RSM-based weight model has a relatively low accuracy and high degree of uncertainty, whereas other models achieve satisfactory results. As only a small amount of data is required by the approximation of the reduced-order quartic RSM-based model, the sample containing 9010 points cannot be covered. Accordingly, although ten-fold cross-validation helps keep models accurate by calculating multiple surrogate models of different subsets, the performance of reduced-order RSM-based models is worthless for large sample data. In other words, RBF-based surrogate models are preferred.

Verification
To verify the modeling method, the auxiliary power unit generator (APUG) in a Boeing 787 with known parameters (apparent power P s = 225 kVA, electric loading A = 60,000 A/m, pole pairs p = 2, rotational speed n varies in the range of 10,800-24,000 r/min, nominal speed n N = 12,000 r/min, air-gap induction B δ+ = 0.8 T, pitch ratio β = 0.75, rated voltage E = 230 V, weight W G = 52 kg, power loss P G = 20 kW when n = 24,000 r/min and power factor cosϕ = 0.85) [49] is used as the case study.
According to the method, the multi-level models are established. As some inputs of the multi-physics model, e.g., operating points B p , B ry , B st , B sy , J a , and J f , and the dimension ratio λ, are not given, the uncertainties clearly result in a Pareto set. Therefore, according to ranges of the uncertain variables in Table 1, NSGA-II (non-dominated sorting genetic algorithm II) with 100 populations and 1000 generations is used, and the Pareto plot comprised of 3171 points obtained is illustrated in Figure 8. The design results of simplified models and the actual solution are also shown in Figure 9. Therein, two optimal solutions (M1 and M2) on the Pareto plot of the multi-physics model are given as examples for the convenience of the analysis.

Verification
To verify the modeling method, the auxiliary power unit generator (APUG) in a Boeing 787 with known parameters (apparent power Ps = 225 kVA, electric loading A = 60,000 A/m, pole pairs p = 2, rotational speed n varies in the range of 10,800-24,000 r/min, nominal speed nN = 12,000 r/min, airgap induction Bδ+ = 0.8 T, pitch ratio β = 0.75, rated voltage E = 230 V, weight WG = 52 kg, power loss PG = 20 kW when n = 24,000 r/min and power factor cosφ = 0.85) [49] is used as the case study.
According to the method, the multi-level models are established. As some inputs of the multiphysics model, e.g., operating points Bp, Bry, Bst, Bsy, Ja, and Jf, and the dimension ratio λ, are not given, the uncertainties clearly result in a Pareto set. Therefore, according to ranges of the uncertain variables in Table 1, NSGA-II (non-dominated sorting genetic algorithm II) with 100 populations and 1000 generations is used, and the Pareto plot comprised of 3171 points obtained is illustrated in Figure  8. The design results of simplified models and the actual solution are also shown in Figure 9. Therein, two optimal solutions (M1 and M2) on the Pareto plot of the multi-physics model are given as examples for the convenience of the analysis.

Multi-Physics Model Analysis
The actual solution is nearly located on the Pareto plot of the multi-physics model with the minimum errors (0.006% of the weight and 2.7% of the power loss, compared with solution M2),

Multi-Physics Model Analysis
The actual solution is nearly located on the Pareto plot of the multi-physics model with the minimum errors (0.006% of the weight and 2.7% of the power loss, compared with solution M2), which verifies the multi-physics model. The uncertainties lead the calculated weight to a variation ranging from −3.82% to 4.15% and the power loss to a variation ranging from −2.56% to 16.27% of the actual solution, which suggests that the weight of optimal solutions varies less than the power loss. Therefore, although the decision-making model is not studied here, more importance will be attached to the power loss by designers.

Dimensionality Reduction Analysis
To verify the dimensionality reduction of simplified models, data distributions of design variables of the optimal solutions are analyzed. As the power of the excitation system is relatively low and leads to few fluctuations of the results, their design variables are ignored here.
It is noteworthy that the Pareto set can be divided into two parts with the demarcation point M1 by the dimension ratio λ; see Figure 9. Figure 10 demonstrates the data distribution of the dimension ratio λ in the two parts. It can be concluded that a value of the dimension ratio λ exists in the range from 1.263 to 1.267 as the extreme point of the weight model, whereas it exists in the range from 1.591 to 1.593 as the extreme point of the power loss model. The selection of the dimension ratio is determined by the designer, depending on which index is more important. According to the power loss-preferred design, solutions of the 1728 blue-colored circle points in Figure 9 are probably better choices. which verifies the multi-physics model. The uncertainties lead the calculated weight to a variation ranging from −3.82% to 4.15% and the power loss to a variation ranging from −2.56% to 16.27% of the actual solution, which suggests that the weight of optimal solutions varies less than the power loss. Therefore, although the decision-making model is not studied here, more importance will be attached to the power loss by designers.

Dimensionality Reduction Analysis
To verify the dimensionality reduction of simplified models, data distributions of design variables of the optimal solutions are analyzed. As the power of the excitation system is relatively low and leads to few fluctuations of the results, their design variables are ignored here.
It is noteworthy that the Pareto set can be divided into two parts with the demarcation point M1 by the dimension ratio λ; see Figure 9. Figure 10 demonstrates the data distribution of the dimension ratio λ in the two parts. It can be concluded that a value of the dimension ratio λ exists in the range from 1.263 to 1.267 as the extreme point of the weight model, whereas it exists in the range from 1.591 to 1.593 as the extreme point of the power loss model. The selection of the dimension ratio is determined by the designer, depending on which index is more important. According to the power loss-preferred design, solutions of the 1728 blue-colored circle points in Figure 9 are probably better choices.  machine is well integrated. Therefore, they are local design variables, which can be ignored in the system-level design. In contrast, values of the current densities are spread relatively evenly in the full ranges given in Table 1, compared with other variables; see Figure 11e-f. The maximum changes of weight and power loss resulting from variations of the armature winding current density Ja and the excitation winding current density Jf are 6.88% and 1.91%, respectively, which suggests that current densities can also be ignored. In addition, it is notable that the models are convex functions of the current densities, as few data exist for the current  machine is well integrated. Therefore, they are local design variables, which can be ignored in the system-level design. In contrast, values of the current densities are spread relatively evenly in the full ranges given in Table 1, compared with other variables; see Figure 11e-f. The maximum changes of weight and power loss resulting from variations of the armature winding current density J a and the excitation winding current density J f are 6.88% and 1.91%, respectively, which suggests that current densities can also be ignored. In addition, it is notable that the models are convex functions of the current densities, as few data exist for the current densities of the armature and excitation windings in the ranges of 7.2-7.6 and 5-6.5 A/mm 2 , respectively.   In conclusion, among the undetermined design variables, λ has a relatively large influence on the weight and power loss, but can be ignored with a high-weighted coefficient of power loss considered in the design and evaluation, whereas component operating points (B p , B ry , B st , B sy , J a , and J f ) can also be reduced due to the optimal design results with few variations. Therefore, the dimensionality reduction is proved conclusive.

Simplified Model Analysis
Compared with the actual solution, the errors of the simplified weight and power loss models are acceptably low, at 1.06% and −2.03%, respectively. As the surrogate model is a 'black-box' model, the established models are of low complexity. Compared with the weight and power loss of simplified models, the relative variation ranges of the optimal solutions obtained from the multi-physics model on the right of M1 are −3.574-3.067% and −0.526-5.299% respectively, which indicates that the models output almost the same results. In summary, they are proved to be accurate, simple, and consistent.
In conclusion, the simplification is proved to be reliable, and the multi-physics model and simplified models are of high accuracy and consistency.

Comparisons
In order to demonstrate the superior performance of the proposed modeling method in an application of aircraft power systems, the proposed modeling method is compared with the baseline model of the permanent magnet synchronous motor (PMSG) in [29] and the SL-based model of PMSM in [30].
As the design of the whole aircraft power system is a multi-disciplinary problem resulting from strongly coupled parameters, to achieve the maximum benefit, global optimization solutions are the targets of designers, which indicates that all design variables should be considered at the system level. However, the large amount of on-board equipment leads system-level variables to a high-dimensional design space. Therefore, in addition to the precision, the dimensionality of design space is selected as the compared index to evaluate the performance of the model in system-level applications. As the calculation speed cannot be assessed in the same condition, the complexities of models are compared simply in theory.

Comparison with the Baseline Model
To model the PMSG, we replace the electrical excitation modules (silicon steel pole modules and excitation winding modules) by NdFeB magnet poles, and obtain the simplified PMSG model. Design variables are deduced and obtained combining parameters in [29] and equations in this paper, e.g., rated power P = 7.753 kW, pole pairs p = 3, rotational speed n = 735.2 r/min, electric loading A = 3.69 × 10 4 A/m, air-gap flux density B δ = 0.856T, and pitch ratio β = 2/3.
According to the improvements of the model mentioned in Section 3.8, the influences of the structure change on the accuracy can be ignored due to the modular modeling structure, whereas those of the integration of design codes and extended input variables remain to be considered, e.g., impacts of different empirical formulas for the calculation of air-gap lengths. As detailed information is not given, only the total model behaviors are compared; see Table 4. It can be noted that, compared with the prototype, there is a remarkable agreement between the results of the proposed simplified model and the corresponding models in [29].
The input number of the weight model in [29] is eight, whereas that of the simplified weight model in this paper is six, which indicates that a lower dimensionality of design space is needed to be considered in the system-level studies using the proposed method, leading to easier designs. Therefore, with the help of the proposed model, lower computational cost and faster development with considerable accuracy can be achieved.
Moreover, the 'black-box' characteristic of the proposed model indicates that its calculation speed is clearly superior to that of the baseline model. To model the PMSM, we change the radial positions of the pole module and slot/armature winding modules inside the PMSG model. As products in [50,51] are of the same series, one reference is selected to identify the variables of the built model, e.g., the electric loading A = 5.62 × 10 4 A/m, air-gap flux density B δ = 0.85 T, and pitch ratio β = 2/3; these values can be used for the whole series.
The relationship between the rated torque and weight is demonstrated in Figure 12. It can be concluded that, although the actual weight includes external accessories besides the active mass, e.g., the resolver and cable, both models achieve results with relatively high accuracy. To model the PMSM, we change the radial positions of the pole module and slot/armature winding modules inside the PMSG model. As products in [50,51] are of the same series, one reference is selected to identify the variables of the built model, e.g., the electric loading A = 5.62 × 10 4 A/m, airgap flux density Bδ = 0.85 T, and pitch ratio β = 2/3; these values can be used for the whole series.
The relationship between the rated torque and weight is demonstrated in Figure 12. It can be concluded that, although the actual weight includes external accessories besides the active mass, e.g., the resolver and cable, both models achieve results with relatively high accuracy. The SL-based method requires three inputs for the weight calculation, i.e., rated torque and weight of the reference machine, and rated torque of the machine to be calculated, which is less than the input number of our model. The power loss estimation requires 10 inputs, i.e., values of power and efficiency of the reference in two different conditions of loads (torque and speed) to calculate the coefficients of the Joule loss and iron loss, and torque and speed of the machine to be calculated. The total required input number of the SL-based method is greater than that of the proposed model. The SL-based method requires three inputs for the weight calculation, i.e., rated torque and weight of the reference machine, and rated torque of the machine to be calculated, which is less than the input number of our model. The power loss estimation requires 10 inputs, i.e., values of power and efficiency of the reference in two different conditions of loads (torque and speed) to calculate the coefficients of the Joule loss and iron loss, and torque and speed of the machine to be calculated. The total required input number of the SL-based method is greater than that of the proposed model.
The SL-based model is comprised of several simple polynomials and some data, whereas the proposed model is a relatively complex data network. Therefore, although they have a similar calculation speed, the occupied memory of the proposed model is higher than that of the SL-based model in the optimization process, which can be ignored compared with the generated solutions.
Moreover, the SL-based model cannot reflect the couplings of multi-physics, which cannot satisfy the study requirements of transport electrification.
Therefore, compared with other models used for the integrated design of aircraft power systems, the proposed method is demonstrated to be of similar accuracy, lower computational cost, and faster calculation speed, and promises to be a powerful tool in the multi-level design of aircraft power systems.

Conclusions
The objective of this paper was to develop an efficient modeling tool for machine design tasks in different design phases of aircraft power systems. For this purpose, the paper improved the existing multi-physics model and proposes a novel simplification process based on MDO. The major features of the proposed method are: 1.
The proposed modular modeling structure of the multi-physics model divides the machine into several replaceable sub-models, enabling the evaluation of different machine types by simply changing sub-models and their integration forms. Design codes are embedded in the model to reduce professional requirements for system integrators and to realize the combination of the principle and practical technologies. Moreover, this modeling method can also be applied to other equipment in the aircraft power system and enables efficient and quick designs with optimization algorithms.

2.
With the knowledge analysis of the improved multi-physics model, sensitive variables are selected from the multi-physics model as the inputs of the simplified models. The similar results of two-level models indicate that the nesting designs of different levels can be achieved with consistent results. Therefore, models of different levels are associated and consistent in aspects of interface and results, and the proposed modeling framework offers a systematic evaluation tool for designs of different levels.

3.
Compared with existing models for aircraft power system integration, the proposed simplified models of electrical machines have the advantages of lower-dimensional design space and black-box characteristics without reducing the precision. Both fast and accurate designs can be achieved by the method.
Future studies will focus on the simulation-based DOE to extract more detailed performance of machines for holistic evaluations. In addition, power allocations of multiple energy systems, the optimal operating voltage level of power systems, and couplings among More-electric systems (i.e., the gearbox and output frequency, fuel pumps, and cooling, etc.) are planned to be studied for future applications.