Hybrid vertical take-off and landing aircraft energy management method based on model predictive control algorithm

Abstract

The invention relates to a hybrid vertical take-off and landing aircraft energy management method based on a model predictive control algorithm, which belongs to the technical field of aircraft energy management, and comprises the following steps: establishing a size model and a performance model of a multi-rotor configuration based on a top demand, a hybrid architecture and a propulsion configuration by using an iterative scheme of size design optimization; forming constraint of size iteration, and optimizing required overall size parameters of the aircraft through an fmincon function; a component-level turboshaft engine model and a battery internal resistance model are provided, a theoretical model is established based on gas-heat dynamic characteristics of each component of the turboshaft engine, the battery internal resistance model is established based on battery testing, and a battery peak current is obtained; and an energy management strategy based on an MPC algorithm, DP optimization information, supplementary CS and a charge strategy.

Description

Technical Field

[0001] The present invention belongs to the field of aircraft energy management technology, and relates to a hybrid vertical take-off and landing aircraft energy management method based on a model predictive control algorithm. Background Art

[0002] For years, global dependence on fossil fuels has significantly exacerbated energy shortages and pollutant emissions. To achieve efficient energy use and zero emissions, technological research and development is underway at an increasingly intensive pace worldwide. Electrification of the automotive and aviation industries is considered an effective path to achieving these goals. Many vehicle and aircraft manufacturers have already adopted new propulsion technologies, including hybrid electric propulsion systems and pure electric propulsion systems. HEPS, by combining the advantages of two different energy sources, is an effective alternative propulsion system until battery technology matures. HEPS can be powered by a variety of combinations, including batteries and engines, batteries and fuel cells, and solar cells and batteries. HEPS offers significant advantages in improving fuel economy and reducing emissions. It also offers advantages in improving safety, reducing noise, and reducing the weight of aircraft batteries.

[0003] To achieve the aforementioned performance of HEPS, energy management plays a crucial role. HEPS energy management is defined as determining the most efficient way to distribute power / torque requirements among the various power sources. Generally speaking, the goal of HEPS energy management is to maximize powertrain efficiency and minimize fuel consumption. The algorithm used to achieve this goal is called an energy management system (EMS).

[0004] In recent years, driven by the demand for faster and more convenient travel, urban air mobility (UAM) scenarios have captured the attention of many visionaries worldwide and are considered a real opportunity for achieving sustainable transportation. Technological advances in propulsion and mechanical manufacturing have accelerated the gradual realization of VTOLs, making UAM scenarios possible. The introduction of VTOLs has further enriched the composition of three-dimensional transportation networks. Similarly, EMS plays a crucial role in maintaining the efficient operation of hybrid VTOLs powered by HEPS. Current EMS research for hybrid VTOLs is still in its early stages, and two factors complicate this issue. First, the driving conditions of hybrid VTOLs are highly complex, encompassing multiple operating modes, such as climb, cruise, hover, and vertical takeoff and landing. The energy demand characteristics of these different operating modes vary significantly, making accurate energy demand prediction particularly important. Second, the architecture of the hybrid VTOL propulsion system dictates the complexity of its energy flow distribution. In the HEPS of a hybrid VTOL, the power requirements for these different operating modes are supplied by the same battery pack and engine generator set. This power mismatch results in higher fuel consumption and a more limited flight range. Therefore, energy flow distribution for hybrid VTOLs remains a significant challenge, requiring the exploration of guiding solutions. It can be seen that both hybrid electric vehicles and hybrid electric aircraft can provide references for hybrid VTOL design. Regarding takeoff and landing configurations, fixed-wing and rotary-wing aircraft can provide a reference for the design of lift and propulsion systems for hybrid VTOLs. In the area of ​​energy management, EMSs from the automotive and aerospace fields can also be applied to hybrid VTOLs. For example, when developing an EMS for a hybrid electric aircraft, it is necessary to comprehensively consider both high-power phases, such as takeoff, and low-power phases, such as cruising, to allocate energy flow. Furthermore, energy flow control schemes for emergency failure scenarios are also necessary. Hybrid VTOLs also have significant variations in power requirements during operation, necessitating safeguard mechanisms. Therefore, the EMS of hybrid electric aircraft provides valuable insights for developing an EMS for hybrid VTOLs. Summary of the Invention

[0005] In view of this, the purpose of the present invention is to provide a hybrid vertical take-off and landing aircraft energy management method based on a model predictive control algorithm, which effectively avoids the problems of SOC exceeding the limit and frequent start and stop of APU while improving the real-time control capability, thereby ensuring fuel economy.

[0006] In order to achieve the above object, the present invention provides the following technical solutions:

[0007] A hybrid vertical take-off and landing aircraft energy management method based on a model predictive control algorithm comprises the following steps:

[0008] S1: Based on the top-level requirements of the aircraft, a multidisciplinary optimization method is used to build a performance and size model for the hybrid multi-rotor VTOL. On this basis, a hybrid VTOL power system-level design framework is established to optimize the optimal size parameters.

[0009] S2: Based on the optimization results, further complete the selection of the hybrid system. By combining mechanistic modeling with experimental modeling, a turboshaft engine component-level model and a power battery system model based on an equivalent circuit are established, providing data support for the formulation and optimization of energy management strategies.

[0010] S3: Based on the constructed hybrid system model and established test conditions, a rule-based energy management strategy was designed. The dynamic programming (DP-EMS) global optimization method was introduced to construct an optimization model with SOC as the state variable and fuel economy as the target. A complete energy management strategy based on the DP algorithm was designed, and a diagram of the maximum probable output power of the APU was obtained through simulation.

[0011] S4: A prediction model is constructed with SOC as the state variable and fuel economy as the optimization target. Combined with the constrained discretization method, the global optimization characteristics and online adaptive capabilities of DP are integrated. CS and charge strategies are introduced to provide solutions for energy management under complex working conditions, and the working condition verification strategy is updated.

[0012] The beneficial effects of the present invention are: It designs and verifies a hybrid VTOL energy management strategy based on a turboshaft engine system and a power battery system, taking into account the power battery's capacity limits, the turboshaft engine's performance, and the real-time nature of the strategy, ensuring fuel economy and real-time control capabilities for flight missions. Specifically, the advantages include the following:

[0013] (1) By integrating the size model and performance model, the optimal size parameters of the hybrid multi-rotor VTOL are optimized and iterated.

[0014] (2) A refined component-level turboshaft engine model was established. The PSFC diagram obtained by simulation had a small error with the Gasturb software. The power battery model was relatively simple, but a highly reliable power battery peak current curve was obtained through multiple rate tests, providing sufficient data support for the formulation and optimization of energy management strategies.

[0015] (3) A rule-based energy management strategy (RB-EMS) was designed to achieve power allocation between the battery and the range extender system by setting the SOC threshold. At the same time, a dynamic programming (DP-EMS) global optimization method was introduced to construct an optimization model with SOC as the state variable and fuel economy as the target, taking into account both power demand and system constraints.

[0016] (4) An MPC-based energy management strategy was designed. A prediction model was established with SOC as the state variable and APU power as the control input. Global optimal power allocation within the local time domain was achieved through dynamic programming and rolling optimization mechanisms. By extracting the global optimization characteristics of the DP strategy and introducing strategies similar to CS and Charge in the RB strategy, real-time control capabilities were significantly improved. While ensuring fuel economy, the improved strategy effectively avoided the problems of SOC exceeding the limit and frequent APU starts and stops.

[0017] Other advantages, objects, and features of the present invention will be described in part in the following description and, in part, will be apparent to those skilled in the art upon examination of the following description or may be learned from practice of the present invention. The objects and other advantages of the present invention may be realized and obtained through the following description. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] In order to make the purpose, technical solutions and advantages of the present invention more clear, the present invention will be described in detail below with reference to the accompanying drawings, in which:

[0019] Figure 1 A simplified logic diagram for the hybrid VTOL energy management strategy of the present invention;

[0020] Figure 2 Designed for route conditions;

[0021] Figure 3 The appearance of the selected model;

[0022] Figure 4 The hybrid system structure diagram;

[0023] Figure 5 Iterate the flow chart for rotor weight;

[0024] Figure 6 It is an overall schematic diagram of the aircraft selection process;

[0025] Figure 7 is the compressor characteristic diagram;

[0026] Figure 8 is the turbine characteristic diagram;

[0027] Figure 9 This is the PSFC diagram of the turboshaft engine;

[0028] Figure 10 This is the battery internal resistance model diagram;

[0029] Figure 11 This is the peak current fitting result diagram;

[0030] Figure 12 Design diagrams for flight sub-mission profiles;

[0031] Figure 13 Power diagram required for test conditions;

[0032] Figure 14 is the APU maximum probability power map;

[0033] Figure 15 Required power diagram for new test conditions;

[0034] Figure 16 The SOC change diagram under three strategies;

[0035] Figure 17 The cumulative fuel consumption diagram under three strategies. DETAILED DESCRIPTION

[0036] The following describes the embodiments of the present invention by means of specific examples, and those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that the illustrations provided in the following embodiments are only schematic illustrations of the basic concept of the present invention, and the following embodiments and features in the embodiments can be combined with each other without conflict.

[0037] It should be noted that the illustrations provided in the following embodiments are merely schematic illustrations of the basic concept of the present invention. Therefore, the illustrations only show components related to the present invention and are not drawn according to the number, shape, and size of components in actual implementation. In actual implementation, the type, quantity, and proportion of each component may be changed arbitrarily, and the component layout may also be more complex.

[0038] In the following description, numerous details are discussed to provide a more thorough explanation of the embodiments of the present invention. However, it will be apparent to those skilled in the art that the embodiments of the present invention may be practiced without these specific details. In other embodiments, well-known structures and devices are shown in block diagram form rather than in detail to avoid obscuring the embodiments of the present invention.

[0039] Example 1:

[0040] like Figure 1 As shown, the present invention provides a hybrid vertical take-off and landing aircraft energy management method based on a model predictive control algorithm, comprising the following steps:

[0041] S101: Based on the top-level requirements of the aircraft, a multidisciplinary optimization method was used to construct a performance and size model for a hybrid multi-rotor VTOL. Based on this model, a hybrid VTOL powertrain system-level design framework was established to optimize the optimal size parameters.

[0042] S102: Based on the optimization results, further complete the selection of the hybrid system. By combining mechanism modeling with experimental modeling, a turboshaft engine component-level model and a power battery system model based on an equivalent circuit are established, providing data support for the formulation and optimization of energy management strategies.

[0043] S103: Establish test conditions. Design a rule-based energy management strategy. Then, introduce the dynamic programming (DP-EMS) global optimization method to build an optimization model with SOC as the state variable and fuel economy as the target. Establish an energy management strategy based on the DP algorithm and obtain the APU maximum probability output power relationship diagram.

[0044] S104: A prediction model is constructed with SOC as the state variable and fuel economy as the optimization target. The constrained discretization method is used to integrate the global optimization characteristics and online adaptive capabilities of DP. Strategies similar to CS and Charge in the RB strategy are introduced. An MPC algorithm-based whole-machine energy management strategy is designed to provide a solution with greater engineering practice potential for energy management under complex working conditions and update the working condition verification.

