Water-air amphibious unmanned aerial vehicle compound control method based on motion state separation and multi-target distribution

Abstract

The invention provides a water-air amphibious unmanned aerial vehicle compound control method based on motion state separation and multi-target distribution. The method comprises the steps of constructing a motion mathematical model and a disturbance model of an unmanned aerial vehicle, and designing an expansion observer based on state separation. Designing a nonlinear mapping rule to optimize the target inclination angles of different modes of the unmanned aerial vehicle; and designing a multi-objective cost function, solving an expected thrust and an expected inclination angle by using an improved genetic algorithm, and calculating an output thrust and a blade rotation angular velocity of each motor. The extended state observer based on state separation accurately estimates slowly changing lumped disturbance which is composed of persistent force / torque such as water flow and wind power and really needs to be compensated by a controller by augmenting an oscillatory motion state; through a multi-target cost function, the problem of single energy consumption minimization is controlled and distributed, so that the inclination angle tracking precision and the thrust change smoothness are achieved; and a high-quality global optimal solution is found at a higher probability by improving the genetic algorithm.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of unmanned aerial vehicle control, and particularly relates to a water-air amphibious unmanned aerial vehicle composite control method based on motion state separation and multi-target allocation. BACKGROUND

[0002] In recent years, with the increasing demand for water operations, water-air amphibious unmanned aerial vehicles as a new emerging high-efficiency platform have been widely concerned. However, the existing technology still has many limitations and objective defects in the design and control of water-air amphibious unmanned aerial vehicles, especially in dealing with complex and variable cross-domain environments.

[0003] Firstly, there is a fundamental defect in the state estimation and disturbance processing of the unmanned aerial vehicle. The existing control method, such as the standard extended state observer, takes the core idea of considering all unmodeled dynamics and external disturbances as a single, slowly varying lumped disturbance for estimation. This method is effective in dealing with slowly varying disturbances such as air resistance or deep water flow. However, when the unmanned aerial vehicle is sailing on the water surface or near the water surface, the motion measured by the sensor is the superposition of the response of the unmanned aerial vehicle to the slowly varying persistent force and the response to the highly oscillatory force. This directly destroys the core assumption of the slowly varying disturbance on which the extended observer relies, resulting in the observer producing sharp jitter to track the rapidly changing wave force, and ignoring the compensation for the slowly varying disturbance force that really needs to be compensated by the controller, greatly affecting the stability of the system.

[0004] Secondly, in terms of control allocation for water-air variable structure unmanned aerial vehicles, the existing technology also has obvious shortcomings. Simple methods such as pseudo-inverse method cannot handle the motor saturation problem. While the more advanced standard quadratic programming method, the optimization objective is usually one-sided, i.e. only considers minimizing motor thrust energy consumption. This design completely ignores the huge dynamic cost brought by the driving tilt structure, including the high energy consumption of the rudder itself, the mechanical wear caused by rapid reciprocating motion, and the system instability introduced by the dramatic change in configuration.

[0005] Therefore, the solution obtained by the existing allocation method is often a theoretically optimal solution, but in actual application, it may lead to reduced flight quality and shortened component life, and is not a real system-level optimal solution. SUMMARY

[0006] In view of the deficiencies of the prior art, the present application provides a water-air amphibious unmanned aerial vehicle composite control method based on motion state separation and multi-target allocation to solve the problems in the background art.

[0007] The technical scheme of the present application is as follows: a water-air amphibious unmanned aerial vehicle composite control method based on motion state separation and multi-target allocation, comprising the following steps:

[0008] S1), constructing a motion mathematical model of the unmanned aerial vehicle;

[0009] S2), constructing a disturbance model for simulating the oscillation disturbance of irregular wave force on the unmanned aerial vehicle near the water surface to generate oscillation displacement and oscillation Euler angle influence on the motion state of the unmanned aerial vehicle;

[0010] S3), designing a position and attitude extended observer based on state separation;

[0011] S4), designing a nonlinear mapping rule of roll motion and pitch angle, and optimizing the target pitch angle of different modes of the unmanned aerial vehicle;

[0012] S5), constructing a multi-objective cost function, and using an improved genetic algorithm to solve the expected thrust and expected pitch angle, and calculating the output thrust of each motor and the blade rotation angular velocity.

[0013] As preferred, in step S1), the motion mathematical model of the unmanned aerial vehicle is constructed, specifically including the following steps:

[0014] S11), establishing a body coordinate system and a space coordinate system for the unmanned aerial vehicle;

[0015] S12), establishing a position dynamics model of the unmanned aerial vehicle;

[0016] S13), establishing an attitude dynamics model of the unmanned aerial vehicle.

[0017] As preferred, in step S2), the oscillation disturbance model suffered by the unmanned aerial vehicle is:

[0018]

[0019] Y = CH; (10)

[0020] In the formula, is the first derivative of the oscillation motion state vector H; A is the system state matrix, K is the system input matrix; C is the system output matrix; u is the unit white noise input; Y is the complete oscillation motion state output vector.

[0021] As preferred, in step S3), the position extended observer based on state separation is represented as:

[0022]

[0023] In the formula, e p is a three-dimensional error vector of the estimated value and the true value of the position; is a three-dimensional estimated value vector of the position state of the unmanned aerial vehicle; y p represents the position state of the unmanned aerial vehicle; z p0 ,z p1 ,z p2The position and state vectors x p0 ,x p1 ,x p2 The estimated values ​​are the estimates of the oscillating displacement, three-dimensional position, and three-dimensional velocity of the UAV. z p0 ,z p1 ,z p2 The first derivative; z p3 It is an estimate of all positional disturbances experienced by the UAV; γ p0 ,γ p1 ,γ p2 ,γ p3 All are three-dimensional diagonal gain matrices of the position observer; A p =diag(A1,A2,A3); A1,A2,A3 are the state matrices of degrees of freedom; ζ i The damping ratio associated with the i-th degree of freedom is used to match the sharpness of the wave spectrum peaks; ω p represents the peak energy frequency of the wave spectrum of the target water area; diag(·) denotes a diagonal matrix with its elements as the diagonal; For H p The first derivative; η is the oscillating displacement vector; i Let ξ be the i-th degree of freedom in the oscillating motion state; R represents the rotation matrix from the body to the spatial coordinate system; F represents the thrust of the UAV in the body coordinate system; C p = [C1C2C3] is the output matrix for each degree of freedom, satisfying C i =[1 0].

[0024] Preferably, in step S3), the attitude expansion observer based on state separation is represented as:

[0025]

