Research on Cooperative Formation Path Planning for Unmanned Aerial Vehicles Based on Cell-like Membrane Computing

1. Introduction

Unmanned aerial vehicles (UAVs), characterized by their flexible structural design, high maneuverability, compact size, and cost-effectiveness, have found widespread applications across diverse operational domains. In recent years, with the rapid advancement of the low-altitude economy, UAV technology has progressed remarkably, transitioning from simple single-vehicle remote control to sophisticated multi-vehicle cooperative operations. A single UAV, constrained by its payload capacity, endurance limits, and sensor field of view, often struggles to efficiently complete large-scale, highly dynamic, or heterogeneous composite missions. By contrast, multi-UAV systems, through coordinated teamwork, can achieve rapid coverage of mission areas, distributed load sharing, and complementary functional enhancement, thereby significantly improving task execution efficiency, robustness, and success rates.

Key application scenarios of multi-UAV systems are summarized in Table 1.

Typical Applications of Multi-UAV Systems

The research landscape of UAVs has evolved into a highly interdisciplinary field integrating artificial intelligence, communication networks, and advanced sensor technologies. The core research directions include autonomous navigation and obstacle avoidance, formation keeping and transformation control, and cooperative task allocation. Among these, cooperative formation path planning is the foundational technology that directly determines the safety and mission effectiveness of UAV swarm operations. The essence of cooperative formation path planning is to plan feasible or optimal paths for each UAV within a formation under multiple constraints, including kinematic limits, energy consumption, obstacle avoidance, and timing sensitivity, while maintaining a predefined geometric formation through stable cooperative flight. This challenging problem encompasses two interlinked aspects: the generation of globally or locally optimal reference paths satisfying mission objectives and safety constraints, and the formation control that maintains the desired formation shape amid external disturbances.

Existing path planning algorithms for UAVs often falter in highly complex environments with dense obstacles and dynamic threats. Traditional optimization algorithms like particle swarm optimization and genetic algorithms are prone to premature convergence and exhibit slow convergence rates. While fast heuristic methods cannot guarantee global optimality, machine learning-based approaches require extensive labeled trajectory samples, which are expensive and risky to obtain in real-world flight scenarios. In formation planning, the inherent coupling between formation keeping and obstacle avoidance creates a fundamental conflict: over-emphasizing formation geometry may lead to delayed collision avoidance, while prioritizing avoidance may compromise formation integrity or induce inter-UAV collisions.

Remarkably, the biological information processing mechanisms observed in nature offer an alternative computation paradigm. Membrane computing, also known as P systems, was pioneered by Păun in 1998 as a branch of natural computing inspired by the structure and function of biological cells. The hierarchical compartmentalization and parallel chemical reactions within cells motivate the design of novel computational models capable of capturing the dynamic regulatory mechanisms inherent in living systems. Membrane systems are broadly classified into three fundamental models: cell-like P systems, tissue-like P systems, and neural-like P systems. Among these, cell-like P systems abstract the internal biological structure of an individual cell, with a hierarchical tree-like topology of nested membranes. The integration of the spiking mechanism from spiking neural P systems into cell-like P systems yields a hybrid model that combines layered nesting structure with pulse-based communication, offering an ideal computational framework for modeling the parallel working scenarios of multi-UAV systems.

This thesis systematically investigates cooperative formation path planning for unmanned aerial vehicles based on cell-like membrane computing. The organization of the work comprises two main research thrusts: first, the proposal of a hybrid algorithm integrating a spiking cell-like membrane system with an improved grey wolf optimizer for multi-UAV cooperative path planning in three-dimensional environments; second, the design of a cell-like membrane system-based formation path planning method that combines a virtual leader architecture with an improved artificial potential field approach.

2. Theoretical Foundations

2.1 Coordinate Systems and Transformations

