Quadrotor unmanned aerial vehicle multi-physical domain wind disturbance control method and system based on digital twinning

Abstract

The present application relates to the technical field of unmanned aerial vehicle flight control, in particular to a quadrotor unmanned aerial vehicle multi-physical domain wind disturbance control method and system based on digital twinning, comprising the following steps: based on the corresponding multi-physical domain digital twin of the quadrotor unmanned aerial vehicle, synchronously solving the wind field aerodynamic force, the body motion and the motor-to-rotor thrust establishment process according to the calibrated wind field parameters and the current feedback state of the unmanned aerial vehicle; determining the response sequence relationship between the wind load change and the rotor thrust establishment, obtaining the wind load execution response state for representing the coordination degree of wind disturbance effect and actual execution mechanism response. The present application synchronously solves the wind field aerodynamic force, the body motion and the motor-to-rotor thrust establishment process, identifies the sequence relationship between the wind load response time and the execution response time, and combines the wind load change degree, the rotor thrust following degree and the position and attitude deviation to form the wind load execution response state, providing a cross-physical domain dynamic basis for wind resistance control.

Description

Technical Field

[0001] This invention relates to the field of unmanned aerial vehicle (UAV) flight control technology, and in particular to a multi-physics domain wind disturbance control method and system for quadrotor UAVs based on digital twins. Background Technology

[0002] Quadcopter drones are characterized by their simple structure, high maneuverability, and strong vertical takeoff and landing capabilities, and are widely used in inspection, surveying, and low-altitude operations. In actual flight environments, wind fields typically exhibit characteristics such as average incoming flow, random turbulence, and sudden wind disturbances. Especially in areas with buildings and wind turbine wakes, wind speed fluctuations can cause drone positional shifts and attitude disturbances, placing high demands on trajectory tracking stability.

[0003] Existing quadrotor control methods mostly employ PID or fuzzy PID to adjust the control input based on position error, attitude error, and error rate of change. However, they typically focus on the airframe's motion response and fail to fully consider the dynamic differences across multiple physical domains, including wind farm aerodynamic changes, motor electrical dynamics, rotor speed establishment, and thrust generation. Especially when wind load has already acted on the airframe but rotor thrust has not yet been established, if the controller still only adjusts based on instantaneous control error enhancement parameters, integral accumulation, response lag, and subsequent overcompensation problems are likely to occur. Furthermore, some simulation models directly use ideal speed sources, making it difficult to realistically reflect the dynamic process of the control signal passing through electrical drive, motor torque, and finally rotor thrust.

[0004] Therefore, there is an urgent need to design a quadcopter UAV control method that can combine multi-physical domain digital twin state recognition to identify the relationship between wind load and actuator response sequence, and coordinate wind-resistant control parameters accordingly. Summary of the Invention

[0005] This invention provides a multi-physics domain wind disturbance control method and system for quadrotor unmanned aerial vehicles based on digital twins.

[0006] On the one hand, a multi-physics domain wind disturbance resistance control method for quadrotor UAVs based on digital twins includes the following steps:

[0007] Based on the multi-physics domain digital twin corresponding to the quadcopter drone, the wind field aerodynamics, body motion, and the thrust establishment process from the motor to the rotor are solved simultaneously according to the calibrated wind field parameters and the current feedback state of the drone.

[0008] Determine the sequential relationship between wind load changes and rotor thrust establishment to obtain the wind load execution response state, which characterizes the degree of coordination between wind disturbance and the actual response of the actuator.

[0009] Based on the control error and error change rate corresponding to the wind load execution response state and the target trajectory, the adaptive fuzzy PID controller is input, and the correction intensity of the proportional, integral and derivative parameters is adjusted according to the response sequence relationship between wind load change and rotor thrust establishment.

[0010] When wind load changes precede rotor thrust establishment, the rapid suppression and damping effects are enhanced and the integral accumulation is limited. When rotor thrust completes to follow and control error enters a convergent state, steady-state correction is restored.

[0011] Based on the smooth update of the adjusted control parameters, the phase-coordinated wind-resistant control quantity is obtained.

[0012] Based on the phase coordination wind resistance control quantity, it is converted into independent control signals for the four rotors according to the rotor layout of the quadcopter UAV, driving the corresponding actuators to complete position and attitude correction.

[0013] The corrected UAV position, attitude, and rotor speed are obtained as closed-loop feedback states, and the closed-loop feedback states are sent back to the multi-physics domain digital twin to update the wind load execution response state for the next control cycle.

[0014] On the other hand, the multi-physics domain wind disturbance resistance control system for quadrotor UAVs based on digital twins includes the following modules:

[0015] A multi-physical domain digital twin module is used to simultaneously solve the wind field aerodynamics, airframe motion, and rotor thrust establishment process based on the calibrated wind field parameters and the current feedback state of the UAV, so as to obtain the wind load execution response state.

[0016] An adaptive fuzzy PID control module is used to adjust the proportional parameter, integral parameter and derivative parameter according to the wind load execution response state, control error and error change rate to obtain the phase-coordinated wind-resistant control quantity.

[0017] A control allocation module is used to convert the phase-coordinated wind-resistant control quantity into independent control signals for the four rotors according to the X-shaped rotor layout of the quadcopter UAV.

[0018] An actuator module is used to adjust the thrust of each rotor according to the independent control signal through an electrical drive, motor torque establishment and rotor speed establishment process, so as to complete the position and attitude correction of the UAV;

[0019] The state perception and feedback module is used to acquire the corrected UAV position, attitude and rotor speed and form a closed-loop feedback state. The closed-loop feedback state is sent back to the multi-physical domain digital twin module to update the wind load execution response state for the next control cycle.

[0020] The beneficial effects of this invention are:

[0021] This invention, based on multi-physical domain digital twin wind load execution response state discrimination, identifies the sequential relationship between wind load response time and execution response time by simultaneously solving the wind field aerodynamic force, airframe motion, and the thrust establishment process from motor to rotor. It also combines the degree of wind load change, rotor thrust following degree, and position and attitude deviation to form the wind load execution response state, providing a cross-physical domain dynamic basis for wind resistance control.

[0022] This invention, based on adaptive fuzzy PID coordinated control of response phase relationship, introduces wind load execution response state on the basis of traditional error and error change rate parameter tuning. When the wind load is in advance, the proportional and derivative corrections are enhanced and the integral accumulation is suppressed. After the rotor thrust completes the follow-up, the steady-state parameters are gradually restored. Then, the quadcopter is driven by smooth update and X-type control distribution, so as to achieve coordination between the disturbance rejection response and the actual dynamic of the actuator. Attached Figure Description

[0023] To more clearly illustrate the technical solutions in this invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only for this invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0024] Figure 1 This is a schematic diagram of the method flow of Embodiment 1 of the present invention;

[0025] Figure 2 This is an architectural diagram of the multi-physical domain digital twin of Embodiment 1 of the present invention;

[0026] Figure 3 This is a diagram of the adaptive fuzzy PID structure of Embodiment 1 of the present invention;

[0027] Figure 4 This is a system module diagram of Embodiment 2 of the present invention. Detailed Implementation

[0028] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. For some well-known technologies, those skilled in the art may also use other alternative methods to implement the invention. Moreover, the accompanying drawings are only for more specific description of the embodiments and are not intended to specifically limit the present invention.

[0029] Example 1

[0030] like Figures 1-3 As shown, the multi-physics domain wind disturbance resistance control method for quadrotor UAVs based on digital twins includes the following steps:

[0031] S1 runs a multi-physics domain digital twin corresponding to the quadcopter UAV. Based on the calibrated wind field parameters and the current feedback state of the UAV, it synchronously solves the wind field aerodynamics, airframe motion, and the thrust establishment process from the motor to the rotor. It determines the response sequence between wind load changes and rotor thrust establishment, and obtains the wind load execution response state used to characterize the degree of coordination between wind disturbance and the actual response of the actuator.

[0032] S11, the construction method of the multi-physics domain digital twin is as follows: Based on the Modelica unified modeling environment, an environment module, a body dynamics module, an actuator module, and a state perception and feedback module are established respectively, and coupled according to the actual energy transfer and state feedback relationship of the quadcopter UAV; wherein, the environment module generates average wind, random turbulence, and corresponding wind field aerodynamic forces and aerodynamic moments based on calibrated wind field parameters, and loads them to the body dynamics module through a multi-body physics interface; the body dynamics module solves for the position and attitude changes of the UAV based on the UAV's mass, inertial parameters, rotor installation position, wind field aerodynamic forces, and the thrust of each rotor; the actuator... The structural module establishes models of the battery, inverter, permanent magnet DC motor, speed / current dual closed loop, and rotor according to the transmission relationship of control signals, electrical drive, motor current, electromagnetic torque, rotor speed, and rotor thrust. The four rotor thrusts are applied to the corresponding installation positions of the airframe through independent physical interfaces. The state perception and feedback module is connected to the airframe dynamics module and the actuator module to obtain the UAV's position, attitude, and rotor speed, and feeds them back to the control module and the digital twin solution process of the next control cycle, thus forming a multi-physics domain digital twin with continuous coupling of environmental load, airframe motion, actuator response, state perception, and closed-loop feedback.