[0045] Furthermore, the step S101 specifically includes the following steps:

[0046] S1011: Clarify top-level requirements, aircraft types, flight conditions, hybrid architecture and other related parameters;

[0047] S1012: Construct dimensional models of components including fuselage, rotors, turboshaft engines, thermal management systems, landing gear, and power distribution lines;

[0048] S1013: Based on momentum theory and empirical formulas, establish power calculation models for multi-rotor configurations in hovering, vertical take-off and landing, and cruise phases;

[0049] S1014: Based on the previous dimensional and performance models, an overall design framework based on multidisciplinary design optimization is proposed to obtain the optimal dimensional parameters of a hybrid multirotor with a design range of 100 km and a hybrid factor of 0.3.

[0050] Further, the flight mission profile of step S1011 is as follows: Figure 2 As shown, the appearance of the model is as follows Figure 3 As shown, the hybrid architecture is Figure 4 The flight conditions are shown in Table 1:

[0051] Table 1

[0052]

[0053] The selected parameters that affect the size of the aircraft are shown in Table 2:

[0054] Table 2

[0055]

[0056]

[0057] The selected parameters that affect the size of the aircraft are shown in Table 3:

[0058] Table 3

[0059]

[0060] Furthermore, the steps for building the fuselage size model in step S1012 are as follows:

[0061] (1) Skin weight estimation

[0062] The fuselage skin is the main structure of the fuselage exterior, and its weight is directly related to its surface area. This paper uses the ellipsoidal surface area approximation formula to calculate the outer surface area S of the skin wet :

[0063]

[0064] Where, L f ,W f and H f are the length, width and height of the fuselage respectively. The skin weight is calculated by the outer surface area S wet and unit area weight ρ areal The product of gets:

[0065] m skin =S wet ·ρ areal

[0066] Where, ρ areal is the weight per unit area, which is determined by the minimum thickness t of the skin material bid , core material thickness t core and coating thickness t paint and its density calculated.

[0067] (2) Bulkhead weight estimation

[0068] The bulkhead is a transverse support component inside the fuselage used to enhance the structural strength. This paper assumes that the bulkhead has a circular cross-section. Based on the bulkhead's geometric shape and its unit area weight, its weight m bulkhead Estimated by the following formula:

[0069]

[0070] (3) Canopy weight estimation

[0071] The weight of the canopy as part of the fuselage is calculated by multiplying one-eighth of the outer surface area by the canopy's thickness and its density:

[0072]

[0073] Where, t canopy is the canopy thickness; ρ canopy is the canopy density.

[0074] (4) Keel weight estimation

[0075] The keel is the main load-bearing structure of the fuselage, used to withstand bending, torsion and landing loads during flight. The calculation of keel weight is also divided into the keel weight under the three requirements of bending load, torsion load and landing load:

[0076] m keel =m bend +m torsion +m land

[0077] Keel weight m under bending load requirement bend , assuming that the maximum bending moment M b From the maximum lift F L produce:

[0078] F L =n g W tow ·sf

[0079]

[0080] Where n g is the maximum overload factor; W tow is the weight of the aircraft; sf is the safety factor. Keel cross-sectional area A b Calculation via bending moment and material stress:

[0081]

[0082] Where h beam is the keel height; σ uni is the material stress. The weight of the keel under bending load is obtained by multiplying the cross-sectional area, length and material density:

[0083] m bend =A b ·L f·ρ uni

[0084] Keel weight m under torsional load requirements torsion , assuming that the keel is subjected to a torque M t By the lift F L The relationship with span is:

[0085]

[0086] Keel cross-sectional area A t Calculated from the keel height and keel width, the keel height is the fuselage height h beam One tenth of the keel width w beam Take one third of the fuselage width and thickness t t Through the torque M t and material shear strength τ bid calculate:

[0087] A t =h beam w beam

[0088]

[0089] The keel weight required for torsional load can be obtained by multiplying the cross-section perimeter, thickness and material density:

[0090] m torsion =2·(h beam +w beam )·t t ·ρ bid

[0091] For landing load, the keel weight m land , assuming the lateral force F during landing land The aircraft weight and landing load factor n l produce:

[0092]

[0093] Keel cross-sectional area A under landing load requirements l and thickness t land By the lateral force F land Calculate the volume V of the reinforcement material pad By thickness t land It turns out that:

[0094]

[0095]

[0096] Where d is the bolt diameter; σbearing is the material bearing strength. The weight of the keel is obtained by multiplying the reinforcement volume and the material density:

[0097] m land =4·V pad ·ρ bid

[0098] In summary, the total weight of the fuselage is m fuselage It can be expressed as the sum of the weight of the skin, bulkhead, canopy and keel:

[0099] m fuselage =m skin +m bulkhead +m canopy +m keel

[0100] The steps for building the rotor blade size model in step S1012 are as follows:

[0101] (1) Geometric modeling

[0102] The geometric parameters of the rotor blades are defined as follows: the blade radius is R, the chord length is c = 0.1R, and the number of blades is n b =3. Using NACA4 series airfoil, its thickness distribution is described by the formula:

[0103]

[0104] Where, t avg is the average relative thickness; x is the chord position coordinate. The airfoil coordinate system is based on the shear center x shear =0.25c is the origin, the leading edge shear web is arranged in the interval [0.25c, 0.35c], and the length of the root assembly section is L root =0.1R.

[0105] (2) Load distribution model

[0106] The thrust is distributed quadratically along the blade extension direction, which is mathematically described as:

[0107]

[0108] Where s f The safety factor is taken as 1.5; T is the maximum thrust. This distribution form is derived from the blade element momentum theory and reflects the radial variation characteristics of the induced velocity at the propeller disk.

[0109] The centrifugal load is generated by the angular velocity of rotation:

[0110]

[0111] Where M tipis the blade tip Mach number; a is the speed of sound.

[0112] (3) Structural strength analysis

[0113] The entire blade adopts a composite laminated structure design, and the main components include:

[0114] ① Skin: bears torsional load, made of bidirectional carbon fiber.

[0115] ② Leading edge reinforcement: bears bending load, made of unidirectional carbon fiber.

[0116] ③ Shear web: bears shear load, made of the same material as the skin.

[0117] ④Core material: polyurethane foam.

[0118] The key strength check formula is as follows:

[0119] Torsional stiffness is provided by the closed cell skin, and the required thickness is calculated using Bredt's formula:

[0120]

[0121] Where A e is the cross-sectional area of ​​the closed chamber; τ allow is the allowable shear stress of bidirectional carbon fiber. When the calculated value is lower than the process limit t min = 0.5mm to enable the minimum thickness constraint.

[0122] The bending stress is borne by the leading edge hat beam, and its thickness design must meet the following requirements:

[0123]

[0124] The first item M f Corresponding to the aerodynamic bending moment; the second term F c Represents centrifugal force. c is the moment of inertia of the hat beam; L c is the hat beam length. Accurate calculations are performed using the discrete integration method to ensure that the stress assessment error is less than 1%.

[0125] The shear strength is borne by the shear web, and the thickness is designed by the shear force V generated by the lifting z calculate:

[0126]

[0127] (4) Weight iterative calculation

[0128] Based on the above formula model, an iterative convergence algorithm is used to calculate the rotor blade weight that meets the requirements. Since the weight distribution affects the centrifugal load, the following method can be used: Figure 5The iterative process is shown in Figure 2. Numerical verification shows that the algorithm converges within five iterations, achieving computational efficiency that meets engineering requirements. The final weight incorporates a 20% process correction factor to account for minor factors such as adhesive layers and corrosion protection coatings. By decomposing the load path and establishing an explicit stress-weight mapping, the model achieves rapid and accurate weight estimation.

[0129] The turboshaft engine, thermal management system, and landing gear size models in step S1012 are as follows:

[0130] The size model and fuel consumption model of the turboshaft engine can be expressed by the following formula:

[0131] M=0.625·(P MC +200) 0.8

[0132]

[0133] The weight of the thermal management system is mainly composed of three parts: the weight of the battery water cooling system, the weight of the water cooling system for the motor and other electrical components, and the weight of the cold-end forced convection cooling rotor. First, the weight of the battery water cooling system is approximately proportional to the maximum output power of the power battery system. Referring to the statistical analysis of the water cooling system of electric vehicle power batteries, this article uses 0.75kW / kg as the power-to-weight ratio. Secondly, the weight of the water cooling system for the motor and other electrical components is related to the required heat dissipation power. This article takes the heat dissipation power-to-weight ratio as 0.68kW / kg. Finally, the weight of the cooling rotor is related to the maximum electrical power consumed during the flight mission. Referring to the technical parameters of some cooling rotor products, this article uses 0.2kW / kg as the power consumption-to-weight ratio.

[0134]

[0135] Where m TMS is the selected weight of the thermal management system; P Bs , Q es 、P TMSprop They are the maximum output power of the battery, the heat dissipation power of other electrical components, and the maximum power consumption of the cooling rotor; f BTMS 、Qwratio eTMS 、f TMSprop They are the ratio of the maximum output power of the battery to the weight of the water cooling system, the ratio of the heat dissipation power of other electrical components to the weight of the water cooling system, and the ratio of the power consumption of the cooling rotor to its weight.

[0136] Based on the relevant literature on helicopter landing gear, this patent assumes that the weight of the landing gear for both configurations is 2% of the maximum takeoff weight:

[0137] m lg =0.02 MTOW

[0138] Another major component that cannot be ignored is the weight of the power distribution cables. This patent, referring to relevant literature, established a power distribution cable model based on the total power flowing to each motor and the motor's position within the aircraft, providing a rough estimate of the distribution cable weight. Furthermore, the anti-collision seat weighs 15 kg. The electric actuator weight is assumed to be 0.65 kg, with a typical multirotor configuration having eight. Finally, an additional 10% margin is added to the total weight to account for accessories and miscellaneous hardware.

[0139] Furthermore, the performance model building steps of step S1013 are as follows:

[0140] Hovering power P of multi-rotor configuration MR,hover It can be expressed as:

[0141] P MR,hover =P i,hover +P profile,hover

[0142] Where, P MR,hover Indicates the hovering power of the multi-rotor configuration; P i,hover Indicates hovering induced power; P profile,hover Indicates the resistive power. The induced power can be expressed as follows:

[0143]

[0144] In the formula, the load on the action disk is simplified to the gravity W of the takeoff weight of the aircraft tow ; A disk is the effective propeller disc area. The cruising speed in the hovering phase is zero, and the drag power can be expressed as follows:

[0145]

[0146] During the flight of a multi-rotor aircraft, the change in the rotor angle of attack will cause the aircraft to adjust its attitude and change its motion state. Considering that the maximum angle of attack of this design is controlled within 12°, the aerodynamic interference effect it produces can be ignored. The shaft power P in the cruise state MR,cruise It consists of four core elements: the induced power component P i,cruise , resistance loss power P profile,cruise , parasitic resistance power P parasite,cruise and vertical maneuvering power P w,climb The vertical maneuvering power represents the energy demand when the altitude changes, and this parameter automatically returns to zero under constant altitude cruise conditions.

