Self-adaptive penalty parameter adjustment criterion based on spectral gradient method

Abstract

The invention discloses a penalty parameter adaptive adjustment criterion based on a spectral gradient method, and the method comprises the steps: setting an adaptive penalty parameter for each dual variable in the process of solving an optimal control model through building the optimal control model, and obtaining the spatial-temporal trajectory of an unmanned plane cluster. According to the method disclosed by the invention, under different task scenes, the number of iterations of algorithm solution can be reduced, and the efficiency is improved.

Description

Technical Field

[0001] The invention relates to a penalty parameter adaptive adjustment criterion based on a spectral gradient method, and belongs to the technical field of aircraft control. Background Art

[0002] Trajectory optimization plays a vital role in enhancing the autonomy of the swarm. Trajectory optimization can be formulated as a multi-agent optimal control (MAOC) problem, which is to optimize the control history of each agent to minimize the global general performance indicator while complying with multiple nonlinear constraints.

[0003] Most MAOC problems are handled using the alternating direction method of multipliers (ADMM). Due to the strong non-convexity of the drone cluster trajectory optimization problem, coupled with the influence of dynamic constraints, collision constraints, communication distance maintenance and other factors, ADMM is difficult to directly apply to actual tasks. Therefore, in the prior art, drone cluster trajectory optimization generally uses ADMM-based distributed computing and DDP-based local planning framework to solve the distributed spatiotemporal trajectory optimization problem of drone clusters. Compared with the sequential convex programming method and the quadratic programming method, the fully distributed characteristics of this method show more excellent scalability.

[0004] However, this method requires setting the terminal time for each UAV in advance, and also needs to consider two additional issues: time optimization and convergence acceleration. Among them, the convergence performance depends largely on the penalty parameter. The existing penalty parameter generally adopts residual balance (RB), which is based on the idea of ​​balancing the original residual and the dual residual to make them approximately equal, and improves the convergence by increasing or decreasing the penalty parameter in the iteration. However, the performance of the RB method varies greatly with the task, and the convergence speed is slow.

[0005] There is also a method in the prior art of accelerating the convergence speed by using the Nesterov acceleration method, but this method is not effective in most cases.

[0006] Therefore, it is necessary to further study the existing UAV cluster trajectory optimization methods to solve the above problems. Summary of the invention

[0007] In order to overcome the above problems, the inventors conducted in-depth research and proposed a penalty parameter adaptive adjustment criterion based on the spectral gradient method, including:

[0008] Establish optimal control model;

[0009] Solve the optimal control model to obtain the spatiotemporal trajectory of the drone cluster;

[0010] The UAV flies according to the obtained space-time trajectory.

[0011] In a preferred embodiment, the optimal control model is expressed as:

[0012]

[0013] Among them, i represents the number of different drones in the drone cluster, N represents the total time, represents the minimum value of the control input of UAV i at all times, represents the final flight time of drone i, t i,N represents the flight time of drone i, X i Represents the stack of state vectors of drone i at all times, U i represents the stack of control inputs of drone i at all times, J i (X i ,U i ) represents the local objective function of the i-th UAV, st represents the constraint condition, x i,k represents the state vector of the i-th UAV at time k, u i,k represents the control vector of the ith UAV at time k, f i () represents the discrete system dynamics with nonlinear characteristics, For x i,k A safe copy variable, a i,k () is a nonlinear function used to describe path constraints, for u i,k A safe copy variable, b i,k () represents the control torque constraint function, t i,N A safe copy variable, c i,N () represents a nonlinear function, represents the set of state vectors of all neighboring drones of drone i at time k, represents the set of state constants between UAV i and its neighboring UAVs at time k, represents the set of time variables of all neighboring drones of drone i at time k, represents the set of time constraints between UAV i and its neighboring UAVs at time k, represents the set of state variables that have reached consensus, A collection of variables representing the time at which consensus is reached.

[0014] In a preferred embodiment, the spatiotemporal trajectory acquisition process of the drone cluster includes the following steps:

[0015] S21. At time k, each drone obtains its own state vector x i,k , control vector u i,kand flight time t i,n ;

[0016] S22. At time k, each drone obtains its own secure copy variable

[0017] S23, each drone communicates with all its neighbor drones to obtain the state variables of drone i observed by neighbor drone j and the flight time of drone i observed by neighboring drone j

[0018] S24. Each drone calculates the consensus state variable z based on its own variables and the variables obtained from its neighbors. i,k and the time variable s to reach consensus i,N ;

