Nonlinear model predictive control method for collaborative transportation and control of suspension load of multi-rotor unmanned aerial vehicle

Abstract

The invention discloses a nonlinear model prediction control method for collaborative transportation and control of a suspension load by a multi-rotor unmanned aerial vehicle, and the method comprises the steps: firstly building a system dynamics model comprising a load, a cable and the rotor unmanned aerial vehicle, and the model is based on SE (3) * S * SO (3) manifold; subsequently, a lightweight NMPC approach is proposed that manipulates the load by optimizing the cable tension distribution while taking into account load dynamics and actuator constraints. The key innovation point is that secondary task constraints, such as minimum distance constraints and obstacle avoidance constraints among unmanned aerial vehicles, are designed in an NMPC framework by utilizing system redundancy, and additional tasks are realized while load control is ensured. In order to reduce the calculation complexity and improve the real-time performance, a state vector parameterization method only containing the six-degree-of-freedom state of the load is adopted.

Description

Technical Field

[0001] The present invention relates to the field of multi-rotor UAV clusters, and in particular to a control method for transporting and manipulating suspended loads in an UAV cluster system. Background Art

[0002] With the rapid development of drone technology, drone swarms have shown great application potential in the coordinated transportation and manipulation of flexible loads. However, there are still many technical challenges in using multiple drones to coordinately transport rigid loads through hanging cables. Traditional control methods often have difficulty in effectively dealing with problems such as indirect drive between loads and drones, nonlinear configuration space, and highly coupled system dynamics. In particular, in the design of control strategies, it is necessary to take into account the precise manipulation of the load and the motion restrictions of the drone, while ensuring the stability and safety of the system.

[0003] In the prior art, some methods regard the load as a disturbance to the UAV formation and compensate it through a disturbance estimator, but this method can usually only achieve basic transportation of the load in three-dimensional space, and cannot achieve precise manipulation of six degrees of freedom (6DoF). Other methods try to model the coupled system dynamics of the load and the UAV and propose corresponding controllers, but these methods often have high computational complexity and are difficult to apply in real time. Summary of the invention

[0004] In order to solve the above technical problems existing in the prior art, the present invention aims to provide a nonlinear model predictive control method for cooperative transportation and manipulation of suspended loads by a multi-rotor UAV. The specific technical solution is as follows:

[0005] A nonlinear model predictive control method for cooperatively transporting and manipulating a suspended load of a multi-rotor UAV comprises the following steps:

[0006] Step 1: System modeling: define the system coordinate system, including the inertial coordinate system, the load coordinate system and the coordinate system of each drone; establish the dynamic model of the load, considering the position, velocity, acceleration and attitude of the load, the attitude includes direction and angular velocity; analyze the effect of cable tension on the load, establish the relationship between cable tension and the force and torque of the load; establish the dynamic model of each drone;

[0007] Step 2: System redundancy analysis: When the number of drones is greater than or equal to 3, analyze the redundancy in the system; determine the cable tension distribution matrix P and analyze its null space to perform additional tasks using redundancy;

[0008] Step 3: Nonlinear Model Predictive Control NMPC Design: Define the state vector and input vector of NMPC. The state vector includes the six-degree-of-freedom state of the load, and the input vector includes the expected force and torque of the load and the redundant coefficient vector; design the objective function of NMPC, including the weighted sum of state error and input error, to minimize tracking error and control energy; establish the discrete form of the system dynamics equation; define the constraints of the optimization problem, including load dynamics constraints, actuator constraints, minimum distance constraints between drones, and obstacle avoidance constraints;

[0009] Step 4: Solving the optimization problem: Using the sequential quadratic programming (SQP) method to solve the NMPC optimization problem in real-time iterative RTI; using the solver to solve the quadratic programming subproblem;

[0010] Step 5: UAV controller design: Calculate the required cable tension for each UAV based on the optimal control input output by the NMPC; Design the UAV controller to convert the cable tension into thrust and torque commands for the UAV; Implement attitude and position control of the UAV to execute the required thrust and torque commands.

[0011] The nonlinear model predictive control method for cooperatively transporting and manipulating a suspended load of a multi-rotor UAV of the present invention has the following advantages:

[0012] 1) Six degrees of freedom precision control and versatility

[0013] This inventive method enables multiple UAVs to precisely manipulate rigid loads within six degrees of freedom through cables, while utilizing system redundancy to perform additional tasks such as obstacle avoidance and separation between UAVs, significantly improving the versatility and flexibility of the system.

[0014] 2) Real-time and scalability optimization

[0015] The method of the present invention reduces the computational complexity and improves the real-time performance by adopting lightweight state vector parameterization. The state vector only contains the six-degree-of-freedom state of the payload, which simplifies trajectory planning and makes the system easy to expand to larger-scale drone clusters.

[0016] 3) Advantages of nonlinear model predictive control

[0017] The NMPC method is applied to handle complex nonlinear system dynamics and constraints, optimizing the coordination of load manipulation and UAV motion. Compared with linearized MPC, it can show better performance in highly dynamic tasks.

[0018] 4) Wide application potential

[0019] This invention method is not only applicable to the fields of construction, logistics, etc., but can also be extended to multiple fields such as disaster relief and environmental monitoring. It has broad application prospects and huge market potential. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 It is a general flow chart of the nonlinear model predictive control method for cooperative transportation and manipulation of suspended loads of a multi-rotor UAV of the present invention;

[0021] Figure 2 It is the system convention definition diagram;

[0022] Figure 3 Power source wrench W L The cable tension vector in the null space of the cable tension distribution matrix P in the payload reference frame L is and tension vector Visualization of the relationship between

[0023] Figure 4 is the control block diagram of this method;

[0024] Figure 5 This is a schematic diagram of the system transportation of multiple UAVs under the control of this method;

[0025] Figure 6 This is a schematic diagram of multiple UAVs coordinating transportation, adjusting posture and avoiding obstacles under this method. DETAILED DESCRIPTION

[0026] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0027] like Figure 1-6 As shown, the nonlinear model predictive control method for cooperatively transporting and manipulating a suspended load of a multi-rotor UAV of the present invention comprises the following steps:

[0028] Step 1: System modeling. Define the system coordinate system, including the inertial coordinate system, the load coordinate system, and the coordinate system of each drone; establish the dynamic model of the load, considering the position, velocity, acceleration, and attitude (direction and angular velocity) of the load; analyze the effect of cable tension on the load, and establish the relationship between cable tension and the force and torque of the load; establish the dynamic model of each drone, considering the position, velocity, acceleration, and attitude control of the drone;

[0029] Step 2: System redundancy analysis. When the number of drones is greater than or equal to 3, analyze the redundancy in the system; determine the cable tension distribution matrix P and analyze its null space to utilize the redundancy to perform additional tasks;

[0030] Step 3: Design of nonlinear model predictive control (NMPC) method. Define the state vector and input vector of NMPC. The state vector includes the six-degree-of-freedom state of the load, and the input vector includes the expected force and torque of the load and the redundant coefficient vector; design the objective function of NMPC, including the weighted sum of state error and input error to minimize tracking error and control energy; establish the discrete form of the system dynamics equation for the optimization problem of NMPC; define the constraints of the optimization problem, including load dynamics constraints, actuator constraints, minimum distance constraints between drones, and obstacle avoidance constraints;

[0031] Step 4: Solve the optimization problem. Use the sequential quadratic programming (SQP) method to solve the NMPC optimization problem in a real-time iteration (RTI) scheme; use an appropriate solver (such as HPIPM) to solve the quadratic programming subproblem;

[0032] Step 5: UAV controller design. Calculate the cable tension required for each UAV based on the optimal control input output by the NMPC; design a UAV controller to convert the cable tension into thrust and torque commands for the UAV; implement attitude and position control of the UAV to execute the required thrust and torque commands.

[0033] In step 1, the UAV dynamics model in system modeling must meet the following assumptions:

[0034] ①The load and drag of the quadrotor can be ignored;

[0035] ②The massless cable is attached to the center of mass of the robot;

[0036] ③ Ignore the aerodynamic effect between the robot and the payload.

[0037] The load state vector, load velocity vector and their corresponding derivatives are defined as follows:

[0038]