[0147] P MR,cruise =P i,cruise +P profile,cruise +P parasite,cruise +Pw

[0148] The induced power during the cruise process can be approximately calculated based on the hovering induced power, which is expressed as:

[0149]

[0150] Where v is the magnitude of the cruising speed vector. Forward flight increases the air resistance experienced by the rotating blades. The drag power during cruising can be calculated using the following empirical formula:

[0151]

[0152] The magnitude of parasitic power is related to the forward flight equivalent flat plate area F x Related:

[0153]

[0154] The calculation of climb / descent power is relatively simple. In this case, the speed is substituted into the vertical speed during cruise:

[0155] P w,climb =W tow v climb

[0156] P w,descent =W tow v descent

[0157] The vertical takeoff power of a multirotor aircraft is composed of induced power, climb power and drag power, which can be expressed as:

[0158] P MR,takeoff =P i,takeoff +P w,takeoff +P profile,takeoff

[0159] Where, P i,takeoff represents the vertical takeoff induced power, which can be expressed as:

[0160]

[0161] Where, v takeoff Indicates vertical takeoff speed, vertical takeoff induced speed v i,takeoff The calculation basis is as follows:

[0162]

[0163] Where, v c is the climbing speed, v i,hover is the hovering induced speed. Climbing power P w,takeoff It can be expressed as follows:

[0164] P w,takeoff =W tow v takeoff

[0165] Similarly, the vertical landing power consumption can be expressed as:

[0166] P MR,landing =P i,land +P w,land +P profile,land

[0167]

[0168] P w,land =W tow v landing

[0169] Where, v landing Indicates vertical descent speed, and its value is negative.

[0170] During the vertical takeoff and landing phase, the aircraft's flight speed is relatively low, and the drag power can be considered the same as during hovering:

[0171] P profile,takeoff =P profile,land =P profile,hover

[0172] Furthermore, the aircraft design optimization iteration process of step S1014 is as follows: Figure 6 shown

[0173] The optimization objective is defined as the minimization of MTOW and can be expressed as the following function:

[0174] f=min(MTOW)

[0175] The optimization variable X vector is represented as:

[0176] X=[r prop m battery m motors MTOW m fuel ]

[0177] Where r prop is the rotor radius of the aircraft; m battery is the weight of the battery pack; m motors is the weight of the motor group; MTOW is the maximum take-off weight of the aircraft; m fuel is the total weight of the fuel.

[0178] The optimization problem also involves a series of nonlinear constraints, which are implemented through the conun function. The constraints mainly reflect the physical and engineering limitations that the aircraft must meet during design and operation, such as the minimum value of the rotor radius and the upper limit of the battery weight. The constraint function can be expressed as:

[0179] c(X)≤0(0.1)

[0180] Where c(X) is a nonlinear inequality constraint. Specific constraints include:

[0181] 1) Energy Constraint

[0182] To ensure that the total energy storage of the aircraft hybrid system covers the required range, the following constraints are imposed on the battery weight and fuel weight while considering different HF values:

[0183] c1=HF·E reserve -m battery ·E Battery DoD η b ≤0 (0.2)

[0184] c2=(1-HF)·E reserve -m fuel ·E fuel ·η e ≤0 (0.3)

[0185] Where, E battery is the battery energy density; E fuel is the fuel energy density, which is 12000Wh / kg here, and DoD is the depth of discharge, which is 0.95; η b is the total efficiency of battery energy transmission; η e is the total efficiency of fuel energy transmission, which is determined by the assumed parameters in the previous article; E reserve To consider the total energy required for the total task in the reserve phase, it can be expressed as follows:

[0186] E reserve =P hover ·t hover +P cruise ·t cruise +P loiter ·t loiter (0.4)

[0187] Where, P hover and P hover It can be obtained from the performance model formula above. The input variable of the performance model is r prop and MTOW;t cruiseCalculated based on the cruise speed and design range given above; the reserve phase power is set equal to the cruise power, and the reserve phase mission duration t loiter Set to 20 minutes.

[0188] 2) Battery power constraints

[0189] To ensure flight safety, the aircraft's required shaft power should be less than the maximum motor power, and the required battery power should also be less than the maximum power the battery can provide:

[0190] c3=P max ·HF-m battery ·η b ·P battery ≤0 (0.5)

[0191] c4=P max -m motors ·η em ·P motor ≤0 (0.6)

[0192] Where, P battery is the battery power density; P motor is the motor power density; η em is the motor efficiency, which is given by the parameters assumed above. For multi-rotor configuration, P max is the maximum power value required by the aircraft in each stage. For the lift + cruise configuration, it can be expressed by the following formula:

[0193]

[0194] T LC,max =T LC,hover ·R maxT / W (0.8)

[0195] Where R maxT / W In order to design the maximum thrust-to-weight ratio, this paper takes the value of 1.6.

[0196] 3) Maximum takeoff weight constraints

[0197] Since the margin factor has been fully considered when constructing the size and weight model in the previous article, the maximum takeoff weight constraint can be expressed as follows:

[0198] c5=m mass -MTOW≤0 (0.9)

[0199] Where m mass The total weight of the entire aircraft obtained by adding up all the dimensional models established above is multiplied by the margin factor; MTOW is the maximum take-off weight that meets the performance model constraints.

[0200] 4) Rotor kinetic energy constraints

[0201] In the design process of multi-rotor aircraft, rotor kinetic energy constraints must be considered, which is one of the important conditions for meeting design specifications. By introducing rotor kinetic energy constraints, it can ensure that the design of multi-rotor aircraft meets relevant safety and performance standards, and improve its reliability and applicability. In the optimization problem, the rotor kinetic energy constraint can be expressed as follows:

[0202]

[0203] Where m rotor is the weight of the rotor; v tip is the speed of the rotor tip; v auto The spin-down speed can be obtained from the total weight W tow and rotor radius r prop Find:

[0204]

[0205] In order to solve the above optimization problem, this study adopted the fmincon function based on MATLAB, combined with multiple random restart strategies to ensure that the optimization process can effectively avoid local optimal solutions and thus find the global optimal solution.

[0206] In this study, HF, or the hybrid factor, is defined as the ratio of the available storage energy of the power battery system to the storage energy required for the entire mission profile considering the reserve phase:

[0207]

[0208] From the optimization results, a multirotor aircraft with a range of 100 km and a mixing factor of 0.3 was selected for energy management strategy design. The optimization result data of the aircraft are shown in Table 4.

[0209] Table 4

[0210]

[0211] Furthermore, the step S102 specifically includes the following steps:

[0212] S1021: Establish component-level model of turboshaft engine;

[0213] S1022: Establish a battery internal resistance model;

[0214] Furthermore, the modeling method of step S1021 is as follows:

[0215] In the process of building a turboshaft engine model, aerodynamic thermodynamic calculations play a vital role. Before starting the modeling work of a single-shaft turbojet gas turbine, the properties of the working fluid in its internal gas path must be clarified. The working fluid is mainly air in the inlet and compressor, and after mixing with fuel and igniting in the combustion chamber, it is converted into combustion gas, which then expands to generate work. Theoretically, the definition of air and combustion gas in gas turbine modeling mainly involves three key parameters, namely, specific heat C, p , enthalpy h and entropy S. The specific calculation formulas of these parameters are as follows:

[0216] Air:

[0217]

[0218] Gas:

[0219]

[0220]

[0221] In the above formula, θ cp ,θ h and θ s can be regarded as correction coefficients, and are all functions of temperature T. i 、b i 、c i d i and e i are the known polynomial fitting coefficients. These functions allow for a detailed description of the gas thermodynamics of a turboshaft engine. For example, during the compression or expansion of engine gas, the enthalpy and entropy differences between the inlet and outlet cross-sections of a component can be calculated to determine the performance parameters of that cross-section.

[0222] ①Intake duct model:

[0223] The goal of the inlet model is to calculate the temperature, pressure and flow rate of the external air to ensure the flow rate gain and ensure that the inlet system meets the flight conditions. The input parameters of the inlet are the flight altitude H and the Mach number Ma inlet , the output is the temperature, pressure, flow rate and static pressure of the air inlet.

[0224] (1) Atmospheric conditions

[0225] First, assume that the altitude H of the aircraft has a significant impact on the atmospheric temperature and pressure. According to the International Standard Atmosphere Model, the temperature and pressure calculation formulas at different altitudes are as follows:

[0226] For flight altitudes less than 11 km, the temperature and pressure calculation formulas are:

[0227] T S0=T air0 -6.5H

[0228]

[0229] Where, T air0 =288.15K is the sea level temperature, and H is the altitude (km).

[0230] (2) Inlet temperature and pressure

[0231] Based on the temperature and pressure of the atmosphere, the temperature and pressure at the air inlet can be calculated using the following formula:

[0232] The calculation formula for the air inlet temperature T1 is:

[0233]

[0234] Where, k = 1.402 is the adiabatic index of air; Ma inlet is the Mach number.

[0235] The formula for calculating the pressure at the air inlet is:

[0236]

[0237] (3) Air inlet flow rate

[0238] The flow rate through the inlet is calculated using the Mach number and the standard gas constant using the formula:

[0239]

[0240] Where R air =287 is the gas constant of air.

[0241] ② Compressor model

[0242] (1) Calculation of speed coefficient

[0243] Speed ​​coefficient n LC,HPC The calculation formula is as follows:

[0244]

[0245] Where N H is the input speed, rpm; T2 is the input temperature, K; N cons is the speed at the design point.

[0246] (2) Obtaining compressor characteristic data

[0247] In this model, the pressure ratio and efficiency corresponding to different speeds and flows are characteristic data derived from the Gasturb tool. The compressor characteristic diagram is shown in Figure 7In the compressor characteristic diagram, the black solid line represents the constant speed line, the black dashed line represents the efficiency contour line, and the red dashed line is the surge line. The data in the diagram is derived and interpolated to obtain the compressor performance under specific conditions.

[0248] (3) Interpolation calculation

[0249] In this model, two-dimensional interpolation is used to estimate compressor performance. The independent variables for interpolation include the pressure ratio and the speed factor. Given characteristic data such as pressure ratio, flow rate, and efficiency at different speeds, these data are interpolated to determine actual performance indicators, such as flow rate and efficiency, under actual operating conditions.

[0250] During the interpolation process, the corresponding performance data point is first found based on the given pressure ratio and speed coefficient. Then, the flow rate and efficiency are calculated using bilinear interpolation, and the interpolation coefficient ckzx is calculated. This coefficient is used to represent the position of the current pressure ratio between two adjacent data points:

[0251]

[0252] Where x j is the current pressure ratio; x a and x b are adjacent data points. Then use the linear interpolation formula to calculate the corresponding flow and efficiency:

[0253] G LC =G LCa +ckzx·(G LCb -G LCa )

[0254] η HPC =η HPCa +ckzx·(η HPCb -η HPCa )