[0033] The environment module is used to generate the wind field input and corresponding wind field aerodynamic forces acting on the quadcopter UAV; the airframe dynamics module is used to solve the position and attitude changes of the UAV under the combined action of wind field aerodynamic forces, rotor thrust and gravity; the actuator module is used to establish the dynamic transmission process of control signals through electrical drive, motor torque, rotor speed and finally rotor thrust; the state perception and feedback module is used to obtain the position, attitude and rotor speed of the UAV and use them as the state feedback quantity for the next control cycle.

[0034] First, the calibration wind field parameters are determined; for the first... A hovering wind field calibration data segment, assuming that it has There are [number] valid sampling points, and the horizontal wind speed at each sampling point is [value]. The average wind speed of this hovering data segment is then expressed as:

[0035] ;

[0036] The standard deviation of wind speed in this hovering data segment is:

[0037] ;

[0038] The corresponding equivalent cable flow intensity is:

[0039] ;

[0040] in, For the first The number of effective wind speed samples per hovering data segment; For the first The first hover data segment Horizontal wind speed at each sampling point; This represents the average wind speed corresponding to the hovering data segment; This represents the standard deviation of wind speed for the corresponding hovering data segment. This represents the equivalent cable flow intensity corresponding to the hovering data segment.

[0041] When wind field data at different hovering heights need to be uniformly converted to a reference height, the reference wind speed is converted according to the logarithmic mean wind profile and expressed as follows:

[0042] ;

[0043] ;

[0044] in, For the first The equivalent reference wind speed after conversion of each hovering data segment; The preset reference height is 20ft, which is approximately 6.096m. For the first The average hovering height of each hovering data segment; The minimum effective height is set to avoid numerical anomalies in logarithmic operations at low altitudes; This represents the length of the surface roughness.

[0045] The surface type is determined based on the actual mission area of ​​the quadcopter UAV. The surface roughness length is then determined by combining existing meteorological data, on-site wind measurement data, or multi-altitude hovering wind measurement results for that area. Specifically, under the same wind field conditions, multiple stable hovering segments at different altitudes are selected to obtain the corresponding average wind speed at each hovering altitude. The altitude-average wind speed data is fitted using a logarithmic wind profile, and the roughness parameter corresponding to the smallest fitting error is taken as the surface roughness length. When on-site multi-altitude wind measurement data is insufficient, the corresponding engineering experience roughness length is selected based on the surface cover type of the mission area, and corrected using measured wind speeds during subsequent wind field calibration. After obtaining the surface roughness length, the minimum effective altitude is determined based on the actual minimum allowable flight altitude of the UAV and the numerical stability requirements of the logarithmic wind profile. This minimum effective altitude is greater than the surface roughness length, and digital twin simulation verification is performed starting from the actual minimum allowable flight altitude of the UAV. The candidate altitudes are gradually reduced, and when the average wind speed calculated from the logarithmic wind profile does not show abnormal sudden increases, numerical divergence, or unreasonable jumps in wind speed at adjacent altitudes, the lowest candidate altitude that meets the above conditions is determined as the minimum effective altitude.

[0046] The wind direction angle is represented by a circular average as follows:

[0047] ;

[0048] in, For the first The first hover data segment One wind direction angle; The average wind direction angle is obtained.

[0049] The obtained equivalent reference wind speed Mean wind angle and equivalent flux intensity Together they serve as the calibration wind field parameters.

[0050] In the measured data of X500, the equivalent reference wind speed in the incoming flow area is approximately 4.43 m / s, the average wind direction is approximately 190°, and the equivalent turbulence intensity is approximately 0.08; the equivalent reference wind speed in the wake area is approximately 2.66 m / s, the average wind direction is approximately 191°, and the equivalent turbulence intensity is approximately 0.20.

[0051] In each current control cycle Initially, the current feedback state of the UAV is mapped to a multi-physics domain digital twin; preferably, the initial state of the twin solution is represented as:

[0052] ;

[0053] in: The initial state is obtained by solving the twin solution for the current control cycle; This is the current location of the drone; This represents the current translation speed of the drone; The quaternion represents the current attitude of the drone; The angular velocity of the drone body; The current rotational speeds of the four rotors; This represents the current armature current of the four motors.

[0054] Among them, position, attitude and rotor speed are preferentially determined by the feedback values ​​obtained by the state perception and feedback module; if the translation speed or body angular velocity is not directly measured, the position and attitude changes of adjacent control cycles are used for differential determination; if the UAV is not equipped with motor current feedback, the motor current state corresponding to the end of the multi-physical domain digital twin of the previous control cycle is used as the initial current state of the current cycle to maintain the continuity of the internal state of the actuator between consecutive control cycles.

[0055] S12, the average wind speed based on the real-time flight altitude of the drone is expressed as:

[0056] ;

[0057] ;

[0058] in, The average wind speed at the current flight altitude; This is the current flight altitude of the drone; For a moment The corresponding effective flight altitude.

[0059] Based on the average wind direction angle, the average wind can be decomposed into an inertial coordinate system as follows:

[0060] ;

[0061] Random turbulence is represented as:

[0062] ;

[0063] in, For random turbulent flow velocity; The average wind direction angle; , , These represent the standard deviations of longitudinal, lateral, and vertical wind speeds, respectively. , , This represents the zero-mean random perturbation after filtering.

[0064] Stable hovering wind measurement data is selected in the mission area corresponding to the quadcopter UAV. The three-dimensional wind speed obtained at each sampling time is decomposed according to the wind direction coordinate system established by the average wind direction. The wind speed component along the average wind direction is taken as the longitudinal wind speed, the wind speed component perpendicular to the average wind direction in the horizontal plane is taken as the lateral wind speed, and the wind speed component in the vertical direction is taken as the vertical wind speed. The average value of the wind speed in each direction in the corresponding hovering data segment is calculated, and the fluctuation of each sampled wind speed relative to the average value in that direction is obtained. Then, the standard deviation is calculated based on all the fluctuations in each direction to obtain the standard deviation of longitudinal wind speed, lateral wind speed, and vertical wind speed.

[0065] The wind speed of the actual input digital twin, formed by superimposing the mean wind and random turbulence and applying a first-order filter, is expressed as:

[0066] ;

[0067] ;

[0068] in, For the first Wind speed input at each twin solution moment; Solve for the step size inside the digital twin; The wind field filtering time constant; The wind field filtering coefficient; It is a natural constant.

[0069] Candidate solution step sizes are determined based on the fastest-changing dynamic elements in the multi-physics domain digital twin. The fastest-changing dynamic elements are preferably the motor current establishment, motor speed establishment, or high-pulsation wind field change process in the actuator module. Under the condition of maintaining the same initial state, control command, and wind field input, the digital twin is run with different candidate solution step sizes from large to small in sequence, and the position, attitude, rotor speed, and wind field aerodynamic response under two adjacent solution step sizes are compared. When the peak value, response time, and steady-state value changes of the above key state quantities are all lower than the preset solution convergence error after further reducing the solution step size, the largest candidate step size that meets this condition is determined as the internal solution step size of the digital twin.

[0070] The mean wind component is removed from the X500 stable hovering measured wind speed sequence corresponding to the current calibrated wind field conditions to obtain a wind speed fluctuation sequence. Autocorrelation analysis is performed on this sequence to determine the time delay at which the correlation decays to approximately 1 / e of the initial correlation level. This time delay is then used as the initial value of the wind field filtering time constant. Subsequently, this initial value is written into the environmental module, and a digital twin is run under the same reference wind speed, mean wind direction, and turbulence intensity conditions as the measured wind field. The standard deviation, main fluctuation period, or power spectrum distribution of the filtered wind speed sequence are compared with those of the measured wind speed sequence, and adjustments are made within a preset candidate range. The candidate value that makes the fluctuation amplitude and main frequency characteristics of the digital twin wind speed closest to the measured wind speed will be determined as the final wind field filtering time constant.

[0071] Next, calculate the relative velocity of the incoming flow; let the velocity in the UAV's inertial coordinate system be... The relative wind speed, when converted to the body coordinate system, is expressed as:

[0072] ;

[0073] in, To be based on attitude quaternions The obtained rotation matrix from the body coordinate system to the inertial coordinate system; The relative wind speed is in the body coordinate system.

[0074] The corresponding wind field aerodynamic forces are expressed as follows:

[0075] ;

[0076] And converted to inertial coordinates as follows: ;

[0077] The aerodynamic moment of the wind field is expressed as: ;

[0078] in, air density; , The drag coefficients of the UAV along its three body axes; This represents the equivalent windward area in the three directions; The vector pointing from the center of mass of the UAV to the point of application of the equivalent aerodynamic force. and These refer to the aerodynamic force and aerodynamic torque of the wind field, respectively.

[0079] The air drag coefficient and equivalent frontal area are determined based on the manufacturer's aerodynamic parameters and are kept consistent within the digital twin of the same UAV, and are not used as parameters that change arbitrarily during the control process.

[0080] The airframe dynamics module solves for positional motion based on wind aerodynamic forces, total rotor thrust, and gravity, as follows:

[0081] ;

[0082] ;

[0083] ;

[0084] in, For drone quality; For the translational acceleration of the drone; Let be the translational velocity of the UAV in the inertial coordinate system; It is the acceleration due to gravity; The total thrust generated by the four rotors; For the first The thrust generated by each rotor; Let be the first derivative of the position vector with respect to time.

[0085] The attitude motion solution is expressed as:

[0086] ;

[0087] in, The inertial matrix of the UAV body; The angular velocity of the machine body; The roll, pitch, and yaw control moments are generated by the differential thrust of the four rotors; It is the angular acceleration of the machine body; The wind field aerodynamic torque is calculated by the environmental module.