[0039] where x L , The position, velocity and acceleration of the payload relative to the inertial reference frame I; q L is the quaternion representation of the direction of the payload relative to I; Ω L , is the angular velocity and angular acceleration of the payload relative to the load coordinate system L. L By taking the derivative of both sides of the equation, the kinematic model of the payload can be obtained as follows:

[0040]

[0041] According to the classical rigid body dynamics and load dynamics equations, it can be concluded that:

[0042]

[0043] In the above formula, F is the sum of the tension of all cables in I, and M is the total torque generated by the cable tension in I on the load, which can be further expressed in matrix form as follows:

[0044]

[0045] In the above formula, μ k is the cable tension of the kth UAV, ρ k is the location of the kth UAV connection point in I, R L is the direction of the load reference system relative to the inertial reference system. The motion equation of the load is expressed in matrix form as follows:

[0046]

[0047]

[0048] Where W∈R 6 represents the power source wrench on the payload, G∈R 6 A vector representing additional physical terms such as gravity and angular momentum effects, represents the general inertia matrix of the system.

[0049] The k-th drone dynamics are modeled, and the dynamics are constrained on a sphere centered at the corresponding attachment point k. Assuming that the cable is always taut, the position of the k-th drone corresponding to the attachment point k is equivalent to the position of I using the vector p k ∈R 3 Indicates that l k is the length of the cable attached to the kth drone, x k is the position of the kth UAV, which can be derived from the geometric position constraint:

[0050] x k =p k -l k ξ k

[0051] The dynamic model of the drone with a hanging load can be derived by the Newton-Euler method or the Lagrangian dynamic equation, and the final form of the two is equivalent. The Newton-Euler method mainly uses Newtonian mechanics for analysis, and the Lagrangian function L defines the kinetic energy E of a mechanical system. k and potential energy E qThe difference between the two is defined as the total external force on the mechanical system is F, which satisfies the following formula:

[0052]

[0053] In the above formula, q can be defined as a displacement variable or an angle variable, and can be defined as a single drone system. The difference between its kinetic energy and potential energy defines the Lagrangian function L, and q is defined as the position x of the drone. k , the total external force is the total thrust of the drone minus the tension of the rope and the gravity of the drone itself. The derivation of Lagrangian dynamics analysis is more complicated, so the Newton-Euler method is used. For the previous geometric constraints, it can be expanded to:

[0054] x k =x L +R L ρ k -l k ξ k

[0055] l k is the length of the rigid body cable connecting the kth UAV, ξ k is the unit vector from the kth drone to its rigid body connection point in I.

[0056] For a straight rigid cable, the two ends are connected to a rigid body, and the connection points at the two ends are a and b respectively. The displacement relationship, velocity relationship, and acceleration relationship are as follows:

[0057] x a =x b +lξ

[0058]

[0059]

[0060] Among them, x a and x b represents the spatial position of a and b, l and ξ represent the cable length and cable unit direction vector respectively. For this UAV system, the following equation is satisfied:

[0061] x k =x L +R L ρ k -l k ξ k

[0062]

[0063]

[0064] One of them satisfies the following equation:

[0065]

[0066]

[0067] in Expressed as Ω L The antisymmetric matrix or Lie algebra element of .

[0068] Secondly, the following equation is satisfied:

[0069] ξ k ·ξ k =1

[0070]

[0071]

[0072]

[0073] The unit direction vector angular velocity is w k , we can get the following formula:

[0074]

[0075]

[0076] where u k is the input force of the kth UAV, a k is the acceleration of the attachment point k equivalent to I.

[0077] Step 2 System redundancy analysis shows that the system has redundancy when the number of UAVs n is greater than or equal to 3. This method provides how to use the extra redundancy in controller design to distribute secondary tasks (such as separation and obstacle avoidance between UAVs).

[0078] In order to facilitate subsequent analysis and controller design, all forces are converted to the load reference system L for representation. The total force F on the load and the tensile force μ k Transferred to the load coordinate system L, they are expressed as F L and

[0079]

[0080] Substituting the above formula into the matrix form of the motion equation of the load, the expression in the load reference system can be obtained as follows:

[0081]