[0019] S25. Each drone communicates with all its neighbor drones to obtain the state variable z calculated by the neighbor drone j of the drone for reaching consensus. j,k and the time variable s used to reach consensus j,N ;

[0020] S26, each drone uses the state variable z that it has obtained to reach consensus i,k and the time variable s to reach consensus i,N and the consensus state variable z of neighbor drone j j,k and the time variable s to reach consensus j,N Based on the basic data, the ADMM method is used to solve the optimal control variables at time k, obtain the dual variables at time k, and obtain the trajectory of the UAV at time k by minimizing the residual;

[0021] S27. Repeat steps S21-S26 to obtain the trajectory of each drone at the next moment in turn.

[0022] In a preferred embodiment, during the process of acquiring the spatiotemporal trajectory of the drone cluster, a penalty parameter is set for each dual variable to accelerate the convergence speed of the solution.

[0023] In a preferred embodiment, the penalty parameter is used to accelerate the convergence process by the following method:

[0024] Update the dual variables;

[0025] Get the iterative spectrum step size;

[0026] Get the penalty parameter for the next iteration based on the spectrum step size.

[0027] In a preferred embodiment, the updating process of the dual variables is expressed as:

[0028]

[0029]

[0030]

[0031]

[0032]

[0033] The superscript n represents the number of iterations. represents the dual variable of the ith UAV at the nth iteration, is the dual variable at the n+1th iteration, is the dual variable ζ i,k , i,k 、v i,N ,y i,k , η i,N The penalty parameter corresponding to the nth iteration.

[0034] In a preferred embodiment, the iterative spectrum step length acquisition method comprises the following steps:

[0035] S31, obtaining the optimal penalty parameter of iteration;

[0036] S32, obtaining an estimate of iterative curvature;

[0037] S33, obtaining the iterative spectrum step size.

[0038] In a preferred embodiment, in S31, x i,k and Constructed in dual form:

[0039]

[0040] in,

[0041]

[0042]

[0043]

[0044] By minimizing The residual at , obtains the optimal penalty parameter of the iteration, expressed as:

[0045]

[0046] in, For custom functions, is the indicator function, H * () represents the Fenchel conjugate of H(), G * () represents the Fenchel conjugate of G(), Indicates H * The sub-gradient of (), express The sub-gradient of is the intermediate variable, is a constant vector, is the dual function The first-order linear approximation of is the dual function The first-order linear approximation of Indicates the spectrum step size of the nth iteration, the indicator function In the subscript f i 、a i,k Represents a set, and the indicator function is defined as: if the function variable is in the set, the indicator function is 0; if the function variable is not in the set, the indicator function is positive infinity.

[0047] In a preferred embodiment, in S32, the estimation of the iterative curvature is expressed as:

[0048]

[0049]

[0050]

[0051]

[0052]

[0053]

[0054] in,

[0055] Indicates the curvature at the nth iteration The estimated value of Indicates the curvature at the nth iteration The superscript SD represents the steepest descent of the parameter, MG represents the minimum gradient of the parameter, and the superscript n l Represents the result of the last iteration of the parameter.

[0056] In a preferred embodiment, in S33, the iterative spectrum step size is expressed as:

[0057]

[0058]

[0059]

[0060] in, represents the spectrum step size of the n+1th iteration, ∈ cor is the accuracy threshold of the curvature estimate, is the dual function The cosine of the angle between the change in the dual subgradient and the change in the dual variable, is the dual function The cosine of the angle between the change in the dual subgradient and the change in the dual variable.

[0061] The beneficial effects of the present invention include:

[0062] (1) Reduce the number of algorithm iterations and improve efficiency;

[0063] (2) Able to handle different task scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 A schematic diagram of a penalty parameter adaptive adjustment criterion flow chart based on the spectral gradient method according to a preferred embodiment of the present invention is shown. DETAILED DESCRIPTION

[0065] The present invention will be further described in detail below through the accompanying drawings and embodiments. Through these descriptions, the characteristics and advantages of the present invention will become more clear and distinct.

[0066] The word "exemplary" is used exclusively herein to mean "serving as an example, embodiment, or illustration." Any embodiment described herein as "exemplary" is not necessarily to be construed as preferred or advantageous over other embodiments. Although various aspects of the embodiments are shown in the drawings, the drawings are not necessarily drawn to scale unless otherwise noted.

[0067] According to the present invention, a penalty parameter adaptive adjustment criterion based on the spectral gradient method is provided, such as Figure 1 As shown, the following steps are included:

[0068] S1. Establish the optimal control model;

