Electric hovercar battery charging optimization method considering operation safety and economical efficiency

Abstract

The invention relates to an electric hovercar battery charging optimization method considering operation safety and economical efficiency, and belongs to the technical field of batteries. The method comprises the following steps: establishing a distributed electrothermal coupling model and a semi-empirical performance attenuation model driven by physical-data fusion of a power battery, and coupling the model to construct a battery fast digital twinning model oriented to full life cycle simulation; considering battery attenuation equivalent cost and charging electricity price cost, and establishing a battery charging economy model; based on a battery fast digital twinborn model and a charging economy model, multiple optimization sub-targets are considered to construct a multi-target charging optimization problem, and an advanced optimization algorithm is utilized to realize optimal solution; and analyzing influence factors of the optimal charging scheme based on the Pareto frontier theory, and determining a multi-stage constant-current optimal charging scheme for the urban air traffic operation scene. The invention can ensure that the battery system of the electric flying car can supplement energy safely and efficiently and improve the operation economy of the whole car.

Description

Technical Field

[0001] The invention belongs to the technical field of batteries and relates to a battery charging optimization method for an electric flying car taking into account operational safety and economy. Background Art

[0002] The battery specific energy of the current mainstream electrochemical system (graphite anode + liquid electrolyte + transition metal oxide cathode) has not yet achieved an effective breakthrough. The development of battery fast charging is an effective measure to alleviate range anxiety and improve user satisfaction. For electric flying cars, high-safety fast charging technology can effectively improve the company's operating income and energy efficiency. Therefore, it is particularly important to study battery fast charging control in urban air traffic scenarios. In order to achieve high-safety fast charging of batteries, various charging schemes have been reported in the literature, among which typical charging schemes include constant current-constant voltage CC-CV charging, multi-stage constant current-constant voltage MCC-CV charging, constant power-constant voltage CP-CV charging, pulse charging PC, boost charging and variable current charging VCC.

[0003] Charging optimization is essentially to optimally plan the battery charging current based on various operating constraints to avoid battery overcharging, overheating, lithium deposition and even thermal runaway. At present, many studies on lithium-ion battery fast charging optimization for electric vehicles have been reported in the literature. For example, the battery fast charging strategy based on the model predictive control MPC algorithm considers SOC and temperature as control constraints, and optimizes the charging time and maximum temperature rise at the same time. Based on linear time-varying MPC and control-oriented electrochemical-thermal coupling model, a fast charging strategy that takes into account battery health is developed to effectively constrain the electrochemical reaction state or parameters inside the battery, such as inhibiting battery lithium deposition, to avoid battery safety risks and extend battery life. Or establish a control-oriented electrochemical-thermal-aging multi-field coupling model to achieve the best compromise between battery charging time and aging through multi-objective optimization. Economic charging usually considers the optimal energy management of the power grid but ignores the cost of use for consumers. The power battery of electric flying cars needs to meet the charging problems of safety, efficiency and economy at the same time. Summary of the invention

[0004] In view of this, the purpose of the present invention is to provide an electric flying car battery charging optimization method taking into account operational safety and economy, which is based on the battery's electric-thermal-aging multi-field coupling model, and proposes an optimal charging strategy that considers charging cost (including electricity cost and battery aging cost), charging time and average temperature rise, wherein the adopted MCC charging protocol is optimized by the NSGA-II algorithm.

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

[0006] A method for optimizing battery charging of an electric flying car taking into account operational safety and economy, the method comprising the following steps:

[0007] S1: Establish a distributed electric-thermal coupling model and a semi-empirical performance degradation model driven by physical-data fusion for power batteries, and couple the two to build a fast digital twin model of batteries for full life cycle simulation;

[0008] S2: Considering the equivalent cost of battery degradation and the cost of charging electricity, a battery charging economic model is established;

[0009] S3: Based on the battery fast digital twin model and charging economy model, multiple optimization sub-objectives are considered to construct a multi-objective charging optimization problem and use advanced optimization algorithms to achieve the optimal solution;

[0010] S4: Analyze the influencing factors of the optimal charging scheme based on the Pareto frontier theory, and determine the multi-stage constant current optimal charging scheme for urban air traffic operation scenarios.

[0011] Further, in step S1, the distributed electrothermal coupling model driven by physical-data fusion includes a distributed thermal model driven by physical-data fusion and a first-order RC equivalent circuit model, the distributed thermal model driven by physical-data fusion is composed of a physics-based two-dimensional distributed thermal model and a temperature calibration model based on a gated recurrent neural network, and the physics-based two-dimensional distributed thermal model is composed of a battery body spatial distribution thermal model and a positive and negative pole ear concentrated mass thermal model.

[0012] In the battery fast digital twin model for full life cycle simulation, which is constructed by coupling the distributed electro-thermal coupling model driven by physical-data fusion and the semi-empirical performance degradation model, the input vector u of the battery digital twin model is dt , state vector x dt and the output vector y dt They are respectively expressed as follows:

[0013]

[0014] Among them, the subscript dt represents the digital twin DT, k represents the time step; the input vector contains the current and ambient temperature at time k, and the state vector contains the state of charge SOC and polarization voltage U at time k. p , the time coefficient of temperature χ and the state of health SOH, the output vector contains the terminal voltage U at time k t , each key temperature T and capacity loss rate Q loss .

[0015] Further, in step S1, the process of establishing a distributed electrothermal coupling model driven by physical-data fusion is as follows:

[0016] The spatial distribution thermal model of the battery body based on the physical two-dimensional distributed thermal model obeys the transient heat conduction equation that satisfies the Robin boundary condition. The expression is as follows:

[0017] ρC h D t Tk x D xx Tk y D yy T=q

[0018] Subject to the following boundary conditions:

[0019]

[0020] Among them, T(x,y,t) and q(x,y,t) represent the temporal and spatial distribution temperature and volumetric heat generation rate of the battery, ρ, C h and k are the battery volume average density, specific heat capacity and anisotropic thermal conductivity; and x∈[0,w], y∈[0,l], where w and l represent the width and length of the battery; subscripts r, l, t, and b refer to the right, left, top, and bottom boundaries of the battery, and h x ∈[h r ,-h l ] and h y ∈[h t ,-h b ] represents the equivalent heat dissipation coefficient of the battery in the x and y directions; the values ​​of the heat dissipation coefficient are all positive, and “±” is used to distinguish the heat dissipation direction in the Cartesian coordinate system; T ∞,x =[T r,∞ ,T l,∞ ] and T ∞,y =[T t,∞ ,T b,∞ ] represents the external cooling environment temperature of the battery in both directions, assuming that the cooling environment around the battery in the incubator is consistent.

[0021] The Chebyshev-Galerkin approximation method and time-space separation technology are introduced to reduce the model order. After coordinate change, boundary condition homogenization, auxiliary temperature function solution, and arbitrary temperature function solution, the PDE type control equation is transformed into multiple ODE type control equations for solution. The deformed state space expression is as follows:

[0022]

[0023] Among them, E, A, B and C are system matrices, T e is any time-invariant temperature function related to the spatial position, the system state x = (χ 00 ,x 10 ,…,χN0 ,x 01 ,x 11 ,…,χ N1 ,x 0N ,x 1N ,…,χ NN ) T ; System input u = [q 1] T , q is the heat generation rate per unit volume; the system output y is defined according to the actual control system requirements; extract T2~T4 and the volume average temperature T b As the output temperature of the thermal model, that is:

[0024]

[0025] The mathematical expressions of the positive and negative lug concentrated mass thermal models based on the physical two-dimensional distributed thermal model are as follows:

[0026]

[0027] Among them, T t and T ∞ Refers to the lug temperature and ambient temperature respectively; m t , C pt 、h t and A t are the mass, specific heat capacity, convection heat transfer coefficient and convection heat transfer area of ​​the lug respectively; q t and q ct They represent the heat generation rate of the tab and the heat flux between the tab and the body respectively; for q ct , the flow into the pole ear is negative, and the flow out of the pole ear is positive; the heat generation of the pole ear q t Ohmic heat only, calculated by the following formula:

[0028] q t =I 2 R st =I 2 (R ot +R ct )

[0029] Among them, R ot and R ct They are the ohmic internal resistance of the tab and the contact internal resistance between the tab and the wire.

[0030] An empirical formula is introduced to calculate the heat exchange between the battery tab and the body:

[0031]

[0032] Among them, T cb is the temperature of the discrete battery unit in direct contact with the tab, which is approximately equal to the tab temperature; T mis the temperature of the measuring point on the same side as the lug, h cb , k y and Δx are the contact coefficient between the battery tab and the body that hinders heat transfer, the thermal conductivity of the battery body in the y direction, and the temperature is determined by T cb Change to T m The length traveled, A ct It is the cross-sectional area for actual heat conduction between the tab and the body.

[0033] Heat flow q between the tab and the body ct The distribution on the battery body is as follows:

[0034]

[0035] Where i,j=1,…,N g -1, k = 1, ..., M, M is the number of discrete units of the battery body, take M = (N g -1)×(N g -1), N g is the number of discrete points in the x and y directions of the battery body; the number of discrete units in the x and y directions of the battery body is equal, and this value can be determined according to the computing and storage capabilities of the controller in the actual application scenario; q ct,k is the heat flow q ct Assigned to the battery body discrete unit (x i ,y j ) of heat, β ct is the adjustment factor, indicating q ct,k The heat that actually contributes to the temperature change of the unit is in q ct,k The proportion in the range is [0,1], x t and t It is the center coordinate of the connection between the battery tab and the body.

[0036] The forward propagation mechanism of the gated recurrent neural network model is represented by the following set of equations:

[0037]

[0038] Among them, σ and tanh represent gated activation functions, σ refers to the sigmoid function, expressed as σ = 1 / (1 + exp(-x)), tanh refers to the hyperbolic tangent function, expressed as tanh = (exp(x) - exp(-x)) / (exp(x) + exp(-x)), represents a linear activation function, ⊙ and [x,x] represent element-wise multiplication and matrix concatenation, respectively, and W r , W z , W h and W0 are the weight matrices of the reset gate, update gate, short-term output and output layer respectively; br 、b z 、b h and b0 is the corresponding bias matrix, r t and z t Refers to the outputs of the reset gate and update gate respectively, and h t Denote candidate hidden states and hidden states respectively, y t Represents the output of GRU-NN.

[0039] During the training phase, the average mean square error (MSE) of all outputs is selected as the loss function of GRU-NN:

[0040]

[0041] Among them, M and N represent the predicted output length and sequence data length of GRU-NN at a certain time step, respectively, and y i,j and are the temperature T at time step j respectively. i The estimated value and reference value of GRU-NN use the loss function gradient to continuously update the weights and biases in the network through back propagation.

[0042] The first-order RC equivalent circuit model is selected as its electrical sub-model. The time domain discrete form of the model is expressed as:

[0043]

[0044] U t,k =U oc,k -U p,k -R s,k I k

[0045] Among them, U t , U oc , R s and I represent the terminal voltage, open circuit voltage, ohmic internal resistance and current of the battery respectively. p and C p Represent polarization resistance and polarization capacitance respectively, C b is the current available capacity of the battery, Δt is the sampling interval, and the subscript k represents the time k.

[0046] Further, in step S1, the semi-empirical performance degradation model is expressed as:

[0047]

[0048] Among them, Q loss Indicates the battery capacity loss rate in %. ezp , β rateand z are respectively the exponential prefactor, the battery aging acceleration factor related to the current rate, and the power function factor related to the ampere-hour throughput; E a is the battery aging activation energy, in J·mol -1 ;I rate is the battery loading current, in C; R g is the ideal gas constant; T b is the average temperature of the battery volume, in °C, calculated by a two-dimensional distributed thermal model; Ah is the ampere-hour throughput of the battery during charge and discharge. As the number of cycles increases, Ah=2×num_cycle×DOD×C bat , where num_cycle represents the number of charge and discharge cycles of the battery.

[0049] Further, in step S2, the battery attenuation equivalent cost model considers the use cost caused by the performance degradation of a single charge and discharge cycle, which is related to the price of a new battery, the resale price of the battery when retired, and the number of equivalent cycles of the battery in its entire life cycle. The residual value of the battery when retired is calculated as follows:

[0050]

[0051] Among them, RV fresh and RV EOL Represent the sales price or residual value of new batteries and retired batteries, respectively, life Refers to the service life of the battery for electric flying cars; dep represents the annual depreciation rate of the battery; dr represents the discount rate; based on this, combined with the battery performance attenuation model, the battery ampere-hour throughput Ah in the entire life cycle tot It is expressed as:

[0052]

[0053] Among them, Q loss Indicates the mandatory retirement threshold of the power battery capacity attenuation of electric flying cars;

[0054] The equivalent number of cycles of the battery in its entire life cycle is N cycle It is expressed as:

[0055]

[0056] The cost of battery performance degradation during a single charge and discharge cycle is expressed as follows:

[0057]

[0058] The amount of electricity that the charging pile needs to provide during the entire charging process is:

[0059]