[0026] In the formula, e Θ It is the three-dimensional error vector between the estimated and true attitude values; It is a three-dimensional estimated vector of the UAV's attitude state; y Θ Indicates the attitude state of the drone; z Θ0 ,z Θ1 ,z Θ2 These are the attitude state vectors x and x. Θ0 ,x Θ1 ,x Θ2 The estimated values ​​are the estimated values ​​of the oscillating Euler angles, three-dimensional Euler angles, and three-dimensional angular velocity of the UAV. z Θ0 ,z Θ1 ,z Θ2the first derivative of H Θ3 is the estimation of all the attitude disturbance moments experienced by the UAV; γ Θ0 ,γ Θ1 ,γ Θ2 ,γ Θ3 are three-dimensional diagonal gain matrices of the attitude observer; A Θ = diag(A4, A5, A6); is the first derivative of H Θ ; is the oscillatory Euler angle vector; C Θ = [C4C5C6]; J is the inertia matrix; τ is the moment of force generated by the propellers on the body axes.

[0027] As preferred, in step S4), when the UAV is navigating in water, it is divided into two modes as follows:

[0028] Symmetrical mode of left and right roll angles: applied when the UAV is in motion, at this time the roll angles satisfy β l + β r = 0;

[0029] Asymmetrical mode of left and right roll angles: applied when the UAV is in hover in water, at this time the roll angles should satisfy β l ≠ β r , and only when β l = β r = 0, the roll angles are equal; in this state, the roll angle of one side of the UAV rotates, and the roll angle of the other side changes according to the actual situation.

[0030] As preferred, in step S4), the nonlinear mapping rule of roll motion and roll angle is designed according to the roll angle and roll angle velocity, that is:

[0031] β d = β max sat(h(e φ ,p), 1, -1); (34)

[0032]

[0033] In the formula, β d represents the roll angle designed according to the roll motion; β max represents the maximum value of the preset roll angle; h(·) is the core part of the entire nonlinear mapping rule design; tanh(·) is the hyperbolic tangent function; sat(·) is the saturation function; ω nl , ω l are the weights of the nonlinear part and the linear part, respectively; k nl , k l are the gain coefficients of the nonlinear and linear inputs, respectively; f(eφ p) represents a linear combination term of the current roll angle error e φ and the current roll angle velocity p; φ is the current roll angle; φ d is the desired roll angle of the upper controller; k φ is the proportional gain of the roll angle error, used to adjust the response strength of the upper controller to static or slowly varying attitude deviation; k p is the differential gain of the roll angle velocity, used to provide system damping to suppress rapid oscillation of the attitude and ensure smoothness of dynamic process.

[0034] As preferred, in step S4), when the UAV is in the left-right roll angle symmetric mode, the target roll angle β target is:

[0035]

[0036] When the UAV is in the left-right roll angle asymmetric mode, the target roll angle β target is:

[0037]

[0038] wherein β ltar , β rtar represent the target roll angle of the left roll angle β l and the right roll angle β r respectively.

[0039] As preferred, in step S5), the multi-objective cost function designed according to the control allocation relationship and the target roll angle β target is:

[0040]

[0041] wherein J(β) is the multi-objective cost function; β represents the rudder roll angle vector currently optimized; w uj represents the weight coefficient of the single motor thrust in the control energy consumption penalty term; u fj represents the output thrust of the jth motor; Δβ tra represents the tracking error of the roll angle; W β represents the weight matrix of the roll angle tracking penalty term; Δu var represents the rate of change of the motor thrust; W △u represents the weight matrix of the motor thrust rate change penalty term; Δu var represents the rate of change of the motor thrust; Δβ var represents the rate of change of the roll angle; W △β represents the weight matrix of the roll angle rate change penalty term; the superscript T represents the transpose operation of the vector; u prerepresents the motor thrust vector of the last control cycle; β pre represents the rudder angle vector of the last control cycle; β max , β min respectively represent the physical upper and lower limits of the rudder angle.

[0042] The beneficial effects of the present application are:

[0043] 1. The state separation-based extended state observer of the present application can accurately separate the reciprocating motion caused by waves from the total measurement of the sensor by augmenting the oscillatory motion state, so that the extended state of the observer can no longer be disturbed by the wave oscillation effect, thereby accurately estimating the slowly varying disturbance composed of water flow, wind force and other persistent forces / torques that truly need to be compensated by the controller.

[0044] 2. The present application dynamically associates the gain with the real-time depth of the unmanned aerial vehicle and the wave environment, so that the observer can fully separate the oscillatory state when near the water surface, and automatically and smoothly transform into a linear extended observer that focuses on estimating the slowly varying disturbance when submerged to a certain depth or flying in the air.

[0045] 3. The present application constructs a nonlinear mapping rule between the roll attitude and the rudder angle, so that the unmanned aerial vehicle can no longer passively respond to the torque command, but can actively predict and utilize its unique variable structure advantage to quickly and smoothly resist external disturbances, significantly improving the attitude stability in complex water conditions.

[0046] 4. The multi-objective cost function of the present application improves the control allocation from a single energy minimization problem to four core indicators of motor energy consumption, rudder angle tracking accuracy, thrust change smoothness, and rudder motion loss; the solution is a globally optimal execution strategy that achieves the best balance between energy efficiency, stability, smoothness, and mechanical life.

[0047] 5. The improved genetic algorithm of the present application uses a low-dimensional coding centered on the rudder angle and introduces heuristic population initialization based on the target rudder angle, effectively avoiding blind search of the algorithm and greatly speeding up the convergence speed; in addition, an adaptive adjustment mechanism is introduced in the crossover and mutation operators, so that the algorithm can perform extensive global exploration in the early evolution and fine local optimization in the later period, effectively balancing exploration and convergence, and can find a high-quality global optimal solution with a higher probability. BRIEF DESCRIPTION OF DRAWINGS

[0048] Figure 1 is a flowchart of the method of the present application;

[0049] Figure 2 is a schematic diagram of the definition of the coordinate system of the unmanned aerial vehicle of the present application;

[0050] Figure 3 This is a schematic diagram of the left-right tilt angle symmetrical mode of the UAV of the present invention;

[0051] Figure 4 This is a schematic diagram of the asymmetrical tilt mode of the UAV of the present invention.

[0052] Figure 5 This is a flowchart illustrating the genetic algorithm of the present invention. Detailed Implementation

[0053] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings:

[0054] like Figure 1 As shown, this embodiment provides a composite control method for amphibious unmanned aerial vehicles based on motion state separation and multi-target allocation, including the following steps:

[0055] S1) Constructing a motion mathematical model for unmanned vehicles; specifically including the following steps:

[0056] S11) Establish a volume coordinate system and a spatial coordinate system for the UAV, such as Figure 2 As shown;

[0057] S12) Establish the position dynamics model of the UAV, namely:

[0058]

[0059] In the formula, P represents the three-dimensional position of the UAV in the spatial coordinate system; Let P be the second derivative, representing the acceleration of the UAV in the spatial coordinate system; R represents the rotation matrix from the airframe to the spatial coordinate system; F represents the thrust of the UAV in the airframe coordinate system; G represents the effect of gravity or buoyancy on the UAV in the air or water; ρ is the density of water or air; V is the displacement of the UAV; g is the acceleration due to gravity; F D The unknown perturbation force experienced by the drone is represented by m; m is the mass of the drone.