[0082] where P∈R 6×3nIt means mapping the power source wrench vectors of all n-grid quadrotor drones to the load reference system L, and the specific form of the matrix of P can be derived as follows:

[0083]

[0084] Observe that all its sub-blocks are constant matrices, because I 3×3 is the identity matrix, is the vector ρ in the load reference system L k are partially symmetric matrices, all of which are fixed vectors in L. Second, P exhibits a constant number of rows, independent of the number of drones in the system. For n ≥ 3 drones, the number of columns in P exceeds the number of rows, which means that the dimension of the domain of P exceeds the dimension of its image, so according to the theory of linear algebra, the matrix P must have a non-trivial null space. According to the Rank-Nullity Theorem, for any m×n matrix A, its rank (rank(A)) plus the dimension of its null space (nullity(A)) equals the number of columns of A (n). Therefore, there is an extra zero in the matrix P that the system can use to perform secondary tasks, such as obstacle avoidance or separation between robots. The dimension of the matrix P is 6×3n, meaning it has 6 rows and 3n columns. In this case, the number of rows in P (6) is much smaller than the number of columns (3n). According to the Rank-Nullity Theorem, the rank of P defaults to 6 (the general case) and the dimension of its null space is 3n-6.

[0085] set up is an arbitrary vector in the null space of the matrix P, defining the basis matrix N∈R of the null space of the matrix P 3n ×(3n-6) , define Λ∈R 3n-6 is an arbitrary coordinate vector about the basis matrix P. Since the cable mechanism's attachment to the matrix not only determines that the matrix P remains unchanged, we can obtain:

[0086]

[0087] Substituting the above formula into the matrix form of the load motion equation, we can get:

[0088] W L =P(μ L +NΛ)

[0089] By increasing the tension vector μ in the null space of the matrix P with an additional vector L , different tension distributions can be obtained, resulting in equivalent control wrench forces on the payload. This highlights the potential to achieve the desired power source wrench W on the payload by appropriately designing μ and Λ and modifying the shape of the UAV handling payload to meet secondary mission requirements, enabling additional tasks such as UAV separation and obstacle avoidance.

[0090] In step 3, NMPC is used to manipulate the payload’s attitude. NMPC finds a sequence of states {X0, X1, … X N} and input {U0,U1,…U N-1}, while its optimization objective function respects nonlinear constraints and system dynamics. The objective function is composed of the operating cost h(X,U) and the terminal cost h N (X) is composed of and represented by:

[0091]

[0092] satisfy:

[0093] X0=X(t0),

[0094] g(X i ,U i )≤0,

[0095] In the above formula, X i+1 =f(X i ,U i ) is the nonlinear system dynamics in discrete form, g(X,U) is the additional state and input constraints. The optimization takes place under the initial condition X0.

[0096] According to step 1, the state vector and input vector of NMPC are defined as follows:

[0097]

[0098] The goal is to manipulate the load to track the desired trajectory in SE (3). By analogy with the linear quadratic function (LQR) in optimal control theory, the following objective function is selected:

[0099]

[0100] The above formula The state and power source wrench error are determined by comparing the predicted state and trajectory reference inputs as follows:

[0101]

[0102] According to step 1, once the desired payload speed and acceleration are obtained, the desired power source wrench F can be derived. des ,M des As shown below:

[0103]

[0104] The variable dimension of NMPC in this method is 13+3n, which can utilize the null space in step 2 to solve secondary tasks such as obstacle avoidance and spatial separation between drones while satisfying nonlinear system dynamics and actuator constraints.

[0105] The dynamic equation of state X can be obtained as follows:

[0106]

[0107] In order to incorporate the above dynamic equations into the discrete time formulation, the fourth-order Runge-Kutta method is applied. i and input U i In the case of , the sampling time dt is numerically integrated as follows:

[0108] X i+1 =f RK4 (X i ,U i ,dt)

[0109] For the constraint problem, according to the content described in step 2, the required power source wrench can be obtained by appropriately designing μ and Λ to meet other potential tasks. This method proposes the following design:

[0110]