[0255] Where G LCa and G LCb is the flow rate at the interpolation point; η HPCa and η HPCb is the efficiency of the interpolation point.

[0256] (4) Compressor output flow

[0257] Compressor output flow G HPC It can be calculated by the following formula:

[0258]

[0259] Where P2 is the input pressure; T2 is the input temperature; 101325 is the standard atmospheric pressure (unit: Pa)

[0260] (5) Isentropic compressor outlet temperature

[0261] Based on the assumption of an isentropic process, the compressor outlet temperature T3 can be calculated by the following steps:

[0262] Calculate the input enthalpy and the outlet enthalpy for an isentropic process:

[0263] H2=f t2hc (T2,0)

[0264] H c,id =f t2hc (T c,id ,0)

[0265] According to the efficiency of the compressor, the ideal outlet enthalpy is corrected to obtain the actual outlet enthalpy H3:

[0266]

[0267] Calculate the actual outlet temperature T3:

[0268] T3=f h2tc (H3,0)

[0269] Where, f h2tc 、f t2hc They are functions for calculating the corresponding temperature based on enthalpy value and the corresponding enthalpy value based on temperature, respectively.

[0270] (6) Calculation of power and flow output

[0271] To maintain compressor stability, the model also considers cooling and bleed air flows. Cooling flow primarily cools the vanes and blades, typically accounting for 5% of the total flow. Bleed air is used for sealing and other purposes and typically accounts for 1% of the total flow. The compressor's power and flow output can be calculated using the following formulas:

[0272] W HPC =(H3-H2)·G HPC (1-0.05-0.01)

[0273] G3=G HPC (1-0.05-0.05-0.01)

[0274] ③Combustion chamber model

[0275] (1) Mass flow and pressure loss in the combustion chamber

[0276] Total output flow of the combustion chamber G out is the sum of the air flow and fuel flow entering the combustion chamber:

[0277] G out =G in +Gfuel

[0278] At the same time, the pressure drop caused by factors such as friction and heat loss of the air flow in the combustion chamber is taken into account. The outlet pressure P of the combustion chamber out There is usually a certain amount of loss, and the model assumes that the outlet pressure is 97% of the inlet pressure:

[0279] P out =P in 0.97

[0280] (2) Calculation of enthalpy of combustion chamber inlet gas

[0281] Enthalpy of the gas at the combustion chamber inlet H in Through its temperature T in Calculated using the function between enthalpy and temperature:

[0282] H in =f t2hc (T in ,0)

[0283] Where 0 means air is an ideal gas.

[0284] (3) Calculation of air-fuel ratio

[0285] The fuel-air ratio is defined as the ratio of fuel flow to air flow:

[0286]

[0287] (4) Calculation of enthalpy of mixed gas

[0288] According to the enthalpy of the intake gas, the lower calorific value of the fuel and the combustion efficiency, the enthalpy of the mixed gas at the combustion chamber outlet is H out The calculation is as follows:

[0289]

[0290] Where H fuel =8.6×10 5 J / kg is the calorific value of fuel, which indicates the total energy per unit mass of fuel. 7 It is the lower calorific value of fuel, which indicates the heat released when unit mass of fuel is completely burned.

[0291] (5) Calculation of combustion chamber outlet temperature

[0292] According to the total enthalpy H of the mixed gas out And the fuel air-fuel ratio FAR4, calculate the temperature T at the combustion chamber outlet out :

[0293] T out=f h2tc (H out ,FAR4)

[0294] ④ Hybrid cavity model

[0295] (1) Calculation of enthalpy of input gas

[0296] The modeling of the mixing chamber first requires the calculation of the enthalpy of the input gas. The fuel gas H4 is calculated from its temperature:

[0297] H4=f t2hc (T4,FAR4)

[0298] Where T4 is the temperature of the fuel gas; FAR4 is the fuel air-fuel ratio.

[0299] (2) Calculation of output gas flow rate and mixing enthalpy

[0300] Total output flow of the mixing chamber G out The sum of the air flow rate and the fuel gas flow rate G4, plus the cooling air flow rate G cooling :

[0301] G out =G4+G cooling

[0302] Air flow G air It is the total flow rate of the mixed gas minus the flow rate of the fuel gas:

[0303] G air =G out -G fuel

[0304] Output fuel air-fuel ratio FAR out Indicates the amount of fuel contained in a unit mass of air, which is the ratio of the fuel gas flow rate to the air flow rate:

[0305]

[0306] Enthalpy of the gas mixture H mix is the weighted average of air and fuel gas taking into account their respective mass flow rates:

[0307]

[0308] (3) Calculation of mixed gas temperature

[0309] Temperature of the mixed gas T mix It can be expressed by its total enthalpy H mix and fuel air-fuel ratio FAR out Calculation yields:

[0310] T mix=f h2tc (H mix ,FAR out )

[0311] ⑤Turbine model

[0312] (1) Calculation of turbine expansion ratio and speed coefficient

[0313] First, the expansion ratio Pr is defined according to the expansion process of the turbine HPT , and its calculation formula is:

[0314]

[0315] Where, P 41 is the turbine inlet pressure; P 45 is the turbine outlet pressure.

[0316] Next, define the reduced speed coefficient n LC，HPT , which takes into account the relationship between intake temperature and turbine speed. The calculation formula of the reduced speed coefficient is as follows:

[0317]

[0318] Where N H is the actual turbine speed; T 41 is the turbine inlet temperature; δ is the reduced temperature coefficient; N cons represents the design point speed.

[0319] (2) Turbine characteristic data and interpolation

[0320] In this model, these characteristic data are derived from Gasturb, including the turbine performance under different speed and pressure ratio conditions. Figure 8 As shown in the turbine characteristic diagram, the black solid line represents the constant speed line, and the black dashed line represents the efficiency contour line. Through interpolation calculation, the flow rate and efficiency under specific operating conditions can be estimated.

[0321] Similar to the compressor method, the flow rate G is calculated based on the interpolation coefficient ckzx LC and efficiency η HPT :

[0322] G LC =G LCa +ckzx·(G LCb -G LCa )

[0323] η HPT =η HPTa +ckzx·(η HPTb -η HPTa )

[0324] Where G LCa and G LCb is the flow rate at the interpolation point; η HPTa and η HPTb is the efficiency of the interpolation point.

[0325] The interpolation result is flow G LC and efficiency η HPT , according to the compressor outlet pressure P 41 and intake air temperature T 41 , the actual turbine flow G can be further calculated HPT :

[0326]

[0327] Where, P 41 is the turbine inlet pressure; T 41 is the turbine inlet temperature.

[0328] (3) Power output and efficiency calculation

[0329] Turbine output enthalpy H T,out It is the enthalpy value corrected according to the actual efficiency. First calculate the turbine inlet enthalpy H T,0 and ideal outlet enthalpy H T,id , and then correct the enthalpy value H according to the efficiency of the turbine T,out :

[0330] H T,0 =f t2hc (T 41 ,FAR out )

[0331] H T,id =f t2hc (T T,id ,FAR out )

[0332] H T,out =(H T,id -H T,0 )·η HPT +H T,0

[0333] Where η HPT is the efficiency of the turbine.

[0334] Through the actual outlet enthalpy H T,out , the turbine output temperature T can be calculated 45 :

[0335] T 45 =f h2tc (H T,out ,FAR out )

[0336] The turbine's power output W HPT It is determined by the difference between the outlet and inlet enthalpies of the turbine and the flow rate and can be calculated using the following formula:

[0337] W HPT =-(H T,out -H T,0 )·G HPT

[0338] ⑥PSFC solution

[0339] In this study, the turboshaft engine specific fuel consumption (PSFC) calculation method is based on the flow and energy balance of each engine component. First, the turboshaft engine output power is calculated using the known output torque and engine speed. The formula is:

[0340]

[0341] Where, T out is the output torque, in Nm; N H is the engine speed in rpm. This formula converts torque and speed into power output to calculate specific fuel consumption. Specific fuel consumption (PSFC) is a key indicator of engine performance, indicating the amount of fuel consumed per unit of power. The PSFC calculation formula is as follows:

[0342]

[0343] Where G fuel is the fuel flow rate, in kg / h; P out is the engine's output power in kW. The formula calculates the specific fuel consumption per unit power by the ratio of fuel flow to power output.

[0344] In order to ensure the matching of the performance of each engine component, error balance calculations were performed. The errors mainly involved the following two aspects:

[0345] (1) Power error

[0346] The power error reflects the deviation between the power output from the turbine and the theoretically calculated power. The calculation formula is:

[0347]

[0348] Where W HPT is the output power of the high-pressure turbine; W HPC is the outlet power of the high-pressure compressor; η rotor is the turbine rotor efficiency; P outis the calculated engine output power. This error measures the energy conversion efficiency between the turbine and compressor.

[0349] (2) Flow error

[0350] The flow error reflects the balance of the engine airflow, especially the flow balance between the mixing chamber and the cooling airflow. Its calculation formula is:

[0351]

[0352] Where G HPC is the flow rate at the high-pressure turbine outlet; 0.1G HPC is the cooling air flow rate; G fuel is the fuel flow rate; G3 is the flow rate at the high-pressure compressor outlet. This error is used to ensure the balance between the cooling airflow, fuel flow rate, and compressor outlet flow rate. Through the above calculation, the final error balance consists of two parts:

[0353]

[0354] These two errors reflect the balance between power and flow and help ensure that the engine components work together under design conditions.

[0355] Through the above calculations, this study established a simulation calculation method based on a turboshaft engine. Based on the thermodynamic model of the turboshaft engine, the PSFC characteristic diagram of the turboshaft engine can be obtained by setting the speed range and torque range. Figure 9 As shown in the figure, the PSFC value of the design point is 331.8g / kw*h. Compared with the data obtained by Gasturb, the error is smaller, which proves that the thermodynamic model can better reflect the fuel consumption characteristics of the turboshaft engine.

[0356] Furthermore, the modeling method of step S1022 is as follows:

[0357] Establish an equivalent internal resistance model such as Figure 10 shown.

[0358] Output power P of the power battery b It is expressed as follows:

[0359]

[0360] Where V oc is the open circuit voltage of the power battery; R int is the internal resistance of the battery; I b is the charge and discharge current of the battery; I b It can be expressed as follows:

[0361]

[0362] The battery SOC can be determined by the integration method, which is specifically expressed as follows:

[0363]

[0364] Where SOC0 represents the initial SOC value of the battery; Q b Represents the battery capacity. Differentiating this equation yields:

[0365]

[0366] When P b When it is positive, the power battery is discharging, otherwise, it is charging.

[0367] Next, parameter identification is performed through static capacity test, open circuit voltage test and HPPC test to obtain the peak current corresponding to different SOC and discharge rate.

[0368] The peak current can be expressed as follows:

[0369]

[0370] The peak current surface is obtained by fitting the fourth-order polynomial fitting method. Figure 11 shown.

[0371] Furthermore, the step S103 specifically includes the following steps:

[0372] S1031: Establish test conditions;