The precise description of UAV spatial positions and attitudes requires well-defined coordinate systems. This study adopts the Earth inertial coordinate system (NED) for global position representation and the body-fixed coordinate system for attitude description. The rotation matrix from the body coordinate system to the Earth inertial coordinate system is expressed in Eq. (1), using the Z-Y-X Euler angle convention composed of roll \(\varphi\), pitch \(\theta\), and yaw \(\psi\).

\[
R_{p} = \begin{bmatrix}
\cos\theta \cos\psi & \sin\varphi \sin\theta \cos\psi – \cos\varphi \sin\psi & \cos\varphi \sin\theta \cos\psi + \sin\varphi \sin\psi \\
\cos\theta \sin\psi & \sin\varphi \sin\theta \sin\psi + \cos\varphi \cos\psi & \cos\varphi \sin\theta \sin\psi – \sin\varphi \cos\psi \\
-\sin\theta & \sin\varphi \cos\theta & \cos\varphi \cos\theta
\end{bmatrix}
\tag{1}
\]

2.2 Quadrotor UAV Dynamic Model

Under the assumptions of ideal rigid body, symmetric structure, uniform mass distribution, and negligible air resistance, the translational dynamics of a quadrotor UAV can be described by:

\[
\ddot{x} = \frac{1}{m} F_{t} (\cos\varphi \sin\theta \cos\psi + \sin\varphi \sin\psi)
\tag{2}
\]

\[
\ddot{y} = \frac{1}{m} F_{t} (\cos\varphi \sin\theta \sin\psi – \sin\varphi \cos\psi)
\tag{3}
\]

\[
\ddot{z} = \frac{1}{m} F_{t} \cos\varphi \cos\theta – g
\tag{4}
\]

The rotational dynamics are governed by:

\[
\dot{p} = \frac{I_{y} – I_{z}}{I_{x}} qr + \frac{l b}{I_{x}} (\omega_{4}^{2} – \omega_{2}^{2})
\tag{5}
\]

\[
\dot{q} = \frac{I_{z} – I_{x}}{I_{y}} pr + \frac{l b}{I_{y}} (\omega_{3}^{2} – \omega_{1}^{2})
\tag{6}
\]

\[
\dot{r} = \frac{I_{x} – I_{y}}{I_{z}} pq + \frac{d}{I_{z}} (\omega_{1}^{2} – \omega_{2}^{2} + \omega_{3}^{2} – \omega_{4}^{2})
\tag{7}
\]

2.3 Path Planning Constraints and Evaluation Functions

The path planning problem for unmanned aerial vehicles requires the simultaneous satisfaction of multiple constraints. The key constraints considered in this thesis are summarized in Table 2.

Key Constraints in Multi-UAV Path Planning

| Constraint Type | Mathematical Expression | Description |
|—|—|—|
| Maximum flight distance | \(\sum_{i=1}^{n} l_{i} \leq L_{max}\) | Energy limitation |
| Flight altitude | \(H_{min} \leq h_{i} \leq H_{max}\) | Safety and performance |
| Horizontal turning angle | \(\alpha_{i} \leq \alpha_{max}\) | Maneuverability |
| Pitch angle | \(\psi_{i} \leq \psi_{max}\) | Maneuverability |
| Threat avoidance | \(d_{i} \geq R_{j}\) | Obstacle clearance |
| Inter-UAV collision | \(d_{ab} \geq d_{min}\) | Safety distance |
| Time coordination | \(\bigcap_{i} [t_{i,start}, t_{i,end}]
eq \emptyset\) | Temporal synchronization |

The comprehensive cost function \(J\) combining flight distance, altitude, turning angle, threat, collision, and time coordination costs is formulated as:

\[
J = w_{1}F_{1} + w_{2}F_{2} + w_{3}F_{3} + w_{4}F_{4} + w_{5}F_{5} + w_{6}F_{6}
\tag{8}
\]

where \(\sum_{i=1}^{6} w_{i} = 1\). The individual cost components are defined as:

\[
F_{1} = k_{1} \sum_{i=1}^{n} l_{i}
\tag{9}
\]

\[
F_{2} = \begin{cases}
k_{2}(H_{min} – h_{i}), & h_{i} < H_{min} \\
0, & H_{min} \leq h_{i} \leq H_{max} \\
k_{2}(h_{i} – H_{max}), & h_{i} > H_{max}
\end{cases}
\tag{10}
\]

\[
F_{3} = k_{3}\sum_{i=1}^{n} (f_{\alpha_{i}} + f_{\psi_{i}})
\tag{11}
\]

\[
F_{4} = k_{4}\sum_{j=1}^{m} p_{d_{j}} (d_{j} – R_{j}), \quad d_{j} \leq R_{j}
\tag{12}
\]

2.4 Cell-like Membrane Computing Foundations

A cell-like P system is formally defined as the tuple:

\[
\Pi = (O, \mu, \omega_{1}, \dots, \omega_{m}, R_{1}, \dots, R_{m}, i_{o})
\tag{13}
\]

where \(O\) is a non-empty finite alphabet, \(\mu\) is the membrane structure of degree \(m\), \(\omega_{i}\) are multisets of objects in membrane \(i\), \(R_{i}\) are finite sets of evolution rules of the form \(u \to v\), and \(i_{o}\) is the output region. The spiking cell-like P system variant is formally expressed as:

\[
\Pi = (O, H, \mu, \sigma_{1}, \dots, \sigma_{m}, R_{1}, \dots, R_{m}, syn, i_{o})
\tag{14}
\]

where each neuron \(\sigma_{i}\) contains a finite set of spikes and rules of the form \(E / a^{c} \to a\), with \(E\) being a regular expression over \(a\), and \(syn \subseteq \{1,2,\dots,m\} \times \{1,2,\dots,m\}\) defines the synaptic connections between membranes.

The inherent parallelism of membrane systems provides a natural computational framework for multi-UAV cooperative control, where each membrane can be mapped to an independent UAV agent, and communication rules between membranes emulate inter-UAV information exchange.

3. Multi-UAV Cooperative Path Planning Using CSP-GWO

3.1 Standard Grey Wolf Optimizer

The grey wolf optimizer (GWO) is a meta-heuristic algorithm inspired by the hunting behavior and social hierarchy of grey wolf packs. The social hierarchy consists of four levels: \(\alpha\), \(\beta\), \(\delta\), and \(\omega\), corresponding to the leader, sub-leader, third-level, and remaining wolves, respectively. The mathematical model of GWO involves three main phases:

Encircling prey:

\[
\vec{D} = |\vec{C} \cdot \vec{x}_{p}(t) – \vec{x}(t)|
\tag{15}
\]

\[
\vec{x}(t+1) = \vec{x}_{p}(t) – \vec{A} \cdot \vec{D}
\tag{16}
\]

where \(\vec{A} = 2a \cdot \vec{r}_{1} – a\) and \(\vec{C} = 2\vec{r}_{2}\). The parameter \(a\) decreases linearly from 2 to 0 during iterations.

Hunting: The positions of the \(\alpha\), \(\beta\), and \(\delta\) wolves are used to guide the remaining wolves:

\[
\vec{x}(t+1) = \frac{1}{3} (\vec{x}_{1} + \vec{x}_{2} + \vec{x}_{3})
\tag{17}
\]

where \(\vec{x}_{1}\), \(\vec{x}_{2}\), and \(\vec{x}_{3}\) denote the adjusted positions based on the \(\alpha\), \(\beta\), and \(\delta\) wolves, respectively.

3.2 Tent Map Initialization

Standard GWO uses random initialization, which may lead to uneven distribution of the initial population and compromise the exploitation of the solution space. To address this limitation, this thesis introduces the Tent chaotic map for population initialization. The Tent map is defined as:

\[
x_{n+1} = \begin{cases}
r \cdot x_{n}, & x_{n} \leq 0.5 \\
r \cdot (1 – x_{n}), & x_{n} > 0.5
\end{cases}
\tag{18}
\]

Through systematic experiments with \(r \in \{1.1, 1.2, \dots, 2.0\}\), the value \(r = 1.5\) achieves the optimal balance between ergodicity and uniformity, outperforming random initialization and the Logistic map.

3.3 Dynamic Fuzzy Weighting Strategy

The traditional GWO uses the average of \(\alpha\), \(\beta\), and \(\delta\) wolves’ positions, which fails to reflect the social hierarchy properly. A Mamdani-type fuzzy inference system with five triangular membership functions is designed to generate dynamic weights \(W_{1}\), \(W_{2}\), and \(W_{3}\) for the \(\alpha\), \(\beta\), and \(\delta\) wolves, respectively. The position update equation becomes:

\[
\vec{x}(t+1) = \frac{W_{1}\vec{x}_{1} + W_{2}\vec{x}_{2} + W_{3}\vec{x}_{3}}{W_{1} + W_{2} + W_{3}}
\tag{19}
\]

The fuzzy rules are designed as follows:

| Rule | Iteration Phase | \(W_{1}\) | \(W_{2}\) | \(W_{3}\) |
|—|—|—|—|—|
| R1 | Early | Medium | Medium | Medium |
| R2 | Early-Mid | High | Medium | Low-Medium |
| R3 | Middle | High | Medium | Low |
| R4 | Mid-Late | Higher | Low-Medium | Low |
| R5 | Late | Highest | Low | Lowest |

3.4 CSP-GWO Framework Design

The CSP-GWO framework integrates cell-like membrane computing with the improved GWO through a hierarchical parallel architecture. Each elementary membrane within the membrane system hosts an independent GWO subpopulation. Tent map initialization is applied in parallel across all elementary membranes, preserving population diversity and preventing premature convergence. After each iteration, optimal solutions \(g_{best,i}\) from each membrane are transmitted via spiking rules to the outer membrane, where global competition identifies the current \(\alpha\) wolf. The global best solution is then copied back into each elementary membrane to guide subsequent evolution.

The complete algorithm procedure is outlined in Table 3.

CSP-GWO Algorithm Procedure

| Step | Description |
|—|—|
| S1 | Initialize the number of elementary membranes \(m\), maximum iterations, and dimension \(dim\) |
| S2 | Apply Tent map to initialize grey wolf populations in parallel within each membrane |
| S3 | Sort populations by fitness, identify the global best individual, and feedback to all membranes |
| S4 | Calculate position updates using fuzzy dynamic weights |
| S5| Transmit optimal populations via spiking rules to other membranes and update the system global best |
| S6 | Check termination criteria; if not met, return to S3; otherwise output the optimal solution |

3.5 Experimental Results on Benchmark Functions

To evaluate the performance of CSP-GWO, comprehensive experiments were conducted using 20 benchmark functions including unimodal functions (F1–F7), multimodal functions (F8–F12), and fixed-dimension multimodal functions (F13–F20). CSP-GWO was compared against standard GWO, SOGWO, MP-GWO, particle swarm optimization (PSO), whale optimization algorithm (WOA), ant lion optimizer (ALO), and sine cosine algorithm (SCA). All algorithms used a population size of 50, 500 iterations, and 25 independent runs. The comparative results for selected functions are presented in Table 4.

Comparison of Algorithm Performance on Benchmark Functions