[0069] S2. Solve the optimal control model and obtain the spatiotemporal trajectory of the UAV cluster.

[0070] According to the present invention, after the space-time trajectory of the drone cluster is obtained, the drones fly according to the obtained space-time trajectory, and the flight control of the drones is achieved by tracking the space-time trajectory.

[0071] In S1, the optimal control model is expressed as a MACO problem that can be distributed optimized, which is used to optimize the trajectory of the drone cluster.

[0072] The optimal control model is expressed as:

[0073]

[0074] Among them, i represents the number of different drones in the drone cluster, N represents the total time, represents the minimum value of the control input of UAV i at all times, represents the final flight time of drone i, t i,N represents the flight time of drone i, X i Represents the stack of state vectors of drone i at all times, U i represents the stack of control inputs of drone i at all times, J i (X i ,U i ) represents the local objective function of the i-th UAV, st represents the constraint condition, x i,k represents the state vector of the i-th UAV at time k, u i,k represents the control vector of the ith UAV at time k, f i () represents the discrete system dynamics with nonlinear characteristics, For x i,k A safe copy variable, a i,k () is a nonlinear function used to describe path constraints, for u i,k A safe copy variable, b i,k () represents the control torque constraint function, t i,N A safe copy variable, c i,N () represents a nonlinear function, represents the set of state vectors of all neighboring drones of drone i at time k, represents the set of state constants between drone i and its neighboring drones, represents the set of time variables of all neighboring drones of drone i at time k, represents the set of time constraints between UAV i and its neighboring UAVs, represents the set of state variables that have reached consensus, A collection of variables representing the time at which consensus is reached.

[0075] a i,k (), b i,k (), c i,NThe specific setting form of () is not limited in the present invention, and technical personnel in this field can freely set it according to actual needs and experience.

[0076] are all local constraints of a single UAV. is the state constant between UAVs, expressed as:

[0077]

[0078]

[0079] in, represents the state variable of neighboring drone j observed from drone i, d ij,k () represents the state constraint between drone i and its neighbor drone j, represents the flight time of neighboring drone j observed from drone i, e ij,N () represents the time constraint between UAV i and its neighbor UAV j, represents the neighbor set of UAV i excluding UAV i itself.

[0080] According to the present invention,

[0081] Among them, z j,k The state variable that represents the consensus, s j,N represents the time variable for reaching consensus, represents the neighbor set of drone i.

[0082] According to the present invention, if drone j is within the effective communication range of drone i, then drone j is the neighbor of drone i, and all neighbors of drone i and drone i together constitute the neighbor set of drone i.

[0083] Preferably, the local objective function of the i-th UAV is set as:

[0084]

[0085] Among them, l i () represents a scalar function, φ i () represents terminal cost.

[0086] Further,

[0087]

[0088] Right now

[0089]