[0088] The actuator module then establishes electromechanical dynamics for each motor; One motor satisfies:

[0089] ;

[0090] Electromagnetic torque satisfies: ;

[0091] The mechanical motion of the motor and rotor satisfies: ;

[0092] The corresponding rotor thrust is expressed as: ;

[0093] in, For the first The equivalent drive voltage of each motor; For the armature inductance of the motor; For the first The rate of change of armature current of an individual motor with respect to time; For the first The armature current of each motor at the current solution time; This refers to the armature resistance of the motor. It is the back electromotive force coefficient; For the first The electromagnetic torque generated by each motor based on the current armature current; This is the electromagnetic torque coefficient; This is the equivalent rotational inertia of the motor and rotor. This is the mechanical damping coefficient; For the first The aerodynamic load torque corresponding to each rotor; For the first Rotor speed; This is the rotor lift coefficient; This corresponds to the rotor thrust.

[0094] The back EMF coefficient, electromagnetic torque coefficient, mechanical damping coefficient, and rotor lift coefficient are obtained as follows: For the permanent magnet DC motor and corresponding rotor used, the rated parameters or test parameters provided by the motor and rotor manufacturers are first read, and calibration is performed through bench testing. Specifically, the back EMF coefficient is obtained by disconnecting the motor's active drive or by placing the motor in a low-load stable rotation state, measuring the motor terminal back EMF and corresponding angular velocity at multiple stable speed points, linearly fitting the back EMF and angular velocity data, and determining the slope of the fitted line as the back EMF coefficient. The electromagnetic torque coefficient is obtained by simultaneously measuring the armature current and output electromagnetic torque at multiple stable operating points using a motor dynamometer, linearly fitting the electromagnetic torque and armature current data after deducting no-load mechanical losses, and determining the slope of the fitted line as the back EMF coefficient. The slope is determined as the electromagnetic torque coefficient; the mechanical damping coefficient is obtained through free deceleration test or no-load steady-state test of the motor rotor assembly. Specifically, the motor is first driven to multiple preset speeds, and after the driving force is cut off, the decay process of the rotor speed over time is recorded. Combined with the known equivalent rotational inertia of the motor and rotor, the speed decay process is fitted, and the damping parameter that minimizes the error between the digital twin speed decay curve and the measured speed decay curve is determined as the mechanical damping coefficient; the rotor lift coefficient is obtained through rotor thrust bench test. The corresponding actual thrust is measured at multiple stable rotor speeds, and after converting each measured speed into angular velocity, the square data of rotor thrust angular velocity is fitted, and the proportional coefficient in the fitted relationship is determined as the rotor lift coefficient.

[0095] In a control cycle Internally, the step size is solved using a digital twin. The simultaneous solution for the aforementioned wind field, structure, and actuators is expressed as follows:

[0096] ;

[0097] ;

[0098] in: That is, the wind-borne airframe executes a synchronous response sequence corresponding to the current control cycle; Number the control cycle; For the first time in the current control cycle The digital twin solution time; This represents the solution sampling sequence number within the current control cycle, with values ​​ranging from 0 to... ; For a moment The aerodynamic vector of the wind field acting on the UAV body is obtained by the environment module; For a moment The corresponding wind field aerodynamic moment vector; For a moment The position vector of the UAV in the inertial coordinate system; For a moment The attitude quaternion of the drone; , , and They are time points The actual thrust generated by the first to fourth rotors; The number of solution intervals obtained within a control cycle according to the solution step size inside the digital twin; This is one control cycle of the adaptive fuzzy PID controller.

[0099] S13, based on the wind-borne airframe's synchronous response sequence, determine the wind load response time when the wind field aerodynamic force changes effectively and the execution response time when the rotor thrust changes effectively. According to the order and interval of the two, divide the current control cycle into wind load-first state, wind load and execution synchronization state, or execution completion follow-up state, and obtain the wind load execution timing discrimination result.

[0100] To avoid directly comparing aerodynamic forces, aerodynamic torques, and rotor thrust with different dimensions, we first establish dimensionless variations.

[0101] The degree of wind load variation is expressed as follows:

[0102] ;

[0103] in, The degree of dimensionless variation of wind load; This is the start time of the current control cycle; The characteristic length of the UAV is preferably the average distance from the center of mass to the center of the rotor.

[0104] The degree of change in actuator thrust is expressed as follows:

[0105] ;

[0106] ;

[0107] in, The degree of dimensionless change in the thrust of the actuator; This represents the baseline thrust for a single rotor in an ideal hovering state.

[0108] Set effective wind load change thresholds respectively and effective change threshold of the implementing agency Represented as:

[0109] ;

[0110] ;

[0111] in, , Within the reference operating range The mean and standard deviation; , Within the reference operating range The mean and standard deviation.

[0112] The reason for using the mean plus three times the standard deviation is to exclude sensor noise, small fluctuations in the solver, and normal motor speed ripple from the effective variation.

[0113] The wind load response time is expressed as: ;

[0114] The execution response time is represented as: ;

[0115] To prevent misjudgment caused by a single spike in random turbulence, only when the corresponding conditions are continuously maintained... The corresponding response time is only confirmed after a twin solution step.

[0116] The Based on the solution step size setting, it is preferable to express the corresponding duration as follows: ;

[0117] in, The continuous holding time for determining effective change; To maintain the number of sampling points required for continuous decision-making; The sampling step size is used to solve for the internal structure of the digital twin.

[0118] The sampling time should not exceed 10% of a control cycle and should not be less than 3 consecutive solution sampling points; this can filter out isolated random spikes without significantly increasing the time series discrimination delay.

[0119] Then, the wind load execution response time difference is determined to be expressed as: ;

[0120] in, The time difference between the effective thrust response of the actuator and the effective change in wind load; This refers to the wind load response time. This is the time to execute the response.

[0121] Further configure response synchronization tolerance The response synchronization tolerance The baseline response of the UAV actuator was obtained through benchmark tests. Under windless or stable wind conditions, multiple small-amplitude step control signals were applied to each motor. The normal time difference between the wind load solution channel and the actuator solution channel in the digital twin was statistically analyzed, and the 95th percentile of the absolute value of the normal time difference was used as... ;at the same time Not less than a digital twin solution step size .

[0122] Stable period of implementing agency The settling time was determined by the rotor thrust step response, specifically when the rotor thrust entered and remained within ±2% of the target thrust. After statistical analysis of four rotors and multiple test conditions, the 95th percentile value was selected as the optimal value. This avoids incorrect judgments of completing the follow-up due to a single, accidental, rapid response.

[0123] Based on the above time relationship, the wind load execution timing discrimination result is expressed as follows:

[0124] when ;

[0125] Furthermore, the implementing agency has not yet undergone a period of stabilization. At that time, it was determined that the wind load was in the leading position.

[0126] when At that time, it was determined that the wind load and execution were in a synchronized state;

[0127] When it has been detected ,and At that time, it is determined that the actuator has completed the thrust build-up process corresponding to this wind load change, and the current control cycle is set as the execution completion follow-up state, wherein: This refers to the current moment.

[0128] S14, the wind load execution timing judgment result is associated with the degree of wind field aerodynamic change, rotor thrust following degree and UAV position and attitude deviation degree within the corresponding control cycle to form the wind load execution response state, which is used to characterize the current wind disturbance has been applied to the body but the actuator has not yet completed the response, the actuator is following the wind load change or the actuator has completed the corresponding adjustment, and outputs it to the adaptive fuzzy PID controller.

[0129] First, based on the obtained dimensionless variation of wind load... The degree of change in wind field aerodynamics is expressed as follows:

[0130] ;

[0131] in: This represents the degree of change in wind field aerodynamics during the current control cycle.

[0132] when When this occurs, it indicates that the current wind load change has not yet exceeded the effective change threshold; when This indicates that an effective wind load change has occurred that requires a controller response.

[0133] The four target rotor thrusts corresponding to the current control quantity are expressed as follows:

[0134] ;

[0135] The actual rotor thrust is expressed as:

[0136] ;

[0137] The dimensionless thrust following residual is defined as follows:

[0138] ;

[0139] The degree of rotor thrust following is defined as follows:

[0140] ;

[0141] in, The combined thrust of the four rotors follows the residual. The degree of rotor thrust following is the range of values. ; The threshold value for thrust following residual in normal actuators.

[0142] The The values ​​were obtained from the actuator benchmark following test. Under stable wind conditions, different rotor speed commands were given, and the residuals between the steady-state actual thrust and the target thrust of each rotor were recorded. The mean of the residuals of all normal operating data plus three times the standard deviation was used as the threshold for the normal actuator thrust following residual. .

[0143] The degree of drone position offset is further defined as follows: ;

[0144] The degree of attitude deviation is expressed as: ;

[0145] in, The target location corresponding to the target trajectory; This is the current location of the drone; , , These are the current roll angle, pitch angle, and yaw angle of the drone, respectively. , , The corresponding target attitude angle; This is the allowable position deviation threshold; This is the allowable attitude deviation threshold.

[0146] The position allowable deviation threshold and attitude tolerance threshold The maximum permissible position error and maximum permissible attitude error are determined according to the current UAV flight mission. When the specific flight mission does not predetermine the corresponding values, the position error and attitude error are counted separately during the windless C0 reference trajectory tracking process, and their 99.7th percentile values ​​are used as the initial permissible position deviation threshold and attitude deviation threshold.