[0111] In the above formula, P ↓ is the pseudo-inverse of P, the cable tension distribution P ↓ W L It is the minimum norm solution. By minimizing the square norm of Λ, the total norm of the cable rope tension vector can be minimized, saving energy for the drone. Since the cable rope tension is along the cable direction, it can be derived that:

[0112]

[0113] in, is the cable direction vector, the position of the kth UAV in the load reference system L, and the robot position constraint can be constructed as follows:

[0114]

[0115] where d r is the minimum distance between selectable drones.

[0116] For obstacle avoidance constraints, let the position of the obstacle in I be p O The following constraints can be defined in NMPC to ensure that the drone group and the payload maintain a safe distance d from obstacles in I Or ,d OL .

[0117]

[0118]

[0119] For actuator constraints, the tension norm can be further constrained to provide some boundaries for the actuator:

[0120]

[0121] In step 4, real-time iteration (RTI) solves the NMPC problem proposed in step 3 through sequential quadratic programming (SQP), and HPIPM is selected as the quadratic programming solver of SQP.

[0122] Step 5: For the optimization problem described in step 3, select the initial moment input To obtain the required tension for the kth drone in the input time series:

[0123]

[0124] The above formula can be expressed as follows: once the desired tension μ is obtained des , then the tension input of each cable of a single quadrotor drone is selected as the projection of the required tension on the corresponding cable rope:

[0125]

[0126] Required direction k,des and the angular velocity w required for the cable rope of the kth UAV k,des It can be expressed as:

[0127]

[0128] in Yes k,des The derivative of the direction vector represents the rate of change of the direction vector over time, which provides information about how the direction of the cable changes over time. Physically, this can be understood as the rate and direction of rotation of the cable in space. The required angular velocity w k,des In order to change the direction of the cable from the current state to the desired state, the quadcopter needs to generate an angular velocity. According to the content in step 1, the thrust f acting on the kth drone is k and moment M k It can be expressed as:

[0129]

[0130]

[0131] Among them, K R and K Ωis a diagonal positive constant matrix, and are the orientation and angular velocity errors defined as follows:

[0132]

[0133]

[0134] It can be seen that the input force vector of the drone is decomposed into the thrust component (z axis) and horizontal moment component (xy plane) It is also clearly stated that only the thrust component needs to be considered when calculating the total thrust. This decomposition is crucial for controlling the flight of a UAV because it allows the thrust and attitude of the UAV to be controlled independently. k and moment M k Expression analysis, its attitude and angular velocity control, K R e Rk and K Ω e Ωk They are used to control the attitude and angular velocity of the drone to make them close to the expected values; Ω k ×J k Ω k Compensates for the force generated by the drone's own rotation; when the drone is rotating, in order to track the desired angular velocity, in addition to directly generating a torque related to the desired angular velocity difference, it is also necessary to consider the additional torque generated by the current rotation. This is to compensate for the additional torque generated by the current rotation. The transformation of the desired angular acceleration represents the required torque calculated based on the desired angular acceleration. Both the transformation of the desired torque and the transformation of the desired angular acceleration are used to convert the desired attitude and angular velocity trajectory into the control torque in the current attitude coordinate system so that the drone can accurately track the trajectory.

[0135] According to step 1, for input and For example, its specific form is as follows:

[0136]

[0137]

[0138] Represents the input force component used to control the drone's attitude. This component is related to the drone's attitude error and angular velocity error, and also needs to take into account cable tension and other dynamic effects. The force and torque generated by angular velocity is the acceleration generated by the projection of the angular velocity in the direction perpendicular to the cable; the acceleration a at the connection point k,c It will also generate forces and moments on the drone. In particular, when the acceleration direction is not parallel to the cable, a force component perpendicular to the cable will be generated. When the desired angular velocity changes, a force component perpendicular to the cable is also generated due to the inertia of the cable. This force can be expressed as Expressed as feedback control terms for cable direction and angular velocity errors.

[0139] is the thrust component that the drone needs to apply along the cable direction. μ k is the direct tension on the cable, which must be offset by the drone in the direction of the cable to keep the system balanced; when the cable rotates around its connection point, a centrifugal force is generated. This centrifugal force is related to the length of the cable, the angular velocity, and the mass of the drone, and can be expressed as The force caused by the acceleration at the cable connection point can be expressed as