[0373] S1032: Design an energy management strategy based on the DP algorithm and simulate to obtain a diagram showing the maximum probability output power of the APU.

[0374] Furthermore, the test condition of step S1031 is established as follows;

[0375] The test task condition contains three identical subtasks, and the subtask profiles are as follows: Figure 12 As shown, each subtask is mainly divided into hovering, cruise climb, cruise and cruise descent. The route mileage of the entire flight subtask is 10km. The second and third subtasks are started after the previous stage is completed and the aircraft descends to an altitude of 15m. After completing all subtasks, the aircraft begins to descend vertically and land. The mission power requirements of the entire aircraft are as follows: Figure 13 shown.

[0376] Furthermore, the strategy design of step S1032 is as follows;

[0377] The specific control rules of the rule-based energy management strategy are as follows:

[0378] (1)CD mode

[0379] When SOC>SOC CS When the system is in CD mode, the system is mainly driven by the battery. The change of SOC directly affects the discharge state of the battery. The system determines whether to continue to maintain battery drive mode or switch to CS mode by monitoring SOC.

[0380] (2) CS mode

[0381] When SOC min <SOC≤SOC CS When the system is in flight, it will activate engine drive mode to prevent the battery from being overly discharged. In this mode, the engine provides power to the aircraft, and the battery is used only for auxiliary purposes. This mode can balance the power requirements of the battery and engine, reducing dependence on a single energy source.

[0382] (4) Charge mode

[0383] When SOC≤SOC min In this mode, the power provided by the engine is used to charge the battery to restore the battery's SOC.

[0384] Furthermore, the strategy design of step S1033 is as follows;

[0385] The design of the energy management strategy based on the DP algorithm first determines the state variable as SOC. The SOC value changes in each stage and is affected by the battery charging and discharging power in the current stage. The decision variable in each stage is the APU power generation power P APU By controlling the power to match the aircraft's power requirements, the aircraft's power needs are met, and the battery's SOC value will also change. The choice of power generation needs to balance fuel consumption and SOC changes. That is:

[0386] x k =SOC k

[0387] u k =P APU (k)

[0388] The calculation and update methods of SOC are as described above. The DP cost function is composed of multiple factors, mainly including fuel consumption and SOC deviation. Fuel consumption is related to the power generation of the engine and can be calculated using the following formula:

[0389]

[0390] Where, f(P APU(i)) represents the relationship between APU power generation and fuel consumption. The cost function of the DP algorithm is expressed as follows:

[0391]

[0392] Where N is the mission step size, which is set to 1032 in this paper. The terminal penalty coefficient α=10 6 , which is used to ensure the convergence of the SOC final value. The overall limitations of the system are as follows:

[0393]

[0394] Where, P dmd Indicates power demand; (P APU , P b ) represent the output power of the APU and power battery system, respectively. The SOC is discretized to obtain a state variable matrix. The DP algorithm first traverses each state variable at each time node in reverse order to obtain the optimal control table. It then continues the forward solution to derive the target control sequence. The inverse process is expressed as follows:

[0395] J(k,i)=min{m F (k,i,j)+J(k+1,j)}

[0396] k=0,1,...,N; i,j=0,1,...,M

[0397] Where J(k+1,j) represents the slave node SOC k+1,j To the final node SOC N Minimum fuel consumption; m F (k,i,j) represents the battery from SOC k,i To SOC k+1,j The instantaneous fuel consumption changes; M is the number of SOC discrete grids.

[0398] When the initial SOC and target SOC are determined, the global optimization result of the DP strategy is fixed. From the optimization results, the optimal power allocation at each required power in each starting state and expected ending state can be statistically calculated. If the starting SOC, target ending SOC, required power and current SOC state are used as input to construct the relationship between the current state and the APU output power, this high-dimensional mapping is too redundant and the mapping relationship is difficult to establish. The data volume is large, and the engine output power may fluctuate greatly at similar current SOCs, which is not conducive to engine control. In this paper, the starting SOC, target ending and required power are used as inputs, and the characteristic parameter of the current SOC is discarded. The APU output power with the highest probability of occurring under the required power is used as the target power of the required power, and then the power is used as the target power at each moment in the prediction time domain, thereby calculating the SOC reference trajectory in the prediction time domain. When the initial SOC is set to 0.8 and the target ending is 0.35, the corresponding relationship between the required power and the APU output power is as follows: Figure 14 shown.

[0399] Furthermore, the step S104 specifically includes the following steps:

[0400] S1041: Use the APU maximum probability output power relationship diagram obtained by DP to design an energy management strategy based on the MPC algorithm;

[0401] S1042: Introducing CS and charge modes to improve the energy management strategy;

[0402] Furthermore, the energy management strategy of step S1041 is formulated as follows;

[0403] Based on the MPC principle, the battery SOC is selected as the state variable and the APU output power is selected as the control variable. A model predictive control strategy for the hybrid power system is designed to optimize fuel economy. The specific steps are as follows:

[0404] (1) Establishing a prediction model: Since the operating power of the aircraft is relatively fixed, this paper assumes that the required power remains unchanged within the prediction time domain k~k+p. That is, the required power P(k) at time k is used as the required power at each moment in the prediction time domain [P(k+1), P(k+2),…, P(k+p)]. On this basis, the upper and lower limits of the SOC in the prediction time domain are estimated to determine the feasible domain of the SOC.

[0405] (2) Solving the optimization problem: In the prediction time domain k~k+p, the hybrid system power distribution optimization objective function is established. Under the system constraints, the APU maximum probability output power diagram obtained in the previous step is used to solve the optimal control in the SOC feasible domain and obtain the optimal control sequence [u(k),u(k+1|k),…,u(k+p|k)].

[0406] (3) Optimal control implementation: In the MPC strategy, instead of applying the optimal control sequence obtained at the current moment to the time series k~k+p, only the first element u(k) of the control sequence is applied to the current moment and solved, and then the system state is updated.

[0407] As mentioned above, a hybrid system can be described as:

[0408]

[0409] in:

[0410]

[0411] In this process, the optimization is carried out in the interval k to k+p, so it is necessary to establish the objective function in this interval, as shown in the following formula:

[0412]

[0413] p is the length of the prediction time domain, m f is the engine fuel consumption, SOC(t) and SOC ref is the SOC at time t and the reference value; ω is the penalty coefficient.

[0414] At the same time, the following constraints need to be met:

[0415]

[0416] Where max and min are the upper and lower limits of the corresponding quantity.

[0417] The MPC solution process is roughly as follows:

[0418] (1) Application of DP results in MPC

[0419] The solution process can be expressed as:

[0420] ① Divide the prediction time domain k~k+p into p+1 stages, i.e. [k, k+1,…, k+p];

[0421] ② Discretize the feasible domain of SOC prediction within k~k+p;

[0422] ③ Use the APU maximum probability output power map obtained in the previous step to optimize the objective function and obtain [P APU (k),P APU (k+1|k),P APU (k+2|k),…,P APU (k+p|k)], and apply the first element to the controlled system to update the state.

[0423] (2) SOC range determination and discretization

[0424] When performing optimization solutions, the SOC range needs to be limited. On the one hand, this is to avoid excessive use of the battery to reduce fuel consumption, resulting in a too-low SOC at the end of the trip; on the other hand, a reasonable SOC range can also reduce the amount of MPC calculations. This requires calculating the SOC range in the prediction time domain based on the current SOC and required power, so as to further narrow the SOC range.

[0425] At any time, the following power constraints need to be met:

[0426] P(k)=P APU (k)+P b (k)

[0427] Then the maximum output power of the battery is:

[0428] P b,max (k)=min(P(k),P b,dismax (SOC(k)))

[0429] The minimum output power of the battery is:

[0430] P b,min (k)=max(P(k)-P APU,max ,P b,chgmax (SOC(k)))

[0431] Taking SOC(k) as the starting point and P b,max (k) is the SOC lower boundary power, P b,min (k) is the upper boundary power of SOC, and the upper and lower boundaries of SOC in the entire prediction time domain can be calculated.

[0432] Where, P b,dismax (SOC(k)) represents the maximum discharge power of the battery under SOC(k); and P b,chgmax (SOC(k)) represents the maximum charging power of the battery under SOC(k); P b,max (k) and P b,min (k) is the maximum and minimum power of the battery under the required power P(k).

[0433] (3) Determination of control variable range

[0434] The previous step describes how to reduce the amount of MPC calculations by using the feasible region of SOC. In addition, the amount of MPC calculations can also be reduced by limiting the range of control variables. Under the condition of required power P(k), the corresponding P APU scope of work.

[0435] The maximum output power of the APU is:

[0436] P APU,max (k)=min(P(k),P APU,max )

[0437] The minimum output power is:

[0438] P APU,min (k)=max(P(k)-P b,dismax (SOC(k)),0)

[0439] Where, P APU,max (k) and P APU,min (k) is the maximum and minimum power under the required power P(k). On this basis, the discrete sequence of the control variable can be obtained according to a certain discrete interval.

[0440] Furthermore, the energy management strategy of step S1042 is improved as follows;

[0441] Due to the uncertainty of the operating conditions, under the above SOC trajectory constraints, the SOC may not drop to the target SOC at the end of the operating condition, and may hit the bottom earlier. Therefore, after the SOC hits the bottom, a strategy similar to CS and Charge is added to ensure that the SOC is close to the target end SOC and that frequent engine starts and stops do not occur. The details are as follows:

[0442] (1) When SOC < SOC target When -0.02, the APU outputs at maximum data power;

[0443] (2) When SOC is reduced to SOC target After -0.02, it returned to SOC target +0.02 to enable MPC optimization again;

[0444] (3) When SOC target -0.01<SOC<SOC target When the value is +0.01, the APU output power is determined by the required power. That is, when the required power P is less than 60kW, the APU does not work and the battery provides all the energy. When the required power P is greater than 60kW, the APU outputs power at its maximum.

[0445] (4) Not within the control range of the above SOC, that is, SOC target -0.02<SOC<SOC target -0.01 and SOC target +0.01<SOC<SOC targetWhen the SOC is +0.02, the control boundary of SOC is involved. Unreasonable distribution will cause SOC to oscillate at the control boundary, resulting in frequent start and stop of APU. Therefore, under this condition, APU is allowed to recharge at maximum power output to avoid frequent start and stop.

[0446] To verify the effectiveness of the MPC strategy, the hovering time was changed to 90s and the cruising altitude was changed to 500m. The new test conditions were as follows: Figure 15 As shown in the figure, the simulation verification was carried out again. Compared with the original working condition, the power duration distribution of this working condition has also changed accordingly. The power duration of hovering, cruise ascent and cruise descent has become longer, and the cruise power time has been reduced accordingly. This can avoid the leakage of working condition information, thereby indirectly improving the reference of the MPC-OPT control effect.

[0447] Furthermore, the energy management strategy described in step S1041 is combined with the improved method and new test conditions in step S1042, and the fuel consumption and SOC comparison chart is finally obtained by simulation. Figure 16 and 17 The fuel consumption comparison is shown in Table 5.