[0147] Finally, the wind load execution timing discrimination results will be obtained. Response time difference The degree of change in wind field aerodynamics Rotor thrust following degree Degree of positional offset and attitude deviation degree The wind load execution response state corresponding to the current control cycle is represented as follows:

[0148] ;

[0149] in, This refers to the wind load execution response status corresponding to the current control cycle; This is the status identifier corresponding to the wind load-first state, the wind load and execution synchronized state, or the execution completed and followed state. The degree of change in wind field aerodynamics; The degree of positional offset; This represents the degree of attitude deviation.

[0150] when In the case of wind load as the primary condition and Larger At a lower level, it indicates that the current wind disturbance has been significantly applied to the drone's body, while the actuators have not yet generated sufficient corresponding thrust; when When the wind load and the operation are synchronized, it indicates that the rotor thrust is building up in response to changes in wind load; when... To complete the follow state and higher and When the temperature has dropped, it indicates that the actuator has completed the corresponding adjustment, and the drone is gradually entering the convergence phase after the wind disturbance.

[0151] S2, the control error and error change rate corresponding to the wind load-execution response state and the target trajectory are input into the adaptive fuzzy PID controller. The correction intensity of the proportional, integral and derivative parameters is adjusted according to the response sequence between wind load change and rotor thrust establishment. When wind load change precedes rotor thrust establishment, the rapid suppression and damping effect is enhanced and the integral accumulation is limited. When rotor thrust completes following and the control error enters the convergence state, the steady-state correction is restored, and the adjusted control parameters are smoothly updated to obtain the phase-coordinated wind-resistant control quantity.

[0152] S21, based on the target position and attitude corresponding to the target trajectory in the current control cycle, and the current feedback state of the UAV obtained by the state perception and feedback module, establish a position control channel and an attitude control channel respectively; for any control channel The control error of the current control cycle is expressed as:

[0153] ;

[0154] in, For the first Control channel in each control cycle Control error; For the target trajectory in the th The target state quantity corresponding to each control cycle; This refers to the actual state quantity corresponding to the target state quantity in the current feedback state of the UAV; To control the channel number, it can correspond to the position direction or the roll, pitch and yaw attitude direction; Number the control cycle.

[0155] The rate of change of control error, determined according to the control error of adjacent control cycles, is expressed as follows:

[0156] ;

[0157] in, For the first The rate of change of control error per control cycle; This represents the control error from the previous control cycle. This refers to the control cycle of the adaptive fuzzy PID controller.

[0158] To ensure that the control errors of different position and attitude channels fall within a unified fuzzy inference range, the control error and the rate of change of control error are normalized respectively:

[0159] ;

[0160] ;

[0161] in, This is the normalized control error; This represents the normalized rate of change of control error. To control the channel The upper limit of error quantization; To control the channel The upper limit of the quantization of the rate of change of error; Indicates the variable It is limited to the range [-1, 1].

[0162] The control channel Upper limit of error quantization The maximum position or attitude error allowed for the corresponding flight mission is determined and verified using a windless baseline, a measured stable inflow condition, and a sudden wind disturbance condition. If a mission-defined maximum allowable error exists during flight, that allowable error is directly used as the initial quantization upper limit. If no clear indicator exists, the controller's error under a preset typical trajectory is statistically analyzed to determine the error value that covers normal control errors without causing a large number of inputs to remain in the saturation range for an extended period. .

[0163] The preset typical trajectory is pre-selected based on the actual flight mission of the quadcopter UAV. It is preferred to set representative trajectories that can cover common motion states such as hovering, straight flight, turning and climbing. The target position, target attitude and duration of each trajectory are determined according to the UAV's allowed flight speed, attitude change range and mission area size. Then, the trajectory is run under windless reference conditions. After confirming that the UAV can track stably, it is used as the preset typical trajectory for control error statistics and parameter calibration.

[0164] The control channel Upper limit of error change rate quantization The determination is based on the statistical results of the control error change rate in the preset typical trajectory changes and wind disturbance tests, specifically for no wind, steady wind, and sudden wind disturbance conditions. Statistical analysis was conducted, and the highest quantization value corresponding to the normal and effective control response was taken as the initial quantization upper limit.

[0165] Will and The language levels are divided into seven fuzzy levels: NB, NM, NS, ZO, PS, PM, and PB. Triangular or trapezoidal membership functions are used for fuzzification. For the normalization interval, ZO corresponds to the area near zero error, and adjacent language levels overlap continuously, thereby avoiding abrupt changes in parameter correction when the error crosses the level boundary.

[0166] Table 1 is a table of preset fuzzy rules.

[0167]

[0168] Using preset fuzzy rules to and Inference is performed to obtain the normalized fuzzy correction outputs corresponding to the proportional, integral, and differential parameters. , and The basic parameter correction amount is then converted and expressed as follows:

[0169] ;

[0170] ;

[0171] ;

[0172] in, , and These are the basic parameter corrections for proportional, integral, and differential parameters, respectively. , and This is the normalized output after fuzzy inference; and These represent the maximum allowable single correction amplitude for the three PID parameters.

[0173] The three maximum single correction amplitudes were obtained through PID parameter perturbation tests. Taking the basic PILD parameters that can stably complete target trajectory tracking under windless conditions as the center, the proportional, integral, and derivative parameters were increased and decreased respectively. Closed-loop tests were conducted under windless, stable wind, and high-pulsating wind conditions to determine the parameter change range that would not cause continuous oscillation, significant overshoot, or continuous saturation of the actuator. The common range that can work stably under all conditions was taken as the allowable correction range, thereby determining the corresponding maximum single correction amplitude.

[0174] From this, we can conclude that:

[0175] ;

[0176] in, For the first Control channel in each control cycle The corresponding basic parameter correction is obtained by the adaptive fuzzy PID controller through fuzzy inference based on the control error and the rate of change of error, and is used as the input for the subsequent parameter coordination process. For the first Control channel in each control cycle The basic correction amount of the proportional parameter is used to adjust the intensity of the proportional control effect; For the first Control channel in each control cycle The basic correction amount of the integral parameter is used to adjust the strength of the steady-state error elimination effect; For the first Control channel in each control cycle The differential parameter is used to adjust the strength of damping and suppression of the trend of control error change; The control channel number corresponds to the position control channel or the roll, pitch, and yaw attitude control channel; This is the control cycle number; the superscript 0 indicates that the parameter correction has not yet been phase-coordinated based on the wind load response status; the superscript... This represents the matrix transpose, used to form a column vector from the three parameter corrections.

[0177] S22, the wind load execution response status of the current control cycle is read as follows:

[0178] ;

[0179] in, This refers to the wind load execution response status corresponding to the current control cycle; This is the status identifier corresponding to the wind load-first state, the wind load and execution synchronized state, or the execution completed and followed state. The degree of change in wind field aerodynamics; The degree of positional offset; This represents the degree of attitude deviation.

[0180] when This indicates that the wind load is currently in a leading state, meaning that the aerodynamic forces of the wind field have changed effectively, while the rotor thrust has not yet been established accordingly. At this time, the basic parameter correction amount is not directly used, but the proportional, differential and integral correction strengths are coordinated secondaryly based on the wind load intensity and response time difference.

[0181] The advance disturbance rejection correction for the proportional parameter is expressed as:

[0182] ;

[0183] in, This is the proportional parameter correction amount under wind-load-first conditions; The degree of change in wind field aerodynamics; This is a sign function used to maintain the original direction of the fundamental parameter correction; This is the maximum phase coordination correction amplitude allowed by the scaling parameter.

[0184] because This has been obtained by normalizing the current degree of change in wind field aerodynamics relative to the effective change threshold of wind load. Therefore, when the wind load just reaches the effective change threshold... The value is close to 1, and the stronger the wind load, the larger the proportional correction amplitude, thus enhancing the rapid suppression effect. The maximum allowable phase coordination correction amplitude for the proportional parameter is... The stable operating range of the proportional parameter obtained from the aforementioned PID stability parameter test is determined so that the sum of the basic proportional parameter and its correction does not exceed the closed-loop stable operating boundary.

[0185] The prior disturbance rejection correction for the differential parameter is expressed as:

[0186] ;

[0187] ;

[0188] in, This is the differential parameter correction amount under the wind load-first condition; To retain only the non-negative time difference between the wind load and the executed response; The time synchronization tolerance used for the aforementioned wind load response and execution response timing discrimination; This represents the maximum phase coordination correction amplitude allowed by the differential parameters.

[0189] When the wind load application time and the rotor thrust response time are basically synchronized, the differential correction remains close to the basic correction amount. As the lag of the executed response relative to the wind load change increases, the intensity of the differential correction gradually increases to improve damping capability and suppress rapid error propagation. The maximum allowable phase coordination correction amplitude of the differential parameters is specified. Similarly, it is determined based on the stable operating range of the PID parameters.

[0190] The integral parameter is then reduced based on the degree of rotor thrust following, and is expressed as follows:

[0191] ;

[0192] in, This is the integral parameter correction amount under the wind-load-first condition; This represents the degree to which the current rotor thrust follows the target thrust; the closer it is to 1, the closer the actual rotor thrust is to the corresponding target thrust. For the first Control channel in each control cycle The basic correction amount of the integral parameter.

[0193] Before the rotor thrust is fully established, When the value is low, the integral correction automatically weakens; as the rotor thrust gradually completes its follow-up, the integral parameters gradually recover.

