Multi-Batch Carrier-Based UAV Formation Rendezvous Method Based on Improved Sequential Convex Programming

Abstract

The limitations of the existing catapults necessitate multiple batches of take-offs for carrier-based unmanned aerial vehicles (UAVs) to form a formation. Because of the differences in takeoff time and location of each batch of UAVs, ensuring the temporal and spatial consistency and rendezvous efficiency of the formation becomes crucial. Concerning the challenges mentioned above, a multi-batch formation rendezvous method based on improved sequential convex programming (SCP) is proposed. A reverse solution approach based on the multi-batch rendezvous process is developed. On this basis, a non-convex optimization problem is formulated considering the following constraints: UAV dynamics, collision avoidance, obstacle avoidance, and formation consistency. An SCP method that makes use of the trust region strategy is introduced to solve the problem efficiently. Due to the spatiotemporal coupling characteristics of the rendezvous process, an inappropriate initial solution for SCP will inevitably reduce the rendezvous efficiency. Thus, an initial solution tolerance mechanism is introduced to improve the SCP. This mechanism follows the idea of simulated annealing, allowing the SCP to search for better reference solutions in a wider space. By utilizing the initial solution tolerance SCP (IST-SCP), the multi-batch formation rendezvous algorithm is developed correspondingly. Simulation results are obtained to verify the effectiveness and adaptability of the proposed method. IST-SCP reduces the rendezvous time from poor initial solutions without significantly increasing the computing time.

1. Introduction

Compared to manned aircraft, unmanned aerial vehicles (UAVs) have gradually become a research hotspot due to their advantages such as their low cost, strong maneuverability, simple structure, and good stealthiness [1,2,3]. In July 2013, the US military’s unmanned flying-wing layout aircraft, X47-B, accomplished a fully autonomous landing, marking the arrival of the era of carrier-based UAVs [4]. However, due to the limited payload capacity of a single UAV, it is challenging to fulfill complex mission requirements using them [5]. To overcome this drawback, multi-UAV collaboration and clustering have been proposed [6,7,8]. In multi-UAV coordination aboard carriers, the problem of multi-UAV formation rendezvous stands as one of the foundational aspects for achieving efficient collaboration.

The takeoff process of carrier-based UAVs relies on aircraft carriers, which face limitations due to their available deck space [9,10]: aircraft carrier decks typically have only 2–4 takeoff runways/catapults. Moreover, after the preceding aircraft takes off, there are delays due to the recharging of catapults and the aircraft wake vortices [11]. Subsequent aircraft using the same runway need to wait for a certain period before taking off. When the number of aircraft needing to be assembled exceeds the number of takeoff runways, the aircraft awaiting assembly cannot take off simultaneously. This results in the issue of multiple batches of formation rendezvous, causing aircraft from different batches to have different takeoff times. Moreover, during the batch takeoff process, due to the carrier’s movement, UAVs from different batches take off from different locations. The carrier-based UAV formation rendezvous problem aims to allow UAVs to take off at various times and locations and rendezvous into a specified formation.

The current research on UAV formation coordination and rendezvous primarily focuses on two types of problems.

The first type of rendezvous problem entails solving for the rendezvous location. Rucco et al. [12] and Wang et al. [13] addressed the rendezvous of small UAVs towards a moving unmanned ground vehicle (UGV), utilizing a rendezvous time predictive method to pre-determine the rendezvous time. Manyam et al. [14] developed a rendezvous method based on Bayesian quadratic curves, presetting the rendezvous time and planning for the rendezvous location and rendezvous path. McDonald et al. [15] simplified the rendezvous process into uniform and uniformly accelerated motions to predict and determine the feasible rendezvous time, thereby solving for the rendezvous locations and trajectories.

The second type of rendezvous problem entails solving for the rendezvous time. This problem was first investigated by McLain et al. [16,17], who aimed to plan the path for multiple UAVs in hostile environments to minimize the length of the rendezvous path, ensuring they all arrive at the target simultaneously for a coordinated rendezvous. Chen et al. [18] preset the rendezvous location and solved for the rendezvous time by developing a decentralized two-layer coordinative framework. Zhao et al. [19] developed a distributed approach that, under the pre-determined rendezvous location, optimized the rendezvous time for multiple UAVs.

The rendezvous problem studied in this paper is more complex than both the first type of problem and the second type of problem due to the inter-batch coordination. The rendezvous positions and times of the preceding batches influence those of the subsequent batches. When aiming to minimize the rendezvous time for the entire formation, a critical issue that demands resolution is determining the optimal time and location for each batch of UAVs to rendezvous into the formation. Thus, the rendezvous problem of carrier-based UAVs necessitates optimizing for both the rendezvous time and the rendezvous location. Research on this type of rendezvous problem can hardly be found in the existing literature. It is worthwhile to explore a batch-coordinated rendezvous method for carrier-based UAVs.

Currently, classical methods for formation rendezvous include velocity–consistency-based approaches [20,21], guidance-based methods [22], and strategies based on multi-agent cooperative path planning [23]. With the rapid advancement of path planning algorithms in recent years, rendezvous methods based on collaborative path planning strategies are receiving more focus. Jouffroy et al. [24] utilized particle swarm optimization to optimize the parameters of Dubins curves, thereby forming trajectories for the UAVs. Yao et al. [25] employed the A* algorithm to acquire the motion trajectories for multiple UAVs to rendezvous, and subsequently smoothed these paths into curves meeting the flight requirements for the UAVs using third-order B-spline curves. Shao et al. [26] utilized an enhanced particle swarm optimization algorithm to optimize the PH path parameters of individual UAVs, achieving the objective of formation rendezvous.

The aforementioned research focused on optimizing the geometric parameters of the UAV’s path, considering spatial consistency during UAV formation rendezvous. Temporal consistency constraints of the problem are usually transformed into spatial consistency constraints by assuming uniform motion of the aircraft. However, in the context of multi-batch carrier-based UAV formation rendezvous, the departures of aircraft from different batches do not occur simultaneously or at the same location. Under the constraint of consistent final velocities, it is impossible to form a formation using the assumption of uniform motion. Thus, the kinematic constraints of aircraft and the consistency constraints in both the temporal and spatial dimensions need to be considered.

When confronting offline trajectory planning problems that consider aircraft kinematic equations and spatiotemporal coupling consistency, the aforementioned planning approaches, such as swarm optimization algorithms, require the discretization of variables for encoding. During the process of space discretization, such methods might encounter issues like dimension explosion [27], high computational complexity, and slow convergence [28]. Hence, exploring methods suitable for handling dynamic constraints, multi-agent collision avoidance, and time-varying consistency constraints in this problem is crucial.

Convex optimization methods, as typical mathematical programming approaches, have a well-developed foundation in optimization theory and make use of numerous efficient polynomial-time solving algorithms. Based on the convex optimization concept, sequential convex programming (SCP) guarantees optimality while maintaining a high computational efficiency [29].

In the field of UAV formation rendezvous, this method is a promising branch [30]. Wang and McDonald [13,15] optimized the rendezvous location of UAVs based on SCP. Chen [18] and Zhao [19] optimized the rendezvous time based on SCP. Therefore, the SCP method is also applicable to the inter-batch coordination issue addressed in this paper.

When the final rendezvous state varies with time, it is difficult to obtain an appropriate initial guess for the SCP. Mao et al. [31] introduced trust region and virtual control techniques into SCP, thereby allowing the initial solution to not necessarily comply with the original problem’s dynamic constraints. Based on this work, they extended the technique to non-convex state and control constraints [32], thereby enabling the initial solution to not necessarily satisfy the original problem’s other non-convex constraints. The successive convexification algorithm (Scvx) they proposed serves as the baseline SCP employed in this paper. These works reduced the difficulty in obtaining an initial solution, yet there remains room for improvement in the sensitivity to initial solutions of the SCP [33]. Banerjee et al. [34] developed a learning-based method to offer a better initial solution for SCP. Instead of optimizing the initial solution, an initial solution tolerance mechanism is developed which decreases the algorithm’s sensitivity for initial values. Inspired by the simulated annealing (SA) algorithm, a novel SCP method that was improved by introducing an initial solution tolerance mechanism is proposed in this article.