| Function | Metric | GWO | CSP-GWO | SOGWO | MP-GWO | PSO |
|—|—|—|—|—|—|—|
| F1 | Mean | 3.35E-33 | 6.54E-55 | 3.86E-24 | 3.86E-24 | 3.86E-24 |
| F3 | Mean | 2.59E-08 | 2.78E-14 | 3.33E-08 | 1.02E-09 | 46.606 |
| F5 | Mean | 26.714 | 25.755 | 26.735 | 28.218 | 83.961 |
| F8 | Mean | 1.552 | 2.27E-15 | 1.390 | 151.524 | 46.459 |
| F10 | Mean | 6.02E-03 | 0 | 2.11E-03 | 4.22E-03 | 9.16E-03 |
| F11 | Mean | 0.026 | 0 | 0.025 | 9.090 | 1.92E-07 |
| F14 | Mean | 1.90E-03 | 3.10E-04 | 1.21E-03 | 6.17E-03 | 9.10E-04 |
| F18 | Mean | -3.86195 | -3.86278 | -3.86208 | -3.86220 | -3.86278 |

The convergence curves demonstrate that CSP-GWO achieves the fastest convergence speed and highest final accuracy on most test functions. The superior performance can be attributed to the synergistic combination of Tent chaotic map initialization, dynamic fuzzy weighting, and the membrane system’s parallel architecture.

3.6 Path Planning Experiments

To validate the practical applicability of CSP-GWO, path planning experiments were conducted in a three-dimensional environment of dimensions 20 km × 20 km × 20 km. The environment contains 16 cylindrical no-fly zones with parameters listed in Table 5.

Threat Zone Information (km)

| Threat Center | Radius | Threat Center | Radius |
|—|—|—|—|
| (4, 11, 12) | 2.5 | (6, 12, 8) | 2.1 |
| (10, 2, 10) | 2.0 | (12, 4, 4) | 1.6 |
| (5, 6, 10) | 1.5 | (3, 7, 5) | 1.3 |
| (16, 16, 0) | 4.0 | (19, 17, 3) | 3.0 |
| (8, 5, 11) | 1.4 | (7, 1, 5) | 1.6 |
| (11, 3, 15) | 2.7 | (11, 13, 9) | 1.2 |
| (20, 6, 9) | 2.3 | (17, 15, 7) | 2.4 |
| (16, 4, 14) | 2.1 | (5, 12, 2) | 1.1 |

B-spline curves are employed for trajectory smoothing to satisfy kinematic constraints:

\[
C(u) = \sum_{i=0}^{n} N_{i,k}(u) \cdot P_{i}
\tag{20}
\]

where \(N_{i,k}(u)\) are the B-spline basis functions defined recursively by the Cox-de Boor formula.

Path planning experiments were conducted for 4, 8, and 16 UAVs. The experimental results are summarized in Tables 6-8.

Path Planning Results for Four UAVs

| Algorithm | Average Cost | Variance | Collisions | Time Violations | Runtime (s) |
|—|—|—|—|—|—|
| GWO | 2.402 | 0.1064 | 0 | 0 | 13.6 |
| SOGWO | 2.365 | 0.075 | 0 | 0 | 12.43 |
| MP-GWO | 2.636 | 0.066 | 0 | 0 | 20.45 |
| PSO | 2.601 | 0.061 | 0 | 0 | 13.12 |
| WOA | 2.760 | 0.062 | 0 | 0 | 16.14 |
| CSP-GWO | 2.098 | 0.041 | 0 | 0 | 11.14 |

Path Planning Results for Eight UAVs

| Algorithm | Average Cost | Variance | Collisions | Time Violations | Runtime (s) |
|—|—|—|—|—|—|
| GWO | 5.08 | 0.10 | 1 | 0 | 26.4 |
| SOGWO | 4.89 | 0.11 | 0 | 0 | 24.6 |
| MP-GWO | 5.11 | 0.13 | 0 | 0 | 38.8 |
| PSO | 5.15 | 0.15 | 1 | 0 | 25.2 |
| WOA | 5.22 | 0.16 | 0 | 0 | 31.4 |
| CSP-GWO | 4.31 | 0.09 | 0 | 0 | 11.6 |

Path Planning Results for Sixteen UAVs