[0194] Simultaneously, the state of integral accumulation is restricted and represented as follows:

[0195] ;

[0196] in, This represents the integral accumulation state under the wind-load-first condition; This represents the integral state of the previous control cycle; The maximum allowable integral accumulation under normal steady-state control conditions was determined through steady-state tracking and actuator saturation tests. Under windless and stable wind conditions, the continuous deviation was gradually increased, and the integral state at which the steady-state error was eliminated without causing the corresponding motor control quantity to reach continuous saturation was recorded. The smaller value of the allowable integral state under each typical operating condition was taken as the maximum allowable integral state. .

[0197] The correction amount of the advance disturbance rejection parameter, consisting of the proportional, integral, and derivative correction results, is expressed as follows:

[0198] ;

[0199] in, For the first Control channel in each control cycle The corresponding correction amount for the advance disturbance rejection parameters.

[0200] The above treatment prioritizes enhancing the proportional rapid correction and differential damping effects when wind load changes have already acted on the airframe but the actual rotor thrust is still in the establishment phase, while reducing integral accumulation. This reduces the possibility of excessive compensation due to electrical, mechanical, and rotor thrust establishment delays in the actuators.

[0201] S23, in subsequent control cycles, the wind load execution response state is obtained from the multi-physical domain digital twin, and the wind load execution timing discrimination result, rotor thrust following degree, control error and control error change rate are jointly judged.

[0202] When the wind load execution timing determination result changes from the wind load-first state to the execution-completed-following state, the further determination is as follows:

[0203] ;

[0204] ;

[0205] ;

[0206] When all of the above conditions are met, the corresponding control channel is determined to have entered the control error convergence state.

[0207] in, The threshold value for the degree of thrust following is set to the preset thrust following requirement; To control the channel The preset error convergence threshold; To control the channel The preset error rate of change convergence threshold.

[0208] Preset thrust follow-up requirements and corresponding thrust follow-up degree threshold The data was obtained through normal rotor thrust tracking under windless baseline and stable inflow conditions; specifically, it was continuously collected after the UAV had entered a stable trajectory tracking state. When the distribution is approximately normal, the mean minus three times the standard deviation is used as the initial threshold; when the distribution deviates significantly from normal, the low quantile of the normal following sample is used as the threshold.

[0209] Preset error convergence threshold The control error is determined by the statistical results of the control error under windless stable tracking conditions and shall not exceed the control error allowed for the corresponding flight mission. When the stable error is approximately normally distributed, the mean and three times the standard deviation of the statistical range of the absolute value of the steady-state control error are used to determine the initial convergence range. When quantile statistics are used, the 99.7th percentile of the normal steady-state error is used as the initial value.

[0210] Preset error rate of change convergence threshold The control error change rate during the stable trajectory tracking phase was determined using statistical results. The setting principle was to cover normal steady-state micro-oscillations while being significantly smaller than the error change rate when wind disturbances occur or the trajectory changes rapidly.

[0211] After confirming that the control error has converged, the PID parameters are not immediately switched to their original steady-state values. Instead, the advance disturbance rejection parameter correction is gradually restored to the basic parameter correction. The restoration coefficient is defined as follows:

[0212] ;

[0213] in, The steady-state recovery coefficient for each control cycle takes a value between 0 and 1; To control the cycle; To determine the steady-state parameter recovery time constant, multiple candidate parameters are set under the same target trajectory and the same wind disturbance conditions. The second overshoot, PID parameter recovery speed, and time required to re-enter the stable error range were recorded after the rotor thrust completed following. Candidate values ​​causing significant second overshoot, parameter abrupt changes, or excessively long recovery times were eliminated. The smallest time constant among the remaining candidate values ​​that allows the proportional and derivative corrections to weaken smoothly, the integral action to gradually recover, and the control error to stabilize quickly was selected as the control error control value. .

[0214] The proportional, integral, and derivative parameters are expressed as follows:

[0215] ;

[0216] in, ; This represents the steady-state recovery parameter correction amount corresponding to the current control cycle; when just entering the control error convergence state... Take the correction amount of the advance disturbance rejection parameter corresponding to the previous control cycle; for subsequent cycles, take the correction amount of the steady-state recovery parameter of the previous cycle. The steady-state recovery coefficient for each control cycle.

[0217] The gradual restoration of the integration limit is expressed as follows:

[0218] ;

[0219] in, The integral limit allowed in the current control cycle; the initial value is the integral limit when the wind load precedence state is exited. This is the limit for the integral in normal steady state; The steady-state recovery coefficient for each control cycle.

[0220] The steady-state recovery parameter correction obtained from the above processing is expressed as follows:

[0221] ;

[0222] in, It is the steady-state recovery parameter correction amount; its function is to gradually weaken the proportional and differential correction effects enhanced by the wind load in the early stage, while removing the integral accumulation limit, so that the integral parameter can once again assume the function of eliminating steady-state deviation. This is the proportional parameter correction amount during the steady-state recovery phase; This is the integral parameter correction amount during the steady-state recovery phase; This represents the differential parameter correction amount during the steady-state recovery phase.

[0223] S24, determine the parameter correction stage corresponding to the current control cycle based on the current wind load execution timing judgment result; when the current cycle is still in the wind load leading state, the aforementioned leading disturbance rejection parameter correction amount is adopted; when the rotor thrust has completed following and has entered the control error convergence state, the aforementioned steady-state recovery parameter correction amount is adopted; for the brief wind load and execution synchronization state transitioning from the wind load leading state to the execution completion following state, maintain the phase coordination parameter correction direction of the previous control cycle, and transition to the current basic parameter correction amount through subsequent smooth updates to avoid PID parameter mutations at the state switching boundary.

[0224] Therefore, the phase coordination parameter correction for the current control cycle is expressed as:

[0225] ;

[0226] in, For the first Control channel in each control cycle The corresponding phase coordination parameter correction amount; This is the correction amount for the phase coordination proportional parameter in the current control cycle; This is the correction amount for the phase coordination integral parameter in the current control cycle; This is the correction amount for the phase coordination differential parameter in the current control cycle.

[0227] Combining this with the basic proportional parameters, basic integral parameters, and basic derivative parameters of the adaptive fuzzy PID controller, the parameters to be updated in the current period are expressed as follows:

[0228] ;

[0229] ;

[0230] ;

[0231] in, , and Control channels The fundamental proportional parameters, fundamental integral parameters, and fundamental differential parameters; , and The parameters to be updated are obtained after adding the phase coordination parameter correction.

[0232] The basic PID parameters are determined by the existing controller calibration method, using response time, overshoot, steady-state error, and control output under windless hovering, trajectory tracking, and stable wind disturbance conditions. This allows the UAV to maintain basic closed-loop stability without phase coordination correction. The phase coordination process then dynamically corrects the response sequence between wind disturbance and the actuator. The proportional action is increased for large errors, the integral action is increased for small errors, and the derivative action is increased when the error changes rapidly.

[0233] To avoid abrupt changes in control parameters across control cycles, the parameters to be updated are represented using first-order smoothing as follows:

[0234] ;

[0235] in, ; To smoothly update the PID control parameters; These are the control parameters from the previous control cycle; The control parameters are expected to be updated this week; The parameter is the smoothing coefficient.

[0236] The smoothing coefficient of the parameter is determined according to: It is confirmed that, among them, The smoothing time constant is used to control parameters.

[0237] The The parameters were obtained through a step update experiment using PID parameters. While keeping the target trajectory and wind field input constant, preset changes were applied to the proportional, integral, and derivative target parameters. Different candidate smoothing time constants were used, and the actual parameter change curves, control output jumps, attitude overshoot, and settling time were compared. Candidate values ​​that avoid significant control output jumps and do not cause significant delays in wind disturbance response were selected as the optimal values. .

[0238] To further reduce the amplification effect of measurement noise and high-frequency errors on the differential term, a first-order filter is applied to the rate of change of the control error, which is expressed as:

[0239] ;

[0240] in, This represents the rate of change of the control error after filtering in the current cycle. This represents the rate of change of control error after filtering in the previous control cycle. These are the differential filter coefficients.

[0241] The differential filter coefficients are expressed as:

[0242] ;

[0243] in, is the time constant of the differential filter.

[0244] Differential filter time constant Based on the high-frequency noise characteristics of the state feedback signal in a stable hovering state, spectral analysis is performed on the position, attitude, or error signals to ensure that the determined differential filter cutoff frequency is lower than the frequency range where noise dominates, while being higher than the main frequency range of the UAV's normal attitude and position closed-loop response. This allows for the suppression of high-frequency measurement noise while preserving the true control error variation.

[0245] Ultimately, for any control channel The control parameters, after smooth update, are calculated and expressed as follows:

[0246] ;

[0247] in, For the first Control channel in each control cycle Output control quantity; For the first The proportional parameters updated after phase coordination and smoothing in each control cycle; This refers to the control error between the target state corresponding to the target trajectory and the actual feedback state of the UAV in the current control cycle. These are the integral parameters after phase coordination and smoothing updates; This represents the integral state of the control error after integral accumulation and limitation during the current control cycle. These are the differential parameters after phase coordination and smoothing updates; This is the filtered error rate obtained after differential filtering of the control error rate of change; and Control channels The minimum and maximum control values ​​allowed for output are determined based on the actual controllable range of the corresponding actuator. For the limiting function, when the calculated control quantity is lower than... Time to take Higher than Time to take When the value is in between, the original calculated value is retained.