The objective of this paper is to propose an efficient rendezvous method for addressing the formation rendezvous problem for carrier-based UAVs that take off in multiple batches. Considering the time intervals and starting/ending location differences during the takeoff process that occur between different carrier-based UAVs, a reverse solution approach is proposed. Based on this, a non-convex optimization problem is derived mathematically, and the optimization model is established to minimize the formation rendezvous time cost. In the existing literature, there is a lack of research on trajectory planning problems for multi-batch UAVs’ rendezvous; therefore, different rendezvous strategies considering different terminal specifications are urgently needed. Subsequently, under the concept of convex optimization, the relaxation technique is introduced to convexify these original non-convex problems. To develop an efficient solver, a rendezvous algorithm is proposed based on SCP that is improved by introducing an initial solution tolerance mechanism. Consequently, the issue that an inappropriate initial solution reduces the rendezvous efficiency is successfully circumvented. Finally, experimental scenarios for the formation of carrier-based UAVs are designed; thereafter, simulation experiments are conducted to validate the effectiveness and adaptability of the algorithm proposed in this paper. Specifically, the initial solution tolerance SCP (IST-SCP) performs better in terms of the rendezvous time compared with the standard SCP, and it has a shorter computing time compared with the initial solution optimization method.

The rest of the paper is organized as follows: Section 2 presents the multi-batch formation rendezvous optimal control problem formulated by the reverse solution approach. Section 3 offers the details of formation rendezvous algorithms based on IST-SCP. The simulation results and analysis are presented in Section 4. Finally, the conclusions are drawn in Section 5.

2. Problem Formulation

In this section, the dynamics model of UAVs is introduced. Then, the process of the multi-batch formation rendezvous of carrier-based UAVs is illustrated and, based on this, a reverse solution approach for formation rendezvous is proposed. Finally, other constraints for the formation rendezvous problem are added, and the non-convex optimal control problem is formulated.

2.1. Dynamics Model of UAVs

The subject of this study is fixed-wing carrier-based UAVs, and this study aims to solve the problem of their rendezvous into a specified formation through multiple batches. The utilized drone dynamics model, employed as the dynamic constraints for the optimal control problem, is as follows.