[0140] Using the derivation of the connection point acceleration described in step 1, we can get the acceleration a k,c The representation is as follows:

[0141]

[0142] in and is a diagonal positive constant matrix, and are the cable direction and cable angular velocity errors respectively, satisfying the following formula:

[0143]

[0144]

[0145] Table 1 Parameter meaning table.

[0146] The specific embodiments are as follows:

[0147] 1) Scenario description:

[0148] On a construction site, a heavy glass panel needs to be transported and precisely manipulated to be installed on the exterior wall of a high-rise building. Due to the heavy weight and special shape of the glass panel, traditional manual or ground transportation methods are difficult to accomplish. Therefore, multiple quadcopters can be used to cooperate to complete this task by hanging cables.

[0149] 2) System preparation:

[0150] Drone and payload configuration: Choose three or more quadcopters, each equipped with a lightweight cable of sufficient length. The cable ends are connected to a special hanging frame that can firmly fix the glass panel.

[0151] State estimation and positioning: Install necessary sensors (such as IMU, GPS and camera) on the drone and glass panel. Figure 4 vicon in the image processing unit (the same as the positioning function in the image processing unit), so as to estimate their position, velocity and attitude in real time.

[0152] Controller deployment: An advanced computing unit (such as an embedded system or a high-performance computer) is deployed on the UAV to run the nonlinear model predictive control (NMPC) algorithm proposed in this method.

[0153] 3) Model building:

[0154] Load dynamics modeling: According to this method, the dynamics of the glass sheet (load) are modeled, including its position, velocity and acceleration in three-dimensional space, as well as rotational motion and moments.

[0155] System redundancy analysis: Since three or more UAVs are used, there is redundancy in the system. According to the analysis in this method, this redundancy can be used to perform additional tasks such as spatial separation between UAVs and obstacle avoidance.

[0156] Drone Dynamics Modeling: Model the dynamics of each drone, considering its motion characteristics under the constraints of the suspension cables.

[0157] 4) NMPC algorithm implementation

[0158] State vector and input vector definition: Define the system's state vector (including the load's position, velocity, posture, etc.) and input vector (including cable tension, etc.).

[0159] Cost Function Design: Design a cost function to minimize the error between the payload position and the desired trajectory, taking into account the system dynamics and constraints.

[0160] System dynamics and constraints: The dynamic equations of the payload and the UAV are incorporated into the NMPC framework, and necessary constraints are added, such as the minimum distance between UAVs, obstacle avoidance requirements, and cable tension limits.

[0161] Optimization solution: The SQP method mentioned in this method is used to solve the NMPC problem under the real-time iteration (RTI) scheme.

[0162] 5) Task Execution

[0163] Trajectory planning: Plan a desired trajectory for the load from the starting point to the end point, which should take into account the specific layout of the construction site and the distribution of obstacles.

[0164] Control command generation: The NMPC algorithm generates the optimal control command (i.e., cable tension distribution) based on the current state and desired trajectory to control the load to move along the desired trajectory.

[0165] Real-time adjustment: During the mission execution, the NMPC algorithm will continuously adjust according to the new state feedback to ensure the precise control of the payload and the safe flight of the drone. The implementation effect is visible Figure 5 , the collaborative obstacle avoidance effect is visible Figure 6 .

[0166] 6) Mission Completed

[0167] Precise positioning of loads: When the load is successfully transported to the designated location on the exterior wall of a high-rise building, the drone team will accurately adjust the load's posture and position to facilitate installation workers to fix and install it.

[0168] Drone Recovery: Once the installation is complete, the drone team releases the cables and returns safely to base or the next mission point.