[0060] P mcci =|U t (τ)I(τ)|+I(τ)(U oc (τ)-U t (τ)

[0061] Among them, E mcc Provides the power required for the entire charging process, in kWh; N is the number of constant current charging segments in the MCC charging protocol, P mcci and t mcci are the power and time required for each constant current charging stage respectively.

[0062] The battery charging electricity cost is expressed as:

[0063] C elec =p e E mcc

[0064] In the formula, p e The time-of-use electricity price is per kilowatt-hour.

[0065] Further, in step S3, the multi-objective charging optimization problem of the battery considers a multi-stage constant current charging mode, and the multi-optimization charging sub-objectives include charging time, charging cost and maximum temperature rise, which are expressed as follows:

[0066] ①Charging time t mcc It is expressed as follows:

[0067]

[0068] Where N is the number of constant current charging stages, t mcc The time taken for the entire charging process, in minutes; t N The duration of the entire charging process in seconds.

[0069] ②Charging cost C mcc It is expressed as:

[0070] C mcc =C age +C elec

[0071] ③The maximum temperature rise is characterized by the highest temperature in the battery plane, expressed as:

[0072] T mcc,max =max([T1 T2 T3 T4 T5])=max(T1)

[0073] Among them, T i Indicates the entire charging process T i Temperature series at the measuring points.

[0074] Each sub-objective function is normalized, and the normalized expression of each sub-objective function is as follows:

[0075]

[0076] Among them, J mcc,t , J mcc,c and J mcc,tr They represent the battery charging time sub-function, charging cost sub-function and charging maximum temperature rise sub-function respectively. There is a mutual exclusive relationship between the three and they cannot reach the optimum at the same time; t min and t max Respectively represent the minimum and maximum charging time of the battery; T min and T max are the minimum and maximum charging temperatures of the battery respectively.

[0077] By introducing the weight factors of each sub-objective function, we can get the complete multi-objective charging function:

[0078]

[0079] α+β+γ=1

[0080] Charging constraints include time and electric and thermal operation constraints, which are described as follows:

[0081]

[0082] Among them, I min and I max Respectively represent the minimum and maximum charging current allowed by the battery; U min and U max Respectively represent the minimum and maximum values ​​of the battery operating voltage; SOC min and SOC max Respectively represent the upper and lower limits of the battery SOC operating range; SOH min and SOH max They represent the upper and lower limits of the battery SOH respectively.

[0083] Further, in step S3, the optimal solution process of the multi-objective charging optimization problem using an advanced optimization algorithm includes:

[0084] The sub-population is divided by the median of the ascending fitness value, and dynamic adjustment is performed in each iteration; the particles with fitness values ​​less than the median are set as "top particles", and the particles with fitness values ​​greater than the median are set as "bottom particles", and equal numbers of particles are extracted from them to form a local PSO model, which are defined as the "dominant group", the "mixed group" and the "inferior group" of the local PSO model, respectively. Among them, the "dominant group" contains all top particles, the "inferior group" contains all bottom particles, and the "mixed group" of the local PSO model includes both top particles and bottom particles.

[0085] For the update of the dominant group, the random step length theory of Levy flight is introduced to continuously adjust the iteration step length. The updating principle of Levy flight step length is as follows:

[0086]

[0087] in, and x g (k) represents the global optimal particle position and the global optimal particle position of the population after Levy flight update at the kth iteration; α s is the step size adjustment factor, and Levy(λ) is the random search path.

[0088] The flight trajectory is simulated using the Mantegna algorithm, and the mathematical description is as follows:

[0089]

[0090] Among them, Levy~μ=k -λ , 1<λ≤3, and μ and υ obey the following normal distribution:

[0091]

[0092] Among them, the standard deviation σ μ and σ υ Determined by the following formulas:

[0093]

[0094] Where Γ(·) is the gamma function.

[0095] For the update of the disadvantaged group, the optimal solution is searched in the entire space by mutation, and whether the particle mutates is determined according to the following formula:

[0096]

[0097] GV(x)=x×(1+N(0,1))

[0098] Among them, r and N(0,1) are random numbers between [0,1], x is the original parameter value, GV(x) is the mutated value; kmax is the maximum number of iterations; when the above inequality holds, the Gaussian mutation operation is performed on the particles; after multiple iterations, the probability of judging that the above inequality holds gradually decreases, and the probability of mutation gradually decreases;

[0099] The concept of hybrid particles is introduced to improve the traditional velocity update formula, and the hybrid particles are recorded as p mix (k), mixed particle p mix The dimension of (k) is a random mixture of the historical optimal values ​​of each particle. Therefore, the position and velocity update method of the disadvantaged group particles is expressed as:

[0100] v id (k+1)=wv id (k)+c1r1(p id (k)-x id (k))+c2r2(p gd (k)-x id (k))+c3r3(p mix (k)-x id (k))

[0101] x id (k+1)=x id (k)+v id (k+1)

[0102] Where c3 is the hybrid learning factor and r3 is a random number between (0,1).

[0103] For the update of the mixed group, cosine and sine functions are introduced into the learning factor to make the self-learning factor monotonically decreasing and the population learning factor monotonically increasing. The speed update formula of the mixed group is expressed as:

[0104]

[0105] For the setting of nonlinear inertia weight during update, refer to the following formula:

[0106]

[0107] The beneficial effects of the present invention are:

[0108] (1) Based on the battery distributed electrothermal coupling model applicable to the short time domain and the semi-empirical performance attenuation model applicable to the long time domain, the present invention constructs a battery fast digital twin model for the entire life cycle, which can perform efficient and accurate full life cycle performance simulation of high-power and large-size electric flying cars, and accelerate the iterative optimization and on-site deployment of the charging optimization algorithm.

[0109] (2) The present invention constructs a multi-objective charging optimization problem, taking into account multiple optimization sub-objectives such as battery charging time, charging cost and maximum temperature rise, and uses an improved dynamic multi-population particle swarm optimization algorithm to achieve the optimal solution, which makes up for the defect that existing studies rarely consider the operating costs and benefits of vehicles in the long-term use process. In addition, based on the Pareto frontier theory, the influencing factors affecting the optimal charging of batteries in urban air traffic operation scenarios are deeply analyzed to ensure the safe, economical and efficient energy replenishment of electric flying cars.

[0110] Other advantages, objectives and features of the present invention will be described in the following description to some extent, and to some extent, will be obvious to those skilled in the art based on the following examination and study, or can be taught from the practice of the present invention. The objectives and other advantages of the present invention can be realized and obtained through the following description. BRIEF DESCRIPTION OF THE DRAWINGS

[0111] 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 in conjunction with the accompanying drawings, wherein:

[0112] Figure 1 This is a flow chart of the battery charging optimization method for an electric flying car of the present invention;

[0113] Figure 2 A schematic diagram of a battery distributed electrothermal model driven by physical-data fusion of the present invention;

[0114] Figure 3 A schematic diagram of a battery distributed thermal model driven by physical-data fusion of the present invention;

[0115] Figure 4 It is a schematic diagram of the arrangement of the battery thermocouple of the present invention;

[0116] Figure 5 It is a schematic diagram of a first-order RC equivalent circuit model of the present invention;

[0117] Figure 6 A schematic diagram of the coupling mechanism of the battery electric-thermal-attenuation model of the present invention;

[0118] Figure 7 This is a multi-objective charging optimization framework diagram of the electric flying car battery of the present invention;

[0119] Figure 8 Flowchart for implementing the advanced optimization algorithm of the present invention. DETAILED DESCRIPTION

[0120] The following describes the embodiments of the present invention by 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 only illustrate the basic concept of the present invention in a schematic manner, and the following embodiments and features in the embodiments can be combined with each other without conflict.

[0121] Among them, the drawings are only used for illustrative explanations, and they only represent schematic diagrams rather than actual pictures, and should not be understood as limitations on the present invention. In order to better illustrate the embodiments of the present invention, some parts of the drawings may be omitted, enlarged or reduced, and do not represent the size of actual products. For those skilled in the art, it is understandable that some well-known structures and their descriptions in the drawings may be omitted.

[0122] The same or similar numbers in the drawings of the embodiments of the present invention correspond to the same or similar parts; in the description of the present invention, it should be understood that if the terms "upper", "lower", "left", "right", "front", "rear", etc. indicate the orientation or position relationship, they are based on the orientation or position relationship shown in the drawings, which is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operate in a specific orientation. Therefore, the terms describing the position relationship in the drawings are only used for illustrative purposes and cannot be understood as limiting the present invention. For ordinary technicians in this field, the specific meanings of the above terms can be understood according to specific circumstances.

[0123] See also Figures 1 to 8 , which is an optimization method for charging batteries of electric flying cars taking into account operational safety and economy.

[0124] Example

[0125] like Figure 1 The overall flow chart of the battery charging optimization method of an electric flying car taking into account operational safety and economy of the present invention is shown, which specifically includes the following steps:

[0126] S1: Establish a distributed electro-thermal coupling model and a semi-empirical performance degradation model driven by physical-data fusion for high-power and large-size lithium-ion power batteries, and couple them to build a fast digital twin model of batteries for full life cycle simulation;

[0127] S2: Considering the equivalent cost of battery degradation and the cost of charging electricity, a battery charging economic model is established;

[0128] S3: Based on the battery fast digital twin model and charging economy model, considering multiple optimization sub-goals such as battery charging time, charging cost and maximum temperature rise, a multi-objective charging optimization problem is constructed and the optimal solution is achieved using advanced optimization algorithms;

[0129] S4: Based on the Pareto frontier theory, the influencing factors of the optimal charging solution are analyzed in detail, and the multi-stage constant current optimal charging solution for urban air traffic operation scenarios is clarified.

[0130] In step S1 of this embodiment, Figure 2 As shown in the figure, the established distributed electrothermal coupling model of high-power lithium-ion power battery driven by physical-data fusion includes a distributed thermal model driven by physical-data fusion and a simplified first-order RC equivalent circuit model. Among them, the distributed thermal model driven by physical-data fusion includes a physics-based two-dimensional distributed thermal model and a temperature calibration model based on a gated recurrent neural network, as shown in Figure 3 More specifically, the physics-based two-dimensional distributed thermal model includes the battery body spatial distribution thermal model and the positive and negative lug concentrated mass thermal model. Specifically, in the process of establishing a fast digital twin model of the battery, the following specific steps are included:

[0131] Preferably, the process of establishing the spatial distribution thermal model of the battery body includes: For large-sized soft-pack lithium-ion batteries, due to the uneven distribution of heat generation and the heat flow exchange between the tabs and the body, there is obvious temperature inconsistency in the battery body. Therefore, at least a two-dimensional distributed thermal model needs to be established to accurately simulate the temperature distribution of the battery body. This type of thermal problem obeys the typical transient heat conduction equation that satisfies the Robin boundary condition, and the expression is as follows:

[0132] ρC h D t Tk x D xx Tk y D yy T=q

[0133] Subject to the following boundary conditions:

[0134]

[0135] Where T(x,y,t) and q(x,y,t) represent the temporal and spatial distribution temperature and volumetric heat generation rate of the battery. h and k are the average volume density, specific heat capacity and anisotropic thermal conductivity of the battery. Considering the lamination manufacturing process of the soft-pack battery, it can be approximately considered that the thermal conductivity is consistent in the x and y directions. and x∈[0,w], y∈[0,l], and w and l represent the width and length of the battery. The subscripts r, l, t, and b refer to the right, left, top, and bottom boundaries of the battery, and h x ∈[h r ,-h l ] and h y ∈[h t ,-h b ] represents the equivalent heat dissipation coefficient of the battery in the x and y directions (taking into account the influence of convection and radiation heat dissipation). In this embodiment, it is stipulated that the values ​​of the heat dissipation coefficient are all positive, and "±" is used to distinguish the heat dissipation direction in the Cartesian coordinate system. T ∞,x =[T r,∞ ,T l,∞ ] and T ∞,y =[T t,∞ ,T b,∞ ] represents the external cooling environment temperature of the battery in both directions, assuming that the cooling environment around the test battery in the incubator is consistent.

[0136] In order to be able to carry out real vehicle applications, the Chebyshev-Galerkin approximation method and time-space separation technology are introduced to reduce the model order. After coordinate change, boundary condition homogenization, auxiliary temperature function solution, arbitrary temperature function solution and other processes, the PDE type control equation is transformed into multiple ODE type control equations for solution. The deformed state space expression is as follows:

[0137]

[0138] Among them, E, A, B and C are system matrices, T e is any time-invariant temperature function related to the spatial position, the system state x = (χ 00 ,x 10 ,…,X N0 ,X 01 ,X 11 ,…,X N1 ,X 0N ,X 1N ,…,χ NN ) T System input u = [q 1] T , q is the heat generation rate per unit volume. The system output y can be defined according to the actual control system requirements. Generally, the thermal management system requires the thermal model to accurately monitor the maximum temperature, maximum temperature gradient and volume average temperature of the battery. In this embodiment, T2~T4 and volume average temperature T b As the output temperature of the thermal model,

[0139]

[0140] Preferably, the process of establishing the concentrated mass thermal model of the positive and negative tabs includes: considering that the tabs are made of metal with good thermal conductivity and have a small area, it is approximately considered that the temperature distribution on the tabs is uniform, and the concentrated mass model can be used to represent its thermodynamics, and its mathematical expression is as follows:

[0141]

[0142] Among them, T t and T ∞ Refers to the lug temperature and ambient temperature respectively. t , C pt 、h t and A t are the mass, specific heat capacity, convection heat transfer coefficient and convection heat transfer area of ​​the lug respectively. t and q ct They represent the heat generation rate of the tab and the heat flux between the tab and the body, respectively. It is easy to see that as long as the tab heat generation and tab temperature are obtained in real time, the heat flux between the tab and the body can be quantified. ct In this embodiment, the inflow to the tab is negative and the outflow to the tab is positive. t It is only ohmic heat, which can be calculated by the following formula:

[0143] q t =I 2 R st =I 2 (R ot +R ct )

[0144] Among them, R ot and R ct They are the ohmic internal resistance of the tab and the contact internal resistance between the tab and the wire.

[0145] The battery body is directly connected to the two tabs and there is a significant temperature gradient, so there must be heat exchange between the two. In order to establish a control-oriented heat transfer model, this embodiment introduces the following empirical formula to calculate the heat exchange between the battery tabs and the body:

[0146]

[0147] Among them, T cb T is the temperature of the discrete cell in direct contact with the tab. Due to the small volume of the discrete cell, its temperature can be approximately considered to be equal to the tab temperature. m The temperature of the measuring point on the same side as the lug (temperature measuring point T2 or T4, for thermocouple arrangement, please refer to Figure 4 ). cb , k y and Δx are the contact coefficient between the battery tab and the body that hinders heat transfer, the thermal conductivity of the battery body in the y direction, and the temperature is determined by Tcb Change to T m The length traveled. ct It is the cross-sectional area for actual heat conduction between the tab and the body.

[0148] Preferably, the process of establishing a physically based two-dimensional distributed thermal model based on the spatially distributed thermal model of the battery body and the concentrated mass thermal model of the positive and negative tabs is as follows: when calculating the heat source term of the two-dimensional distributed thermal model developed based on the Chebyshev-Galerkin approximation method and the time-space separation technology, since the two-dimensional spatial distribution sub-model is essentially still a numerical calculation model, the heat flow exchange between the battery tab and the body needs to be superimposed on the uneven heat generation of each discrete unit of the battery body. The heat flow q between the tab and the body ct The distribution on the battery body is as follows:

[0149]

[0150] Where i,j=1,…,N g -1, k = 1, ..., M, M is the number of discrete units of the battery body, take M = (N g -1)×(N g -1). N g is the number of discrete points in the x and y directions of the battery body. For ease of processing, the number of discrete units in the x and y directions of the battery body is equal. This value can be determined based on the computing and storage capabilities of the controller in the actual application scenario. ct,k is the heat flow q ct Assigned to the battery body discrete unit (x i ,y j ) of heat. ct is the adjustment factor, indicating q ct,k The heat that actually contributes to the temperature change of the unit is in q ct,k The value range is [0,1]. t and t It is the center coordinate of the connection between the battery tab and the body.

[0151] Preferably, the forward propagation mechanism of the gated recurrent neural network model can be expressed by the following set of equations:

[0152]

[0153] Among them, σ and tanh represent gated activation functions. σ refers to the sigmoid function, expressed as σ = 1 / (1 + exp(-x)), and tanh refers to the hyperbolic tangent function, expressed as tanh = (exp(x) - exp(-x)) / (exp(x) + exp(-x)). represents a linear activation function. ⊙ and [x,x] represent element-wise multiplication and matrix concatenation, respectively.r , W z , W h and W0 are the weight matrices of the reset gate, update gate, short-term output, and output layer, respectively. Similarly, b r 、b z 、b h and b0 is the corresponding bias matrix. t and z t Refers to the outputs of the reset gate and update gate respectively. and h t Denote candidate hidden states and hidden states respectively. t Represents the output of GRU-NN. During the training phase, the mean square error (MSE) of all outputs is selected as the loss function of GRU-NN, as shown below:

[0154]

[0155] Where M and N represent the predicted output length and sequence data length of GRU-NN at a certain time step, respectively. i,j and are the temperature T at time step j respectively. i The estimated value and reference value of GRU-NN use the loss function gradient to continuously update the weights and biases in the network through back propagation. This embodiment selects the Adam optimization algorithm to minimize the loss function because Adam, as the most typical gradient-based optimization method, can adaptively adjust the learning rate without requiring too much storage space.

[0156] See also Figure 5 , taking the ternary lithium-ion battery as an example, the first-order RC equivalent circuit model is selected as its electrical sub-model. The time domain discrete form of the model is expressed as:

[0157]

[0158] U t,k =U oc,k -U p,k -R s,k I k

[0159] Among them, U t , U oc , R s and I represent the terminal voltage, open circuit voltage, ohmic internal resistance and current of the battery respectively. p and C p Represent polarization resistance and polarization capacitance respectively, C b is the current available capacity of the battery, Δt is the sampling interval, and the subscript k represents the time k.

[0160] The semi-empirical performance attenuation model is specifically expressed as:

[0161]

[0162] Among them, Q loss Indicates the battery capacity loss rate in %. exp , β rate and z are the exponential prefactor, the battery aging acceleration factor related to the current rate, and the power function factor related to the ampere-hour throughput, respectively. a is the battery aging activation energy, in J·mol -1 I rate is the battery load current, in C. R g is the ideal gas constant, which is 8.314 J·mol -1 ·K -1 . T b is the average temperature of the battery volume, in °C, which can be calculated by a two-dimensional distributed thermal model. Ah is the ampere-hour throughput of the battery during charging and discharging. In this embodiment, since the battery is fully charged and discharged each time, it should be taken as 2C. bat As the number of cycles increases, Ah = 2 × num_cycle × DOD × C bat , where num_cycle represents the number of charge and discharge cycles of the battery.

[0163] The battery fast digital twin model built for full life cycle simulation takes into account the battery's electrical, thermal, and attenuation processes. It only needs to input current and ambient temperature to sense the key internal states and parameters of the battery and output terminal voltage and key temperature information, thus realizing a fast digital twin of the battery entity's full life cycle. dt , state vector x dt and the output vector y dt They are respectively expressed as follows:

[0164]

[0165] The subscript dt represents the digital twin DT, and k represents the time step. The input vector contains the current and ambient temperature at time k, and the state vector contains the state of charge SOC, polarization voltage U at time k. p , the time coefficient of temperature χ and the state of health SOH, the output vector contains the terminal voltage U at time k t , each key temperature T and capacity loss rate Q loss The state and parameter transfer rules between sub-models are as follows: Figure 6 .

[0166] In step S2 of this embodiment, the battery attenuation equivalent cost model considers the use cost caused by the performance degradation of a single charge and discharge cycle, which is related to the price of a new battery, the resale price of the battery when retired, and the number of equivalent cycles of the battery in its entire life cycle. Among them, the residual value of the battery when retired can be calculated as follows:

[0167]

[0168] Among them, RV fre and RV EoL Represents the sales price (or residual value) of new batteries and retired batteries respectively. life Refers to the service life of the battery for an electric flying car. In this embodiment, the value is 2 years, which is a huge difference from the 8-10 year service life of a ground vehicle. dep The annual depreciation rate of the battery is 50%, which is much higher than the 20% depreciation rate of fuel vehicles. Electric vehicles, especially aviation battery packs, have serious performance degradation due to high-frequency use during service, so their value retention rate is much lower. dr represents the discount rate, which is 10%, which is also higher than the 8% of fuel vehicles. For the large-size soft-pack NCM battery cell used in this embodiment, the purchase price on the open market when it is brand new is 160 yuan, that is, RV fresh =160. In this embodiment, it is assumed that when the power battery capacity of the electric flying car decays to 85% of the initial capacity, that is, Q loss =15%, it needs to be retired compulsorily. Based on this, and combined with the battery performance degradation model, the battery ampere-hour throughput Ah in the entire life cycle tot It is expressed as:

[0169]

[0170] The battery’s ampere-hour throughput in a single UAM cycle is 2×DOD×C bat , where DOD = 1. Then the equivalent number of cycles of the battery in the whole life cycle is N cycle It is expressed as:

[0171]

[0172] From the above, we can see that the cost of using the battery due to the degradation of battery performance in a single charge and discharge cycle is expressed as follows:

[0173]

[0174] The battery charging electricity price cost model is affected by the local time electricity price and the power drawn from the charging pile side for a single time. The power drawn from the charging pile side includes the power charged into the battery and the power dissipated in the form of heat energy. Considering that the charging protocol is MCC fast charging, the power that the charging pile side needs to provide during the entire charging process is:

[0175]

[0176] P mcci =|U t (τ)I(τ)|+I(τ)(U oc (τ)-U t (τ)

[0177] Among them, E mcc Provides the power required for the entire charging process, in kWh. N is the number of constant current charging segments in the MCC charging protocol, P mcci and t mcci are the power and time required for each constant current charging stage. Therefore, the battery charging electricity cost can be expressed as:

[0178] C elec =p e E mcc

[0179] In the formula, p e The time-of-use electricity price is per kilowatt-hour.

[0180] In step S3 of this embodiment, Figure 7 As shown in the figure, the multi-objective charging optimization problem of the battery considers the multi-stage constant current charging mode, and the charging objectives consider the charging time, charging cost, maximum temperature rise and other charging sub-objectives, which are expressed as follows:

[0181] ①Charging time t mcc It is expressed as follows:

[0182]

[0183] Where N is the number of constant current charging stages, t mcc The time taken for the entire charging process, in minutes. N The duration of the entire charging process in seconds.

[0184] ②Charging cost C mcc It is expressed as:

[0185] C mcc =C age +C elec

[0186] ③The maximum temperature rise is characterized by the highest temperature in the battery plane, expressed as:

[0187] T mcc,max =max([T1 T2 T3 T4 T5])=max(T1)

[0188] Among them, T i Indicates the entire charging process T i Temperature series at the measuring points.

[0189] Usually, the calculated values ​​of the above three sub-objectives are not in the same order of magnitude. For example, the charging temperature varies from a few minutes to tens of minutes, while the battery charging cost is much less than 1 yuan. In order to eliminate the amplitude difference and improve the rationality of the objective function, each sub-objective function needs to be normalized. The normalized expressions of each sub-objective function are as follows:

[0190]

[0191] Among them, J mcc,t , J mcc,c and J mcc,tr They represent the battery charging time sub-function, charging cost sub-function and charging maximum temperature rise sub-function respectively. There is a mutually exclusive relationship between the three and they cannot reach the optimum at the same time. min and t max Respectively represent the minimum and maximum charging time of the battery. In this embodiment, t min = 0 minutes and t max = 120 minutes. min and T max are the minimum and maximum charging temperatures of the battery respectively, and T min =15℃ and T max =45°C. By introducing the weight factors of each sub-objective function, a complete multi-objective charging function can be obtained:

[0192]

[0193] α+β+γ=1

[0194] Charging constraints include time and electric and thermal operation constraints, which are described as follows:

[0195]

[0196] Among them, I min and I max Respectively represent the minimum and maximum charging current allowed by the battery, which are 0C and 6C respectively. min and U max Represents the minimum and maximum operating voltage of the battery, which are 2.75V and 4.20V respectively. min and SOC max They represent the upper and lower limits of the battery SOC operating range, which are 30% and 95% respectively. The applicable scenario for charging optimization here is the standard urban air traffic operating condition, and about 30% of the power needs to be reserved for the safe landing and emergency transfer of the flying car. The target SOC is taken as 0.95, mainly because the super fast charging under UAM conditions only considers the MCC stage, and does not include the constant voltage charging stage. SOH min and SOHmax They represent the upper and lower limits of the battery SOH, which are set to 0 and 1 respectively. The initial values ​​of each operating constraint are set as follows: SOC0 = SOC min 、SOH0=SOH max , T b0 =T ∞ .

[0197] The improved dynamic multi-population particle swarm optimization algorithm introduces the parallel strategy of the swarm intelligence algorithm, divides the entire population into multiple sub-populations, and realizes parallel calculations of multiple sub-populations without losing the optimality of the solution by designing the information interaction mechanism between sub-populations, greatly improving the algorithm's operating efficiency and global search capabilities, and can effectively avoid falling into local optimal solutions when solving complex multi-peak problems.

[0198] The sub-population is divided by the median of the fitness value in ascending order, and dynamic adjustments are made in each iteration. The particles with smaller fitness values ​​are set as "top particles" and the particles with larger fitness values ​​are set as "bottom particles". An equal number of particles are extracted from each of them to form a local PSO model to ensure that the three populations are of equal size, which are defined as the "dominant group" closer to the optimal solution, the "mixed group" of the local PSO model, and the "inferior group". The particles in the "dominant group" have smaller fitness values, inherit the position information of the optimal solution in the population, reduce the step size and conduct a more refined search. On the contrary, the particles in the "inferior group" should increase the step size, explore new local optimal solutions while approaching the global extreme value to improve the mining capacity of the sub-group. For the "mixed group", it contains both better particles and worse particles, which can not only absorb individual experience, but also share global information, and can dynamically adjust the diversity of the population. The specific implementation process of the algorithm is as follows. Figure 8 .

[0199] The update mechanism of each group is described as follows:

[0200] ① Advantageous group

[0201] The main reason why the convergence speed of the particle swarm algorithm drops significantly in the later stage of optimization is that it is difficult to get rid of the current local extreme value and obtain the global optimum. In order to enhance the ability of the dominant group to jump out of the local optimum, the Lévy flight random step length theory is introduced to continuously adjust the iteration step length, which can effectively enhance the activity and jumping ability of particles, expand the particle search range, and thus improve the particle diversity, avoiding the algorithm from falling into the local optimum. The Lévy flight step length update principle is as follows:

[0202]

[0203] in, and x g (k) respectively represent the global optimal particle position and the population global optimal particle position after Levy flight update at the kth iteration. sis the step length adjustment factor, generally α s = 0.01, search the original low probability exploration area with large and small steps to increase the uniformity of the search area, and Levy (λ) is a random search path. Considering that the Levy distribution is very complex, this embodiment uses the Mantegna algorithm to simulate its flight trajectory, and the mathematical description is as follows:

[0204]

[0205] Among them, Levy~μ=k -λ , 1<λ≤3, and μ and υ obey the following normal distribution:

[0206]

[0207] Among them, the standard deviation σ μ and σ υ Determined by the following formulas:

[0208]

[0209] Where Γ(·) is the gamma function, and the constant χ is generally taken as 1.5.

[0210] ② Disadvantaged group

[0211] The particles in the "disadvantaged group" are far away from the optimal solution of the problem, and fewer particles can be retained in the end. By searching for the optimal solution in the entire space through mutation and finding various possible solution areas, it can effectively avoid premature convergence and increase the diversity of sub-populations. The following formula is used to determine whether a particle has mutated:

[0212]

[0213] GV(x)=x×(1+N(0,1))

[0214] Where r and N(0,1) are random numbers between [0,1], x is the original parameter value, and GV(x) is the mutated value. max is the maximum number of iterations. When the above inequality holds, the Gaussian mutation operation can be performed on the particles. After multiple iterations, the probability of judging that the above inequality holds gradually decreases, and the probability of mutation also gradually decreases. The mutation operation can make particles mutate with a greater probability in the initial stage, expand the search range of particles in the solution space, and ensure the diversity of the particle group. In addition, the concept of mixed particles is introduced to improve the traditional speed update formula, and the mixed particles are recorded as p mix (k), mixed particle p mix The dimension of (k) is a random mixture of the historical optimal values ​​of each particle. Therefore, the position and speed update method of the disadvantaged group particles can be expressed as:

[0215] vid (k+1)=wv id (k)+c1r1(p id (k)-x id (k))+c2r2(p gd (k)-x id (k))+c3r3(p mix (k)-x id (k))

[0216] x id (k+1)=x id (k)+v id (k+1)

[0217] Where c3 is the hybrid learning factor and r3 is a random number between (0,1).

[0218] ③Mixed group

[0219] The mixed group is between the dominant group and the inferior group. Due to the large differences between individuals in the early stage of optimization, the focus is on cognition and communication between groups. In the later stage of optimization, it is necessary to strengthen the guidance of the global extreme value to encourage particles to gather around the optimal solution. Therefore, cosine and sine functions are introduced into the learning factor to make the self-learning factor monotonically decrease and the population learning factor monotonically increase. Based on this, the speed update formula of the mixed group can be expressed as:

[0220]

[0221] ④ Nonlinear inertia weight

[0222] Since the inertia weight factor affects the performance of the PSO algorithm, a larger value is conducive to global search, while a smaller weight is conducive to accelerating the convergence speed of the algorithm and enhancing the local search ability. Therefore, a larger weight factor can be set in the early stage of the search, and a smaller weight factor can be set in the later stage of the search. This nonlinear inertia weight can be used to balance the local and global search capabilities of the entire population. The weight setting is implemented by referring to the following formula:

[0223]

[0224] In step S4 of this embodiment, the influencing factors of the optimal charging solution include charging cut-off voltage, temperature optimization sub-goals, etc. Battery charging is a typical multi-objective optimization problem. Under given resources and conditions, it is impossible to achieve the optimal of each optimization sub-goal at the same time. It is necessary to use the Pareto frontier theory to scientifically and reasonably determine a balanced charging solution that meets the battery operation goals. The optimal charging solution for urban air traffic operation scenarios needs to balance the weights among the three optimization sub-goals of charging time, charging cost and maximum temperature rise to achieve the safe, efficient and economical charging goals of electric flying cars. Through the analysis of the Pareto frontier curve and the factors affecting the optimal charging, the optimal charging solution for specific application scenarios can be clarified.

[0225] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. 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 solution of the present invention can be modified or replaced by equivalents without departing from the purpose and scope of the technical solution, which should be included in the scope of the claims of the present invention.

Claims

1. A battery charging optimization method for electric flying vehicles taking into account operational safety and economy, characterized in that: The method comprises the following steps: S1: Establish a distributed electric-thermal coupling model and a semi-empirical performance degradation model driven by physical-data fusion for power batteries, and couple the two to build a fast digital twin model of batteries for full life cycle simulation; S2: Considering the equivalent cost of battery degradation and the cost of charging electricity, a battery charging economic model is established; S3: Based on the battery fast digital twin model and charging economy model, multiple optimization sub-objectives are considered to construct a multi-objective charging optimization problem and use advanced optimization algorithms to achieve the optimal solution; S4: Analyze the influencing factors of the optimal charging scheme based on the Pareto frontier theory, and determine the multi-stage constant current optimal charging scheme for urban air traffic operation scenarios.

2. The battery charging optimization method for electric flying vehicles taking into account operational safety and economy according to claim 1, characterized in that: In step S1, the distributed electrothermal coupling model driven by physics-data fusion includes a distributed thermal model driven by physics-data fusion and a first-order RC equivalent circuit model, the distributed thermal model driven by physics-data fusion consists of a physics-based two-dimensional distributed thermal model and a temperature calibration model based on a gated recurrent neural network, and the physics-based two-dimensional distributed thermal model consists of a battery body spatial distribution thermal model and positive and negative lug concentrated mass thermal models; In the battery fast digital twin model for full life cycle simulation, which is constructed by coupling the distributed electro-thermal coupling model driven by physical-data fusion and the semi-empirical performance degradation model, the input vector u of the battery digital twin model is dt , state vector x dt and the output vector y dt They are respectively expressed as follows: Among them, the subscript dt represents the digital twin DT, k represents the time step; the input vector contains the current and ambient temperature at time k, and the state vector contains the state of charge SOC and polarization voltage U at time k. p , the time coefficient of temperature χ and the state of health SOH, the output vector contains the terminal voltage U at time k t , each key temperature T and capacity loss rate Q loss .

3. The battery charging optimization method for electric flying vehicles taking into account operational safety and economy according to claim 2, characterized in that: In step S1, the process of establishing a distributed electrothermal coupling model driven by physical-data fusion is as follows: The spatial distribution thermal model of the battery body based on the physical two-dimensional distributed thermal model obeys the transient heat conduction equation that satisfies the Robin boundary condition. The expression is as follows: ρC h D t T-k x D xx T-k y D yy T=q Subject to the following boundary conditions: Where T(x, y, t) and q(x, y, t) represent the temporal and spatial distribution temperature and volumetric heat generation rate of the battery, ρ, C h and k are the battery volume average density, specific heat capacity and anisotropic thermal conductivity; and x∈[0,w],y∈[0,l], where w and l represent the width and length of the battery; subscripts r, l, t and b refer to the right, left, top and bottom boundaries of the battery, h x ∈[h r , -h l ] and h y ∈[h t , -h b ] represents the equivalent heat dissipation coefficient of the battery in the x and y directions; the values ​​of the heat dissipation coefficient are all positive, and "±" is used to distinguish the heat dissipation direction in the Cartesian coordinate system; T ∞,x =[T r，∞ , T l，∞ ] and T ∞,y =[T t，∞ , T b，∞ ] represents the external cooling environment temperature of the battery in both directions, assuming that the cooling environment around the battery in the temperature chamber is consistent; The Chebyshev-Galerkin approximation method and time-space separation technology are introduced to reduce the model order. After coordinate change, boundary condition homogenization, auxiliary temperature function solution, and arbitrary temperature function solution, the PDE type control equation is transformed into multiple ODE type control equations for solution. The deformed state space expression is as follows: Where E, A, B and C are system matrices, Te is an arbitrary time-invariant temperature function related to the spatial position, and the system state x = (χ 00 , χ 10 ,…,χ N0 , χ 01 , χ 11 ,…,χ N1 , χ 0N , χ 1N ,…,χ NN ) T ; System input u = [q 1] T , q is the heat generation rate per unit volume; the system output y is defined according to the actual control system requirements; extract T2~T4 and the volume average temperature T b As the output temperature of the thermal model, that is: The mathematical expressions of the positive and negative lug concentrated mass thermal models based on the physical two-dimensional distributed thermal model are as follows: Among them, T t and T ∞ Refers to the lug temperature and ambient temperature respectively; m t , C pt 、h t and A t are the mass, specific heat capacity, convection heat transfer coefficient and convection heat transfer area of ​​the lug respectively; q t and q ct They represent the heat generation rate of the tab and the heat flux between the tab and the body respectively; for q ct , the flow into the pole ear is negative, and the flow out of the pole ear is positive; the heat generation of the pole ear q t Ohmic heat only, calculated by the following formula: q t =I 2 R st =I 2 (R ot +R ct ) Among them, R ot and R ct They are the ohmic internal resistance of the tab and the contact internal resistance between the tab and the wire; An empirical formula is introduced to calculate the heat exchange between the battery tab and the body: Among them, T cb is the temperature of the discrete battery unit in direct contact with the tab, which is approximately equal to the tab temperature; T m is the temperature of the measuring point on the same side as the lug, h cb , k y and Δx are the contact coefficient between the battery tab and the body that hinders heat transfer, the thermal conductivity of the battery body in the y direction, and the temperature is determined by T cb Change to T m The length traveled, A ct It is the cross-sectional area of ​​actual heat conduction between the tab and the body; Heat flow qc between the tab and the body t The distribution on the battery body is as follows: Where i, j = 1, ..., N g -1, k = 1, ..., M, M is the number of discrete units of the battery body, and M = (N g -1)×(N g -1), N g is the number of discrete points in the x and y directions of the battery body; the number of discrete units in the x and y directions of the battery body is equal, and this value can be determined according to the computing and storage capabilities of the controller in the actual application scenario; q ct,k is the heat flow q ct Assigned to the battery body discrete unit (x i ,y j ) of heat, β ct is the adjustment factor, indicating q ct,k The heat that actually contributes to the temperature change of the unit is in q ct,k The proportion in the range is [0,1], x t and t is the center coordinate of the connection between the battery tab and the body; The forward propagation mechanism of the gated recurrent neural network model is represented by the following set of equations: Among them, σ and tanh represent gated activation functions, σ refers to the sigmoid function, expressed as σ = 1 / (1 + exp(--x)), tanh refers to the hyperbolic tangent function, expressed as tanh = (exp(x) - exp(-x)) / (exp(x) + exp(-x)), represents a linear activation function, ⊙ and [x, x] represent element-wise multiplication and matrix concatenation, respectively, and W r , W z , W h and W0 are the weight matrices of the reset gate, update gate, short-term output and output layer respectively; b r 、b z 、b h and b0 is the corresponding bias matrix, r t and z t Refers to the outputs of the reset gate and update gate respectively, and h t Denote candidate hidden states and hidden states respectively, y t Represents the output of GRU-NN; During the training phase, the average mean square error (MSE) of all outputs is selected as the loss function of GRU-NN: Among them, M and N represent the predicted output length and sequence data length of GRU-NN at a certain time step, respectively, and y i，j and are the temperature T at time step j respectively. i The estimated value and reference value of GRU-NN use the loss function gradient to continuously update the weights and biases in the network through back propagation; The first-order RC equivalent circuit model is selected as its electrical sub-model. The time domain discrete form of the model is expressed as: IN t，k =U oc，k -IN p，k -R s，k AND k Among them, U t , U oc , R s and I represent the terminal voltage, open circuit voltage, ohmic internal resistance and current of the battery respectively. p and C p Represent polarization resistance and polarization capacitance respectively, C b is the current available capacity of the battery, Δt is the sampling interval, and the subscript k represents the time k.

4. The battery charging optimization method for electric flying vehicles taking into account operational safety and economy according to claim 2, characterized in that: In step S1, the semi-empirical performance degradation model is expressed as: Among them, Q loss Indicates the battery capacity loss rate in %. exp , β rate and z are respectively the exponential prefactor, the battery aging acceleration factor related to the current rate, and the power function factor related to the ampere-hour throughput; E a is the battery aging activation energy, in J·mol -1 ;I rate is the battery loading current, in C; R g is the ideal gas constant; T b is the average temperature of the battery volume, in °C, calculated by a two-dimensional distributed thermal model; Ah is the ampere-hour throughput of the battery during charge and discharge. As the number of cycles increases, Ah = 2 × num_cycle × DOD × C bat , where num_cycle represents the number of charge and discharge cycles of the battery.

5. The battery charging optimization method for electric flying vehicles taking into account operational safety and economy according to claim 2, characterized in that: In step S2, the battery attenuation equivalent cost model considers the use cost caused by the performance degradation of a single charge and discharge cycle, which is related to the price of a new battery, the resale price of the battery when retired, and the number of equivalent cycles of the battery in its entire life cycle. The residual value of the battery when retired is calculated as follows: Among them, RV fresh and RV EOL Represent the sales price or residual value of new batteries and retired batteries, respectively, life Refers to the service life of the battery for electric flying cars; dep represents the annual depreciation rate of the battery; dr represents the discount rate; based on this, combined with the battery performance attenuation model, the battery ampere-hour throughput Ah in the entire life cycle tot It is expressed as: Among them, Q lo ss represents the mandatory retirement threshold of the power battery capacity attenuation of the electric flying car; The equivalent number of cycles of the battery in its entire life cycle is N cycle It is expressed as: The cost of battery performance degradation during a single charge and discharge cycle is expressed as follows: The amount of electricity that the charging pile needs to provide during the entire charging process is: P mcci =|U t (τ)I(τ)|+I(τ)(U oc (t)-U t (t)) Among them, E mcc Provides the power required for the entire charging process, in kWh; N is the number of constant current charging segments in the MCC charging protocol, P mcci and t mcci They are the power and time required for each constant current charging stage; The battery charging electricity cost is expressed as: C elec =p e E mcc In the formula, p e The time-of-use electricity price is per kilowatt-hour.

6. The battery charging optimization method for electric flying vehicles taking into account operational safety and economy according to claim 5, characterized in that: In step S3, the multi-objective charging optimization problem of the battery considers a multi-stage constant current charging mode, and the multi-optimization charging sub-objectives include charging time, charging cost and maximum temperature rise, which are expressed as follows: ①Charging time t mcc It is expressed as follows: Where N is the number of constant current charging stages, t mcc The time taken for the entire charging process, in minutes; t N The duration of the entire charging process, in seconds; ②Charging cost C m cc is represented by: C mcc =C age +C elec ③The maximum temperature rise is characterized by the highest temperature in the battery plane, expressed as: T mcc，max =max([T1 T2 T3T4 T5])=max(T1) Among them, T i Indicates the entire charging process T i Temperature series at the measuring points; Each sub-objective function is normalized, and the normalized expression of each sub-objective function is as follows: Among them, J mcc,t , J mcc,c and J mcc,tr They represent the battery charging time sub-function, charging cost sub-function and charging maximum temperature rise sub-function respectively. There is a mutual exclusive relationship between the three and they cannot reach the optimum at the same time; t min and t max Respectively represent the minimum and maximum charging time of the battery; T min and T max are the minimum and maximum charging temperatures of the battery respectively; By introducing the weight factors of each sub-objective function, we can get the complete multi-objective charging function: α+β+γ=1 Charging constraints include time and electric and thermal operation constraints, which are described as follows: Among them, I min and I max Respectively represent the minimum and maximum charging current allowed by the battery; U min and U max Respectively represent the minimum and maximum values ​​of the battery operating voltage; SOC min and SOC max Respectively represent the upper and lower limits of the battery SOC operating range; SOH min and SOH max They represent the upper and lower limits of the battery SOH respectively.

7. The battery charging optimization method for electric flying vehicles taking into account operational safety and economy according to claim 6, characterized in that: In step S3, the optimal solution process of the multi-objective charging optimization problem using an advanced optimization algorithm includes: The sub-population is divided by the median of the fitness value in ascending order, and dynamic adjustment is performed in each iteration; the particles with fitness values ​​less than the median are set as "top particles", and the particles with fitness values ​​greater than the median are set as "bottom particles", and the same number of particles are extracted from them to form a local PSO model, which are defined as the "dominant group", the "mixed group" of the local PSO model, and the "inferior group", among which the "dominant group" contains all top particles, the "inferior group" contains all bottom particles, and the "mixed group" of the local PSO model includes both top particles and bottom particles; For the update of the dominant group, the random step length theory of Levy flight is introduced to continuously adjust the iteration step length. The updating principle of Levy flight step length is as follows: in, and x g (k) represents the global optimal particle position and the global optimal particle position of the population after Levy flight update at the kth iteration; α s is the step size adjustment factor, Levy(λ) is the random search path; The flight trajectory is simulated using the Mantegna algorithm, and the mathematical description is as follows: Among them, Levy~μ=k -λ , 1<λ≤3, and μ and v obey the following normal distribution: Among them, the standard deviation σ μ and σ v Determined by the following formulas: Where Γ(·) is the gamma function; For the update of the disadvantaged group, the optimal solution is searched in the entire space by mutation, and whether the particle mutates is determined according to the following formula: GV(x)=x×(1+N(0,1)) Where r and N(0,1) are random numbers between [0,1], x is the original parameter value, and GV(x) is the mutated value; k max is the maximum number of iterations; when the above inequality holds, the Gaussian mutation operation is performed on the particles; after multiple iterations, the probability of judging that the above inequality holds gradually decreases, and the probability of mutation gradually decreases; The concept of hybrid particles is introduced to improve the traditional velocity update formula, and the hybrid particles are denoted as p mix (k), mixed particle p mix The dimension of (k) is a random mixture of the historical optimal values ​​of each particle. Therefore, the position and velocity update method of the disadvantaged group particles is expressed as: v id (k+1)=wv id (k)+c1r1(p id (k)-x id (k))+c2r2(p gd (k)-x id (k))+c3r3(p mix (k)-x id (k)) x id (k+1)=x id (k)+v id (k+1) In the formula, c3 is the hybrid learning factor, r3 is a random number between (0, 1); For the update of the mixed group, cosine and sine functions are introduced into the learning factor to make the self-learning factor monotonically decreasing and the population learning factor monotonically increasing. The speed update formula of the mixed group is expressed as: For the setting of nonlinear inertia weight during update, refer to the following formula:

Images