[0060] S13) Establish the attitude dynamics model of the UAV.

[0061]

[0062] In the formula, Θ represents the Euler angle of the UAV; Let θ be the second derivative of the Euler angle Θ, representing the three-axis angular acceleration of the UAV; Let be the first derivative of Θ, representing the three-axis angular velocities of the UAV; J is the inertial matrix; G is the first derivative of Θ. a The gyroscopic torque acting on the drone; τ D τ is the unknown disturbance torque experienced by the UAV; τ is the torque generated by the propeller on the airframe shaft.

[0063] S2), construct a disturbance model for simulating the oscillation disturbance of irregular wave force on the unmanned aerial vehicle near the water surface to generate oscillation displacement and oscillation Euler angle influence on the motion state of the unmanned aerial vehicle; specifically as follows:

[0064] S21), set the oscillation motion state ξ caused by the wave force on the unmanned aerial vehicle six degrees of freedom as:

[0065]

[0066] In the formula, is the three-dimensional oscillation displacement of the unmanned aerial vehicle in the spatial coordinate system x, y, z axis; is the three-dimensional oscillation Euler angle of the unmanned aerial vehicle in x, y, z axis;

[0067] S22), for the i-th degree of freedom of the oscillation motion state ξ i (i = 1, 2, 3, 4, 5, 6), the simulation expression is:

[0068]

[0069] In the formula, and are the second-order and first-order derivatives of ξ i , respectively, representing the velocity and acceleration of the oscillation motion; ζ i is the damping ratio related to the i-th degree of freedom, used to match the peak sharpness of the wave spectrum; ω p is the energy peak frequency of the wave spectrum of the target water area; k i is the input gain of the i-th degree of freedom; u is a unit white noise input, used to represent the noise input;

[0070] S23), let the state variable η 1i = ξ i , then and substitute it into equation (5) to obtain:

[0071]

[0072] Rewrite equation (7) into matrix form:

[0073]

[0074] At the same time, the oscillation motion state output equation of this degree of freedom is:

[0075]

[0076] In the formula, are the state variables η 1i , η2i the first derivative of the state vector H; y i an output representing the state of the oscillation motion;

[0077] S24), record the state vector Stack the state vectors of its six degrees of freedom, and there is a complete oscillation motion state vector Complete oscillation motion state output vector Thus, the perturbation model is obtained as:

[0078]

[0079] Y = CH; (10)

[0080] The above formula satisfies:

[0081]

[0082] Wherein: the first derivative of the state vector H; A is the system state matrix, and is a full rank matrix; diag(·) represents a diagonal matrix with elements as the diagonal line; Each degree of freedom state matrix A1,…,A6 satisfies K is the system input matrix; Each degree of freedom input matrix K1,…,K6 satisfies C is the system output matrix; Each degree of freedom output matrix C1,…,C6 satisfies C i =

[10] , (i = 1,…,6).

[0083] S3), design a position and attitude extended observer based on state separation; Specifically, the following steps are included:

[0084] S31), obtain the oscillation displacement vector from the oscillation motion state equation in formula (9) and (10) Oscillation displacement output vector Then the oscillation displacement state equation is:

[0085]

[0086] Y p = C p H p ; (13)

[0087] Wherein, the first derivative of the state vector H; A is the system state matrix, and is a full rank matrix; diag(·) represents a diagonal matrix with elements as the diagonal line; Each degree of freedom state matrix A1,…,A6 satisfies p C p = [C1C2C3];

[0088] S32), according to the oscillation Euler angle vector Oscillation Euler angle output vector​ Then the oscillatory bit Euler angle state equation is:

[0089]

[0090] Y Θ =C Θ H Θ ; (15)

[0091] where, is the first derivative of H Θ ; A Θ =diag(A4,A5,A6); C Θ =[C4 C5 C6];

[0092] S33), according to the dynamic model and the disturbance model of the unmanned aerial vehicle according to formulas (1), (2) and (3), a three-dimensional position state vector x p0 =H p , x p1 =P, Then the state space equation of the position dynamic model of the unmanned aerial vehicle is:

[0093]

[0094] In the formula, respectively represent the first derivatives of the position state vectors x p0 , x p1 , x p2 , x p3 ; y p represents the position state of the unmanned aerial vehicle; δ p is the position measurement accuracy error; f pD (·) represents all the position disturbance force vectors suffered by the unmanned aerial vehicle, and satisfies f pD (·) = G - F D ;

[0095] S34), a three-dimensional attitude state vector x Θ0 =H Θ , x Θ1 =Θ, Then the state space equation of the attitude dynamic model of the unmanned aerial vehicle is

[0096]

[0097] In the formula, respectively represent the first derivatives of the attitude error vectors x Θ0 , x Θ1 , x Θ2 ; y Θ represents the attitude state of the unmanned aerial vehicle; δΘ is the attitude measurement accuracy error; f ΘD (·) represents all the attitude disturbance moments experienced by the UAV, satisfying

[0098] S35), the optimized position extended observer is expressed according to formula (16):

[0099]

[0100] In the formula, e p is a three-dimensional difference vector of the estimated value and the true value of the position; is a three-dimensional estimated value vector of the position state of the UAV; z p0 ,z p1 ,z p2 are the estimated values of the position state vector x p0 ,x p1 ,x p2 , that is, the estimated values of the oscillation displacement, three-dimensional position and three-dimensional velocity of the UAV; are the first-order derivatives of z p0 ,z p1 ,z p2 ; z p3 is the estimated value of all the position disturbance forces experienced by the UAV;

[0101] γ p0 ,γ p1 ,γ p2 ,γ p3 are three-dimensional diagonal gain matrices of the position observer;

[0102] The optimized attitude extended observer is expressed according to formula (17):

[0103]

[0104] In the formula, e Θ is a three-dimensional error vector of the estimated value and the true value of the attitude; is a three-dimensional estimated value vector of the attitude state of the UAV; z Θ0 ,z Θ1 ,z Θ2 are the estimated values of the attitude state vector x Θ0 ,x Θ1 ,x Θ2 , that is, the estimated values of the oscillation Euler angle, three-dimensional Euler angle and three-dimensional angular velocity of the UAV; are the first-order derivatives of z Θ0 ,z Θ1 ,z Θ2 ; z Θ3 is the estimated value of all the attitude disturbance moments experienced by the UAV; γΘ0 ,γ Θ1 ,γ Θ2 ,γ Θ3 are three-dimensional diagonal gain matrices of the attitude observer;

[0105] The slowly-varying aggregate disturbance force z p3 and the slowly-varying aggregate disturbance moment z Θ3 After the output of the upper controller is compensated, the control command T c is obtained, which is used for subsequent control allocation.

[0106] In this embodiment, according to the extended observer of the position and the attitude, the error of each degree of freedom in the position and the attitude and the observer gain are stacked, and the observation error is expressed as:

[0107]

[0108] wherein, The observation error vectors e1, e2, e3 satisfy The first derivative of the error vectors e1, e2, e3 is represented as; the gain matrices γ1, γ2, γ3 satisfy γ j = diag(γ pj ,γ Θj ), (j = 1, 2, 3);

[0109] For the i-th (i = 1, …, 6) degree of freedom in the error derivative vector , there is:

[0110]

[0111] In the formula, γ 1i ,γ 2i ,γ 3i respectively represent the i-th diagonal element of the gain matrices γ1, γ2, γ3; e bi represents the i-th element of the output error vector e b ; D i represents the i-th element of the disturbance vector D;

[0112] From the above formula, the characteristic polynomial f(λ i ) of this degree of freedom is:

[0113] f(λ i ) = λ i 3 + γ 1i λ i 2 + γ 2i λ i + γ​3i ; (22)

[0114] where λ i represents the eigenvalue of the degree of freedom;

[0115] According to the ideal characteristic polynomial f * (λ i )

[0116] f * (λ i ) = (λ i + ω ob,i ) 3 ; (23)

[0117] If the observer error converges, by comparing equation (22) and (23), we have:

[0118]

[0119] where ω ob,i represents the observer bandwidth of the i-th degree of freedom that satisfies the convergence condition.

[0120] In this embodiment, the oscillation state gain γ p0 , γ Θ0 of the observer are set; the core target of the oscillation state gain mainly has two points; one is to consider the influence of water depth on wave force, accurately evaluate the influence of wave on the unmanned aerial vehicle under water; the second is to ensure that the estimation error of the oscillation state can be stably converged at a faster speed than the dynamic of the wave itself, so as to realize the separation of the influence of the oscillation disturbance; accordingly, the oscillation state gain γ p0 , γ Θ0 of the observer are represented as:

[0121]

[0122] wherein,

[0123] where h w represents the real-time diving depth of the unmanned aerial vehicle; f(h w ) represents the scheduling function of water depth; the elements γ pΘ , γ vω in the matrixes Γ p and Γ Θ are the reference oscillation gain of the observer, γ pΘ,i , γ pΘ,i , γ vω,i , γ vω,i(i = 1, 2, 3) is the reference gain of the oscillation angular velocity;

[0124] The water depth scheduling function f(h) is mentioned above. w A smooth logic function is used, with the value range designed to be [0,1], to determine the real-time diving depth h of the UAV. w Dynamically adjust the observer gain γ p0 γ Θ0 The strength;

[0125] The water depth scheduling function f(h) w The expression for ) is:

[0126]

[0127] In the formula, e is the natural index; λ p The dominant wavelength of the wave can be determined using the relationship between it and the peak frequency. Calculate c0, where c0 is the transition steepness, controlling the width of the transition zone. A smaller c0 value results in a smoother, more gradual transition. c1 is a dimensionless design coefficient, selectable based on the required attenuation sensitivity. Since the depth of wave influence underwater is directly related to its wavelength, when the water depth reaches half a wavelength, the amplitude of underwater particle motion has attenuated to approximately 4% of the surface amplitude. It is generally considered that the wave influence has essentially disappeared, so c1 is typically set to a value of [value missing].

[0128] When the drone is traveling on the water surface or underwater, i.e. h w When ≥0,

[0129] When drones operate on the water surface or in very shallow water, the depth h w It is close to zero. At this point, c0(h) on the exponent term is close to zero. w -c1λ p If ) is a large negative number, then the exponent term is such that It approaches 0. Therefore,

[0130] f(h w The value of γ is approximately equal to 1, at which point the dynamic gain is approximately equal to its reference gain, i.e., γ. p0 ≈Γ p γ Θ0 ≈Γ Θ ;

[0131] And as depth h w The increase of the exponential term continuously increasing, f(h) w The dynamic gain γ0 also decreases as the input signal γ0 decreases. According to the proportional characteristics of a linear system, this is equivalent to the input signal γ0 decreasing. p0 e p γΘ0 e Θ Decay, Also decay. Until the depth h w Increases to infinity, at which time f(h w ) is also infinitely close to 0, the dynamic gain γ p0 , γ Θ0 is approximately equal to 0; thus we know:

[0132] For the optimized extended observer, the and in equations (18) and (19) are expressed as:

[0133]

[0134] According to equations (27) and (28), this means that the UAV is not affected by the oscillatory disturbance, and the extended observer is converted into a linear extended observer to estimate the slowly varying disturbance in water.

[0135] When the UAV is flying above the water surface, i.e., h w <0, at this time f(h w ) = 0, then γ p0 , γ Θ0 are zero matrices; the extended observer is naturally converted into a linear extended observer to estimate the disturbance in air.

[0136] In the setting process of the observer reference oscillation gain, let f(h w ) = 1, then γ p0 = Γ p , γ Θ0 = Γ Θ ; at this time, for the oscillation position and attitude state error e p = z p0 - x0, e Θ = z Θ0 - x Θ0 , the simplified error derivative equation is obtained:

[0137]

[0138] In the equation, and are the first derivatives of the oscillation position state error e p and the attitude state error e Θ , respectively.

[0139] According to equations (12), (13), (14), and (15), for the six-degree-of-freedom error dynamic matrix The i-th (i = 1, …, 6) degree-of-freedom A err,iis expressed as follows:

[0140]

[0141] wherein A err,i is the error dynamic matrix of the ith degree of freedom;

[0142] The characteristic polynomial P(s) = det(sI-A err,i of the error system of the degree of freedom is as follows:

[0143]

[0144] wherein s is the Laplace operator.

[0145] For this, the desired dynamic behavior of the observer error system is defined as a standard stable second-order system, and the desired ideal characteristic equation P di (s) is as follows:

[0146]

[0147] wherein ζ gobs,i is the damping ratio of the ideal second-order system; and ω gobs is the bandwidth of the ideal second-order system.

[0148] By comparing the coefficients of P i (s) in equation (32) and P di (s) in equation (32), it can be obtained that γ pΘ,i , γ vω,i satisfy the following equation:

[0149]

[0150] As can be seen from the above equation, by adjusting the parameters ζ gobs,i and ω gobs,i to obtain the ideal second-order system, the overshoot, oscillation characteristics of the estimation process and the convergence speed of the estimation error can be optimized, so that the error dynamics can be converged to zero most quickly and most stably.

[0151] S4), a nonlinear mapping rule of the roll motion and the pitch angle is designed, and the target pitch angle of the unmanned aerial vehicle in different modes is optimized; specifically including the following steps:

[0152] S41), the tilt angle of the rudder tilt structure in the positive direction of the y-axis is set as the left side pitch angle β l , and the tilt angle of the rudder tilt structure in the negative direction of the y-axis is set as the right side pitch angle β r ; at the same time, it is stipulated that when parallel to the positive direction of the x-axis, the left and right side pitch angles β l = β r= 0; the positive direction of y-axis rotation is positive direction; when the UAV is in air flight mode, the left and right rudder angles are set to 90°, i.e. parallel to the z-axis;