| Algorithm | Average Cost | Variance | Collisions | Time Violations | Runtime (s) |
|—|—|—|—|—|—|
| GWO | 9.42 | 0.25 | 2 | 2 | 49.3 |
| SOGWO | 9.38 | 0.23 | 2 | 0 | 48.2 |
| MP-GWO | 9.84 | 0.21 | 1 | 2 | 73.2 |
| PSO | 10.17 | 0.21 | 1 | 1 | 49.0 |
| WOA | 10.21 | 0.26 | 1 | 2 | 60.7 |
| CSP-GWO | 8.32 | 0.10 | 0 | 0 | 12.8 |

The experimental results clearly demonstrate the superiority of CSP-GWO. For the sixteen-UAV scenario, CSP-GWO reduced the average path cost by 11.3% compared to SOGWO, maintained a zero collision rate and zero time coordination violations, and reduced the average runtime by 73.4% compared to SOGWO. The near-constant runtime of CSP-GWO across increasing UAV scales highlights the benefits of the membrane system’s parallel computing architecture, demonstrating excellent scalability for large-scale UAV swarm applications.

4. Multi-UAV Formation Path Planning Based on Cell-like Membrane Systems

4.1 Virtual Leader Formation Model

The virtual leader method abstracts the reference benchmark as a mathematical entity rather than a physical UAV, fundamentally eliminating the single-point failure risk inherent in the classical leader-follower approach. The virtual leader defines the desired formation configuration, while each UAV tracks its expected position relative to this virtual point.

The desired position of the \(i\)-th UAV is given by:

\[
\vec{p}_{i}^{des} = \vec{p}_{V} + R(\psi_{V}) \cdot \vec{r}_{i}, \quad i = 1, 2, \dots, N
\tag{21}
\]

where \(\vec{p}_{V}\) is the virtual leader position, \(\psi_{V}\) is its heading angle, \(\vec{r}_{i}\) is the expected relative position vector, and \(R(\psi_{V})\) is the rotation matrix expressed as:

\[
R(\psi_{V}) = \begin{bmatrix}
\cos\psi_{V} & -\sin\psi_{V} & 0 \\
\sin\psi_{V} & \cos\psi_{V} & 0 \\
0 & 0 & 1
\end{bmatrix}
\tag{22}
\]

The expected velocity including the angular velocity compensation is:

\[
\vec{v}_{i}^{des} = \vec{v}_{V} + \vec{\omega}_{V} \times (R(\psi_{V}) \cdot \vec{r}_{i})
\tag{23}
\]

4.2 Cell-like Membrane System Formation Model

The cell-like membrane system formation model maps the virtual leader logic onto a parallel membrane computing architecture. The formation membrane system is formally defined as:

\[
\Pi_{UAV} = (O, \mu, \omega_{1}, \dots, \omega_{n}, R_{1}, \dots, R_{i}, i_{o})
\tag{24}
\]

In this mapping, the surface membrane corresponds to the virtual leader and its operating environment, storing global mission objective objects, obstacle objects, and the virtual leader’s state objects. Each elementary membrane represents an individual UAV, encapsulating its state objects and local error evolution logic. The communication rules between membranes correspond to the broadcasting of virtual leader state information. The system operates under the principle of maximally parallel execution: in each computation step, all elementary membranes simultaneously trigger their internal control law evolution rules.

The formation membrane system executes five mission phases through different evolution rules, as summarized in Table 9.

Multi-Stage Formation Control Rules

The Bézier curve smoothing for the takeoff-to-cruise transition is given by:

\[
\vec{p}_{L}(\tau) = (1-\tau)^{3}\vec{P}_{0} + 3(1-\tau)^{2}\tau\vec{P}_{1} + 3(1-\tau)\tau^{2}\vec{P}_{2} + \tau^{3}\vec{P}_{3}
\tag{25}
\]

4.3 Improved Artificial Potential Field Method

The traditional artificial potential field (APF) method suffers from two critical drawbacks: local minima problems and the unreachable target problem (GNRON). To overcome these limitations, this thesis introduces an improved repulsive potential function that incorporates a distance adjustment factor:

\[
U_{rep}(q) = \begin{cases}
\frac{1}{2}\eta \left(\frac{1}{\rho(q, q_{obs})} – \frac{1}{\rho_{0}}\right)^{2} \rho^{n}(q, q_{d}), & \rho(q, q_{obs}) \leq \rho_{0} \\
0, & \rho(q, q_{obs}) > \rho_{0}
\end{cases}
\tag{26}
\]

The distance adjustment factor \(\rho^{n}(q, q_{d})\) ensures that the repulsive potential smoothly decreases to zero as the UAV approaches the target point, thereby eliminating the GNRON problem.

Furthermore, a distance-weighted dynamic control allocation strategy is introduced to resolve the priority conflict between formation keeping and autonomous obstacle avoidance:

\[
\vec{u}_{form} = k_{p}(\vec{q}_{L} + \vec{r}_{i} – \vec{q}_{i}) + k_{v}(\vec{v}_{L} – \vec{v}_{i})
\tag{27}
\]

The weight factor \(\alpha(t)\) is designed using a cubic spline to achieve smooth transition between formation and obstacle avoidance priorities:

\[
\alpha(t) = \begin{cases}
1, & d_{m} > \rho_{0} \\
\alpha_{min} + (1-\alpha_{min})(3t^{2} – 2t^{3}), & d_{min} \leq d_{m} \leq \rho_{0} \\
\alpha_{min}, & d_{m} < d_{min}
\end{cases}
\tag{28}
\]

The final synthesized control law is:

\[
\vec{u}_{i} = Sat(\alpha \cdot \vec{u}_{form} + (1-\alpha)\vec{u}_{safe} – \beta_{global}\vec{v} – \beta_{local}\vec{v}_{\perp})
\tag{29}
\]

where \(Sat(\cdot)\) is the acceleration limiting function.

4.4 Stability and Convergence Analysis

The stability analysis of the proposed formation control law is conducted using Lyapunov stability theory. For the takeoff assembly and mission maneuver phases (S1/S2), the Lyapunov function is constructed as:

\[
V_{1} = \sum_{i=1}^{N} \left(\frac{1}{2}k_{p}\vec{e}_{pi}^{T}\vec{e}_{pi} + \frac{1}{2}\vec{e}_{vi}^{T}\vec{e}_{vi}\right)
\tag{30}
\]

Taking the derivative of \(V_{1}\) along the closed-loop trajectory yields:

\[
\dot{V}_{1} = \sum_{i=1}^{N} \left(k_{p}\vec{e}_{pi}^{T}\vec{e}_{vi} + \vec{e}_{vi}^{T}(-k_{p}\vec{e}_{pi} – k_{v}\vec{e}_{vi})\right) = -\sum_{i=1}^{N} k_{v}\|\vec{e}_{vi}\|^{2} \leq 0
\tag{31}
\]

By LaSalle’s invariance principle, the system states converge to the largest invariant set where \(\dot{V}_{1} = 0\), which implies \(\vec{e}_{vi} = 0\) and consequently \(\vec{e}_{pi} = 0\). Therefore, the formation system is globally asymptotically stable at the equilibrium point.

For the obstacle avoidance phase, the input-to-state stability (ISS) of the system is established. The total potential field function serves as the Lyapunov function:

\[
V_{2} = \frac{1}{2}\vec{e}_{v}^{T}\vec{e}_{v} + U_{att}(\rho) + U_{rep}
\tag{32}
\]

The derivative satisfies:

\[
\dot{V}_{2} = -(\beta_{g} + \beta_{l})\|\vec{v}\|^{2} \leq 0
\tag{33}
\]

which confirms that the improved APF-based controller guarantees ISS stability under disturbance inputs.

4.5 Full-Process Simulation Results