$$\left\{\begin{array}{c}{\dot{x}}_g=Vcos\gamma cos\chi \\ {\dot{y}}_g=Vcos\gamma sin\chi \\ {\dot{z}}_g=Vsin\gamma \\ \dot{V}=g(n_x-sin\gamma )\\ {\dot{\chi }}_{sin}={\displaystyle \frac{g{n}_y}{Vcos\gamma }}{\chi }_{cos}\\ {\dot{\chi }}_{cos}=-{\displaystyle \frac{g{n}_y}{Vcos\gamma }}{\chi }_{sin}\\ \dot{\gamma }={\displaystyle \frac{g(n_z-cos\gamma )}{V}}\end{array}\right.$$

(1)

where $x_g,y_g,z_g$ represent the ground-based coordinates; $V,\chi ,\gamma$ stand for the UAV’s speed, heading angle, and flight path angle, respectively; $\mathit{u}={\left[n_x,n_y,n_z\right]}^T$ represent UAV control variables, indicating the overload in the $xyz$ directions within the flight path coordinated frame. Since the heading angle $\chi$ changes periodically, the cosine and sine values are utilized in the state space instead of $\chi$ itself, as ${\chi }_{sin}=sin\chi , {\chi }_{cos}=cos\chi$.

2.2. Formation Rendezvous Process in Multiple Batches

For analytical convenience, let us assume that each batch of UAVs arrives at the designated formation position simultaneously. For $N$ UAVs launched in $q$ batches, each batch consisting of $\left(p_1,p_2,...,p_q\right)$ UAVs, these UAVs need to assemble into a specific formation. Assuming that UAVs within the same batch enter the formation simultaneously, this problem can be expressed as follows:

$$\begin{gathered} t_{if}=t_f^j \\ \sum _{m=1}^{j-1}{p}_m<i⩽\sum _{m=1}^j{p}_m,1⩽j⩽q \end{gathered}$$

(2)

where $i$ represents the ith UAV, $j$ represents the jth batch of UAVs, $t_{if}$ denotes the time at which the ith UAV enters the formation, and $t_f^j$ represents the time at which the jth batch of UAVs enters the formation.

Due to the time intervals between batches during takeoff and the inability of fixed-wing aircraft to hover, this paper adopts a rendezvous mode for UAV formations where the UAVs circle and rendezvous at a specified location and altitude, as illustrated in Figure 1.

In Figure 1, the process of forming a hexagonal formation using three batches of UAVs as an example is described: each batch of UAVs is launched from a path-known carrier at a specified time at a certain interval. They take off and fly towards a circling orbit to enter the formation. The position of each UAV within the formation is predetermined. The first batch of UAVs to take off enters the inner circling orbit, and these UAVs circle within the orbit to await the later UAVs. Upon entering the orbit, they directly assume their designated formation positions, forming a partial formation. The coordinates of the formation centroid when the jth batch of aircraft enters the formation are $\left(x_c^j,y_c^j\right)$. Each UAV uses the radius $R_c$ of the orbit where the formation centroid is located as the reference radius, calculates the radius $R_i$ of the circling orbit, and calculates the angular position relative to the formation centroid based on the geometric shape of the formation to attain its designated formation position. The angular velocity ${\omega }_c$ of the UAV formation remains constant throughout this process, and the formation conducts uniform circular motion at the given orbit radius $R_c$. When the last UAV assumes its formation position, the formation rendezvous task is completed.

In terms of the rendezvous sequence, a Gantt chart illustrating the circling rendezvous process is depicted in Figure 2.

In this problem, the optimization objective is the completion time $t_f$ of the rendezvous. Because the UAVs’ velocity is greater than that of the carrier, the rendezvous time of the front batch is shorter than that of the following batch when the circling orbit is designated, which means that $∆t_f^j⩾0$. Thus, optimizing the completion time $t_f$ is equivalent to optimizing the time when the last batch of aircraft reaches the designated formation position $t_f^q$.

2.3. Reverse Solution Approach

To address the challenges of spatiotemporal coupling during the rendezvous process while maintaining the optimality of the rendezvous time, a reverse solution approach is proposed.

As depicted in Figure 3, this approach allows us to calculate the shortest time for the final batch of UAVs to enter the circling orbit. Using this time as a reference, the entry time and trajectories for the other batches of UAVs are reverse solved. This adjustment enables these batches to synchronize their entry into the circling orbit according to the entry time of the final batch of UAVs, meeting the consistency constraints.

The objective function for each batch is represented as follows:

$$J=w_t{t}_f^j+w_z{g}_h\left(z_g\right)$$

(3)

where $g_h\left(z_g\right)=‖z_c-z_g‖$ represents the altitude penalty term, driving the UAVs’ climb to the designated altitude $z_c$.

For the final batch of UAVs, this objective function serves as the optimization goal for this problem. However, for the preceding q − 1 batches of UAVs, utilizing this objective can facilitate these batches’ swift entry into the formation, avoiding any impact on subsequent batches.

In this process, the circling orbit radius and circling angular velocity of the UAV formation are set. Therefore, the final constraint for the ith UAV can be represented in the following form:

$${\mathit{x}}_f^i=\left[\begin{array}{c}{x}_H+R_{i}cos\left({\omega }_c\left(t_f^j-t_f^q\right)-\left({\displaystyle \frac{\pi }{2}}-{\chi }_i\right)\right)\\ y_H+R_{i}sin\left({\omega }_c\left(t_f^j-t_f^q\right)-\left({\displaystyle \frac{\pi }{2}}-{\chi }_i\right)\right)\\ z_c\\ V_i\\ sin\left({\omega }_c\left(t_f^j-t_f^q\right)+{\chi }_i\right)\\ cos\left({\omega }_c\left(t_f^j-t_f^q\right)+{\chi }_i\right)\\ 0\end{array}\right]$$

(4)

where $\left(x_H,y_H\right)$ represents the coordinates of the circling orbit’s center, and $t_f^q$ denotes the time when the final batch of aircraft enters the formation. The formation’s geometric parameters only affect $R_i,{\chi }_i$, avoiding the influence of the formation shape on the problem-solving process.

2.4. Other Constraints and Non-Convex Optimization Problem Formulation

Besides meeting the dynamic constraints and the consistency constraints mentioned earlier, the rendezvous process also needs to satisfy the following constraints.

Considering the spatial constraints of the formation rendezvous task and the limitations of UAVs’ flight performance, we impose constraints on the state and control variables mentioned in (1) as follows:

$$\begin{gathered} {\mathit{x}}_{\mathrm{m}\mathrm{i}\mathrm{n}}⩽\mathit{x}⩽{\mathit{x}}_{\mathrm{m}\mathrm{a}\mathrm{x}} \\ {\mathit{u}}_{\mathrm{m}\mathrm{i}\mathrm{n}}⩽\mathit{u}⩽{\mathit{u}}_{\mathrm{m}\mathrm{a}\mathrm{x}} \end{gathered}$$

(5)

The UAVs need to avoid the no-fly zones along the path while heading to the formation area. The no-fly zone avoidance constraint is as follows:

$$\begin{array}{c}‖\mathit{C}{\mathit{x}}_i-{\mathit{x}}_{ob}^l‖\ge r+R_{ob}^l\\ i=1, 2, \dots , N\\ l=1, 2, \dots , L\end{array}$$

(6)

where $\mathit{C}$ represents a matrix that extracts the UAV’s state information $\left(x_i^{\tau },y_i^{\tau }\right)$ from the state vector ${\mathit{x}}_i$; ${\mathit{x}}_{ob}^l$ represents the center coordinates of the lth no-fly zone; $r$ denotes the safety radius of the aircraft, and $R_{ob}^l$ represents the radius of the lth no-fly zone.

To prevent collisions among the UAVs within the formation, the collision avoidance constraint is as follows:

$$\begin{array}{c}‖\mathit{C}\left({\mathit{x}}_i-{\mathit{x}}_l\right)‖\ge 2r\\ i=1, 2, \dots , N\\ l=1, 2, \dots , N\\ i\ne l\end{array}$$

(7)

where ${\mathit{x}}_i,{\mathit{x}}_l$ represent the states of the ith and lth UAVs, respectively.

The initial state constraint is as follows:

$${\mathit{x}}_0^i={\mathit{c}}^i$$

(8)

where ${\mathit{c}}^i$ is a constant, and the state at which each carrier-based UAV takes off is predetermined.

Thus, the non-convex optimal control problem, the MFRP (multi-batch formation rendezvous problem), is formulated as follows:

$$\begin{array}{cc}min& J=w_t{t}_f^j+w_z{g}_h\left(z_g\right)\\ s.t.& \left(1\right)\left(2\right)\left(4\right)\left(5\right)\left(6\right)\left(7\right)\left(8\right)\end{array}$$

(9)

Among these, constraints (1), (4), (6), and (7) are non-convex and require convexification.

3. Multi-Batch Formation Rendezvous Algorithm

In this section, convexification and discretization techniques are introduced first. Thereafter, an improved SCP algorithm incorporating an initial solution tolerance mechanism is derived. Subsequently, a rendezvous optimization method is developed with the improved SCP as its core for the formation rendezvous problem.

3.1. Constraints Convexification and Discretization

As mentioned previously, the MFRP contains non-convexity. Therefore, it is necessary to convexify the non-convex constraints around a certain reference solution $\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$. In this process, introducing a trust region constraint is crucial to ensuring the effectiveness of the convexification approximation. However, the solution space near the reference solution might not always contain a solution after convexification. Hence, the relaxation factors are added to the constraints to prevent situations with no feasible solution.

The linear convexification of (1), (4), (6), (7) is shown in (10–13).

$${\dot{\mathit{x}}}_{\mathit{i}}={\mathit{A}}_{\mathit{i}}{\mathit{x}}_{\mathit{i}}+{\mathit{B}}_{\mathit{i}}{\mathit{u}}_{\mathit{i}}+{\mathit{D}}_{\mathit{i}}$$

(10)

$$‖\mathit{C}{\overline{\mathit{x}}}_i-{\mathit{x}}_{ob}^l‖+{\displaystyle \frac{(\mathit{C}{\overline{\mathit{x}}}_i-{\mathit{x}}_{ob}^l{)}^T}{‖\mathit{C}{\overline{\mathit{x}}}_i-{\mathit{x}}_{ob}^l‖}}(\mathit{C}{\mathit{x}}_i-\mathit{C}{\overline{\mathit{x}}}_i)\ge r+R_{ob}^l$$

(11)

$${\displaystyle \frac{({\overline{\mathit{x}}}_i-{\overline{\mathit{x}}}_l{)}^T{\mathit{C}}^T}{‖\mathit{C}({\overline{\mathit{x}}}_i-{\overline{\mathit{x}}}_l)‖}}\mathit{C}({\mathit{x}}_i-{\mathit{x}}_l)\ge 2r$$

(12)

$${\mathit{x}}_f^i=\left[\begin{array}{c}{x}_H+R_{i}cos\left({\overline{\chi }}_i^{t_f}-{\displaystyle \frac{\pi }{2}}\right)-R_{i}sin\left({\overline{\chi }}_i^{t_f}-{\displaystyle \frac{\pi }{2}}\right){\omega }_c∆t_f^j\\ y_H+R_{i}sin\left({\overline{\chi }}_i^{t_f}-{\displaystyle \frac{\pi }{2}}\right)+R_{i}cos\left({\overline{\chi }}_i^{t_f}-{\displaystyle \frac{\pi }{2}}\right){\omega }_c∆t_f^j\\ z_c\\ V_i\\ cos{\overline{\chi }}_i^{t_f}-{\omega }_{c}sin{\overline{\chi }}_i^{t_f}∆t_f^j\\ sin{\overline{\chi }}_i^{t_f}+{\omega }_{c}cos{\overline{\chi }}_i^{t_f}∆t_f^j\\ 0\end{array}\right]$$

(13)

The coefficient matrices A, B, and D in (10) are as shown in (14).

$$\begin{gathered} {\mathit{A}}_{\mathit{i}}=\left[\begin{array}{ccccccc}0& 0& 0& \mathrm{c}\mathrm{o}\mathrm{s}\overline{\chi }\mathrm{c}\mathrm{o}\mathrm{s}\overline{\gamma }& \overline{V}\mathrm{c}\mathrm{o}\mathrm{s}\overline{\gamma }& 0& -\overline{V}\mathrm{c}\mathrm{o}\mathrm{s}\overline{\chi }\mathrm{s}\mathrm{i}\mathrm{n}\overline{\gamma }\\ 0& 0& 0& \mathrm{s}\mathrm{i}\mathrm{n}\overline{\chi }\mathrm{c}\mathrm{o}\mathrm{s}\overline{\gamma }& 0& \overline{V}\mathrm{c}\mathrm{o}\mathrm{s}\overline{\gamma }& -\overline{V}\mathrm{s}\mathrm{i}\mathrm{n}\overline{\chi }\mathrm{s}\mathrm{i}\mathrm{n}\overline{\gamma }\\ 0& 0& 0& \mathrm{s}\mathrm{i}\mathrm{n}\overline{\gamma }& 0& 0& \overline{V}\mathrm{c}\mathrm{o}\mathrm{s}\overline{\gamma }\\ 0& 0& 0& 0& 0& 0& -g\mathrm{c}\mathrm{o}\mathrm{s}\overline{\gamma }\\ 0& 0& 0& {\displaystyle \frac{g{\overline{n}}_y}{{\overline{V}}^{2}cos\overline{\gamma }}}\mathrm{s}\mathrm{i}\mathrm{n}\overline{\chi }& 0& -{\displaystyle \frac{g{\overline{n}}_y}{\overline{V}cos\overline{\gamma }}}& -{\displaystyle \frac{g{\overline{n}}_y}{\overline{V}co{s}^2\overline{\gamma }}}\mathrm{s}\mathrm{i}\mathrm{n}\overline{\chi }sin\overline{\gamma }\\ 0& 0& 0& -{\displaystyle \frac{g{\overline{n}}_y}{{\overline{V}}^{2}cos\overline{\gamma }}}\mathrm{c}\mathrm{o}\mathrm{s}\overline{\chi }& {\displaystyle \frac{g{\overline{n}}_y}{\overline{V}cos\overline{\gamma }}}& 0& {\displaystyle \frac{g{\overline{n}}_y}{\overline{V}co{s}^2\overline{\gamma }}}\mathrm{c}\mathrm{o}\mathrm{s}\overline{\chi }sin\overline{\gamma }\\ 0& 0& 0& -{\displaystyle \frac{g({\overline{n}}_z-cos\overline{\gamma })}{{\overline{V}}^2}}& 0& 0& {\displaystyle \frac{gsin\overline{\gamma }}{\overline{V}}}\end{array}\right] \\ {\mathit{B}}_{\mathit{i}}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ g& 0& 0\\ 0& -{\displaystyle \frac{g\mathrm{s}\mathrm{i}\mathrm{n}\overline{\chi }}{\overline{V}cos\overline{\gamma }}}& 0\\ 0& {\displaystyle \frac{g\mathrm{c}\mathrm{o}\mathrm{s}\overline{\chi }}{\overline{V}cos\overline{\gamma }}}& 0\\ 0& 0& {\displaystyle \frac{g}{\overline{V}}}\end{array}\right] \\ {\mathit{D}}_{\mathit{i}}={\dot{\overline{\mathit{x}}}}_{\mathit{i}}-{\mathit{A}}_{\mathit{i}}{\overline{\mathit{x}}}_{\mathit{i}}-{\mathit{B}}_{\mathit{i}}{\overline{\mathit{u}}}_{\mathit{i}} \end{gathered}$$

(14)

Then, introduce the trust region constraint as shown in (15). Note that the radius of the trust region is updated following the solving of each sequential subproblem using the update strategy proposed in Section 3.2.