[0169] Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A nonlinear model predictive control method for cooperative transportation and manipulation of suspended loads by a multi-rotor UAV, comprising the following steps: Step 1: System modeling: define the system coordinate system, including the inertial coordinate system, the load coordinate system and the coordinate system of each drone; establish the dynamic model of the load, considering the position, velocity, acceleration and attitude of the load, the attitude includes direction and angular velocity; analyze the effect of cable tension on the load, establish the relationship between cable tension and the force and torque of the load; establish the dynamic model of each drone; Step 2: System redundancy analysis: When the number of drones is greater than or equal to 3, analyze the redundancy in the system; Determine the cable tension distribution matrix P and analyze its null space to exploit redundancy to perform additional tasks; Step 3: Nonlinear model predictive control NMPC design: Define the state vector and input vector of NMPC, where the state vector includes the six-degree-of-freedom state of the load, and the input vector includes the expected force and torque of the load and the redundant coefficient vector; design the objective function of NMPC, including the weighted sum of state error and input error, to minimize tracking error and control energy; Establish the discrete form of the system dynamics equation; define the constraints of the optimization problem, including load dynamics constraints, actuator constraints, minimum distance constraints between drones, and obstacle avoidance constraints; Step 4: Solving the optimization problem: Using the sequential quadratic programming (SQP) method to solve the NMPC optimization problem in real-time iterative RTI; Use solvers to solve quadratic programming subproblems; Step 5: UAV controller design: Calculate the required cable tension for each UAV based on the optimal control input output by NMPC; Design a UAV controller to convert cable tension into thrust and torque commands for the UAV; implement attitude and position control of the UAV to execute the required thrust and torque commands.

2. The nonlinear model predictive control method for cooperative transportation and manipulation of suspended loads by a multi-rotor UAV according to claim 1 is characterized in that: The UAV dynamics model in step 1 must meet the following assumptions: ①The load and drag of the quadrotor can be ignored; ②The massless cable is attached to the center of mass of the robot; ③ Ignore the aerodynamic effect between the robot and the payload. The load state vector, load velocity vector and their corresponding derivatives are defined as follows: where x L , The position, velocity and acceleration of the payload relative to I; q L is the quaternion representation of the direction of the payload relative to I; Ω L , is the angular velocity and angular acceleration of the payload relative to the load coordinate system L; L Taking the derivative of both sides of the equation, the kinematic model of the payload is as follows: According to the classical rigid body dynamics and load dynamics equations, it can be concluded that: Where F is the sum of the tension of all cables in I, and M is the total torque generated by the cable tension in I on the load, which can be further expressed in matrix form as follows: Among them, μ k is the cable tension of the kth UAV; the motion equation of the load is expressed in matrix form as follows: Where W∈R 6 represents the power source wrench on the payload, G∈R 6 A vector representing additional physical terms such as gravity and angular momentum effects, represents the general inertia matrix of the system; The k-th UAV dynamics are modeled, and the dynamics are constrained on a sphere centered on the corresponding attachment point k; the position of the k-th UAV corresponding to the attachment point k is equivalent to the position of I using the vector p k ∈R 3 Indicates that l k is the length of the cable tied to the kth UAV, which can be derived from the geometric position constraint: x k =p k -l k x k By deriving the dynamic model of the UAV with hanging load, the following equation is satisfied: In the above formula, q can be defined as a displacement variable or an angle variable. To define a single UAV system, the difference between its kinetic energy and potential energy defines the Lagrangian function L, and q is defined as the position x of the UAV. k , the total external force is the total thrust of the drone minus the tension of the rope and the gravity of the drone itself; for the previous geometric constraints, it can be expanded to: x k =x L +R L r k -l k x k After derivation, it satisfies the following formula: Which satisfies the following formula: in Expressed as Ω L The antisymmetric matrix or Lie algebra element of ; Secondly, the following equation is satisfied: The unit direction vector angular velocity is w k , we can get the following formula: where u k is the input force of the kth UAV, a k is the acceleration of the attachment point k equivalent to I.

3. The nonlinear model predictive control method for cooperative transportation and manipulation of suspended loads by a multi-rotor UAV according to claim 2 is characterized in that: In step 2, all forces are converted to the load reference system L for expression; the total force F and the tensile force μ on the load are k Transferred to the load coordinate system L, they are expressed as F L and Substituting the above formula into the matrix form of the motion equation of the load, the expression in the load reference system can be obtained as follows: where P∈R 6×3n It means mapping the power source wrench vectors of all n-grid quadrotor drones to the load reference system L, and the specific form of the matrix of P can be derived as follows: All sub-blocks of P are constant matrices, I 3×3 is the identity matrix, is the vector ρ in the load reference system L k The partially symmetric matrices are all fixed vectors in L; set up is an arbitrary vector in the null space of the matrix P, defining the basis matrix N∈R of the null space of the matrix P 3n×(3n-6) , define Λ∈R 3n-6 For any coordinate vector about the basis matrix P, we can get: Substituting the above formula into the matrix form of the load motion equation, we can get: IN L =P(μ L +NΛ) By increasing the tension vector μ in the null space of the matrix P with an additional vector L , different tension distributions are obtained, resulting in equivalent controlled wrench power on the payload.