To validate the effectiveness of the proposed formation control method, full-process simulation experiments were conducted in a three-dimensional environment of dimensions 260 m × 160 m × 50 m, with four UAVs executing a five-phase mission profile. The simulation parameters are listed in Table 10.

Simulation Parameter Configuration

| Symbol | Description | Value |
|—|—|—|
| \(k_{p}, k_{v}\) | Formation position/velocity gains | 11.5 s⁻², 21.2 s⁻¹ |
| \(\eta\) | Repulsive gain | 20 m² |
| \(\rho_{0}\) | Influence distance | 5 m |
| \(a_{max}\) | Maximum acceleration | 15 m/s² |
| \(v_{max}\) | Maximum horizontal velocity | 20.0 m/s |
| \(d_{safe}\) | Obstacle buffer threshold | 10.0 m |
| \(d_{min}\) | Inter-UAV minimum distance | 2.0 m |
| \(\beta_{g}\) | Global damping coefficient | 0.5 s⁻¹ |

Experimental results demonstrate that:

(1) In the takeoff phase (S1), the Bézier curve smoothing effectively eliminates the sharp trajectory discontinuity at the vertical-to-horizontal transition point. UAVs using smoothed trajectories exhibit significantly reduced attitude oscillations and smooth acceleration profiles.

(2) In the cruise maneuver phase (S2), the velocity feedforward compensation enables UAVs on the outer turn trajectory to automatically compensate with higher linear velocity while inner UAVs execute tighter turns, achieving high-precision heading coordination.

(3) In the obstacle avoidance phase (S3), the improved APF algorithm successfully overcomes the local minima and GNRON problems. Compared to standard APF, which causes attitude oscillations of approximately ±20 degrees, the improved algorithm limits attitude variations to within 5 degrees while maintaining smooth velocity and acceleration profiles.

(4) In the formation reconfiguration phase (S4), the linear interpolation rule enables smooth transitions between diamond and line formations. The minimum inter-UAV distance remains above the safety threshold during the entire reconfiguration process, and all parameters converge rapidly after each transformation.

(5) In the descent phase (S5), the radius-decreasing spiral strategy achieves precise convergence to the landing center while maintaining smooth parameter variations throughout the descent.

5. Conclusions and Future Work

This thesis systematically investigated cooperative formation path planning for unmanned aerial vehicles based on cell-like membrane computing. The key contributions are summarized as follows:

First, this thesis established comprehensive mathematical models for quadrotor UAVs, including coordinate systems, translational and rotational dynamics, and systematically formalized the multiple constraints in path planning and formation control.

Second, a hybrid optimization algorithm (CSP-GWO) integrating spiking cell-like membrane systems with an improved grey wolf optimizer was proposed. The algorithm utilizes Tent chaotic mapping for population initialization, fuzzy logic-based dynamic weighting, and membrane system parallel architecture to enhance global exploration capability and convergence speed. Benchmark tests confirmed its superior optimization performance, and path planning experiments with 4, 8, and 16 UAVs demonstrated its excellent scalability, robustness, and computational efficiency.

Third, a cell-like membrane system-based formation path planning method was developed, integrating the virtual leader architecture with an improved artificial potential field approach. The stability of the proposed method was rigorously proven using Lyapunov stability theory and input-to-state stability analysis. Full-process simulations validated the effectiveness of the method across five mission phases, demonstrating smooth formation assembly, synchronized maneuvering, stable obstacle avoidance, flexible formation reconfiguration, and precision descent.

Future research directions include:
(1) Extending the path planning approach to dynamic environments with moving obstacles, since this thesis focuses on offline planning for known static environments.
(2) Conducting physical flight experiments on real UAV platforms to validate the practical performance of the proposed algorithms.
(3) Scaling up the experimental validation to 32, 64, or more UAVs to further assess the scalability of the proposed methods in high-density, large-scale cooperative scenarios.
(4) Exploring the integration of deep reinforcement learning with membrane computing for real-time adaptive path planning.