[0153] S42), when the UAV is in water navigation, it is divided into two modes as follows, as shown in Figure 3 、 4

[0154] Left and right angle symmetry mode: applied to the UAV in motion state, at this time the angle satisfies β l + β r = 0;

[0155] Left and right angle asymmetry mode: applied to the UAV in water hovering state, at this time the angle should satisfy β l ≠ β r , and only when β l = β r = 0, the angles are equal; in this state, the angle of one side of the UAV rotates, and the angle of the other side changes according to the actual situation;

[0156] S42), the nonlinear mapping rule of roll motion and angle is designed according to the roll angle and roll angle velocity, i.e.

[0157] β d = β max sat(h(e φ ,p),1,-1); (34)

[0158]

[0159] In the formula, β d represents the angle designed according to the roll motion; β max represents the maximum value of the preset angle; h(·) is the core part of the entire nonlinear mapping rule design; tanh(·) is the hyperbolic tangent function; sat(·) is the saturation function; ω nl , ω l are the weights of the nonlinear part and the linear part respectively; k nl , k l are the gain coefficients of the nonlinear and linear inputs respectively; f(e φ ,p) represents the linear combination term combining the roll angle error e φ and the current roll angle velocity p; φ is the current roll angle; φ d is the roll angle expected by the upper controller; k φ is the proportional gain of the roll angle error, used to adjust the response strength of the upper controller to static or slowly changing attitude deviation; k p ​is the derivative gain of roll angle velocity, which is used to provide system damping, suppress the rapid oscillation of attitude and ensure the smoothness of dynamic process.

[0160] When input f(e φ ,p) is small, formula (35) is mainly dominated by the nonlinear term tanh(k nl f(e φ ,p)); this means that when the error is small, the system will correct the deviation with the greatest force, trying to make the error quickly return to zero. When the input f(e φ ,p) is large, formula (35) is mainly dominated by the nonlinear term k l f(e φ ,p); in the face of huge errors or shocks, the response of the upper controller is no longer unlimitedly increased, but grows with a more moderate and predictable linear slope. This can effectively prevent control from being oversaturated, system oscillation and severe impact on the actuator.

[0161] S43), when the UAV is in the left-right angle symmetry mode, the target angle β target is:

[0162]

[0163] When the UAV is in the left-right angle asymmetry mode, the target angle β target is:

[0164]

[0165] In the formula, β ltar ,β rtar represent the target angle of the left angle β l and the right angle β r , respectively.

[0166] S5), a multi-objective cost function is constructed, and an improved genetic algorithm is used to solve the expected thrust and expected angle, and the output thrust and blade rotation speed of each motor are calculated; specifically including the following steps:

[0167] S51), when the UAV moves underwater, the motor of the UAV generates torque in the x, y, z axes of the body coordinate system and thrust in the x, z axes of the body coordinate system, and its expression is:

[0168] A(β)u f = T c ; (38)

[0169] In the formula, T c = [τ x ,τ y ,τ z , Tx ,T z ] T is the total thrust control command for the upper controller output after disturbance compensation; u f = [f1, …, f8] T is the thrust generated by each motor; A(β) is the control allocation matrix;

[0170]

[0171] where l is the wheelbase of the motor from the center of the UAV; k is the anti-torque coefficient; β l is the left roll angle; β r is the right roll angle;

[0172] S52), in order to further decouple the roll angle and thrust bi-variables when solving the input u f of the actuator, a multi-objective cost function is designed according to the control allocation relationship and the target roll angle β target , that is:

[0173]

[0174]

[0175] where J(β) is the multi-objective cost function; β represents the rudder roll angle vector of the current optimization solution; w uj represents the weight coefficient of the thrust of a single motor in the control energy consumption penalty term; u fj represents the output thrust of the jth motor; Δβ tra represents the tracking error of the roll angle; W β represents the weight matrix of the roll angle tracking penalty term; Δu var represents the rate of change of the motor thrust; W △u represents the weight matrix of the motor thrust rate penalty term; Δu var represents the rate of change of the motor thrust; Δβ var represents the rate of change of the roll angle; W △β represents the weight matrix of the roll angle rate penalty term; the superscript T represents the transpose operation of the vector; u pre represents the motor thrust vector of the previous control period; β pre represents the rudder roll angle vector of the previous control period; β max , β min respectively represent the physical upper and lower limits of the rudder roll angle.

[0176] is the control energy consumption penalty term; Δβ tra T W β Δβ traFor the tilt angle tracking penalty term, by minimizing this term, both the main control task and the solution that makes the tilt angle as close as possible to the smooth, stable and adaptive target value of the flight state are completed; Δy var T W △u Δu var and Δβ var T W △β Δβ var For the actuator rate penalty term, this term ensures the smoothness of the control signal, and the penalty on the motor thrust rate can reduce current impact and mechanical vibration, and the penalty on the tilt angle rate directly limits the movement speed of the rudder, prevents it from swinging at high speed, thereby reducing mechanical wear and avoiding exciting the high-frequency dynamics of the system;

[0177] S53), solve the multi-objective cost function J(β) using the improved genetic algorithm, and decompose the thrust and tilt angle synchronization optimization problem into a low-dimensional nonlinear optimization problem centered on the tilt angle, as shown in Figure 5 , which specifically includes the following steps:

[0178] S531), gene coding and population initialization;

[0179] The core task of the genetic algorithm is to find the optimal left and right rudder tilt angle, therefore, each chromosome is defined as a vector β i = [β l β r ] T ; the coding method uses real number coding, which directly corresponds to the physical angle of the rudder, which is more intuitive and more efficient in calculation.

[0180] The population is initialized using a hybrid initialization strategy of random seeding and heuristic seeding, wherein random seeding accounts for 70% of the population size and is generated by uniformly random sampling within the physically feasible range of the rudder tilt angle [β min , β max ]; heuristic seeding accounts for 30% of the population size and is generated by the target tilt angle β target Designing the generation method as:

[0181]

[0182] In the formula, is a normally distributed random number with a mean of 0 and a standard deviation of σ init ;

[0183] In addition to being seeded near β target in a Gaussian distribution, the optimal solution of the previous control cycle is also implanted as an individual.

[0184] S532), fitness function design and solution;

[0185] The fitness function is the core standard to evaluate the merits of each individual, which can be directly derived from the cost function designed above. Since genetic algorithm usually solves the maximum problem, the fitness function is as follows:

[0186]

[0187] In the formula, β ind is the inclination vector decoded from the indth chromosome, J(β ind ) is the total value of the cost function calculated; ε is a very small positive number, which is used to prevent division by zero error when J(β ind ) tends to zero;

[0188] For each individual β ind in the population, the calculation process of its fitness is as follows:

[0189] The control distribution matrix A(β ind )u f =T c , A(β ind )u f =T c is solved by using the pseudo-inverse method, and the minimum two-norm solution is obtained, that is:

[0190]

[0191] In the formula, u f,ibd indicates the thrust vector corresponding to the indth inclination, which meets the minimum two-norm solution; A + (β ind ) indicates the pseudo-inverse matrix corresponding to the indth inclination;

[0192] Then, u f,ind , β ind , and the multi-objective cost function in formula (40) are substituted into the fitness function in formula (43) to solve the fitness of the indth chromosome;

[0193] S533), for selection, in each generation, a number of individuals are randomly selected from the current population for comparison by using the tournament selection combined with elitism, and the individual with the highest fitness is selected to enter the mating pool of the next generation;

[0194] S534), for crossover, the adaptive simulated binary crossover operator is used to select two parent individuals P1 and P2 from the parent population to generate two child individuals C1 and C2; for each gene β ind of the chromosome β l and β rThe generation of offspring follows the following simulated binary crossover operator rules:

[0195]

[0196] where ρ c is the expansion factor; μ is a random number with value in (0, 1); ζ is the crossover expansion exponent; ζ is the core parameter that controls the difference between offspring and parents; the larger the value of ζ, the more closely the value of ρ c will be clustered around 1, resulting in offspring that are very similar to parents and have a narrower distribution range; otherwise, the value of ρ c has a greater probability of deviating from 1, resulting in offspring that are very different from parents and have a wider distribution range;

[0197] The function ζ(b) is designed to make ζ change smoothly with the running of the algorithm, which satisfies the following expression:

[0198]

[0199] where b is the current evolution generation; G max is the preset maximum total evolution generation; ζ max and ζ min are the upper and lower limit values of the distribution exponent;

[0200] From the above expression, it can be seen that in the early evolution stage, b is small, so ζ is close to its lower limit ζ min , and a smaller ζ value will produce an expansion factor ρ that has a wider distribution. This means that the generated offspring individuals have a high probability of falling outside their two parent individuals, greatly encouraging the algorithm to explore the entire solution space extensively and in a large range, thereby effectively avoiding the algorithm from prematurely converging to a local optimal trap.

[0201] In the later evolution stage, b is large, so ζ is close to its lower limit ζ max , and a larger ζ value will produce an expansion factor ρ that is very close to 1. This means that the generated offspring individuals will be very closely distributed between the two parent individuals, and at this time the algorithm has found the region containing the optimal solution through the early exploration, and then the algorithm can perform fine and small-range search and fine-tuning around the already found high-quality solution, thereby achieving accurate convergence to the optimal solution.

[0202] For mutation, a small number of individuals are selected for mutation, and a target-guided adaptive mutation operator is used. When a parent individual P F is selected for mutation, the new individual Cp generated is no longer a simple random disturbance of P F , but a result of weighted combination of three vectors:

[0203]

[0204] where P F is the parent individual; V g is the guidance vector, which is a vector pointing directly to P g ; P g is the current best individual; β best is the target inclination angle; w target and w b are the setting weights of P t and β g , respectively, so that the guidance vector tends to P best or β target ; V e is the exploration vector, which is used to find a better solution that may exist in the vicinity of the parent individual P F but has not been evaluated; P M (P, η m ) is the standard polynomial mutation, which is random in direction and generates a new individual around the parent individual P F , with the parameter η m being the distribution index of the mutation, which is used to control the strength of the polynomial mutation; w g and w e are the adaptive weight matrices, which will be updated according to the fitness standard deviation σ fit of the current population;

[0205] where:

[0206]

[0207] where e is the natural constant; d σ is the normalized population diversity index calculated according to σ fit ; σ g is the standard deviation of the Gaussian function; d g controls the speed of the weight transition; B is the peak value of the Gaussian function; d g is the center point of the Gaussian function, which should be set to be smaller than d low ; k is the steepness coefficient; d m is the set value of the low population diversity; g max (·) is the diversity scale function; clip(·) is the clipping function; σ min and σ fit are the reference upper and lower limits set for the fitness standard deviation σ fit , respectively;

[0208] In the early stage of evolution, the population is randomly generated, so the population diversity is rich, and σ g and d are large, so w eAll are small, the overall impact of variation is small, the algorithm mainly relies on crossover operator to carry out large-scale population recombination. With the gradual reduction of population diversity, σ fit , d gradually reduces, w g , w e increases, the guidance and exploration are enhanced; but because d low < d g , w g will increase sharply and reach the peak, become the dominant force in this stage. Until d < d low , w g starts to decrease, the guidance is gradually weakened; while w e is still increasing, the exploration is gradually enhanced. With the further reduction of σ fit , d in the late evolution, w g tends to 0, w e tends to 1, the exploration becomes dominant, weakening the convergence trend, which is conducive to prevent the algorithm from falling into local optimal solution.

[0209] S536), constantly selection, crossover, mutation and fitness calculation, the evolution process will be terminated when any of the following conditions is met, and then return the individual with the highest fitness according to the fitness function β(ind), i.e. the optimal tilt angle β out and the corresponding motor expected thrust solution u fout found so far; the termination conditions are:

[0210] a. The evolution generation of the algorithm reaches the preset upper limit; b. The fitness of the optimal individual in the population does not appear meaningful improvement within a specified number of consecutive generations; c. The execution time of the algorithm exceeds the main loop frequency of the flight control computer.

[0211] In this embodiment, after the optimal tilt angle β out and the corresponding motor expected thrust solution u fout , the following steps are taken: 1) if the thrust solution u fout derived by the genetic algorithm satisfies the physical limitations of the motor, i.e. it satisfies the minimum and maximum thrust limitations that the motor can generate, u min ≤ u fout ≤ u max ; then let u out = u fout ; 2) if the thrust solution derived by the genetic algorithm does not satisfy the physical limitations of the motor, then fix the optimal tilt angle β out , convert the control allocation matrix A(β) into a constant matrix A out = A(β out ); then use the null space pseudo-inverse solution method or the adaptive weighted pseudo-inverse method to solve u foutsolving to obtain a set of solutions that meet the physical limits, or constructing a new cost function s.t.A out u=T c , and the QP is solved to obtain the desired thrust vector u out ; 3) the desired thrust vector u out is obtained, and the required rotational speed of each propeller can be calculated by the following equation:

[0212]

[0213] where c is the power coefficient; u out,i is the desired thrust vector u out ; ω i is the rotational speed of the i-th propeller.

[0214] The above embodiments and descriptions are only illustrative of the principles and the best mode of the present application, and various changes and modifications can be made without departing from the spirit and scope of the present application, and such changes and modifications are intended to fall within the scope of the present application.

Claims