[0448] Table 5

[0449]

[0450] Example 2:

[0451] An electronic device comprising a memory and a processor;

[0452] The memory is used to store computer programs;

[0453] The processor is configured to implement the method described in Example 1 when executing the computer program.

[0454] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not limiting. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention can be modified or replaced by equivalents without departing from the purpose and scope of the technical solutions, which should all be included in the scope of the claims of the present invention.

Claims

1. A hybrid vertical take-off and landing aircraft energy management method based on a model predictive control algorithm, characterized by: The following steps are involved: S1: Based on the top-level requirements of the aircraft, a multidisciplinary optimization method is used to build a performance and size model for the hybrid multi-rotor VTOL. On this basis, a hybrid VTOL power system-level design framework is established to optimize the optimal size parameters. S2: Based on the optimization results, further complete the selection of the hybrid system. By combining mechanistic modeling with experimental modeling, a turboshaft engine component-level model and a power battery system model based on an equivalent circuit are established, providing data support for the formulation and optimization of energy management strategies. S3: Based on the constructed hybrid system model and established test conditions, a rule-based energy management strategy was designed. The dynamic programming (DP-EMS) global optimization method was introduced to construct an optimization model with SOC as the state variable and fuel economy as the target. A complete energy management strategy based on the DP algorithm was designed, and a diagram of the maximum probable output power of the APU was obtained through simulation. S4: A prediction model is constructed with SOC as the state variable and fuel economy as the optimization target. Combined with the constrained discretization method, the global optimization characteristics and online adaptive capabilities of DP are integrated. CS and charge strategies are introduced to provide solutions for energy management under complex working conditions, and the working condition verification strategy is updated.

2. The energy management method for a hybrid vertical take-off and landing aircraft based on a model predictive control algorithm according to claim 1, characterized in that: The step S1 specifically includes the following steps: S11: Clarify top-level requirements, aircraft types, flight conditions, and hybrid architecture parameters; S12: Build a dimensional model of the fuselage, rotor blades, turboshaft engine, thermal management system, landing gear, and power distribution lines; S13: Based on momentum theory and empirical formulas, build performance models, including power calculation models for hovering, vertical take-off and landing, and cruise phases; S14: Based on the size model and performance model, an overall design framework based on multidisciplinary design optimization is proposed to obtain the optimal size parameters of the hybrid multi-rotor.

3. The energy management method for a hybrid vertical take-off and landing aircraft based on a model predictive control algorithm according to claim 2, characterized in that: The steps for constructing the fuselage size model in step S12 are as follows: (1) Skin weight estimation The outer surface area S of the skin is calculated using the ellipsoid surface area approximation formula wet : Where, L f ,W f and H f are the length, width and height of the fuselage respectively; the skin weight is calculated by the outer surface area S wet and unit area weight ρ areal The product of gets: m skin =S wet ·r areal Where, ρ areal is the weight per unit area, which is determined by the minimum thickness t of the skin material bid , core material thickness t core and coating thickness t paint and its density is calculated; (2) Bulkhead weight estimation Assuming the bulkhead has a circular cross-section, based on the bulkhead's geometry and its weight per unit area, its weight m bulkhead Estimated by the following formula: (3) Canopy weight estimation The canopy weight is calculated by multiplying one-eighth of the outer surface area by the canopy thickness and its density: Where, t canopy is the canopy thickness; ρ canopy is the canopy density; (4) Keel weight estimation The calculation of keel weight is divided into three parts: bending load, torsional load and landing load. m keel =m bend +m torsion +m land Keel weight m under bending load requirement bend , assuming that the maximum bending moment M b From the maximum lift F L produce: F L =n g ·W tow ·sf Where n g is the maximum overload factor; W tow is the weight of the aircraft; sf is the safety factor; keel cross-sectional area A b Calculation via bending moment and material stress: Where h beam is the keel height; σ uni is the material stress; the keel weight under bending load is obtained by multiplying the cross-sectional area, length and material density: m bend =A b ·L f ·r uni Keel weight m under torsional load requirements torsion , assuming that the keel is subjected to a torque M t By the lift F L The relationship with span is: Keel cross-sectional area A t Calculated from the keel height and keel width, the keel height is the fuselage height h beam One tenth of the keel width w beam Take one third of the fuselage width and thickness t t Through the torque M t and material shear strength τ bid calculate: A t =h beam w beam The keel weight required for torsional loading is obtained by multiplying the section perimeter, thickness, and material density: m torsion =2·(h beam +w beam )·t t ·ρ bid For landing load, the keel weight m land , assuming the lateral force F during landing land The aircraft weight and landing load factor n l produce: Keel cross-sectional area A under landing load requirements l and thickness t land By the lateral force F land Calculate the volume V of the reinforcement material pad By thickness t land It turns out that: Where d is the bolt diameter; σ bearing is the material bearing strength; the keel weight is obtained by multiplying the reinforcement volume and the material density: m land =4 V pad ·r bid In summary, the total weight of the fuselage is m fuselage Expressed as the sum of the weight of the skin, bulkhead, canopy and keel: m fuselage =m skin +m bulkhead +m canopy +m keel 。 4. The energy management method for a hybrid vertical take-off and landing aircraft based on a model predictive control algorithm according to claim 2, characterized in that: The steps for constructing the rotor blade size model are as follows: (1) Geometric modeling The geometric parameters of the rotor blades are defined as follows: the blade radius is R, the chord length is c = 0.1R, and the number of blades is n b =3; the rotor thickness distribution is described by the formula: Where, t avg is the average relative thickness; x is the chord position coordinate; the airfoil coordinate system is based on the shear center x shear =0.25c is the origin, the leading edge shear web is arranged in the interval [0.25c, 0.35c], and the length of the root assembly section is L root =0.1R; (2) Load distribution model The thrust is distributed quadratically along the blade extension direction, which is mathematically described as: Where s f is the safety factor; T is the maximum thrust; The centrifugal load is generated by the angular velocity of rotation: Where M tip is the blade tip Mach number; a is the speed of sound; (3) Structural strength analysis The entire blade adopts a composite laminate structure design, and the components include: ① Skin: bears torsional load, made of bidirectional carbon fiber; ② Leading edge reinforcement: bears bending load, made of unidirectional carbon fiber; ③ Shear web: bears shear load, made of bidirectional carbon fiber; ④ Core material: polyurethane foam; Key strength check formulas include: Torsional stiffness is provided by the closed cell skin, and the required thickness is calculated using Bredt's formula: Where A e is the cross-sectional area of ​​the closed chamber; τ allow is the allowable shear stress of bidirectional carbon fiber; when the calculated value is lower than the process limit t min =0.5mm, the minimum thickness constraint is enabled; The bending stress is borne by the leading edge hat beam, and its thickness design must meet the following requirements: The first item M f Corresponding to the aerodynamic bending moment; the second term F c I represents centrifugal force; c is the moment of inertia of the hat beam; L c is the hat beam length; The shear strength is borne by the shear web, and the thickness is designed by the shear force V generated by the lifting z calculate: (4) Weight iterative calculation An iterative convergence algorithm is used to calculate the rotor blade weight that meets the requirements, including the following steps: Initialize weight distribution m (0) ; Calculating centrifugal force Updated structural dimensions and weight m (k+1) ; Check convergence conditions |m (k+1) -m (k) |<10 -8 Is it satisfied? If not, continue to update the structural dimensions and weight. If slow, output the rotor weight m; The final weight includes a 20% process correction factor.

5. The energy management method for a hybrid vertical take-off and landing aircraft based on a model predictive control algorithm according to claim 2, characterized in that: In step S12, the size models of the turboshaft engine, thermal management system, landing gear, and power distribution line are constructed as follows: The size model and fuel consumption model of the turboshaft engine are expressed by the following formula: M=0.625·(P MC +200) 0.8 The weight of the thermal management system consists of three parts: the weight of the battery water cooling system, the weight of the water cooling system for the motor and other electrical components, and the weight of the cold-end forced convection cooling rotor. The weight of the battery water cooling system is approximately proportional to the maximum output power of the power battery system; the weight of the water cooling system for the motor and other electrical components is related to the required heat dissipation power; and the weight of the cooling rotor is related to the maximum power consumption during the flight mission. Where m TMS is the selected weight of the thermal management system; P Bs , Q es 、P TMSprop They are the maximum output power of the battery, the heat dissipation power of other electrical components, and the maximum power consumption of the cooling rotor; f BTMS 、Qwratio eTMS 、f TMSprop They are the ratio of the maximum output power of the battery to the weight of the water cooling system, the ratio of the heat dissipation power of other electrical components to the weight of the water cooling system, and the ratio of the power consumption of the cooling rotor to its weight; Assuming the landing gear weight for both configurations to be 2% of the maximum takeoff weight: m lg =0.02·MTOW Distribution line weight: Build a distribution line model based on the total power flowing to each motor and the motor's location in the aircraft to estimate the distribution line weight; It also includes the weight of the anti-collision seat and the weight of the electric actuator; Add extra allowance when calculating total weight to account for accessories and miscellaneous hardware.