$$\begin{gathered} ‖\mathit{x}-\overline{\mathit{x}}‖⩽{\delta }_x \\ ‖\mathit{u}-\overline{\mathit{u}}‖⩽{\delta }_u \\ \left|t-\overline{t}\right|⩽{\delta }_t \end{gathered}$$

(15)

Hence, the convexified MFRCS (multi-batch formation rendezvous convex subproblem) is depicted below.

$$\begin{array}{cc}min& J=w_t{t}_f^j+w_z{g}_h\left(z_g\right)\\ s.t.& \left(10\right)\left(2\right)\left(11\right)\left(5\right)\left(12\right)\left(13\right)\left(8\right)\left(15\right)\end{array}$$

(16)

For the MFRCS, when the trust region is inappropriate, it might lead to conflicts with other constraints. To ensure stability in the solving process, relaxation factors are added to each constraint and penalized in the objective function. This results in the MFRCSR (multi-batch formation rendezvous convex subproblem with relaxation factors), shown as follows:

$$\begin{array}{cc}\mathrm{m}\mathrm{i}\mathrm{n}& J=w_t{t}_f+w_z{g}_h\left(z_g\right)+w_p{p}_u\left({\mathit{v}}_s,{\mathit{v}}_o,{\mathit{v}}_c,{\mathit{v}}_{ic},{\mathit{v}}_{tc}\right)\end{array}$$

(17a)

$$\begin{array}{cc}s.t.& {\dot{\mathit{x}}}_i=\mathit{A}{\mathit{x}}_i+\mathit{B}{\mathit{u}}_i+{\mathit{D}}_i+{\mathit{v}}_s^i\end{array}$$

(17b)

$$‖\mathit{C}{\overline{\mathit{x}}}_i-{\mathit{x}}_{ob}^l‖+{\displaystyle \frac{(\mathit{C}{\overline{\mathit{x}}}_i-{\mathit{x}}_{ob}^l{)}^T}{‖\mathit{C}{\overline{\mathit{x}}}_i-{\mathit{x}}_{ob}^l‖}}(\mathit{C}{\mathit{x}}_i-\mathit{C}{\overline{\mathit{x}}}_i)\ge r+R_{ob}^l+v_o^i$$

(17c)

$${\displaystyle \frac{({\overline{\mathit{x}}}_i-{\overline{\mathit{x}}}_l{)}^T{\mathit{C}}^T}{‖\mathit{C}({\overline{\mathit{x}}}_i-{\overline{\mathit{x}}}_l)‖}}\mathit{C}({\mathit{x}}_i-{\mathit{x}}_l)\ge 2r+v_c^i$$

(17d)

$${\mathit{x}}_0^i={\mathit{c}}^i+{\mathit{v}}_{ic}^i$$

(17e)

$${\mathit{x}}_f^i=g_f\left(\mathit{x},\mathit{u},t\right)+{\mathit{v}}_{tc}^i$$

(17f)

$$g_c\left(t\right)=0$$

(17g)

$$tru\left(\mathit{x},\mathit{u},t\right)⩽\mathit{\delta }$$

(17h)

where ${\mathit{v}}_s,{\mathit{v}}_o,{\mathit{v}}_c,{\mathit{v}}_{ic},{\mathit{v}}_{tc}$ refer to the relaxation factors for the motion constraints, obstacle avoidance constraints, collision avoidance constraints, initial state constraints, and final state constraints, respectively. In (17a), $w_p$ is the penalty weight, and the function $p_u\left({\mathit{v}}_s,{\mathit{v}}_o,{\mathit{v}}_c,{\mathit{v}}_{ic},{\mathit{v}}_{tc}\right)$ represents the penalty function for the relaxation factor to eventually converge to 0; as shown in (18), it adopts the sum of the 1-norms of the relaxation factors. By incorporating these relaxation factors, the algorithm’s initial solution and sequential solutions are no longer required to satisfy the original problem’s constraints [33]. Equations (17g,h) correspond to the time consistency constraint and the trust region constraint, respectively.

$$p_u\left({\mathit{v}}_s,{\mathit{v}}_o,{\mathit{v}}_c,{\mathit{v}}_{ic},{\mathit{v}}_{tc}\right)={‖{\mathit{v}}_s‖}_1+{‖{\mathit{v}}_o‖}_1+{‖{\mathit{v}}_c‖}_1+{‖{\mathit{v}}_{ic}‖}_1+‖{\mathit{v}}_{tc}‖$$

(18)

The MFRCSR still contains time-continuous equations with an infinite-dimensional variable. To simplify the computations, the dynamic equations in (17b) are approximated using linearized rectangular integration. Discretizing the MFRCSR results in the MFRDCS (multi-batch formation rendezvous discrete time convex subproblem), as represented by the following equation.

$$\begin{array}{cc}\mathrm{m}\mathrm{i}\mathrm{n}& J=w_t{t}_f+w_z{g}_h\left(z_g\right)+w_p{p}_u\left({\mathit{v}}_s,{\mathit{v}}_o,{\mathit{v}}_c,{\mathit{v}}_{ic},{\mathit{v}}_{tc}\right)\end{array}$$

(19a)

$$\begin{array}{c} \begin{array}{cc}s.t.& \mathit{x}{\left[k+1\right]}_i=\mathit{x}{\left[k\right]}_i+∆\overline{t}\left(\mathit{A}∆\mathit{x}{\left[k\right]}_i+\mathit{B}∆\mathit{u}{\left[k\right]}_i+{\mathit{D}}_i\right)\end{array}\\ +∆t\left(\mathit{A}∆\overline{\mathit{x}}{\left[k\right]}_i+\mathit{B}∆\overline{\mathit{u}}{\left[k\right]}_i+{\mathit{D}}_i\right)+\mathit{v}{\left[k\right]}_s^i\end{array}$$

(19b)

$$\begin{array}{c}{\displaystyle \frac{(\mathit{C}\overline{\mathit{x}}{\left[k\right]}_i-\mathit{x}{\left[k\right]}_{ob}^l{)}^T}{‖\mathit{C}\overline{\mathit{x}}{\left[k\right]}_i-\mathit{x}{\left[k\right]}_{ob}^l‖}}(\mathit{C}\mathit{x}{\left[k\right]}_i-\mathit{C}\overline{\mathit{x}}{\left[k\right]}_i)\\ +‖\mathit{C}\overline{\mathit{x}}{\left[k\right]}_i-\mathit{x}{\left[k\right]}_{ob}^l‖\ge r_i+R_{ob}^l+v{\left[k\right]}_o^i\end{array}$$

(19c)

$$\begin{array}{c}{\displaystyle \frac{(\overline{\mathit{x}}{\left[k\right]}_i-\overline{\mathit{x}}{\left[k\right]}_l{)}^T{\mathit{C}}^T}{‖\mathit{C}(\overline{\mathit{x}}{\left[k\right]}_i-\overline{\mathit{x}}{\left[k\right]}_l)‖}}\mathit{C}(\mathit{x}{\left[k\right]}_i-\mathit{x}{\left[k\right]}_l)\\ \ge 2r+v{\left[k\right]}_c^i\end{array}$$

(19d)

$$\mathit{x}{\left[1\right]}_i={\mathit{c}}^i+{\mathit{v}}_{ic}^i$$

(19e)

$$\mathit{x}{\left[K\right]}_i=g_5\left(\mathit{x}{\left[k\right]}_i,\mathit{u}{\left[k\right]}_i,t_i\right)+{\mathit{v}}_{tc}^i$$

(19f)

$$g_c\left(t\right)=0$$

(19g)

$$tru\left(\mathit{x}{\left[k\right]}_i,\mathit{u}{\left[k\right]}_i,t_i\right)⩽\mathit{\delta }$$

(19h)

It should be noted that, due to the time intervals between the takeoffs of the different batches of carrier-based UAVs, the $k$ representing time does not correspond to the same moment across different batches. The avoidance constraint mentioned above applies only to UAVs in the same batch. To handle avoidance between UAVs from different batches, it is necessary to compute the position of each UAV at the current time to ensure collision avoidance.

At this point, we have obtained the MFRDCS, which is suitable for convex optimization solvers.

3.2. Initial Solution Tolerance Sequential Convex Programming

When solving the convexified subproblems, convex optimization solvers may be misguided by these convexified constraints. As shown in Figure 4, when the reference solution is on the right side of the obstacle, the penalty generated by the convexified obstacle constraint can only approximately reflect the true penalty situation on the right side. Therefore, the variation in the sequence solution always stays on one side of the obstacle, and is only determined by the initial guess path’s position relative to the obstacle.

For the problem, in which the final state does not change over time, this limitation of the SCP method seems insignificant. However, when the final rendezvous state is closely related to the time, this limitation can lead the SCP method into a quite unsatisfactory local optimal solution.