1. A composite control method for amphibious unmanned aerial vehicles based on motion state separation and multi-target allocation, characterized in that, Includes the following steps: S1) Construct a motion mathematical model for the UAV; S2) Construct a disturbance model to simulate the effect of irregular wave force oscillation disturbance on the motion state of UAVs near the water surface, resulting in oscillation displacement and oscillation Euler angle; S3) Design a position and attitude extension observer based on state separation; S4) Design the nonlinear mapping law between roll motion and tilt angle, and optimize the target tilt angle for different modes of UAV. S5) Construct a multi-objective cost function, and use an improved genetic algorithm to solve for the desired thrust and desired tilt angle, and calculate the output thrust and blade rotational angular velocity of each motor.

2. The amphibious unmanned aerial vehicle (UAV) composite control method based on motion state separation and multi-target allocation according to claim 1, characterized in that: In step S2), the oscillation disturbance model experienced by the UAV is as follows: Y = CH; (10) The above equation satisfies: In the formula, Let H be the first derivative of the oscillating motion state vector H; A is the system state matrix, and it is a full-rank matrix; diag(·) denotes a diagonal matrix with its elements as diagonals; the state matrices A1,…,A6 of each degree of freedom satisfy the following conditions: ξ i The damping ratio associated with the i-th degree of freedom is used to match the sharpness of the peaks in the wave spectrum; ω p K represents the peak energy frequency of the wave spectrum in the target water area; K is the system input matrix; the input matrices K1,…,K6 for each degree of freedom satisfy… k i Let C be the input gain of the i-th degree of freedom; C is the system output matrix; the output matrices C1,…,C6 of each degree of freedom satisfy C i = [1 0], (i = 1, ..., 6); y is the unit white noise input; Y is the complete oscillating motion state output vector.

3. The amphibious unmanned aerial vehicle (UAV) composite control method based on motion state separation and multi-target allocation according to claim 2, characterized in that: In step S3), the position extension observer based on state separation is represented as: In the formula, e p It is a three-dimensional error vector between the estimated and true location values; It is a three-dimensional estimated vector of the UAV's position and state; y p Indicates the location status of the drone; z p0 ,z p1 ,z p2 The position and state vectors x p0 ,x p1 ,x p2 The estimated values ​​are the estimates of the oscillating displacement, three-dimensional position, and three-dimensional velocity of the UAV. z p0 ,z p1 ,z p2 The first derivative; z p3 It is an estimate of all positional disturbances experienced by the UAV; γ p0 ,γ p1 ,γ p2 ,γ p3 All are three-dimensional diagonal gain matrices of the position observer; A p =diag(A1,A2,A3); A1,A2,A3 are the state matrices of degrees of freedom; ζ i The damping ratio associated with the i-th degree of freedom is used to match the sharpness of the peaks in the wave spectrum; ω p represents the peak energy frequency of the wave spectrum of the target water area; diag(·) denotes a diagonal matrix with its elements as the diagonal; For H p The first derivative; η is the oscillating displacement vector; i Let ξ be the i-th degree of freedom in the oscillating motion state; R represents the rotation matrix from the body to the spatial coordinate system; F represents the thrust of the UAV in the body coordinate system; C p = [C1C2C3] is the output matrix for each degree of freedom, satisfying C i =[1 0].

4. The amphibious unmanned aerial vehicle (UAV) composite control method based on motion state separation and multi-target allocation according to claim 3, characterized in that: In step S3), the attitude extension observer based on state separation is represented as: In the formula, e Θ It is the three-dimensional error vector between the estimated and true attitude values; It is a three-dimensional estimated vector of the UAV's attitude state; y Θ Indicates the attitude state of the drone; z Θ0 ,z Θ1 ,z Θ2 These are the attitude state vectors x and x. Θ0 ,x Θ1 ,x Θ2 The estimated values ​​are the estimated values ​​of the oscillating Euler angles, three-dimensional Euler angles, and three-dimensional angular velocity of the UAV. z Θ0 ,z Θ1 ,z Θ2 The first derivative; z Θ3 It is an estimate of all attitude disturbance moments experienced by the UAV; γ Θ0 ,γ Θ1 ,γ Θ2 ,γ Θ3 All are three-dimensional diagonal gain matrices of the attitude observer; A Θ =diag(A4,A5,A6); For H Θ The first derivative; C is the oscillating Euler angle vector; Θ = [C4C5C6]; J is the inertia matrix; τ is the torque generated by the propeller on the fuselage shaft.

5. The amphibious unmanned aerial vehicle (UAV) composite control method based on motion state separation and multi-target allocation according to claim 4, characterized in that: In step S4), the tilt angle of the servo tilting structure in the positive y-axis direction is set as the left tilt angle β. l The tilt angle of the servo tilting structure in the negative y-axis direction is the right tilt angle β. r When the drone is navigating in water, it is divided into the following two modes: Left and right tilt symmetry mode: Applied when the drone is in motion, in which case the left and right tilt angles satisfy β. l +β r =0; Asymmetric tilt mode: Applied when the drone is hovering in water, the tilt angle should satisfy β. l ≠β r There is one and only β l =β r When the tilt angle is 0, the tilt angles are equal; in this state, the tilt angle on one side of the drone rotates, while the tilt angle on the other side changes according to the actual situation.

6. The amphibious unmanned aerial vehicle (UAV) composite control method based on motion state separation and multi-target allocation according to claim 5, characterized in that: In step S4), the nonlinear mapping law between roll motion and tilt angle is designed based on the roll angle and roll angular velocity, namely: b d =b max sat(h(e φ ,p),1,-1);(34) In the formula, β d Indicates the tilt angle designed based on roll motion; β max The maximum value of the preset tilt angle is represented by h(·), which is the core part of the entire nonlinear mapping law design; tanh(·) represents the hyperbolic tangent function; sat(·) is the saturation function; ω nl ω l These are the weights of the nonlinear and linear components, respectively; k nl k l These are the gain coefficients for nonlinear and linear inputs, respectively; f(e φ ,p) represents the combined roll angle error e φ A linear combination of the current roll rate p; φ is the current roll angle; φ d k is the desired roll angle of the upper-level controller. φ This is the proportional gain for the roll angle error, used to adjust the response strength of the upper-level controller to static or gradually changing attitude deviations; k p The differential gain of the roll angular velocity is used to provide system damping, suppress rapid attitude oscillations, and ensure the smoothness of the dynamic process.

7. The amphibious unmanned aerial vehicle (UAV) composite control method based on motion state separation and multi-target allocation according to claim 6, characterized in that: In step S4), when the UAV is in a symmetrical tilt mode, the target tilt angle β is... target for: When the drone is in an asymmetrical tilt mode, the target tilt angle β is as follows. target for: In the formula, β ltar ,β rtar These represent the left tilt angle β. l and right side tilt angle β r The target inclination angle.