6. The energy management method for a hybrid vertical take-off and landing aircraft based on a model predictive control algorithm according to claim 2, characterized in that: The steps for constructing the power calculation model for the hovering, vertical take-off and landing, and cruise phases in step S13 are as follows: Hovering power P of multi-rotor configuration MR,hover Expressed as: P MR,hover =P i,hover +P profile,hover Where, P MR,hover Indicates the hovering power of the multi-rotor configuration; P i,hover Indicates hovering induced power; P profile,hover The induced power is expressed as follows: In the formula, the load on the action disk is simplified to the gravity W of the takeoff weight of the aircraft tow ; A disk is the effective propeller disc area; the cruising speed in the hovering phase is zero, and the drag power is expressed as follows: Shaft power P in cruising state MR,cruise It consists of four core elements: the induced power component P i,cruise , resistance loss power P profile,cruise , parasitic resistance power P parasite,cruise and vertical maneuvering power P w,climb ; The vertical maneuvering power represents the energy demand when the altitude changes, and this parameter automatically returns to zero under the constant altitude cruise condition; P MR,cruise =P i,cruise +P profile,cruise +P parasite,cruise +P w The induced power during the cruise process is approximately calculated based on the hovering induced power, which is expressed as: Where v is the magnitude of the cruising speed vector; The drag power during cruising is calculated using the following empirical formula: The magnitude of parasitic power is related to the forward flight equivalent flat plate area F x Related: Calculate the vertical speed during the cruise process when climbing or descending: P w,climb =W tow v climb P w,descent =W tow v descent The vertical takeoff power of a multirotor aircraft is composed of induced power, climb power and drag power, which can be expressed as: P MR,takeoff =P i,takeoff +P w,takeoff +P profile,takeoff Where, P i,takeoff represents the vertical takeoff induced power, expressed as: Where, v takeoff Indicates vertical takeoff speed, vertical takeoff induced speed v i,takeoff The calculation basis is as follows: Where, v c is the climbing speed, v i,hover is the hovering induced speed; Climb power P w,takeoff It is expressed as follows: P w,takeoff =W tow v takeoff The vertical landing power consumption is expressed as: P MR,landing =P i,land +P w,land +P profile,land P w,land =W tow v landing Where, v landing Indicates vertical landing speed, its value is negative; The drag power of the aircraft during vertical takeoff and landing is the same as when hovering: P profile,takeoff =P profile,land =P profile,hover 。 7. The energy management method for a hybrid vertical take-off and landing aircraft based on a model predictive control algorithm according to claim 2, characterized in that: The overall design framework based on multidisciplinary design optimization described in step S14 has an optimization goal of minimizing MTOW, which can be expressed as the following function: f=min(MTOW) The optimization variable X vector is represented as: X=[r prop m battery m motors MTOW m fuel ] Where r prop is the rotor radius of the aircraft; m battery is the weight of the battery pack; m motors is the weight of the motor group; MTOW is the maximum take-off weight of the aircraft; m fuel is the total weight of the fuel; The constraints of the optimization problem include: 1) Energy Constraint Considering different HF values, the following constraints are imposed on battery weight and fuel weight: c1=HF·E reserve -m battery ·IN Battery ·DoD·η b ≤0 c2=(1-HF)·E reserve -m fuel ·E fuel ·η e ≤0 Where, E battery is the battery energy density; E fuel is the fuel energy density, DoD is the depth of discharge; η b is the total efficiency of battery energy transmission; η e is the overall efficiency of fuel energy transmission; E reserve The total energy required to consider the total task during the reserve phase is expressed as: E reserve =P hover ·t hover +P cruise ·t cruise +P loiter ·t loiter Where, P hover and P hover Calculated by the performance model formula, the input variable of the performance model is r prop and MTOW;t cruise It is obtained based on the cruising speed and design range; the reserve phase power is set to be equal to the cruising power, and the reserve phase mission duration t is set loiter ; 2) Battery power constraints The required shaft power of the aircraft should be less than the maximum power of the motor, and the required battery power should be less than the maximum power the battery can provide: c3=P max ·HF-m battery ·η b ·P battery ≤0 c4=P max -m motors ·η em ·P motor ≤0 Where, P battery is the battery power density; P motor is the motor power density; η em is the motor efficiency, which is given by the assumed parameters; for the multi-rotor configuration, P max is the maximum power value required by the aircraft in each stage; for the lift + cruise configuration, it is expressed as follows: T LC,max =T LC,hover ·R maxT / W Where R maxT / W is the designed maximum thrust-to-weight ratio; 3) Maximum takeoff weight constraints c5=m mass -MTOW≤0 Where m mass The total weight of all dimension models built for the entire aircraft is multiplied by the margin factor; MTOW is the maximum takeoff weight that meets the performance model constraints; 4) Rotor kinetic energy constraints Where m rotor is the weight of the rotor; v tip is the speed of the rotor tip; v auto The spin-down speed can be obtained from the total weight W tow and rotor radius r prop Find The hybrid factor HF is defined as the ratio of the stored energy available in the power battery system to the stored energy required for the entire mission profile taking into account the reserve phase: Select qualified multi-rotor aircraft from the optimization results for energy management strategy design.

8. The energy management method for a hybrid vertical take-off and landing aircraft based on a model predictive control algorithm according to claim 1, characterized in that: The step S2 specifically includes the following steps: S21: Establish a turboshaft engine component-level model and obtain the turboshaft engine's PSFC characteristic diagram through refined modeling, including: The definition of air and gas in gas turbine modeling involves three key parameters, namely, specific heat C p , enthalpy h and entropy S, the specific calculation formula is as follows: Air: Gas: Where θ cp ,θ h and θ s is the correction coefficient, and both are functions of temperature T; a i 、b i 、c i d i and e i are the known polynomial fitting coefficients; These functions describe in detail the gas thermodynamic calculation process of the turboshaft engine: ①Intake duct model: The goal of the inlet model is to calculate the temperature, pressure and velocity of the external air; the input parameters of the inlet are the flight altitude H and the Mach number Ma inlet , the output is the temperature, pressure, flow rate and static pressure of the air inlet; (1) Atmospheric conditions Assuming that the aircraft's flight altitude H has a significant impact on atmospheric temperature and pressure, according to the International Standard Atmosphere Model, the temperature and pressure calculation formulas for different altitude ranges are as follows: For flight altitudes less than 11 km, the temperature and pressure calculation formulas are: T S0 =T air0 -6.5H Where, T air0 =288.15K is the sea level temperature, H is the altitude; (2) Inlet temperature and pressure According to the temperature and pressure of the atmosphere, the calculation formula of the air inlet temperature T1 is: Where, k = 1.402 is the adiabatic index of air; Ma inlet is the Mach number; The formula for calculating the pressure at the air inlet is: (3) Air inlet flow rate The flow rate through the inlet is calculated using the Mach number and the standard gas constant using the formula: Where R air =287 is the gas constant of air; ② Compressor model (1) Calculation of speed coefficient Speed ​​coefficient n LC,HPC The calculation formula is as follows: Where N H is the input speed, rpm; T2 is the input temperature, K; N cons is the speed at the design point; (2) Obtaining compressor characteristic data The pressure ratio and efficiency corresponding to different speeds and flow rates are characteristic data derived from the Gasturb tool; (3) Interpolation calculation Compressor performance is estimated using a two-dimensional interpolation method. The independent variables for interpolation include pressure ratio and speed coefficient. Given characteristic data, including pressure ratio, flow rate, and efficiency at different speeds, these data are interpolated to determine performance indicators under actual operating conditions. During the interpolation process, the corresponding performance data point is first found based on the given pressure ratio and speed coefficient. Then, the flow rate and efficiency are calculated using bilinear interpolation, and the interpolation coefficient ckzx is calculated. This coefficient is used to represent the position of the current pressure ratio between two adjacent data points: Where x j is the current pressure ratio; x a and x b are adjacent data points respectively; Use the linear interpolation formula to calculate the corresponding flow rate and efficiency: G LC =G LCa +ckzx·(G LCb -G LCa ) or HPC =the HPCa +ckzx·(η HPCb -or HPCa ) Where G LCa and G LCb is the flow rate at the interpolation point; η HPCa and η HPCb is the efficiency of the interpolation point; (4) Compressor output flow Compressor output flow G HPC Calculated by the following formula: Where P2 is the input pressure; T2 is the input temperature; 101325 is the standard atmospheric pressure; (5) Isentropic compressor outlet temperature Based on the assumption of an isentropic process, the compressor outlet temperature T3 is calculated by the following steps: Calculate the input enthalpy and the outlet enthalpy for an isentropic process: H2=f t2hc (T2,0) H c,id =f t2hc (T c,id ,0) According to the efficiency of the compressor, the ideal outlet enthalpy is corrected to obtain the actual outlet enthalpy H3: Calculate the actual outlet temperature T3: <h2 style=";text-align:left;direction:ltr">T3=f<h2 style=";text-align:left;direction:ltr"> h2tc <h2 style=";text-align:left;direction:ltr"> (H3,0) Where, f h2tc 、f t2hc They are functions for calculating the corresponding temperature based on the enthalpy value and the corresponding enthalpy value based on the temperature; (6) Calculation of power and flow output Considering the cooling flow and bleed air flow, the cooling flow comes from cooling the stator blades and rotor blades; the bleed air flow is used for sealing. The power and flow output of the compressor are calculated using the following formula: W HPC =(H3-H2)·G HPC ·(1-0.05-0.01) G3=G HPC ·(1-0.05-0.05-0.01) ③Combustion chamber model (1) Mass flow and pressure loss in the combustion chamber: Total output flow of the combustion chamber G out is the sum of the air flow and fuel flow entering the combustion chamber: G out =G in +G fuel The outlet pressure of the combustion chamber P out 97% of inlet pressure: P out =P in ·0.97 (2) Calculation of enthalpy of combustion chamber inlet gas Enthalpy of the gas at the combustion chamber inlet H in Through its temperature T in Calculated using the function between enthalpy and temperature: H in =f t2hc (T in ,0) Where 0 means air is an ideal gas; (3) Calculation of air-fuel ratio The fuel-air ratio is defined as the ratio of fuel flow to air flow: (4) Calculation of enthalpy of mixed gas According to the enthalpy of the intake gas, the lower calorific value of the fuel and the combustion efficiency, the enthalpy of the mixed gas at the combustion chamber outlet is H out The calculation is as follows: Where H fuel =8.6×10 5 J / kg is the calorific value of fuel, which represents the total energy per unit mass of fuel; 4.3124×10 7 It is the lower calorific value of the fuel, which indicates the amount of heat released when a unit mass of fuel is completely burned; (5) Calculation of combustion chamber outlet temperature According to the total enthalpy H of the mixed gas out And the fuel air-fuel ratio FAR4, calculate the temperature T at the combustion chamber outlet out : T out =f h2tc (H out ,FAR4) ④ Hybrid cavity model (1) Calculation of enthalpy of input gas The modeling of the mixing chamber first requires the calculation of the enthalpy of the input gas; the fuel gas H4 is calculated from its temperature: H4=f t2hc (T4,FAR4) Where T4 is the temperature of the fuel gas; FAR4 is the fuel air-fuel ratio; (2) Calculation of output gas flow rate and mixing enthalpy Total output flow of the mixing chamber G out The sum of the air flow rate and the fuel gas flow rate G4, plus the cooling air flow rate G cooling : G out =G4+G cooling Air flow G air It is the total flow rate of the mixed gas minus the flow rate of the fuel gas: G air =G out -G fuel Output fuel air-fuel ratio FAR out Indicates the amount of fuel contained in a unit mass of air, which is the ratio of the fuel gas flow rate to the air flow rate: Enthalpy of the gas mixture H mix is the weighted average of air and fuel gas taking into account their respective mass flow rates: (3) Calculation of mixed gas temperature Temperature of the mixed gas T mix Through its total enthalpy H mix and fuel air-fuel ratio FAR out Calculation yields: T mix =f h2tc (H mix ,FAR out ) ⑤Turbine model (1) Calculation of turbine expansion ratio and speed coefficient First, the expansion ratio Pr is defined according to the expansion process of the turbine HPT , and its calculation formula is: Where, P 41 is the turbine inlet pressure; P 45 is the turbine outlet pressure; Define the reduced speed coefficient n LC，HPT , considering the relationship between intake temperature and turbine speed, the calculation formula of the equivalent speed coefficient is as follows: Where N H is the actual turbine speed; T 41 is the turbine inlet temperature; δ is the reduced temperature coefficient; N cons represents the design point speed; (2) Turbine characteristic data and interpolation Turbine characteristic data includes turbine performance at different speeds and pressure ratios. Through interpolation calculations, flow and efficiency under specific operating conditions are estimated. Calculate the flow G based on the interpolation coefficient ckzx LC and efficiency η HPT : G LC =G LCa +ckzx·(G LCb -G LCa ) or HPT =the HPTa +ckzx·(η HPTb -or HPTa ) Where G LCa and G LCb is the flow rate at the interpolation point; η HPTa and η HPTb is the efficiency of the interpolation point; The interpolation result is flow G LC and efficiency η HPT , according to the compressor outlet pressure P 41 and intake air temperature T 41 , calculate the actual turbine flow G HPT : Where, P 41 is the turbine inlet pressure; T 41 is the turbine inlet temperature; (3) Power output and efficiency calculation Turbine output enthalpy H T,out It is the enthalpy value corrected according to the actual efficiency; first calculate the turbine inlet enthalpy H T,0 and ideal outlet enthalpy H T,id , and then correct the enthalpy value H according to the efficiency of the turbine T,out : H T,0 =f t2hc (T 41 ,FAR out ) H T,id =f t2hc (T T,id ,FAR out ) H T,out =(H T,id -H T,0 )·η HPT +H T,0 Where η HPT is the efficiency of the turbine; Through the actual outlet enthalpy H T,out , calculate the turbine output temperature T 45 : T 45 =f h2tc (H T,out ,FAR out ) Turbine power output W HPT It is determined by the difference between the outlet and inlet enthalpies of the turbine and the flow rate, and is calculated using the following formula: W HPT =-(H T,out -H T,0 )·G HPT ⑥ Solution of PSFC of turboshaft engine The specific fuel consumption (PSFC) of a turboshaft engine is based on the flow and energy balance of each engine component. The output power of the turboshaft engine is calculated using the known output torque and engine speed. The formula is: Where, T out is the output torque, in Nm; N H is the engine speed; To calculate the specific fuel consumption PSFC, the formula is as follows: Where G fuel is the fuel flow rate; P out is the output power of the engine; The calculation of error balance involves the following two errors: (1) Power error The power error reflects the deviation between the turbine output power and the theoretically calculated power. The calculation formula is: Where W HPT is the output power of the high-pressure turbine; W HPC is the outlet power of the high-pressure compressor; η rotor is the turbine rotor efficiency; P out is the calculated engine output power; this error measures the energy conversion efficiency between the turbine and compressor; (2) Flow error The flow error reflects the balance of the engine airflow, and its calculation formula is: Where G HPC is the flow rate at the high-pressure turbine outlet; 0.1G HPC is the cooling air flow rate; G fuel is the fuel flow rate; G3 is the flow rate at the high-pressure compressor outlet; this error is used to ensure the balance between the cooling airflow, fuel flow rate, and compressor outlet flow rate; through the above calculation, the final error balance consists of two parts: These two errors reflect the balance between power and flow; S22: Establish a power battery internal resistance model to clarify the battery capacity boundary; Establish an equivalent internal resistance model, the output power P of the power battery b It is expressed as follows: Where V oc is the open circuit voltage of the power battery; R int is the internal resistance of the battery; I b is the charge and discharge current of the battery; I b It is expressed as follows: The SOC of the battery is determined by the integration method, which is expressed as follows: Where SOC0 represents the initial SOC value of the battery; Q b Represents the battery capacity, and differentiating the equation yields: When P b When it is positive, the power battery is discharging, otherwise, it is charging; Parameter identification is performed through static capacity test, open circuit voltage test and HPPC test to obtain the peak current corresponding to different SOC and discharge rate; The peak current is expressed as follows: The peak current surface was obtained by fitting using a fourth-order polynomial fitting method.