For the q − 1 batches of UAVs, the final rendezvous positions circle within the orbit. As shown in Figure 5, the initial trajectory guess of type A (green dotted line) is on the left side of the centroid of the no-fly zone, and the initial path guess of type B (yellow dotted line) is on the right side of the centroid of the no-fly zone. Both A and B share one initial rendezvous position guess. As the SCP algorithm optimizes the rendezvous time, the rendezvous position moves in the opposite direction of the formation’s circling angular velocity (green arrow and yellow arrow). Path A (green curve), obstructed by the no-fly zone, is optimized to a local optimum because of the poor approximation of the penalty value on the right side of the no-fly zone. However, Path B (yellow curve) reached the global optimal solution with a better initial path guess.

It is difficult to accurately guess the initial solution of multiple UAVs for the complex environment/collision-free constraints/spatiotemporal relationship of the rendezvous problem. Instead of optimizing a better initial solution guess, an improved SCP method inspired by the simulated annealing (SA) algorithm is proposed.

The simulated annealing (SA) algorithm was first proposed by Metropolis et al. [35] in 1953. The SA algorithm starts from a higher initial temperature and, with a continuous decrease in temperature parameters, uses probability-based climbing to randomly search for the global optimal solution of the objective function in the solution space. That is, it can probabilistically jump away from the local optimal solution and eventually approach the global optimal solution.

The problem that the solution is trapped in the local optimum because of an inaccurate initial guess is mainly because the SCP method cannot handle the original problem of the non-convexity of obstacle constraints. It is unable to search for the optimal solution in the global solution space.

To fix this problem, an initial solution tolerance mechanism is combined with SCP. The initial solution tolerance SCP (IST-SCP) algorithm is shown in Figure 6.

At the beginning, an initial guess is provided by setting the initial rendezvous position and defining a linear interpolation between the takeoff and rendezvous states which is the same as the conventional SCP algorithm. Using this initial guess, the current minimum cost solution is initialized, which will be used to record the minimum cost solution. After convex relaxation around the reference solution, the probability-based search module is introduced. This module generates random disturbances with the obstacle avoidance penalty value ${\mathit{v}}_o$. When a random number between 0 and 1 is greater than the search probability P₁, we set the obstacle avoidance penalty to 0 for the optimizer (when calculating termination and update conditions, ${\mathit{v}}_o$ is calculated as the other penalty values). The purpose of doing so is to eliminate obstacles that hinder the rendezvous optimization, enabling the algorithm to explore spatiotemporal coupled trajectories in a wider space. Since $P_1=e^{-1/T}$, T is the annealing temperature, which decreases as the iteration increases, and P₁ will eventually drop to 0. As the annealing temperature decreases, the algorithm gradually stops the search frequency and ultimately optimizes the solution under the constraints of the MFRDCS.

In determining the termination criteria, two criteria need to be considered:

• Whether the search module is triggered;
• Whether the solution of the current iteration matches the solution of the previous iteration.

For criterion 1, when the search module is triggered, the solution obtained by the solver does not consider the penalty of the obstacle avoidance constraints, and this solution may not be feasible for both the MFRDCS and the original MFRP. Thus, the algorithm can only terminate when the search module is not triggered.

For criterion 2, based on the termination strategy in [33], ${\mathit{d}}_s,{\mathit{d}}_o,{\mathit{d}}_c,{\mathit{d}}_{ic},{\mathit{d}}_{tc}$ can be defined as shown in (20). These represent the real constraint deviation during MFRP solving, while their counterparts, represented by the relaxation factors ${\mathit{v}}_s,{\mathit{v}}_o,{\mathit{v}}_c,{\mathit{v}}_{ic},{\mathit{v}}_{tc}$, indicate the linearly approximated constraint deviation.

$${\mathit{d}}_s=\dot{\mathit{x}}-\left[\begin{array}{c}Vcos\gamma cos\chi \\ Vcos\gamma sin\chi \\ Vsin\gamma \\ g(n_x-sin\gamma )\\ {\displaystyle \frac{g{n}_y}{V}}cos\chi \\ -{\displaystyle \frac{g{n}_y}{V}}sin\chi \\ {\displaystyle \frac{g(n_z-cos\gamma )}{V}}\end{array}\right]$$

(20a)

$${\mathit{d}}_o=‖\mathit{C}{\mathit{x}}_i-{\mathit{x}}_{ob}^l‖-\left(r+R_{ob}^l\right)$$

(20b)

$${\mathit{d}}_c=‖\mathit{C}\left({\mathit{x}}_i-{\mathit{x}}_l\right)‖-2r$$

(20c)

$${\mathit{d}}_{ic}={\mathit{x}}_i^{t_0^i}-{\mathit{c}}^i$$

(20d)

$${\mathit{d}}_{tc}={\mathit{x}}_i^{t_f^i}-g_f\left(\mathit{x},\mathit{u},t\right)$$

(20e)

Criterion (2) can be expressed as follows:

$$∆J_1=w_t\left({\overline{t}}_f-t_f^{*}\right)+w_z\left(g_h\left({\overline{z}}_g\right)-g_h\left(z_g^{*}\right)\right)$$

(21)

where ${\overline{t}}_f,{\overline{z}}_g,\overline{\mathit{d}}$, corresponding to $\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}$, represent the solution obtained from solving the previous iteration’s subproblem, which also serves as the relaxed reference solution for this iteration’s subproblem. $t_f^{\mathrm{*}},z_g^{\mathrm{*}},{\mathit{v}}^{\mathrm{*}}$, corresponding to ${\mathit{x}}^{\mathrm{*}},{\mathit{u}}^{\mathrm{*}},t^{\mathrm{*}}$, represent the solution obtained from solving the subproblem in the current iteration ( superscript * stands for the optimal solution/objective value of the subproblem). $∆J_1$ represents the difference between the real (non-approximate) cost from the previous iteration and the convex approximate cost of the current iteration. When $∆J_1<\epsilon$, it is considered that the convex approximate cost of the current iteration’s solution is close enough to the real cost from the previous iteration, indicating convergence of the algorithm.

Note that $\overline{\mathit{v}}$ is not used here to describe the first part of criterion 2. This is because criterion 2 is mainly used to describe whether the solution of the nthiteration is consistent with that of the n − 1th iteration, indicating whether the algorithm has converged. $\overline{\mathit{v}}$ represents the deviation of the convex approximate constraint from that of the previous iteration. Due to different trust region selections between iterations, the level of convex approximation varies. Using the deviation in the previous iteration’s real constraints helps eliminate judgment errors caused by differences in the convex approximation between iterations.

If the current solution does not meet the termination condition, the trust region, reference solution, and annealing temperature need to be updated for the next iteration.

The trust region is updated with the following constraints: similar to $∆J_1$, define $∆J_2$:

$$∆J_2=w_t\left({\overline{t}}_f-t_f^{*}\right)+w_z\left(g_h\left({\overline{z}}_g\right)-g_h\left(z_g^{*}\right)\right)+w_p\left(p_u\left(\overline{\mathit{d}}\right)-p_u\left({\mathit{d}}^{*}\right)\right)$$

(22)

Relaxation acceptability can be expressed as:

$$\rho ={\displaystyle \frac{∆J_2}{∆J_1}}$$

(23)

The term $∆J_2$ represents the difference between the real cost of the previous iteration and the current iteration. $\rho$ stands for the ratio of the difference between the real cost and the convex approximation cost. A higher value of $\rho$ indicates a smaller convex approximation deviation, signifying that the current relaxation is reasonable. Conversely, a lower value of $\rho$ suggests a larger convex approximation deviation, implying that the current relaxation might not be precise enough, and that it might be necessary to promptly adjust the trust region size.

The strategy for updating the trust region is as follows:

$$\begin{array}{ccc}\rho ⩾{\rho }_{up}:& \delta =\alpha \delta & update{\left(\overline{\mathit{x}},\overline{\mathit{u}},t\right)}_{\mathit{n}\mathit{e}\mathit{x}\mathit{t}}\\ {\rho }_{mid}⩽\rho <{\rho }_{up}:& \delta =\delta & update{\left(\overline{\mathit{x}},\overline{\mathit{u}},t\right)}_{\mathit{n}\mathit{e}\mathit{x}\mathit{t}}\\ {\rho }_{low}⩽\rho <{\rho }_{mid}:& \delta =\delta /\alpha & update{\left(\overline{\mathit{x}},\overline{\mathit{u}},t\right)}_{\mathit{n}\mathit{e}\mathit{x}\mathit{t}}\\ \rho <{\rho }_{low}:& \delta =\delta /\alpha & {\left(\overline{\mathit{x}},\overline{\mathit{u}},t\right)}_{\mathit{n}\mathit{e}\mathit{x}\mathit{t}}=\left(\overline{\mathit{x}},\overline{\mathit{u}},t\right)\end{array}$$