4. The nonlinear model predictive control method for cooperative transportation and manipulation of suspended loads by a multi-rotor UAV according to claim 3 is characterized by: In step 3, NMPC is used to manipulate the payload’s attitude; NMPC finds a sequence of states {X0, X1, …X N } and input {U0,U1,…U N-1 }, and optimizes the objective function, which is composed of the operating cost h(X,U) and the terminal cost h N (X) is composed of and represented by: satisfy: X0=X(t0), g(X i ,U i )≤0, In the above formula, X i+1 =f(X i ,U i ) is the nonlinear system dynamics in discrete form, g(X,U) is the additional state and input constraints, and the optimization occurs under the initial condition X0; According to step 1, the state vector and input vector of NMPC are defined as follows: Select the following objective function: The above formula The state and power source wrench error are determined by comparing the predicted state and trajectory reference inputs as follows: According to step 1, once the desired payload speed and acceleration are obtained, the desired power source wrench F can be derived. des ,M des As shown below: The dynamic equation of state X can be obtained as follows: In a given state X i and input U i In the case of , the sampling time dt is numerically integrated as follows: X i+1 =f RK4 (X i ,U i ,dt) For the constraint problem, according to step 2, the required power source wrench is obtained by appropriately designing μ and Λ to meet other potential tasks; The design is as follows: In the above formula, P ↓ is the pseudo-inverse of P, the cable tension distribution P ↓ W L is the minimum norm solution, which minimizes the total norm of the cable rope tension vector by minimizing the square norm of Λ. Since the cable rope tension is along the cable direction, it is derived that: in, is the cable direction vector, the position of the kth UAV in the load reference system L, and the robot position constraint is constructed as follows: where d r is the minimum distance between selectable drones; For obstacle avoidance constraints, let the position of the obstacle in I be p O ; Define the following constraints in NMPC so that the UAV group and the payload maintain a safe distance d from obstacles in I Or ,d OL : For the actuator constraints, the tension norm is further restricted to provide boundaries for the actuator:

5. The nonlinear model predictive control method for cooperative transportation and manipulation of suspended loads by a multi-rotor UAV according to claim 4 is characterized in that: In step 4, HPIPM is selected as the quadratic programming solver for SQP.

6. The nonlinear model predictive control method for cooperative transportation and manipulation of suspended loads by a multi-rotor UAV according to claim 5, characterized in that: Step 5: For the optimization problem described in step 3, select the initial moment input To obtain the required tension for the kth drone in the input time series: The above formula shows that once the desired tension μ is obtained des , then the tension input of each cable of a single quadrotor drone is selected as the projection of the required tension on the corresponding cable rope: Required direction k,des and the angular velocity w required for the cable rope of the kth UAV k,des It can be expressed as: in Yes k,des The derivative of represents the rate of change of the direction vector over time, which provides information about how the cable direction changes over time; Required angular velocity w k,des represents the angular velocity that the quadcopter needs to generate in order to change the direction of the cable from the current state to the desired state. According to step 1, the thrust f acting on the kth UAV k and moment M k It can be expressed as: Among them, K R and K Ω is a diagonal positive constant matrix, and are the orientation and angular velocity errors defined as follows: According to step 1, for the input and For example, its specific form is as follows: represents the input force component used to control the attitude of the drone, is the thrust component that the drone needs to apply along the cable direction. Using the derivation of the connection point acceleration described in step 1, we can get the acceleration a k,c The representation is as follows: in and is a diagonal positive constant matrix, and are the cable direction and cable angular velocity errors respectively, satisfying the following formula:

Images