8. The amphibious unmanned aerial vehicle (UAV) composite control method based on motion state separation and multi-target allocation according to claim 7, characterized in that: In step S5), based on the control allocation relationship and the target tilt angle β... target The designed multi-objective cost function is as follows: In the formula, J(β) is the multi-objective cost function; β represents the servo tilt angle vector of the current optimization solution; w uj This represents the weighting coefficient of the thrust of a single motor in the energy consumption penalty term; u fj Δβ represents the output thrust of the j-th motor. tra The tracking error represents the tilt angle; W β The weight matrix representing the tilt angle tracking penalty term; Δu var W represents the rate of change of the motor thrust. △u The weight matrix representing the penalty term for the rate of change of motor thrust; Δu var Δβ represents the rate of change of the motor thrust. var W represents the rate of change of the tilt angle. △β The weight matrix represents the penalty term for the rate of change of tilt angle; the superscript T indicates the transpose of the vector; u pre β represents the motor thrust vector of the previous control cycle. pre β represents the servo tilt angle vector of the previous control cycle. max β min These represent the upper and lower physical limits of the servo tilt angle, respectively.

9. The amphibious unmanned aerial vehicle (UAV) composite control method based on motion state separation and multi-target allocation according to claim 8, characterized in that: In step S5), an improved genetic algorithm is used to solve the multi-objective cost function J(β), decomposing the thrust and tilt angle synchronization optimization problem into a low-dimensional nonlinear optimization problem with tilt angle as the core; specifically, it includes the following steps: S531), gene coding and population initialization; Each chromosome is defined as a vector β containing two genes. i =[β l ,β r ] T The encoding method uses real number encoding. A hybrid initialization strategy of random seeding and heuristic seeding was employed to initialize the population, with random seeding accounting for 70% of the population size. This initialization was carried out within the physically feasible region [β] of the servo tilt angle. min ,β max Generate the population through uniform random sampling within the target area; heuristic seeding accounts for 30% of the population size, using the target tilt angle β. target The design generation method is as follows: In the formula, With a mean of 0 and a standard deviation of σ init Normally distributed random numbers; Besides in β target The seeds are sown in a Gaussian distribution in the vicinity, and the optimal solution of the previous control cycle is also implanted as an individual. S532), Fitness Function Design and Solution; The fitness function is as follows: In the formula, β ind The tilt vector decoded from chromosome ind, J(β) ind ) represents the total cost value calculated using the cost function; ε is a very small positive number used to prevent the total cost value from being calculated when J(β) is less than or equal to 0. ind A division-by-zero error occurs when the value approaches zero. For each individual β in the population ind The fitness calculation process is as follows: The control allocation matrix A(β) is solved using the pseudo-inverse method. ind )u f =T c We obtain the least-norm 2 solution, that is: In the formula, u f,ind A represents the thrust vector corresponding to the least-norm solution at the ind-th tilt angle; + (β ind ) represents the pseudo-inverse matrix corresponding to the ind-th tilt angle; Then u f,ind β ind The multi-objective cost function in equation (40) is substituted into the fitness function in equation (43) to solve for the fitness of the ind-th chromosome; (S533) For selection, in each generation, a strategy combining tournament selection and elitism is used to randomly select several individuals from the current population for comparison, and the individual with the highest fitness is selected to enter the mating pool of the next generation. (S534) For crossover, an adaptive simulated binary crossover operator is used to select two individuals P1 and P2 from the parent population to generate two offspring individuals C1 and C2; for chromosome β... ind Each gene β l and β r The generation of offspring follows the rules of the simulated binary crossover operator: In the formula, ρ c ζ is the expansion factor; μ is a random number with a value in the range (0,1); ζ is the cross-expansion index; Design a function ζ(b) such that ζ changes smoothly as the algorithm runs, satisfying the following expression: In the formula, b is the current generation number; G max ζ is the preset maximum total number of generations; max and ζ min These are the upper and lower limits of the distribution index; (S535) For mutation, a small number of individuals are selected for mutation, and a target-guided adaptive mutation operator is used. When a parent individual P... F When selected for mutation, the resulting new individual C p No longer P F It is not a simple random perturbation, but the result of a weighted synthesis of three vectors: In the formula, P F For the parent generation; V g It is a guiding vector, a vector that explicitly points to P. g The vector of P; g For the current optimal individual P best and target inclination angle β target Dynamically blended guiding objectives; w b and w t For P g The weighting of the vector makes the guiding vector more inclined towards P. best Or B target V e It is an exploration vector used to discover the parent individual P. F Possible but unverified better solutions in the vicinity; P M (P,η m This is a standard polynomial mutation, with random direction, producing a parent individual P. F The new individuals around it, whose parameter η m w is the distribution index of the variation, used to control the strength of the polynomial variation. g and w e The adaptive weight matrix will be based on the fitness standard deviation σ of the current population. fit Update; in: In the formula, e is the natural constant; d σ According to σ fit The calculated normalized population diversity index; σ g d is the standard deviation of the Gaussian function; controls the rate of weight transition; B is the peak value of the Gaussian function; d g The center point of the Gaussian function is d; k is the kurtosis coefficient; d low This is a set value for population diversity that is too low; this value should be less than d. g g m (·) is the diversity scaling function; clip(·) is the cutoff function; σ max and σ min These are the fitness standard deviations σ and σ, respectively. fit The set upper and lower reference limits; (S536) The evolutionary process continues with selection, crossover, mutation, and fitness calculations. It will terminate when any of the following conditions are met, and then return the individual with the highest fitness obtained according to the fitness function F(ind), which is the currently found optimal tilt angle β. out and the corresponding expected thrust solution of the motor u fout ; The termination condition is as follows: a) The algorithm reaches the preset upper limit of the number of generations; b) The fitness of the best individual in the population does not show meaningful improvement within a specified number of generations; c) The execution time of the algorithm exceeds the frequency of the main loop of the flight control computer.

10. The amphibious unmanned aerial vehicle (UAV) composite control method based on motion state separation and multi-target allocation according to claim 9, characterized in that: In step S5), after obtaining the optimal tilt angle β... out and the corresponding expected thrust solution of the motor u fout Then, proceed as follows:

1. If the thrust solution u obtained by the genetic algorithm fout The physical limitations of the motor must be met, that is, the minimum and maximum thrust limits that the motor can generate must be met. min ≤u fout ≤u max Then let u out =u fout ; 2. If the thrust solution obtained by the genetic algorithm does not meet the physical constraints of the motor, then fix the optimal tilt angle β. out Transform the control allocation matrix A(β) into a constant matrix A out =A(β) out Then, the null space pseudo-inverse method or the adaptive weighted pseudo-inverse method is used to solve u. fout Solving for a set of solutions that satisfy the physical limits, or constructing a new cost function. stA out u = T c The thrust is optimized using QP (Quantization-Physics) to obtain the desired thrust u of the motor. out ; 3. Obtain the desired thrust vector u of the motor. out Then, the required blade speed for each motor is obtained using the following formula. In the formula, c is the dynamic coefficient; u out,i For the desired thrust u out The desired output thrust of the i-th motor; ω i Let be the rotational angular velocity of the i-th blade.

Images