(24)

where ${\rho }_{low}<{\rho }_{mid}<{\rho }_{up}$ represents the evaluation parameter for relaxation acceptability, and $\alpha$ is a constant greater than 1, indicating the update speed of the trust region. If $\rho ⩾{\rho }_{up}$, it indicates a good convex approximation, and the trust region can be moderately expanded to enhance the solving efficiency; otherwise, if ${\rho }_{mid}⩽\rho <{\rho }_{up}$, it indicates an acceptable range for the convex approximation, and the current trust region is maintained; further, if ${\rho }_{low}⩽\rho <{\rho }_{mid}$, it indicates an inadequate convex approximation, yet the obtained solution is within acceptable limits, so the trust region size is decreased for further solving; finally, if $\rho <{\rho }_{low}$, it indicates a failed convex approximation with an unacceptable solution, meaning that updating the base solution and decreasing the trust region are avoided.

Update reference solution: if the current solution is acceptable, the real cost value of the current solution $J^{\mathrm{*}}$ and the minimum cost value $J_{min}$ are compared. Similar to the SA algorithm, if $J^{\mathrm{*}}<J_{min}$, it means that the current solution is the best, so $J_{min}$ is updated and the current solution is used as the reference solution for the next iteration. Otherwise, the current suboptimal solution is accepted based on P₂ as the reference solution for the next iteration. According to the metropolis criterion:

$$P_2=e^{{\displaystyle \frac{-\left(J^{*}-J_{min}\right)}{T}}}$$

(25)

This means that, at high temperatures, new states with significant energy differences from the current state are acceptable; at low temperatures, only new states with a smaller energy difference from the current state are acceptable. As the temperature gradually decreases, P₂ tends to 0, and the algorithm will use the current best solution as the reference solution for the next iteration.

Update annealing temperature: the annealing temperature gradually decreases with each iteration, which can be expressed as:

$$T=\beta T$$

(26)

where $\beta$ is a constant less than 1.

3.3. The Convergence Analysis of IST-SCP

This section focuses on proving the convergence of IST-SCP. First, the general form of the rendezvous problem is presented. Then, according to the theorems proved in [32], the convergence of the SCP without the initial solution tolerance mechanism is illustrated. Finally, the impact of the initial solution tolerance mechanism on the convergence of the SCP is analyzed.

The original non-convex MFRP can be rewritten in the general form, as Equation (27).

$$\begin{array}{ccc}\mathrm{m}\mathrm{i}\mathrm{n}& J\left(\mathit{x},\mathit{u},t\right)& \end{array}$$

(27a)

$$\begin{array}{ccc}s.t.& f_i\left(\mathit{x},\mathit{u},t\right)=0& \forall i\in {\mathcal{I}}_{eq}\end{array}$$

(27b)

$$\begin{array}{ccc}& f_i\left(\mathit{x},\mathit{u},t\right)⩽0& \forall i\in {\mathcal{I}}_{ineq}\end{array}$$

(27c)

$$\begin{array}{ccc}& h_j\left(\mathit{x},\mathit{u},t\right)=0& \forall j\in {\mathcal{J}}_{eq}\end{array}$$

(27d)

$$\begin{array}{ccc}& h_j\left(\mathit{x},\mathit{u},t\right)⩽0& \forall j\in {\mathcal{J}}_{ineq}\end{array}$$

(27e)

where $J\left(\mathit{x},\mathit{u},t\right)$ in (27a) represents the objective function in (3). $f_i\left(x,u,t\right)$ in (27b) represents the non-convex equality constraints, which are the constraints in (1) and (4). $f_i\left(\mathit{x},\mathit{u},t\right)$, in (27c), represents the non-convex inequality constraints, which are the constraints in (6) and (7). $h_j\left(\mathit{x},\mathit{u},t\right)$, in (27d), represents the convex equality constraints, which are the constraints in (2) and (8). $h_j\left(\mathit{x},\mathit{u},t\right)$, in (27e), represents the convex inequality constraints, which are the constraints in (5). ${\mathcal{I}}_{eq}\cup {\mathcal{I}}_{ineq}$ represents the set of non-convex constraint indices, and ${\mathcal{J}}_{eq}\cup {\mathcal{J}}_{ineq}$ represents the set of convex constraint indices.

By utilizing the L1 penalty function, the original problem in (27) is reformulated as in Equation (28).

$$\begin{array}{cc}\mathrm{m}\mathrm{i}\mathrm{n}& J_{L1}\left(\mathit{x},\mathit{u},t\right)=J\left(\mathit{x},\mathit{u},t\right)+\sum _{i\in {\mathcal{I}}_{eq}}{\lambda }_i\left|f_i\left(\mathit{x},\mathit{u},t\right)\right|+\sum _{i\in {\mathcal{I}}_{ineq}}{\lambda }_i\mathrm{m}\mathrm{a}\mathrm{x}\left(0,f_i\left(\mathit{x},\mathit{u},t\right)\right)\\ & +\sum _{\forall j\in {\mathcal{J}}_{eq}}{\sigma }_j\left|h_j\left(\mathit{x},\mathit{u},t\right)\right|+\sum _{\forall j\in {\mathcal{J}}_{ineq}}{\sigma }_j\mathrm{m}\mathrm{a}\mathrm{x}\left(0,h_j\left(\mathit{x},\mathit{u},t\right)\right)\end{array}$$

(28)

To facilitate the analysis, some fundamental theorems and assumption of convex optimization are provided.

Theorem 1. 

If there exist $\overline{\lambda }⩾0,\overline{\sigma }⩾0$ corresponding to $\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$, such that, for any $\lambda ⩾\overline{\lambda },\sigma ⩾\overline{\sigma }$ within an open neighborhood of $\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$, $J_{L1}\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)⩽J_{L1}\left(\mathit{x},\mathit{u},t\right)$, then $\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ is a local optimum to the problem in Equation (27).

Assumption 1. 

The set of the active inequality constraints indices are defined as:

$$\begin{array}{c}{\mathcal{I}}_{ineq}^{ac}≔\left\{i|f_i\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)=0,i\in {\mathcal{I}}_{ineq}\right\}\\ {\mathcal{J}}_{ineq}^{ac}≔\left\{j|h_j\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)=0,j\in {\mathcal{J}}_{ineq}\right\}\end{array}$$

(29)

The constraints in the MFRP are continuously differentiable. Therefore, define $F_{eq}\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ constructed by $\nabla f_i\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ for $\forall i\in {\mathcal{I}}_{eq}$ as its column, define $H_{eq}\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ constructed by $\nabla h_j\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ for $\forall j\in {\mathcal{J}}_{eq}$ as its column, define $F_{ineq}\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ constructed by $\nabla f_i\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ for $\forall i\in {\mathcal{I}}_{ineq}$ as its column, and define $H_{ineq}\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ constructed by $\nabla h_i\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ for $\forall j\in {\mathcal{J}}_{ineq}$ as its column. Then, the linear independence constraint qualification (LICQ) is satisfied at $\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ if the columns of the following matrix are linearly independent:

$$A_c\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)=\left[\begin{array}{cccc}{F}_{eq}\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)& H_{eq}\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)& F_{ineq}\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)& H_{ineq}\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)\end{array}\right]$$

(30)

Theorem 2. 

If the constraints of the problem in Equation (27) satisfy the LICQ assumption in Assumption 1, and $\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ is a local optimum of the problem in Equation (27), then there exist Lagrange multipliers ${\overline{\mu }}_i$ for $\forall i\in {\mathcal{I}}_{eq}$, ${\overline{\mu }}_i⩾0$ for $\forall i\in {\mathcal{I}}_{ineq}$, ${\overline{\tau }}_j$ for $\forall j\in {\mathcal{J}}_{eq}$, and ${\overline{\tau }}_j⩾0$ for $\forall j\in {\mathcal{J}}_{ineq}$, such that:

$$\begin{array}{c}\nabla J\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)+\sum _{\forall i\in {\mathcal{I}}_{eq}}{\overline{\mu }}_i\nabla f_i\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)+\sum _{\forall i\in {\mathcal{I}}_{ineq}}{\overline{\mu }}_i\nabla f_i\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)\\ +\sum _{\forall j\in {\mathcal{J}}_{eq}}{\overline{\tau }}_j\nabla h_j\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)+\sum _{\forall j\in {\mathcal{J}}_{ineq}}{\overline{\tau }}_j\nabla h_j\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)=0\end{array}$$

(31)

The point that satisfies the above conditions is referred to as the KKT point.