9. The energy management method for a hybrid vertical take-off and landing aircraft based on a model predictive control algorithm according to claim 1, characterized in that: The step S3 specifically includes the following steps: S31: Establish a test condition. The test condition contains three identical subtasks. Each subtask is divided into hovering, cruise climb, cruise, and cruise descent. The second and third subtasks are started after the previous stage is completed and the aircraft descends. After completing all subtasks, the aircraft begins to descend vertically and lands. S32: Design a rule-based energy management strategy for the entire machine. The specific control rules of the rule-based energy management strategy are as follows: (1)CD mode When SOC>SOC CS When the system is in CD mode, the system relies on battery drive, and the change of SOC directly affects the discharge state of the battery. The system determines whether to continue to maintain battery drive mode or switch to CS mode by monitoring SOC. (2) CS mode When SOC min <SOC≤SOC CS When the system starts the engine drive mode, the engine provides power to the aircraft and the battery is only used for auxiliary; (4) Charge mode When SOC≤SOC min When the system switches to charge mode, the power provided by the engine is used to charge the battery to restore the battery's SOC; S33: Design a whole-machine energy management strategy based on the DP algorithm and simulate the APU maximum probability output power relationship diagram; first, the state variable is determined to be SOC. The SOC value will change in each stage and is affected by the battery charging and discharging power in the current stage. The decision variable in each stage is the APU power generation power P APU , by controlling the power P APU The power demand of the aircraft is matched to meet the power demand of the aircraft. At the same time, the SOC value of the battery will also change. The selection of power generation power needs to balance the change of fuel consumption and SOC, that is: x k =SOC k u k =P APU (k) The DP cost function is composed of multiple factors, including fuel consumption and SOC deviation. Fuel consumption is related to the engine's power generation and can be calculated using the following formula: Where, f(P APU (i)) represents the relationship between APU power generation and fuel consumption; the cost function of the DP algorithm is expressed as follows: Where N is the mission step length, and the terminal penalty coefficient α is 10 6 , used to ensure the convergence of the SOC final value; the overall limitations of the system are as follows: Where, P dmd Indicates power demand; (P APU , P b ) represent the output power of the APU and power battery system respectively; the SOC is discretized to obtain the state variable matrix; the DP algorithm first traverses each state variable at each time node in reverse order to obtain the optimal control table; then continues to solve forward to derive the target control sequence. The inverse process is expressed as follows: J(k,i)=min{m F (k,i,j)+J(k+1,j)} k=0,1,…,N; i,j=0,1,…,M Where J(k+1,j) represents the slave node SOC k+1,j To the final node SOC N Minimum fuel consumption; m F (k,i,j) represents the battery from SOC k,i To SOC k+1,j The instantaneous fuel consumption when changing; M is the number of discrete grids of SOC; When the initial SOC and target SOC are determined, the global optimization result of the DP strategy is fixed. The optimal power allocation at each required power in each starting state and expected ending state is statistically calculated from the optimization results. The starting SOC, target ending and required power are used as inputs, and the APU output power with the maximum probability of occurring under the required power is used as the target power of the required power. This power is used as the target power at each moment in the prediction time domain to calculate the SOC reference trajectory in the prediction time domain.

10. The energy management method for a hybrid vertical take-off and landing aircraft based on a model predictive control algorithm according to claim 1, characterized in that: The step S4 of implementing the design and verification of the MPC strategy specifically includes the following steps: S41: Using the APU maximum probability output power relationship diagram obtained by DP, design an energy management strategy based on the MPC algorithm, including: According to the MPC principle, the battery SOC is selected as the state variable and the APU output power is selected as the control variable to design a hybrid power system model predictive control strategy for optimizing fuel economy: (1) Establish a prediction model: Assume that the power demand remains unchanged within the prediction time domain k to k+p. That is, the power demand P(k) at time k is used as the power demand at each time in the prediction time domain [P(k+1), P(k+2), …, P(k+p)]. On this basis, the upper and lower limits of the SOC in the prediction time domain are estimated to determine the feasible domain of the SOC. (2) Solving the optimization problem: In the prediction time domain k to k+p, the hybrid system power allocation optimization objective function is established. Under system constraints, the APU maximum probability output power diagram is used to solve the optimal control in the SOC feasible domain and obtain the optimal control sequence [u(k),u(k+1|k),…,u(k+p|k)]; (3) Optimal control implementation: In the MPC strategy, only the first element u(k) of the control sequence is applied to the current moment and solved, and then the system state is updated. The hybrid system is described as: in: In this process, the optimization is carried out in the interval k~k+p, and the objective function in this interval is established as shown in the following formula: p is the length of the prediction time domain, m f is the engine fuel consumption, SOC(t) and SOC ref is the SOC at time t and the reference value; ω is the penalty coefficient; At the same time, the following constraints are met: In the formula, max and min are the upper and lower limits of the corresponding quantity; The MPC solution process is as follows: (1) Application of DP results in MPC. The solution process is expressed as: ① Divide the prediction time domain k~k+p into p+1 stages, i.e. [k, k+1,…, k+p]; ② Discretize the feasible domain of SOC prediction within k~k+p; ③ Use the obtained APU maximum probability output power map to optimize the objective function and obtain [P APU (k),P APU (k+1|k),P APU (k+2|k),…,P APU (k+p|k)], and apply the first element to the controlled system to update the state; (2) SOC range determination and discretization When optimizing the solution, the SOC range is restricted and the SOC range in the prediction time domain is calculated based on the current SOC and required power, so as to further narrow the SOC range. At any time, the following power constraints are met: P(k)=P APU (k)+P b (k) The maximum output power of the battery is: P b,max (k)=min(P(k),P b,dismax (SOC(k))) The minimum output power of the battery is: P b,min (k)=max(P(k)-P APU,max ,P b,chgmax (SOC(k))) Taking SOC(k) as the starting point and P b,max (k) is the SOC lower boundary power, P b,min (k) is the upper boundary power of SOC, and the upper and lower boundaries of SOC in the entire prediction time domain are calculated; Where, P b,dismax (SOC(k)) represents the maximum discharge power of the battery at SOC(k); P b,chgmax (SOC(k)) represents the maximum charging power of the battery under SOC(k); P b,max (k) and P b,min (k) is the maximum and minimum power of the battery under the required power P(k); (3) Determination of the control variable range: By limiting the range of the control variable, the MPC calculation amount is reduced. Under the condition of the required power P(k), the corresponding P APU Scope of work; The maximum output power of the APU is: P APU,max (k)=min(P(k),P APU,max ) The minimum output power is: P APU,min (k)=max(P(k)-P b,dismax (SOC(k)),0) Where, P APU,max (k) and P APU,min (k) is the maximum and minimum power under the required power P(k); on this basis, the discrete sequence of the control variable is obtained according to a certain discrete interval; S42: Introducing CS and charge modes to improve the overall energy management strategy based on the MPC algorithm, including: After the SOC hits bottom, add strategies similar to CS and Charge to ensure that the SOC is close to the target end SOC without frequent engine starts and stops. The details are as follows: (1) When SOC < SOC target When -0.02, the APU outputs at maximum data power; (2) When SOC is reduced to SOC target After -0.02, it returned to SOC target +0.02 to enable MPC optimization again; (3) When SOC target -0.01<SOC<SOC target When the power requirement is +0.01, the APU output power is determined by the power requirement. That is, when the power requirement P is less than 60kW, the APU does not work and the battery provides all the energy. When the power requirement P is greater than 60kW, the APU outputs power at its maximum. (4) Not within the control range of the above SOC, that is, SOC target -0.02<SOC<SOC target -0.01 and SOC target +0.01<SOC<SOC target When it reaches +0.02, the control boundary of SOC is involved. Under this condition, the APU is allowed to recharge at maximum power output to avoid frequent starting and stopping.

Images