[0248] Control Channel The minimum and maximum control quantities allowed for output are determined based on the actual control capabilities of the UAV actuator. Specifically, they are calculated by combining the allowable voltage and current of the motor, the maximum rotor speed, and the corresponding maximum and minimum achievable thrust of the rotor, so that the output required by the controller does not exceed the control range that the actuator can actually establish. When using a layered structure of position outer loop and attitude inner loop, the position outer loop output first forms the target thrust and target attitude, and then the attitude inner loop forms the roll, pitch and yaw control torques.

[0249] The outputs of each control channel are combined according to the predetermined control structure of the position outer loop and attitude inner loop, and the phase-coordinated wind-resistant control quantity for the current control cycle is expressed as follows:

[0250] ;

[0251] in, The phase-coordinated wind-resistant control quantity for the current control cycle; Total thrust for the target; , and These are the target control torques in the roll, pitch, and yaw directions, respectively.

[0252] S3, the phase-coordinated wind-resistant control quantity is converted into independent control signals for the four rotors according to the rotor layout of the quadcopter UAV, driving the corresponding actuators to complete position and attitude correction, obtaining the corrected UAV position, attitude and rotor speed as closed-loop feedback state, and sending the closed-loop feedback state back to the multi-physical domain digital twin to update the wind load execution response state for the next control cycle.

[0253] S31, Read the phase coordination wind resistance control quantity obtained in the current control cycle, wherein the phase coordination wind resistance control quantity is expressed as:

[0254] ;

[0255] in, For the first Phase-coordinated wind-resistant control quantity for each control cycle; Total thrust for the target; For the target control torque of the roll; For pitch control torque; Control torque for yaw target; Number the control cycle.

[0256] Based on the X-shaped rotor layout of the quadcopter drone, the positions of the first to fourth rotors relative to the center of mass of the fuselage are represented as follows:

[0257] ;

[0258] in, For the first The position vector of each rotor center relative to the center of mass of the UAV body; and The first The coordinates of each rotor along the two horizontal axes of the airframe coordinate system; The rotors are numbered; in an X-shaped configuration, the four rotors are located diagonally opposite the fuselage's center of mass, therefore... and The positive and negative relationship is determined based on the actual rotor number and installation direction.

[0259] The quadcopter does not directly apply the total thrust and total torque to the center of mass. Instead, it applies the thrust of each of the four rotors to the corresponding positions on the fuselage, allowing the roll, pitch, and yaw moments to be naturally generated through the action of each rotor.

[0260] The target thrust of the four rotors is expressed as:

[0261] ;

[0262] The control allocation relationship corresponding to the X-type rotor configuration is represented as follows:

[0263] ;

[0264] in, The target thrust vector for the four rotors; For the first The target thrust of each rotor; The control allocation matrix corresponding to the X-type rotor configuration; For the first Phase-coordinated wind-resistant control quantity for each control cycle.

[0265] The control allocation matrix is ​​represented as follows:

[0266] ;

[0267] The first row is used to form the total target thrust; the second row forms the roll control torque based on the lateral distance of each rotor relative to the center of mass; and the third row forms the pitch control torque based on the longitudinal distance of each rotor relative to the center of mass. For the first The yaw moment conversion coefficient of each rotor is determined by the rotation direction of the corresponding rotor.

[0268] and Obtained directly from the actual frame dimensions of the quadcopter drone; The sign is determined by the clockwise or counterclockwise rotation direction of the rotor, and its absolute value is obtained by bench testing the anti-torque of the corresponding motor-rotor combination at different stable speeds. It is obtained by fitting the relationship between the rotor anti-torque and the corresponding thrust. For four rotors using the same type of motor and blades, the absolute value adopts the same calibration value, but is assigned a positive or negative sign according to the rotation direction; adjacent rotors adopt opposite directions and yaw adjustment is achieved through rotor anti-torque.

[0269] Therefore, we can conclude that: When established based on specific rotor numbers When reversible, the thrust of the four target rotors can be obtained through its inverse matrix.

[0270] Based on the aforementioned rotor lift relationship, the target thrust can be further converted into the target angular velocity of each rotor as follows:

[0271] ;

[0272] in, For the first The first rotor in the... The target angular velocity for each control cycle; For the corresponding target thrust; This represents the rotor lift coefficient.

[0273] To avoid the control allocation result exceeding the actual driving capability of the motor rotor assembly, the target angular velocity is limited and expressed as follows:

[0274] ;

[0275] in, For the first The final independent control signal for each rotor; and The first The minimum and maximum operating angular velocities allowed for each motor rotor assembly.

[0276] Minimum operating angular velocity and maximum working angular velocity The speed limit is determined primarily based on the permissible speed range provided by the motor, inverter, and rotor manufacturers, and verified through continuous operation tests on the motor and rotor test bench. The final limit is the speed range where the motor current, temperature rise, and rotor mechanical condition are all within the permissible range.

[0277] Four independent control signals Together they constitute the independent control signals for the four rotors.

[0278] S32, input the independent control signals of the first to fourth rotors to the corresponding actuators; for the first... Each actuator will control the target rotor angular velocity. Compared with the current actual rotor angular velocity By comparison, the rotor speed error is expressed as:

[0279] ;

[0280] in, For the first Rotor speed error of each actuator; The target rotor angular velocity corresponding to the current control cycle; This represents the actual rotor angular velocity at the current moment; This refers to the digital twin solution time within the current control cycle.

[0281] The target current value generated by the speed closed loop in the actuator based on the rotor speed error is expressed as follows:

[0282] ;

[0283] in, For the first The target armature current of each motor; and These are the proportional coefficient and integral coefficient of the velocity closed loop, respectively.

[0284] The motor drive voltage obtained from the current closed loop is expressed as:

[0285] ;

[0286] in, Apply to the inverter to the first The equivalent drive voltage of a permanent magnet DC motor; This is the actual armature current; and These are the proportional coefficient and integral coefficient of the current closed loop, respectively.

[0287] The , , and First, determine the initial values ​​based on the resistance, inductance, moment of inertia, and rated operating range of the motor used. Then, perform calibration using speed step tests and current step tests respectively. Gradually adjust each coefficient to enable the current loop to follow the target current quickly without continuous oscillation, and to enable the speed loop to follow the target speed within the allowable overshoot range. Finally, select the parameter set that simultaneously satisfies the requirements of response time, overshoot, and stability as the final value.

[0288] Subsequently, the motor current, electromagnetic torque, and rotor angular velocity were continuously solved according to the aforementioned electrical equations, electromagnetic torque equations, and rotor mechanical dynamics equations, and based on... Get the first The actual thrust of each rotor.

[0289] The actual total thrust of the four rotors is expressed as:

[0290] ;

[0291] The roll and pitch moments generated by the actual installation positions of each rotor on the fuselage are uniformly expressed as:

[0292] ;

[0293] in, This is the resultant torque of roll and pitch formed by the actual thrust of the four rotors relative to the center of mass of the fuselage.

[0294] The rotational counter-torque of each rotor forms the actual yaw control torque; the actual total thrust and the actual control torques in the roll, pitch and yaw directions together form the airframe action, which is input into the airframe dynamics module and acts on the UAV airframe together with the current wind field aerodynamic force and aerodynamic torque, continuously solving the UAV's position and attitude changes according to the aforementioned translational dynamics and rotational dynamics equations.

[0295] After the actuator response and body motion solution of the current control cycle, the UAV position, UAV attitude and actual rotation speed of each rotor at the end of the current cycle are obtained, and the body state after the control execution is defined as the corrected flight state.

[0296] S33, at the end of the current control cycle, the state perception and feedback module performs unified sampling of the corrected flight state to obtain the corrected UAV position, UAV attitude, and rotor speed of the four rotors.

[0297] The location of the drone is represented as:

[0298] ;

[0299] in, For the first The location of the drone after each control cycle is completed; , and These are the position coordinates in the three directions of the inertial coordinate system.

[0300] The preferred attitude of the UAV is represented in quaternion form, consistent with the aforementioned airframe dynamics module, as follows:

[0301] ;

[0302] in, The corrected attitude quaternion for the UAV; The real part of the quaternion; , and It consists of three imaginary parts.

[0303] The actual rotational speeds of the four rotors are expressed as follows:

[0304] ;

[0305] in, This is the actual rotational speed vector of the four rotors; , , and These represent the actual angular velocities of the first to fourth rotors at the sampling time of the current control cycle.

[0306] The above states will be identified using a unified time identifier. The association is represented as:

[0307] ;

[0308] in, For the first The closed-loop feedback state corresponding to each control cycle; Used as a unified time identifier for the current control cycle; superscript This indicates the matrix transpose.

[0309] In the Modelica unified solution environment used in this implementation, all state variables are solved under the same model clock. Therefore, the unified time identifier is directly taken from the fixed sampling time of the controller and represented as follows:

[0310] ;

[0311] in, This refers to the control cycle of the aforementioned adaptive fuzzy PID controller.

[0312] S34, the closed-loop feedback state The data is sent back to a multi-physical domain digital twin; among other things, the drone's location is... and drone attitude Write the data into the airframe dynamics module to replace the corresponding position and attitude state from the previous cycle; and set the actual rotational speeds of the four rotors. Write to the actuator module to replace the corresponding rotor speed state of the previous cycle.