Given the consistency in the formulation of the problem, the convergence of the SCP without the initial solution tolerance mechanism applies to the conclusion drawn in [32]. The convergence of the SCP without the initial solution tolerance mechanism is shown in Theorem 3.

Theorem 3. 

Given Assumption 1, the SCP without the initial solution tolerance mechanism always has limit points, and any limit point, $\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$, is a stationary point of the non-convex penalty problem in Equation (28). Furthermore, if $\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ is feasible for the original non-convex problem in Equation (27), then it is a KKT point of the MFRP.

For IST-SCP, when the annealing temperature $T$ tends to zero, both $P_1$ and $P_2$ tend to zero, which means that the function of searching for and accepting a less optimal solution of the IST-SCP is deactivated. At this point, the IST-SCP will degenerate into the conventional SCP with an initial solution $\left(\overline{\mathit{x}},\overline{\mathit{u}},\overline{t}\right)$ generated by the initial solution tolerance mechanism.

In conclusion, the initial tolerance mechanism does not influence the convergence of the SCP. The sequence solutions generated by the IST-SCP will converge to the local optimal of the original MFRP.

3.4. Multi-Batch Formation Rendezvous Algorithm

The pseudo-codes for the IST-SCP-based formation rendezvous algorithm are shown in the below Algorithm 1.

Algorithm 1: Multi_batch formation rendezvous algorithm

$\mathrm{I}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{z}\mathrm{e}{ \mathit{x}}^q,{\mathit{u}}^q,t^q\to {\overline{\mathit{x}}}^q,{\overline{\mathit{u}}}^q,{\overline{t}}^q$

$\mathrm{I}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{z}\mathrm{e} T^q,{\mathit{\delta }}^q$

$\mathrm{w}\mathrm{h}\mathrm{i}\mathrm{l}\mathrm{e} ∆J_1^q>\epsilon :$

$\hspace{1em} \mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e} P_1^q\to v_o^q=0?$

$\hspace{1em} \mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}\mathrm{e} \mathrm{M}\mathrm{F}\mathrm{R}\mathrm{D}\mathrm{C}{\mathrm{S}}^{\mathrm{q}}\to {\mathit{x}}^{q\mathrm{*}},{\mathit{u}}^{q\mathrm{*}},t^{q\mathrm{*}}$

$\hspace{1em} \mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e} ∆J_1^q,∆J_2^q,{\rho }^q,P_2^q$

$\hspace{1em} \mathrm{u}\mathrm{p}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{e} {\overline{\mathit{x}}}^q,{\overline{\mathit{u}}}^q,{\overline{t}}^q,{\mathit{\delta }}^q,T^q$

$\mathrm{e}\mathrm{n}\mathrm{d}$

$\mathrm{w}\mathrm{h}\mathrm{i}\mathrm{l}\mathrm{e} j<q:$

$\hspace{1em} \mathrm{I}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{z}\mathrm{e} {\mathit{x}}^j,{\mathit{u}}^j,t^j\to {\overline{\mathit{x}}}^j,{\overline{\mathit{u}}}^j,{\overline{t}}^j$

$\hspace{1em} \mathrm{I}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{z}\mathrm{e} T^j,{\mathit{\delta }}^j$

$\hspace{1em} \mathrm{w}\mathrm{h}\mathrm{i}\mathrm{l}\mathrm{e} ∆J_1^j>\epsilon :$

$\hspace{1em}\hspace{1em} \mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e} P_1^j\to v_o^j=0?$

$\hspace{1em}\hspace{1em} \mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}\mathrm{e} \mathrm{M}\mathrm{F}\mathrm{R}\mathrm{D}{\mathrm{S}}^{\mathrm{j}}\to {\mathit{x}}^{j\mathrm{*}},{\mathit{u}}^{j\mathrm{*}},t^{j\mathrm{*}}$

$\hspace{1em}\hspace{1em} \mathrm{c}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e} ∆J_1^j,∆J_2^j,{\rho }^j,P_2^j$

$\hspace{1em}\hspace{1em} \mathrm{u}\mathrm{p}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{e} {\overline{\mathit{x}}}^j,{\overline{\mathit{u}}}^j,{\overline{t}}^j,{\mathit{\delta }}^j,T^j$

$\hspace{1em} \mathrm{e}\mathrm{n}\mathrm{d}$

$\mathrm{e}\mathrm{n}\mathrm{d}$

The Algorithm 1 states that the rendezvous problem of the final batch of UAVs, MFRDCS^q, is solved first. Then, utilizing the results as the reference to the consistency constraints, the rendezvous problems of other batches of UAVs are solved, denoted as MFRDCS^j.

4. Simulation Results and Analysis

In this section, a common scene related to carrier-based UAV formation rendezvous is addressed, validating the effectiveness of the multi-batch formation rendezvous algorithm. Compared with the conventional SCP, the effectiveness of the IST-SCP is verified. Compared with the algorithms that optimize the initial solution (GA+SCP and SA+SCP), the efficiency of the IST-SCP is proved. Additionally, the adaptability of the algorithm concerning (1) the number of batches and the number of UAVs per batch, and (2) the geometric shape of the UAV formation is evaluated.

The conventional SCP method used for comparison was developed by Acikmese et al. [33], who introduced the relaxation factor and the trust region strategy into the SCP algorithm. In recent years, this method has demonstrated commendable performance in various domains such as Mars entry [36], atmospheric re-entry [37], and power system optimization [38].

The simulations are performed in MATLAB with 8 GB of memory and an Intel Core i7 2.2 GHz processor. The convex solver utilized in the SCP is CVX. The parameters of the IST-SCP and the conventional SCP are shown in

Table 1. The parameters of the IST-SCP and SCP.

Parameter	IST-SCP	SCP
$\mathrm{Weight} \mathrm{of} \mathrm{the} \mathrm{rendezvous} \mathrm{time} w_t$	1	1
$\mathrm{Weight} \mathrm{of} \mathrm{the} \mathrm{altitude} \mathrm{penalty} w_z$	3	3
$\mathrm{Weight} \mathrm{of} \mathrm{the} \mathrm{penalty} \mathrm{on} \mathrm{relaxation} \mathrm{factors} w_p$	1 × 102	1 × 102
$\mathrm{The} \mathrm{trust} \mathrm{region} \mathrm{update} \mathrm{rate} \alpha$	2	2
$\mathrm{The} \mathrm{trust} \mathrm{region} \mathrm{evaluation} ({\rho }_{low},{\rho }_{mid},{\rho }_{up})$	(0.1, 0.25, 0.75)	(0.1, 0.25, 0.75)
$\mathrm{The} \mathrm{initial} \mathrm{trust} \mathrm{region} \mathrm{radius} ({\delta }_{x0},{\delta }_{u0},{\delta }_{t0})$	(1 × 102, 10, 10)	(1 × 102, 10, 10)
$\mathrm{The} \mathrm{temperature} \mathrm{update} \mathrm{rate} \beta$	0.6
$\mathrm{The} \mathrm{range} \mathrm{of} \mathrm{temperature} (T_{min},T_{max})$	(0.1, 1 × 102)

.

The scenario of carrier-based UAV formation rendezvous is constructed as: a carrier starts from the coordinate origin (0,0) and travels northeast at a constant speed of V = 42 km/h. During this process, every 2 min, two carrier-based UAVs simultaneously take off. In total, six isomorphic UAVs with a wingspan of 18 m take off. They are directed to a specified altitude of 600 m, at a specified circling angular velocity of 0.04 rad/s, with a circling radius of 3 km around the designated formation center. The objective is to form a regular hexagonal formation with a center-to-center distance of 24 m between the aircraft. During the formation rendezvous process of UAVs, it is necessary to avoid cylindrical no-fly zones. The detailed parameters and rendezvous scenario are illustrated in

Table 2. Task setting parameters.

Parameter	Value
Start point of carrier	(0 m,0 m)
Velocity of carrier	42 km/h
Time interval for takeoff between batches	2 min
Number of UAVs per batch	2
Number of batches	3
Rendezvous altitude	600 m
Circling angular velocity	0.04 rad/s
Radius of rendezvous orbit	3 km
Center of rendezvous orbit	(11,880 m, 16,120 m)
Geometric shape of formation	Regular hexagonal
Center-to-center distance in the formation	24 m
Wingspan of UAV	18 m

and Figure 7.

The velocity of the carrier is that used in [11], which is around 22 knots. The time interval for takeoff between the batches is that used in [39]. The center-to-center distance in the formation is that used in [40], and adopts a larger value to leave enough space for circling maneuvers.

The green line in the figure represents the trajectory of the carrier’s navigation. The blue circles denote the predetermined circling orbits of the formation centroid. The green squares represent the locations where the UAVs from three different batches take off. The red circular areas delineate the no-fly zones for the carrier-based UAVs.