[0090] In the formula, represents the expected terminal state vector of the ith UAV, and is the weighted diagonal matrix of the cost function, ([x i ,y i ,θ i ] T is the state variable of a single drone, [x i ,y i ] T represents the coordinates of the i-th UAV, θ i is the flight heading angle.

[0091] S2 includes the following sub-steps:

[0092] S21. At time k, each drone obtains its own state vector x i,k , control vector u i,k and flight time t i,N ;

[0093] S22. At time k, each drone obtains its own secure copy variable

[0094] S23, each drone communicates with all its neighbor drones to obtain the state variables of drone i observed by neighbor drone j and the flight time of drone i observed from neighboring drone j

[0095] S24. Each drone calculates the consensus state variable z based on its own and its neighbors’ variables i,k and the time variable s to reach consensus i,N ; Preferably, the calculation is expressed as:

[0096]

[0097] in, represents the set of drones j that consider drone i as their neighbors, represents the state variable of drone i observed by neighbor drone j, represents the flight time of UAV i observed by neighbor UAV j.

[0098] S25. Each drone communicates with all its neighbor drones to obtain the consensus state variable z calculated by the neighbor drone j of the drone. j,k and the time variable s to reach consensus j,N ;

[0099] S26, each drone uses the state variable z that it has obtained to reach consensusi,k and the time variable s to reach consensus i,N and the consensus state variable z of neighbor drone j j,k and the time variable s to reach consensus j,N Based on the basic data, the ADMM method is used to solve the optimal control variables at time k, obtain the dual variables at time k, and obtain the trajectory of the UAV at time k by minimizing the residual;

[0100] S27. Repeat steps S21-S26 to obtain the trajectory of each drone at the next moment in turn.

[0101] According to the present invention, in S21-S22, the DDP method is used. DDP (differential dynamic programming) is a commonly used method in optimal control methods and will not be described in detail in the present invention.

[0102] In the above process, each UAV adopts the DDP method to solve and realize distributed computing. At the same time, it is combined with the ADMM method, and ADMM is used to perform spatiotemporal constraints between UAVs, thereby realizing the spatiotemporal trajectory optimization of the UAV cluster.

[0103] Preferably, before S21, there is further step S20 of initializing the optimal control model, including:

[0104] S201, each UAV uses the DDP method to obtain its state vector under unconstrained conditions Control vector Flight time

[0105] S202, each drone communicates with all its neighbor drones and receives the state vector of the neighbor drones under unconstrained conditions Unconstrained flight time with neighboring drones

[0106] S203: Setting the initial values ​​of the optimal control model parameters.

[0107] In S203, x i,k The initial value of the parameter is set to The initial value of the parameter is set to t i,N The initial value of the parameter is set to The initial value of the parameter is set to The initial value of the parameter is set to The initial value of the parameter is set to The initial value of the parameter is set to The initial value of the parameter is set to ζ i,k , i,k ,yi,k , ν i,N , η i,N The initial values ​​of the parameters are all set to 0.

[0108] in, represents the state vector of the neighboring drones of drone i under unconstrained conditions, represents the final time of the neighboring drones of drone i under unconstrained conditions, ζ i,k , i,k , ν i,N ,y i,k , η i,N represents the dual variable of the constraints in the augmented Lagrange multiplier formula in the ADMM method.

[0109] Preferably, in S2, a penalty parameter is set for each dual variable to accelerate the convergence speed of the solution.

[0110] Preferably, the initial value of the penalty parameter is set to 0.

[0111] Furthermore, when the solution in S2 does not converge, the penalty parameter is used to accelerate the convergence process by the following method:

[0112] Update the dual variables; obtain the iterative spectrum step size;

[0113] Get the penalty parameter for the next iteration based on the spectrum step size.

[0114] Preferably, the above process is not performed in each iteration, but in each C Freq The iterative step performs the above process once, thereby improving reliability, that is, preferably:

[0115] If the number of iterations n is divided by C Freq The remainder is 1, the dual variable is updated, and the spectral step length of the iteration is obtained. The penalty parameter of the next iteration is obtained according to the spectral step length;

[0116] If the number of iterations n is divided by C Freq If the remainder is not 1, the current penalty parameter is used as the penalty parameter for the next iteration.

[0117] Among them, C Freq is a configurable parameter.

[0118] The updating process of the dual variables is expressed as:

[0119]

[0120]

[0121]

[0122]

[0123]

[0124] The superscript n represents the number of iterations. represents the dual variable of the ith UAV at the nth iteration, is the dual variable at the n+1th iteration, is the dual variable ζ i,k , i,k 、v i,N ,y i,k , η i,N The penalty parameter corresponding to the nth iteration.

[0125] The iterative spectrum step length acquisition method comprises the following steps:

[0126] S31, obtaining the optimal penalty parameter of iteration;

[0127] S32, obtaining an estimate of iterative curvature;

[0128] S33, obtaining the iterative spectrum step size.

[0129] In S31, x i,k and Constructed in dual form:

[0130]

[0131] in,

[0132]

[0133]

[0134]

[0135]

[0136]

[0137]

[0138] By minimizing The residual at , obtains the optimal penalty parameter of the iteration, expressed as:

[0139]

[0140] Among them, H(x i,k ), For custom functions, is the indicator function, subscript fi 、a i,k represents a set, which is defined as follows: if the function variable is in the set, the indicator function is 0; if the function variable is not in the set, the indicator function is positive infinity. * () represents the Fenchel conjugate of H(), G * () represents the Fenchel conjugate of G(), Indicates H * The sub-gradient of (), express The sub-gradient of is the intermediate variable, is a constant vector, is the dual function The first-order linear approximation of is the dual function The first-order linear approximation of Indicates the spectrum step size of the nth iteration.

[0141] In S32, the iterative curvature estimation is expressed as:

[0142]

[0143]

[0144]

[0145]

[0146]

[0147]

[0148] in,

[0149] Indicates the curvature at the nth iteration The estimated value of Indicates the curvature at the nth iteration The superscript SD represents the steepest descent of the parameter, MG represents the minimum gradient of the parameter, and the superscript n l Represents the result of the last iteration of the parameter.

[0150] In the present invention, a hybrid step length calculation rule is adopted. and It is calculated based on the intermediate quantities of different UAVs. In the iteration, this update rule has different results, which makes this method have a significant effect on different UAVs in distributed problems.

[0151] In S33, the iterative spectrum step size is expressed as:

[0152]

[0153]

[0154]

[0155] in, represents the spectrum step size of the n+1th iteration, ∈ cor is the accuracy threshold of the curvature estimate, is the dual function The cosine of the angle between the change in the dual subgradient and the change in the dual variable, is the dual function The cosine of the angle between the change in the dual subgradient and the change in the dual variable.

[0156] The inventors found that the gradient The linear approximation may not be accurate, resulting in an invalid step size. In the present invention, the accuracy of the curvature estimation is evaluated by the cosine value of the angle between the change of the dual subgradient and the change of the dual variable, thereby ensuring the safety of the updated step size. Only when the above evaluation result is lower than a preset threshold, the spectral step size is allowed to be updated.

[0157] More preferably, the spectrum step length is updated and the following restrictions are imposed:

[0158]

[0159] Among them, C cg is the convergence constant.

[0160] The above restrictions ensure the numerical stability of the acceleration.

[0161] According to the present invention, the above-mentioned adaptive penalty parameters are used, which helps to reduce the number of iterations of the algorithm and improve efficiency.

[0162] Various embodiments of the method described above in the present invention can be implemented in digital electronic circuit systems, integrated circuit systems, field programmable gate arrays (FPGAs), application specific integrated circuits (ASICs), application specific standard products (ASSPs), systems on chips (SOCs), load programmable logic devices (CPLDs), computer hardware, firmware, software, and / or combinations thereof. These various embodiments may include: being implemented in one or more computer programs, which may be executed and / or interpreted on a programmable system including at least one programmable processor, which may be a dedicated or general programmable processor, which may receive data and instructions from a storage system, at least one input device, and at least one output device, and transmit data and instructions to the storage system, the at least one input device, and the at least one output device.

[0163] It should be understood that the various forms of processes shown above can be used to reorder, add or delete steps. For example, the steps described in the disclosure of the present invention can be performed in parallel, sequentially or in different orders, as long as the desired results of the technical solution disclosed in the present invention can be achieved, and this document does not limit this.

[0164] Example

[0165] Example 1

[0166] In the simulation experiment, each UAV in the cluster is set to fly at a constant speed in a two-dimensional plane. The kinematic model of each UAV is expressed as:

[0167]

[0168]

[0169]

[0170] [x i ,y i ,θ i ] T is the state variable of a single drone, [x i ,y i ] T represents the coordinates of the i-th UAV, θ i is the flight heading angle. The flight speed of the drone in the horizontal plane is set to a fixed value V. The drone is controlled by adjusting the turning angular velocity of the heading angle. Since the turning angular velocity of the drone is limited by physical constraints, the turning angular velocity ω i The range is set to |ω i |≤ω max , where ω max is a positive constant. In addition,min ≤t N ≤t max Set as the flight time limit for each drone.

[0171] The cluster should approach the desired terminal state while respecting the following different constraints. To this end, the quadratic cost function J for all the drones in the cluster is i (X i ,U i ) is calculated as:

[0172]

[0173] In the formula, represents the expected terminal state vector of the ith UAV, and is the weighted diagonal matrix of the cost function; Path constraint: obstacles in the environment during flight are regarded as circular no-fly zones, that is:

[0174]

[0175] In the formula, [x o ,y o ] T represents the center coordinate of the obstacle area, r o is the radius, and the minimum safe distance o∈O that each agent needs to maintain from the obstacle is denoted by d o , the set of obstacles in the environment is represented as O; constraints between agents: each drone in the cluster needs to maintain a safe distance from neighboring drones to avoid collision. And due to the physical limitations of the communication equipment, it also needs to maintain a certain distance from its neighbors to maintain the smooth communication link. Therefore, these constraints can be described as:

[0176] (x i,k -x j,k ) 2 +(y i,k -y j,k ) 2 ≤d col

[0177] (x i,k -x j,k ) 2 +(y i,k -y j,k ) 2 ≥d con

[0178] Where, d col Indicates the safety distance, d conRepresents the maximum effective communication distance between two drones; Terminal timing constraint: It is assumed that the drones arrive at the target in sequence within the specified time interval. This constraint can be expressed mathematically as:

[0179]

[0180] In the formula, is the time interval vector, the matrix A i Represents the time series relationship between agent i and its neighbors:

[0181]

[0182] In the formula, A i The jth row of represents the temporal relationship between drone i and its neighbor j, corresponding to The j-th row of represents the preset time interval between i and j.

[0183] In the simulation experiment, the UAV adopts the penalty parameter adaptive adjustment criterion based on the spectral gradient method, including:

[0184] Establish optimal control model;

[0185] Solve the optimal control model and obtain the spatiotemporal trajectory of the drone swarm.

[0186] The optimal control model is expressed as:

[0187]

[0188] The spatiotemporal trajectory acquisition process of the drone cluster includes the following steps:

[0189] S21. At time k, each drone obtains its own state vector x i,k , control vector u i,k and flight time t i,N ;

[0190] S22. At time k, each drone obtains its own secure copy variable

[0191] S23, each drone communicates with all its neighbor drones to obtain the state variables of drone i observed by neighbor drone j of neighbor drone and the flight time of drone i observed from neighboring drone j

[0192] S24. Each drone calculates the consensus state variable z based on its own and its neighbors’ variables i,k and the time variable s to reach consensus i,N ;

[0193] S25. Each drone communicates with all its neighbor drones to obtain the consensus state variable z calculated by the neighbor drone j of the drone. j,k and the time variable s to reach consensus j,N ;

[0194] S26, each drone uses the state variable z that it has obtained to reach consensus i,k and the time variable s to reach consensus i,N and the consensus state variable z of neighbor drone j j,k and the time variable s to reach consensus j,N Based on the basic data, the ADMM method is used to solve the optimal control variables at time k, obtain the dual variables at time k, and obtain the trajectory of the UAV at time k by minimizing the residual;

[0195] S27. Repeat steps S21-S26 to obtain the trajectory of each drone at the next moment in turn.

[0196] In the process of acquiring the spatiotemporal trajectory of the drone cluster, a penalty parameter is set for each dual variable to accelerate the convergence speed of the solution. The penalty parameter is used to accelerate the convergence process through the following methods:

[0197] If the number of iterations n is divided by C Freq The remainder is 1, the dual variable is updated, and the spectral step length of the iteration is obtained. The penalty parameter of the next iteration is obtained according to the spectral step length;

[0198] If the number of iterations n is divided by C Freq If the remainder is not 1, the current penalty parameter is used as the penalty parameter for the next iteration, where C Freq Set to 10.

[0199] The updating process of the dual variables is expressed as:

[0200]

[0201]

[0202]

[0203]

[0204]

[0205] The updating process of the dual variables is expressed as:

[0206]

[0207]

[0208]

[0209]

[0210]

[0211] The iterative spectrum step length acquisition method comprises the following steps:

[0212] S31, obtaining the optimal penalty parameter of iteration;

[0213] S32, obtaining an estimate of iterative curvature;

[0214] S33, obtaining the iterative spectrum step size;

[0215] In S31, x i,k and Constructed in dual form:

[0216]

[0217]

[0218] in,

[0219]

[0220]

[0221]

[0222]

[0223]

[0224]

[0225] By minimizing The residual at , obtains the optimal penalty parameter of the iteration, expressed as:

[0226]

[0227] In S32, the iterative curvature estimation is expressed as:

[0228]

[0229]

[0230]

[0231]

[0232]

[0233]

[0234] in,

[0235] In S33, the iterative spectrum step size is expressed as:

[0236]

[0237]

[0238]

[0239] Table 1 gives the conditions for the parameters of the cost function:

[0240] Table 1 Initial conditions and aerodynamics

[0241]

[0242] Four experimental scenarios are set during the simulation, including:

[0243] Experimental scenario 1: Four drones, each of which is flying towards the other two drones, with a circular obstacle in the environment, from which all drones need to keep a safe distance. In addition, the flight time of all drones is free and not restricted by constraints, which means that each drone can optimize its own flight time without meeting specific timing requirements with its neighbors (although each drone still needs to share time variables with its neighbors and reach consensus). The four swarm drones in this scenario are neighbors of each other. Table 2 summarizes the specific state of each drone, the location and radius of the obstacle, and the initial value and boundary of the flight time.

[0244] Table 2 Initial conditions of experimental scenario 1

[0245]

[0246]

[0247] Experimental scenario 2: Considering the scenario of five drones intercepting continuously, the drone cluster intercepts different targets in the same area from different directions, and has position and flight path angle constraints. In addition, the drone cluster needs to intercept the target within a certain time interval, that is, to meet the time series constraint. The initial time of all drones is set to 9.0s, and the collision time of all drones meets the equal interval The constraint is set to 0.01s for all drones. In addition, the neighborhood size of each drone is set to The neighborhood set is selected based on the initial position of the drone using the closest distance criterion. Table 3 summarizes the relevant settings of experimental scenario 2, where the same parameters as those of experimental scenario 1 are not repeated.

[0248] Table 3 Initial conditions of experimental scenario 2

[0249]

[0250]

[0251] Experimental scenario 3: 16 drones are evenly distributed in a circular formation, crossing the obstacle in the center of the formation and flying to another point symmetrical to the center of the circle. In the process of flying to the target point, the 16 cluster drones need to communicate with neighboring drones and maintain a safe distance. The neighborhood size of each drone is set to In addition, all drones are required to complete their flight missions simultaneously, which means that each drone needs to optimize its own flight time while keeping the same time with its neighboring drones. Table 4 summarizes the relevant settings of experimental scenario 3.

[0252] Table 4 Initial conditions of experimental scenario 3

[0253]

[0254] Experimental scenario 4: Let 20 drones be evenly distributed in sequence, pass through several random obstacles, and fly to the other side of the scene. In the process of flying to the target point, each cluster drone needs to communicate with a neighborhood size of The UAVs can communicate and maintain a safe distance. Due to the complex environment, the flight path length of each UAV may be different, so each UAV optimizes its own flight time individually without satisfying constraints.

[0255] Comparative Example 1

[0256] The same simulation experiment as in Example 1 is performed, except that the original fixed penalty parameter is used, that is, no penalty parameter is set during the acquisition of the spatiotemporal trajectory of the drone cluster.

[0257] Comparative Example 2

[0258] The same simulation experiment as in Example 1 is performed, except that the residual balance (RB) scheme and the Nesterov acceleration (NA) scheme are used.

[0259] The residual balance (RB) scheme can be found in Saravanos, A.D., Aoyama, Y., Zhu, H., and Theodorou, E.A., “Distributed Differential Dynamic Programming Architectures for Large-Scale Multiagent Control,” IEEE Transactions on Robotics, 2023, pp. 1–21. https: / / doi.org / 10.1109 / TRO.2023.3319894;

[0260] For the Nesterov acceleration (NA) scheme, see Goldstein, T., O'Donoghue, B., Setzer, S., and Baraniuk, R., "Fast Alternating Direction Optimization Methods," SIAM Journal on Imaging Sciences, Vol. 7, No. 3, 2014, pp. 1588–1623. https: / / doi.org / 10.1137 / 120896219.

[0261] The simulation results of Example 1, Example 2 and Comparative Example 1 are shown in Table 5.

[0262] Table 5. The number of iterations when different methods reach convergence

[0263] plan Comparative Example 1 Comparative Example 2-RB Method Comparative Example 2—NA Method Example 1 Scenario 1 102 64 70 48 Scenario 2 128 100 129 95 Scene 3 89 73 69 68 Scene 4 147 71 130 68

[0264] It can be seen from Table 5 that the method in Example 1 can significantly reduce the number of iterations and improve efficiency in different task scenarios compared to other methods.

[0265] The present invention has been described above in conjunction with preferred embodiments, but these embodiments are only exemplary and serve only as an illustration. On this basis, the present invention may be subjected to a variety of substitutions and improvements, all of which fall within the scope of protection of the present invention.

Claims

1. A penalty parameter adaptive adjustment criterion based on the spectral gradient method, characterized in that: include: Establish optimal control model; Solve the optimal control model to obtain the spatiotemporal trajectory of the drone cluster; The UAV flies according to the obtained space-time trajectory.

2. The penalty parameter adaptive adjustment criterion based on the spectral gradient method according to claim 1, characterized in that: The optimal control model is expressed as: s.t.x i，k+1 =f i (x i，k ,u i，k ) Among them, i represents the number of different drones in the drone cluster, N represents the total time, represents the minimum value of the control input of UAV i at all times, represents the final flight time of drone i, t i，N represents the flight time of drone i, X i Represents the stack of state vectors of drone i at all times, U i represents the stack of control inputs of drone i at all times, J i (X i , U i ) represents the local objective function of the i-th UAV, st represents the constraint condition, x i，k represents the state vector of the i-th UAV at time k, u i，k represents the control vector of the ith UAV at time k, fi() represents the discrete system dynamics with nonlinear characteristics, For x i，k A safe copy variable, a i，k () is a nonlinear function used to describe path constraints, for u i，k A safe copy variable, b i，k () represents the control torque constraint function, t i，N A safe copy variable, c i，N () represents a nonlinear function, represents the set of state vectors of all neighboring drones of drone i at time k, represents the set of state constants between UAV i and its neighboring UAVs at time k, represents the set of time variables of all neighboring drones of drone i at time k, represents the set of time constraints between UAV i and its neighboring UAVs at time k, represents the set of state variables that have reached consensus, A collection of variables representing the time at which consensus is reached.

3. The penalty parameter adaptive adjustment criterion based on the spectral gradient method according to claim 1, characterized in that: The spatiotemporal trajectory acquisition process of the drone cluster includes the following steps: S21. At time k, each drone obtains its own state vector x i，k , control vector u i，k and flight time t i，N ; S22. At time k, each drone obtains its own secure copy variable S23, each drone communicates with all its neighbor drones to obtain the state variables of drone i observed by neighbor drone j and the flight time of drone i observed by neighboring drone j S24. Each drone calculates the consensus state variable z based on its own variables and the variables obtained from its neighbors. i，k and the time variable s to reach consensus i，N ; S25. Each drone communicates with all its neighbor drones to obtain the state variable z calculated by the neighbor drone j of the drone for reaching consensus. j，k and the time variable s used to reach consensus j，N ; S26, each drone uses the state variable z that it has obtained to reach consensus i，k and the time variable s to reach consensus i，N and the consensus state variable z of neighbor drone j j，k and the time variable s to reach consensus j，N Based on the basic data, the ADMM method is used to solve the optimal control variables at time k, obtain the dual variables at time k, and obtain the trajectory of the UAV at time k by minimizing the residual; S27. Repeat steps S21-S26 to obtain the trajectory of each drone at the next moment in turn.

4. The penalty parameter adaptive adjustment criterion based on the spectral gradient method according to claim 3 is characterized in that: In the process of acquiring the spatiotemporal trajectory of the UAV cluster, a penalty parameter is set for each dual variable to accelerate the convergence speed of the solution.

5. The penalty parameter adaptive adjustment criterion based on the spectral gradient method according to claim 4, characterized in that: The penalty parameter is used to speed up the convergence process by: Update the dual variables; Get the iterative spectrum step size; Get the penalty parameter for the next iteration based on the spectrum step size.

6. The penalty parameter adaptive adjustment criterion based on the spectral gradient method according to claim 5, characterized in that: The updating process of the dual variables is expressed as: The superscript n represents the number of iterations. represents the dual variable of the ith UAV at the nth iteration, is the dual variable at the n+1th iteration, is the dual variable ζ i，k , i，k 、v i，N ,y i，k , η i，N The penalty parameter corresponding to the nth iteration.

7. The penalty parameter adaptive adjustment criterion based on the spectral gradient method according to claim 5, characterized in that: The iterative spectrum step length acquisition method comprises the following steps: S 31. Obtaining the optimal penalty parameter for iteration; S 32. Obtaining an estimate of the iterative curvature; S 33. Obtain the iterative spectrum step size.

8. The penalty parameter adaptive adjustment criterion based on the spectral gradient method according to claim 7, characterized in that: In S31, x i，k and Constructed in dual form: in, By minimizing The residual at , obtains the optimal penalty parameter of the iteration, expressed as: Among them, H(x i，k ), For custom functions, is the indicator function, H * () represents the Fenchel conjugate of H(), G * () represents the Fenchel conjugate of G(), Indicates H * The sub-gradient of (), Represents G ** The sub-gradient of (), is the intermediate variable, is a constant vector, is the dual function The first-order linear approximation of is the dual function The first-order linear approximation of Indicates the spectrum step size of the nth iteration, the indicator function In the subscript f i 、a i，k Represents a set, and the indicator function is defined as: if the function variable is in the set, the indicator function is 0; if the function variable is not in the set, the indicator function is positive infinity.

9. The penalty parameter adaptive adjustment criterion based on the spectral gradient method according to claim 7, characterized in that: In S 32, the iterative curvature estimate is expressed as: in, Indicates the curvature at the nth iteration The estimated value of Indicates the curvature at the nth iteration The superscript SD represents the steepest descent of the parameter, MG represents the minimum gradient of the parameter, and the superscript n l Represents the result of the last iteration of the parameter.

10. The penalty parameter adaptive adjustment criterion based on the spectral gradient method according to claim 7, characterized in that: In S33, the iterative spectrum step size is expressed as: in, represents the spectrum step size of the n+1th iteration, ∈ cor is the accuracy threshold of the curvature estimate, is the dual function The cosine of the angle between the change in the dual subgradient and the change in the dual variable, is the dual function The cosine of the angle between the change in the dual subgradient and the change in the dual variable.

Images