[0313] For continuous states such as body translational velocity, body angular velocity, and motor current that are not directly updated by the closed-loop feedback state, the final values ​​of the corresponding states within the multi-physics domain digital twin at the end of the current control cycle are used to form the initial state for the twin solution of the next control cycle, expressed as follows:

[0314] ;

[0315] in, For the first The initial state of the twin solution for each control cycle; The translational velocity of the UAV in the body dynamics module at the end of the current control cycle; This is the angular velocity of the UAV at the end of the current control cycle. This represents the armature current vectors of the four motors at the end of the current control cycle.

[0316] The purpose is to use closed-loop feedback to correct position, attitude and rotor speed, while maintaining the continuity of the previous solution cycle for continuous internal states not directly provided by the feedback state, so as to avoid sudden changes in the digital twin solution state caused by manually resetting speed, angular velocity or motor current between adjacent control cycles.

[0317] The first The calibration wind field parameters corresponding to each control cycle are continued to be written into the environment module, and are used as described above. As a new starting point for the solution, the wind-borne airframe execution synchronization response sequence for the next control cycle is obtained through the synchronous operation of the environment module, the body dynamics module, and the actuator module, and is expressed as follows:

[0318] ;

[0319] in, For the first The wind-borne aircraft executes a synchronous response sequence for each control cycle.

[0320] The wind load response time and execution response time for the next control cycle are determined from the synchronous response sequence of the wind-borne airframe, resulting in a new wind load execution timing discrimination result. This result is further updated based on the degree of change in the associated wind field aerodynamic forces, the degree of rotor thrust following, and the degree of UAV position and attitude deviation.

[0321] ;

[0322] in, This refers to the wind load execution response state for the next control cycle.

[0323] This forms a continuous closed loop consisting of phase-coordinated wind-resistant control quantity, independent control signal, corrected flight state, closed-loop feedback state, and wind load execution response state for the next control cycle.

[0324] Example 2

[0325] like Figure 4 As shown, the multi-physics domain wind disturbance resistance control system for quadcopter UAVs based on digital twins includes the following modules:

[0326] A multi-physical domain digital twin module is used to simultaneously solve the wind field aerodynamics, airframe motion, and rotor thrust establishment process based on the calibrated wind field parameters and the current feedback state of the UAV, so as to obtain the wind load execution response state.

[0327] An adaptive fuzzy PID control module is used to adjust the proportional parameter, integral parameter and derivative parameter according to the wind load execution response state, control error and error change rate to obtain the phase-coordinated wind-resistant control quantity.

[0328] A control allocation module is used to convert the phase-coordinated wind-resistant control quantity into independent control signals for the four rotors according to the X-shaped rotor layout of the quadcopter UAV.

[0329] An actuator module is used to adjust the thrust of each rotor according to the independent control signal through an electrical drive, motor torque establishment and rotor speed establishment process, so as to complete the position and attitude correction of the UAV;

[0330] The state perception and feedback module is used to acquire the corrected UAV position, attitude and rotor speed and form a closed-loop feedback state. The closed-loop feedback state is sent back to the multi-physical domain digital twin module to update the wind load execution response state for the next control cycle.

[0331] This invention encompasses any substitutions, modifications, equivalent methods, and solutions made within the spirit and scope of this invention. To provide the public with a thorough understanding of this invention, specific details are described in detail in the following preferred embodiments; however, those skilled in the art will fully understand the invention even without these details. Furthermore, to avoid unnecessary misunderstanding of the essence of this invention, well-known methods, processes, procedures, components, and circuits are not described in detail.

[0332] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A multi-physics domain wind disturbance resistance control method for quadrotor UAVs based on digital twins, characterized in that, Includes the following steps: Based on the multi-physics domain digital twin corresponding to the quadcopter drone, the wind field aerodynamics, body motion, and the thrust establishment process from the motor to the rotor are solved simultaneously according to the calibrated wind field parameters and the current feedback state of the drone. Determine the sequential relationship between wind load changes and rotor thrust establishment to obtain the wind load execution response state, which characterizes the degree of coordination between wind disturbance and the actual response of the actuator. Based on the control error and error change rate corresponding to the wind load execution response state and the target trajectory, the adaptive fuzzy PID controller is input, and the correction intensity of the proportional, integral and derivative parameters is adjusted according to the response sequence relationship between wind load change and rotor thrust establishment. When wind load changes precede rotor thrust establishment, the rapid suppression and damping effects are enhanced and the integral accumulation is limited. When rotor thrust completes follow-up and control error enters a convergent state, steady-state correction is restored. Based on the smooth update of the adjusted control parameters, the phase-coordinated wind-resistant control quantity is obtained; Based on the phase coordination wind resistance control quantity, it is converted into independent control signals for the four rotors according to the rotor layout of the quadcopter UAV, driving the corresponding actuators to complete position and attitude correction. The corrected UAV position, attitude, and rotor speed are obtained as closed-loop feedback states, and the closed-loop feedback states are sent back to the multi-physics domain digital twin to update the wind load execution response state for the next control cycle.

2. The method for multi-physics domain wind disturbance resistance control of a quadrotor UAV based on digital twin according to claim 1, characterized in that, A multi-physics domain digital twin is established, consisting of an environment module, a body dynamics module, an actuator module, and a state perception and feedback module. Calibrated wind field parameters are written into the environment module, and the current feedback state of the UAV is mapped to the body dynamics module and the actuator module. The wind direction angle is processed using a circular average method and expressed as follows: ； in, For the first The first hover data segment One wind direction angle; To obtain the average wind direction angle, For the first The total number of effective wind direction sampling points in each effective wind field data segment Number the wind direction sampling points within the data segment to form the initial state of the twin solution corresponding to the current control cycle; Based on the initial state of the twin solution, the environment module generates a wind field input after superimposing the average wind and random turbulence according to the calibrated wind field parameters, and obtains the wind field aerodynamic forces and aerodynamic moments acting on the UAV body through the aerodynamic model. The corresponding wind field aerodynamic forces are expressed as follows: ； And converted to inertial coordinates as follows: ; The aerodynamic moment of the wind field is expressed as: ; in, This is the aerodynamic vector of the wind field in the body coordinate system; air density; , , These are the body along the body coordinate system. , , Aerodynamic drag coefficient in the direction; , , These are the equivalent frontal areas of the aircraft in the corresponding three directions; , , Relative wind speeds Components in the three directions of the body coordinate system; , , These represent the absolute values ​​of the corresponding relative wind speed components; This represents the aerodynamic vector of the wind field after transformation to the inertial coordinate system. The quaternion of the current attitude of the drone The direction cosine matrix from the determined body coordinate system to the inertial coordinate system; The quaternion represents the current attitude of the drone; The aerodynamic moment vector formed by the wind field aerodynamic forces relative to the center of mass of the UAV body; The position vector of the equivalent point of action of the wind field aerodynamic force relative to the center of mass of the UAV body is given. The position and attitude response of the UAV are obtained synchronously by the body dynamics module, and the thrust response of each rotor is obtained by the actuator module according to the transmission process of the control signal through electrical drive, motor torque, rotor speed to rotor thrust, forming a wind-borne body execution synchronous response sequence corresponding to the same control cycle.

3. The method for multi-physics domain wind disturbance resistance control of a quadrotor UAV based on digital twin according to claim 2, characterized in that, Based on the synchronous response sequence executed by the wind-borne airframe, the wind load response time when the wind field aerodynamic force changes effectively and the corresponding execution response time when the rotor thrust changes effectively are determined. The degree of wind load change is expressed as follows: ； The degree of change in actuator thrust is expressed as follows: ； ； in, For a moment The dimensionless change in wind load relative to the start of the current control cycle; and They are time points and the start time of the current control cycle The wind field aerodynamic vector; and These are the wind field aerodynamic moment vectors at the corresponding times; Represents the 2-norm of a vector; This refers to the total mass of the quadcopter drone; It is the acceleration due to gravity; The characteristic length of the UAV; This is the start time of the current control cycle; The digital twin solution time within the current control cycle; For a moment The dimensionless change in actuator thrust relative to the start of the current control cycle; For the first Each rotor at time The actual thrust; For the first The actual thrust of each rotor at the start of the current control cycle; Number the rotor blades; The reference thrust for a single rotor in the hovering state of the UAV is determined, and the current control cycle is divided into wind load-first state, wind load and execution synchronization state, or execution completion follow state according to the time sequence and interval between the two, so as to obtain the wind load execution timing discrimination result.

4. The method for multi-physics domain wind disturbance control of a quadrotor UAV based on digital twin according to claim 3, characterized in that, The wind load execution timing judgment result is correlated with the degree of wind field aerodynamic change, rotor thrust following degree, and UAV position and attitude deviation degree within the corresponding control cycle to form the wind load execution response state, which characterizes whether the current wind disturbance has been applied to the airframe but the actuator has not yet completed the response, the actuator is following the wind load change, or the actuator has completed the corresponding adjustment. The dimensionless thrust following residual is defined as: ； The degree of rotor thrust following is defined as follows: ； in, The combined thrust of the four rotors follows the residual. The degree to which rotor thrust follows; For the normal actuator thrust following residual threshold, For the first The target thrust is obtained by controlling and distributing the thrust of each rotor according to the phase-coordinated wind resistance control parameters. For the first The actual thrust generated by each rotor after being electrically driven, having its motor torque established, and its rotor speed established is output to an adaptive fuzzy PID controller.