4.1. Effectiveness Analysis of the Multi-Batch Formation Rendezvous Algorithm

The algorithm 1 proposed in Section 3.4 is applied in the aforementioned task scene. Because of the probability-based search module and metropolis criterion, the algorithm 1 is calculated five times, and takes the shortest rendezvous time as the final result, which is shown in Figure 8.

Figure 8a illustrates the 3D flight trajectory of the six UAVs during the formation rendezvous process. It is noticed that, after sequentially taking off in three batches, the aircraft climbed to the designated altitude because of the altitude penalty term in (3). Then, they avoided the cylindrical no-fly zones. Eventually, they reached the designated formation position and executed the circling flight, forming the intended formation.

Figure 8b displays a projection of the flight trajectories of the six UAVs on the horizontal plane. It indicates that all six UAVs found the optimal paths from their starting positions to the final destinations. Under the constraints of consistency, they ultimately formed the specified regular hexagonal formation.

The effectiveness of the multi-batch formation rendezvous algorithm is well proved.

To highlight the effectiveness of IST-SCP, the optimizer in the multi-batch formation rendezvous algorithm is replaced by conventional SCP. The results are shown in

Table 3. Comparison results of IST-SCP and conventional SCP.

	The Type of Time Node	IST-SCP/s	Conventional SCP/s
1st batch	Takeoff	0
	Rendezvous	203.05	213.80
2nd batch	Takeoff	120
	Rendezvous	315.86	336.99
3rd batch	Takeoff	240
	Rendezvous	405.63	428.46

.

It is obvious that, for each batch of UAVs, the rendezvous time generated by IST-SCP is shorter than conventional SCP.

The rendezvous time of each UAV is related to its path and velocity. The velocities during the rendezvous process of IST-SCP and conventional SCP are shown in Figure 9.

It is shown that the velocity calculated by both algorithms first increases and then decreases due to the formation’s circling velocity, except for that of the climbing process at the beginning. This means that both algorithms can generate the maximum velocity of the UAVs flying on this path. Therefore, the difference in rendezvous time is mainly caused by the path.

As shown in Figure 10, the IST-SCP jumps out of the local optimum caused by inaccurate initial solutions and obtains a path with a shorter rendezvous time.

Taking the numerical results in Figure 10b as an example, the optimization details of the IST-SCP are presented. The initial solution of the second batch is on the left sides of no-fly zone 1 and no-fly zone 3. At this point, the annealing temperature was at its maximum, $T=T_{max}$, and so was $P_1=0.99$. Thus, the subproblem was more likely to disregard the influence of the relaxation factor ${\mathit{v}}_o$ when being solved, searching for the optimal rendezvous time and location within the entire task space. However, the trajectory obtained in this way had a greater obstacle avoidance penalty. Under the influence of the greatest $P_2$, the IST-SCP was encouraged to search even when the obstacle avoidance penalty was high. As the annealing temperature decreased, $P_1$ gradually decreased, and the probability of considering the relaxation factor ${\mathit{v}}_o$ in solving the subproblems increased. Simultaneously, $P_2$ was also decreasing, leading the algorithm to favor solutions with a lesser obstacle avoidance penalty as the new reference solution. This continued until the annealing temperature dropped to its minimum,$T=T_{min}$. At this point, the reference trajectory was on the left side of no-fly zone 2 and was on the right side of no-fly zone 3. The IST-SCP fully degenerated into conventional SCP, optimizing the rendezvous time and location until the algorithm’s stopping criteria were met.

To better demonstrate the reduced sensitivity of the IST-SCP to the initial solutions, computational results for the rendezvous time under various initial solutions for the second batch of UAVs are presented in

Table 4. Comparison of the rendezvous time under various initial rendezvous times.

Initial Rendezvous Time/s	Rendezvous Time/s of IST-SCP	Rendezvous Time/s of SCP
320	314.52	314.37
340	315.86	336.99
360	315.93	323.31
380	318.25	357.68

. It can be observed that the initial solution tolerance mechanism effectively reduces the sensitivity to the initial values.

The rendezvous time generated by the IST-SCP is reduced by up to 11.02%, and on average by 8.40%, in comparison to the SCP. The effectiveness of the multi-batch formation rendezvous algorithm has been well demonstrated.

4.2. Efficiency Analysis of the Multi-Batch Formation Rendezvous Algorithm

To highlight the efficiency of IST-SCP, the SA+SCP and GA+SCP methods, which optimize the initial solution, are compared with IST-SCP in computing time.

The SA+SCP and GA+SCP methods contain two layers. The first layer, using SA/GA, takes the initial rendezvous position as the optimization variable. Then, the second layer, using conventional SCP, calculates the cost value for each initial solution and converts it back to the first layer to find the best initial rendezvous position for the conventional SCP.

The rendezvous time and computing time for from six to eight UAVs were simulated. The results are shown in

Table 5. Comparison of the computing time between IST-SCP, GA+SCP, and SA+SCP.

	Rendezvous Time/s	Computing Time/s
	6 UAVs in 3 batches
IST-SCP	405.63	485.99
GA+SCP	400.34	14,266.16
SA+SCP	400.07	16,948.33
	7 UAVs in 3 batches
IST-SCP	404.85	528.26
GA+SCP	405.27	15,654.69
SA+SCP	406.13	18,595.76
	8 UAVs in 4 batches
IST-SCP	520.27	647.70
GA+SCP	524.49	19,415.76
SA+SCP	522.96	22,597.06

. Compared with the GA+SCP, the computing time of the IST-SCP is reduced by 96.60% on average. Compared with the SA+SCP, the computing time of the IST-SCP is reduced by 97.14% on average. It is obvious that IST-SCP, as an initial solution tolerance algorithm, has higher efficiency than such initial solution optimization algorithms.

4.3. Adaptability Analysis of Multi-Batch Formation Rendezvous Algorithm

This section will prove the strong adaptability of the multi-batch formation rendezvous algorithm concerning the following factors: (1) the number of batches and the number of UAVs per batch, and (2) the geometric shape of the UAV formation.

(1) The number of batches and the number of UAVs per batch:

The solution for (1) four batches, each with two UAVs, totaling eight UAVs, to form a regular octagon, (2) three batches, with one UAV in the first batch, two UAVs in the second batch, and three UAVs in the third batch, to form a regular hexagon, (3) and three batches, with one UAV in the first batch, three UAVs in the second batch, and two UAVs in the third batch, to form a regular hexagon, are shown in Figure 11.

The adaptability to the number of batches and the number of UAVs per batch of the multi-batch formation rendezvous algorithm is well demonstrated.

(2) The geometric shape of the UAV formation:

The solution for (1) three batches, with one UAV in the first batch and two UAVs in the second batch and third batches, to form a diamond-shaped formation, and for (2) three batches, with one UAV in the first batch and three UAVs in the second third batches, to form a herringbone formation are shown in Figure 12.

The adaptability to the geometric shape of the multi-batch formation rendezvous algorithm is well proved.

In conclusion, it is well proved that the multi-batch formation rendezvous algorithm has strong adaptability concerning the number of batches, the number of UAVs per batch, and the geometric shape of the formation.

5. Conclusions

In this paper, the multi-batch formation rendezvous problem for carrier-based UAVs is studied.

(1) An innovative initial solution tolerance SCP (IST-SCP) algorithm is proposed to solve the problem that the conventional SCP algorithm is easily affected by the initial solution and falls into the local optimum. By introducing a probability-based search module and metropolis criterion, the IST-SCP algorithm can accept inappropriate initial solutions and jump out of the local optimum. As the convergency of the IST-SCP is well proved, the IST-SCP algorithm can be applied to a general spatiotemporal coupled trajectory planning problem. Then, utilizing IST-SCP, the multi-batch formation rendezvous method is developed to solve the multi-batch carrier-based UAV formation rendezvous problem.

(2) A new type of UAV rendezvous problem is solved for the first time, considering all kinds of rendezvous limitations of carrier-based UAVs. A reverse solution approach is presented for modeling the optimization problems, which simplifies the spatiotemporal coupling consistency constraints without compromising the optimality of the problem. Then, considering the non-convexity arising from constraints such as the UAVs’ dynamics, obstacle avoidance specifications and convex relaxation techniques are introduced to convexify the original problem.

(3) Simulation experiments are carried out to validate the effectiveness of the multi-batch formation rendezvous algorithm. It is proved that IST-SCP generates better solutions than the conventional SCP and has higher efficiency than the initial solution optimization method. The adaptability of the proposed algorithm is also verified in terms of the following factors: the number of batches, the number of UAVs per batch, and the geometric shape of the formation.

In future studies, we will prioritize considerations of wind field models, the carrier motion, environmental uncertainties, communication delays, optimizations in energy consumption, and computational loads, alongside the challenges inherent in practical deployments.