5. The method for multi-physics domain wind disturbance resistance control of a quadrotor UAV based on digital twin according to claim 1, characterized in that, Based on the target position and attitude corresponding to the target trajectory and the current feedback state of the UAV, the control error and error rate of change are determined. The control error and error rate of change are quantized and fuzzified, and the basic parameter correction quantities of the proportional parameter, integral parameter, and derivative parameter are obtained according to the corresponding fuzzy rules, respectively, and expressed as follows: ； ； ； in, , and These are the basic parameter corrections for proportional, integral, and differential parameters, respectively. , and This is the normalized output after fuzzy inference; and The three PID parameters are respectively the maximum allowable single correction amplitude, and the basic parameter correction amount and the wind load execution response state of the current control cycle are jointly input into the parameter coordination link; Based on the wind load execution timing discrimination result, the degree of wind field aerodynamic change, and the rotor thrust following degree in the wind load execution response state, the proportional correction intensity, integral correction intensity, and differential correction intensity of the basic parameter correction are coordinated. The prior disturbance rejection correction of the proportional parameter is expressed as: ； in, This is the proportional parameter correction amount under wind-load-first conditions; The degree of change in wind field aerodynamics; This is a sign function used to maintain the original direction of the fundamental parameter correction; This is the maximum phase coordination correction amplitude allowed by the scaling parameter; The prior disturbance rejection correction for the differential parameter is expressed as: ； ； in, This is the differential parameter correction amount under the wind load-first condition; To retain only the non-negative time difference between the wind load and the executed response; The time synchronization tolerance used for the aforementioned wind load response and execution response timing discrimination; This is the maximum phase coordination correction magnitude allowed by the differential parameters. Let be the time difference between the wind load response time and the execution response time. When the wind load execution timing judgment result is a wind load-first state, the proportional correction intensity is increased according to the degree of wind field aerodynamic change, and the differential correction intensity is increased according to the time difference between the wind load response time and the execution response time. At the same time, the integral correction intensity is reduced according to the rotor thrust following degree, and the current integral accumulation is limited. The limitation on the integral accumulation state is expressed as follows: ； in, This represents the integral accumulation state under the wind-load-first condition; This represents the integral state of the previous control cycle; This represents the maximum allowable integral accumulation under normal steady-state control conditions. For the first Control error per control cycle For the amplitude limiting function, Given the rotor thrust following degree in the current control cycle, a preliminary disturbance rejection parameter correction is obtained to improve rapid suppression and damping capabilities during the period before rotor thrust has fully established.

6. The method for multi-physics domain wind disturbance resistance control of a quadrotor UAV based on digital twin according to claim 5, characterized in that, The wind load execution response state of subsequent control cycles is continuously used to verify the correction amount of the preliminary disturbance rejection parameters. When the wind load execution timing judgment result changes to the execution completion following state, and the rotor thrust following degree reaches the preset thrust following requirement, and the control error and error change rate enter the preset error convergence range and error change rate convergence range respectively, it is determined that the UAV has entered the control error convergence state. The proportional correction intensity and derivative correction intensity in the preliminary disturbance rejection parameter correction amount are gradually reduced, and the integral accumulation limit is removed, restoring the correction effect of the integral parameter on the steady-state deviation. The proportional, integral, and derivative parameters are expressed as follows: ； in, ; This is the correction amount for the steady-state recovery parameters corresponding to the current control cycle; The steady-state recovery coefficient for each control cycle; The gradual restoration of the integration limit is expressed as follows: ； in, This is the integral limit allowed for the current control cycle; The normal steady-state integral limit is used to obtain the steady-state recovery parameter correction.

7. The method for multi-physics domain wind disturbance control of a quadrotor UAV based on digital twin according to claim 6, characterized in that, Based on the current wind load execution timing judgment result, the phase coordination parameter correction amount corresponding to the current control cycle is determined from the prior disturbance rejection parameter correction amount or the steady-state recovery parameter correction amount. This phase coordination parameter correction amount is then combined with the basic proportional parameter, basic integral parameter, and basic derivative parameter of the adaptive fuzzy PID controller. A smooth transition is performed between the control parameters of the previous control cycle and the control parameters to be updated in the current cycle. This involves limiting the integral term, filtering the derivative term, and limiting the control output. The parameters to be updated in the current cycle are expressed as follows: ； ； ； in, , and Control channels The fundamental proportional parameters, fundamental integral parameters, and fundamental differential parameters; , and The parameters to be updated are obtained after incorporating the phase coordination parameter correction. This is the correction amount for the phase coordination proportional parameter in the current control cycle; This is the correction amount for the phase coordination integral parameter in the current control cycle; The phase coordination differential parameter correction for the current control cycle is expressed as follows: The parameter to be updated is represented by a first-order smoothing: ； in, ; To smoothly update the PID control parameters; These are the control parameters from the previous control cycle; The control parameters are expected to be updated this week; The parameter smoothing coefficient is used to calculate the phase-coordinated wind-resistant control quantity for the current control cycle using the smoothed and updated proportional, integral, and derivative parameters.

8. The method for multi-physics domain wind disturbance resistance control of a quadrotor UAV based on digital twin according to claim 1, characterized in that, Read the phase-coordinated wind-resistant control quantities of the current control cycle. Based on the X-shaped rotor layout of the quadcopter UAV and the position and rotation direction of each rotor relative to the center of mass of the fuselage, control the target total thrust, roll target control torque, pitch target control torque, and yaw target control torque in the phase-coordinated wind-resistant control quantities. Determine the target drive quantities corresponding to the first to fourth rotors respectively, and obtain the independent control signals of the four rotors as follows: ;in, For the first The final independent control signal for each rotor; and The first Minimum and maximum operating angular velocities of each motor rotor assembly For the amplitude limiting function, For the first The first rotor in the... The target angular velocity for each control cycle; The independent control signals of the four rotors are input to their respective actuators. Through electrical drive, motor torque establishment, and rotor speed establishment processes, the actual thrust of each rotor is adjusted. The thrust combination among the four rotors forms the airframe action corresponding to the target total thrust, roll target control torque, pitch target control torque, and yaw target control torque. The speed closed loop in the actuator generates a current target value based on the rotor speed error, expressed as follows: ； in, For the first The target armature current of each motor; and These are the proportional and integral coefficients for the speed closed loop, respectively. For the first The rotational speed error of each rotor This is the integral of the rotational speed error over time; The motor drive voltage obtained from the current closed loop is expressed as: ； in, Apply to the inverter to the first The equivalent drive voltage of a permanent magnet DC motor; This is the actual armature current; and These are the proportional coefficient and integral coefficient of the current closed-loop circuit, respectively. The current tracking error is the integral of time, which drives the quadcopter UAV to complete position and attitude correction, and obtain the corrected flight state.

9. The method for multi-physics domain wind disturbance control of a quadrotor UAV based on digital twin according to claim 8, characterized in that, The corrected flight state is synchronously collected using the state perception and feedback module to obtain the corrected UAV position, UAV attitude, and rotor speeds of the first to fourth rotors. The UAV position is represented as: ； in, For the first The location of the drone after each control cycle is completed; , and These are the position coordinates in the three directions of the inertial coordinate system; The optimal attitude of the UAV is represented in quaternion form as follows: ； in, The corrected attitude quaternion for the UAV; The real part of the quaternion; , and It has three imaginary parts; The actual rotational speeds of the four rotors are expressed as follows: ； in, This is the actual rotational speed vector of the four rotors; , , and These are the actual angular velocities of the first to fourth rotors at the sampling time of the current control cycle, and are associated according to the unified time identifier of the current control cycle to form a closed-loop feedback state that characterizes the motion state of the aircraft and the actual response state of the actuator after control execution. The closed-loop feedback state is sent back to the multi-physics domain digital twin, where the UAV position and attitude are updated to the airframe dynamics module, the rotor speeds of the first to fourth rotors are updated to the actuator module, and the updated state is used as the initial state for the twin solution of the next control cycle. Under the action of the corresponding calibrated wind field parameters, the wind load-airframe-execution synchronization response sequence of the next control cycle is re-solved, and the wind load execution response state of the next control cycle is updated according to the wind load-airframe execution synchronization response sequence.

10. A multi-physical domain wind disturbance resistance control system for a quadrotor UAV based on digital twins, used to implement the multi-physical domain wind disturbance resistance control method for a quadrotor UAV based on digital twins as described in any one of claims 1-9, characterized in that, Includes the following modules: A multi-physical domain digital twin module is used to simultaneously solve the wind field aerodynamics, airframe motion, and rotor thrust establishment process based on the calibrated wind field parameters and the current feedback state of the UAV, so as to obtain the wind load execution response state. An adaptive fuzzy PID control module is used to adjust the proportional parameter, integral parameter and derivative parameter according to the wind load execution response state, control error and error change rate to obtain the phase-coordinated wind-resistant control quantity. A control allocation module is used to convert the phase-coordinated wind-resistant control quantity into independent control signals for the four rotors according to the X-shaped rotor layout of the quadcopter UAV. An actuator module is used to adjust the thrust of each rotor according to the independent control signal through an electrical drive, motor torque establishment and rotor speed establishment process, so as to complete the position and attitude correction of the UAV; The state perception and feedback module is used to acquire the corrected UAV position, attitude and rotor speed and form a closed-loop feedback state. The closed-loop feedback state is sent back to the multi-physical domain digital twin module to update the wind load execution response state for the next control cycle.

Images