2.2.1 Dynamics model

The forces experienced longitudinally by the stratospheric aerostat include buoyancy, gravity, aerodynamic resistance and the dynamic equation can be expressed as (positive upwards):

(1) \begin{align}\left\{ {\begin{array}{*{20}{l}}{\left( {{m_{tot}} + {m_{add}}} \right)\frac{{d{v_h}}}{{dt}} + {{\dot m}_{tot}}{v_h} = B - {m_{tot}}g - {D_h}}\\{\frac{{dh}}{{dt}} = {v_h}}\end{array}} \right.\end{align}

where ${m_{tot}}$ is the total mass of the stratospheric aerostat, ${m_{add}}$ is additional quality, $B$ is total buoyancy, ${D_h}$ is aerodynamic resistance in the altitude direction, $h$ is altitude, ${v_h}$ is vertical velocity.

The buoyancy force and longitudinal aerodynamic resistance of the stratospheric aerostat can be expressed as Equation (2).

(2) \begin{align}\left\{ {\begin{array}{*{20}{l}}{B = {\rho _{air}}{V_{total}}g}\\{{D_h} = \frac{1}{2}{\rho _{air}}{v_h}\left| {{v_h}} \right|{C_d}S}\end{array}} \right.\end{align}

where ${V_{total}}$ is the total volume of the aerostat, ${\rho _{air}}$ is atmospheric density. ${C_D}$ is drag coefficient, $S$ is aerostat reference area.

The altitude adjustment of the stratospheric aerostat is mainly achieved by the adjustment of valves and compressor. When the stratospheric aerostat is in the process of rising, the air valve opens, the air in air bag is discharged, and the quality of the air changes as Equation (3):

(3) \begin{align}\frac{{d{m_{air}}}}{{dt}} = - n{A_{valve}}C\sqrt {2{\rm{\Delta }}{P_{valve}}{\rho _{air}}} \end{align}

where $n$ is the number of valves, ${A_{valve}}$ is the cross-sectional area of a single valve, $C$ is the exhaust coefficient, which is obtained by considering the difference between theoretical flow and actual flow, ${\rm{\Delta }}{P_{valve}}$ is pressure difference between inside and outside the valve.

When the aerostat is in the process of descending, the air outside is compressed by the compressor into the air bag, and the mass flow rate can be expressed as:

(4) \begin{align}\frac{{d{m_{air}}}}{{dt}} = {\rho _{air}}{V_{in}}\left( h \right)\end{align}

where ${V_{in}}$ is the air flow into the air bag per unit time.

The horizontal kinetic equation can be expressed as:

(5) \begin{align}\left\{ {\begin{array}{*{20}{l}}{\left( {{m_{tot}} + {m_{add}}} \right)\frac{{d{v_x}}}{{dt}} = \frac{1}{2}{\rho _{air}}\left( {{v_x} + u} \right)\left| {{v_x} + u} \right|{C_d}S}\\{\frac{{dx}}{{dt}} = {v_x}}\\{\left( {{m_{tot}} + {m_{add}}} \right)\frac{{d{v_y}}}{{dt}} = \frac{1}{2}{\rho _{air}}\left( {{v_y} + v} \right)\left| {{v_y} + v} \right|{C_d}S}\\{\frac{{dy}}{{dt}} = {v_y}}\end{array}} \right.\end{align}

where $u$ and $v$ are the wind speed components in the east-west direction and the north-south direction, ${v_x}$ and ${v_y}$ are aerostat speeds in the horizontal direction, $x$ and y are displacements in the horizontal direction.

2.2.2 Description of the optimal coverage problem in the restricted aera for stratospheric aerostats

Considering that $n$ stratospheric aerostats are randomly distributed in the restricted area $Q \subset {R^2}$ , the position information of the aerostat is $P = \left\{ {{p_1},{p_2}, \ldots, {p_n}} \right\}$ , and the sensitivity function corresponding to any point $q$ in the restricted area is $\phi\! \left( q \right)$ , indicating the importance of the point information to the aerostat [Reference Zuo, Yan and Yan23]. To quantitatively describe the coverage effect of stratospheric aerostat on a restricteded area, the cost function is defined:

(6) \begin{align}J\left( {P,W,\phi } \right) = \mathop \sum \limits_{i = 1}^n \mathop \int \nolimits_{{W_i}} f\!\left( {{p_i},q} \right)\phi\! \left( q \right)dq\end{align}

where, $W = \left\{ {{W_1},{W_2}, \ldots, {W_n}} \right\}$ represents the distribution area of each stratospheric aerostat, and $f\!\left( {{p_i},q} \right)$ represents the measurement cost (distance, time, etc.) of the stratospheric aerostats $i$ from its position ${p_i}$ to any point $q$ within the distribution area.

As can be seen from the above formula, the formation coverage effect of multi-stratospheric aerostats is related to the position information $P$ of the aerostat, the area ${W_i}$ assigned to each aerostat, and the sensitivity function $\phi\! \left( q \right)$ [Reference Ma, He, Wang and Yao24]. In this section, the distribution area of each aerostat is adjusted to minimise the total cost function by the stratospheric aerostat dynamics model and altitude adjustment using the wind field model.

Voronoi partition is an important content in computational geometry. According to the degree of adjacence of points ${p_i}$ in a region, the region $Q$ is divided into several sub-regions, which are called sub-regions ${V_i}$ of the region. The Euclidean distance from the point ${p_i}$ to the sub-region is smaller than the distance from the other points ${p_j}$ to the region.

The uniform coverage of stratospheric aerostats in a restricted area is mainly considered in the study, and the spatial configuration of multiple aerostats can be regarded as a two-dimensional area. The cost function $f\!\left( {{p_i},q} \right)$ can be defined as the Euclidean distance from the point $q$ to the stratospheric aerostat ${p_i}$ , and the above formula can be expressed as:

(7) \begin{align}{V_i} = \{ q \in Q|\|{p_i} - q \|\le \|{p_j} - q\|,\forall i \ne j,\forall i,j \in n\} \end{align}

The cost function can be expressed as:

(8) \begin{align}J\!\left( P \right) = \mathop \sum \limits_{i = 1}^n \mathop \int \nolimits_{{V_i}} \|{p_i} - q\|\phi\! \left( q \right)dq\end{align}

Taking the derivative of the cost function along the trajectory of each stratospheric aerostat yields:

(9) \begin{align}\begin{array}{*{20}{l}}{\frac{{\partial J\left( P \right)}}{{\partial {p_i}}} = \mathop \smallint \nolimits_{{\!V_i}} \frac{\partial }{{\partial {p_i}}}{{({p_i} - q)}^2}\phi\! \left( q \right)dq + \mathop \smallint \nolimits_{{\!V_i}} \frac{{\partial {q_{\partial {V_i}}}\left( P \right)}}{{\partial {p_i}}}{n_{\partial {V_i}}}{{({p_i} - q)}^2}dq}\\[5pt]{ + \mathop \sum \limits_{j \in {N_i}} \mathop \smallint \nolimits_{{\!l_{ji}}} \frac{{\partial {q_{{l_{ij}}}}\left( {{p_i},{p_j}} \right)}}{{\partial {p_i}}}{n_{{l_{ji}}}}{{({p_j} - q)}^2}dq}\end{array}\end{align}

where, $\partial {V_i}$ and ${q_{\partial {V_i}}}$ are respectively the boundary and the point on the boundary of the Voronoi subregion ${V_i}$ , are the normal direction of the boundary of the subregion, ${N_i}$ are the set of neighbours of the Voronoi subregion ${V_i}$ , ${l_{ij}}$ are the boundary between the subregions ${V_i}$ and ${V_j}$ , ${q_{{l_{ij}}}}\left( {{p_i},{p_j}} \right)$ and ${n_{ji}}$ are the any point on the boundary and normal direction of ${l_{ij}}$ .

The Voronoi region boundary of a restricted region can be divided into two parts: the internal boundary of a Voronoi unit and the boundary of a restricted region. By the nature of the Voronoi partition, we know that $\partial {V_i}$ consists of components ${l_{ij}}$ , where the directions between ${l_{ij}}$ and ${l_{ji}}$ are opposite. Therefore, the above formula can be expressed as:

(10) \begin{align}\begin{array}{*{20}{l}}{\frac{{\partial J\left( P \right)}}{{\partial {p_i}}} = \mathop \smallint \nolimits_{{\!V_i}} \frac{\partial }{{\partial {p_i}}}{{({p_i} - q)}^2}\phi\! \left( q \right)dq + \mathop \smallint \nolimits_{\!Q} \frac{{\partial {q_{\partial Q}}\left( P \right)}}{{\partial {p_i}}}{n_{\partial {V_i}}}{{({p_i} - q)}^2}dq}\\[5pt]{ + \mathop \sum \limits_{j \in {N_i}} \mathop \smallint \nolimits_{{\!l_{ji}}} \frac{{\partial {q_{{l_{ij}}}}\left( {{p_i},{p_j}} \right)}}{{\partial {p_i}}}{n_{{l_{ji}}}}[{{({p_i} - q)}^2} - \left( {{p_j} - q{)^2}} \right]{{({p_j} - q)}^2}dq}\end{array}\end{align}

Since the boundary of the restricted region $\partial Q$ is fixed and unchanging, and $\frac{{\partial {q_{\partial Q}}\left( P \right)}}{{\partial {p_i}}} = 0$ , Equation (10) can be simplified to:

(11) \begin{align}\frac{{\partial J\!\left( P \right)}}{{\partial {p_i}}} = \mathop \int \nolimits_{{V_i}} \frac{\partial }{{\partial {p_i}}}{({p_i} - q)^2}\phi\! \left( q \right)dq\end{align}

According to the convex optimisation theory, if $\frac{{\partial J\left( P \right)}}{{\partial {p_i}}} = 0$ , we can obtain:

(12) \begin{align}p_i^{\rm{*}} = \mathop \int \nolimits_{{V_i}} q\phi\! \left( q \right)\!dq/\mathop \int \nolimits_{{V_i}} \phi\! \left( q \right)dq = {L_{{V_i}}}/{M_{{V_i}}} = {C_{{V_i}}}\end{align}

where, ${M_{{V_i}}} = \mathop \smallint \nolimits_{{V_i}} \phi\! \left( q \right)dq$ is the mass of the Voronoi subregion ${V_i}$ , ${L_{{V_i}}} = \mathop \smallint \nolimits_{{V_i}} q\phi\! \left( q \right)dq$ is the static moment, and ${C_{{V_i}}}$ is the centre of mass of the Voronoi subregion [Reference Du, Faber and Gunzburger25, Reference Lloyd26].

In the process of aerostat formation coverage, the Lloyd algorithm is employed to dynamically update the Voronoi segmentation, with the cycle being repeated [Reference Lloyd26]. By adopting this method, each aerostat is equipped with a controller, and information is shared among nearby aerostats, which aligns with the principles of distributed control. When an aerostat enters or exits the restricted area, its influence on other aerostats is relatively minimal, thereby demonstrating the robustness of the control method.

Assuming that the coverage radius of each aerostat is ${r_c} = 100{\rm{km}}$ , and the communication radius ${R_c}$ is $\sqrt 3 {r_c}$ . The connectivity between aerostat is defined as follows [Reference Yang and Xia27, Reference Yao, Liao, Liu and Yang28]:

(13) \begin{align}{p_c}\!\left( {{b_i},{b_j}} \right) = \left\{ {\begin{array}{*{20}{l}}{1{\rm{\;\;\;\;\;\;\;\;}}d\!\left( {{b_i},{b_j}} \right) \le {R_c}}\\{0{\rm{\;\;\;\;\;\;\;\;}}d\!\left( {{b_i},,{b_j}} \right) \gt {R_c}}\end{array}} \right.\end{align}

where, ${p_c} = 1$ means that the network single-hop connectivity is realised between the aerostats $i$ and $j$ , and ${p_c} = 0$ means that the aerostats is not connected. $d\!\left( {{b_i},{b_j}} \right)$ indicates the distance between two aerostats. For aerostats that cannot be directly connected, multi-hop connectivity can be performed through other aerostats connected with single hop, and the aerostat not connected with other aerostats is called isolated system [Reference Yang and Xia27]. To describe coverage performance, define multistratospheric aerostats coverage rate:

(14) \begin{align}Rati{o_{\mathop {{\rm{cov}}}\limits_{} }} = \frac{{\left( {{S_{cover,i}} \cup {S_{cover,j}} \cup \cdots } \right) \cap {S_{rest}}}}{{{S_{rest}}}} \times 100{\rm{\% }}\end{align}

where ${S_{cover,i}}$ represents the area covered by each aerostat. ${S_{rest}}$ represents the restricted area to be covered.

Define the connectivity rate:

(15) \begin{align}Rati{o_{con}} = \frac{{{N_{p\_total}} - {N_{islo}}}}{{{N_{p\_total}}}} \times 100{\rm{\% }}\end{align}

According to the established aerostat dynamic model, the optimal coverage of the aerostats in the restricted area (east-west direction $500{\rm{km}}$ , north-south direction $500{\rm{km}}$ ) is calculated under the condition that the aerostat only relies on altitude adjustment and uses the wind.

In the coverage range, 10 particle positions are randomly generated as the initial position of the aerostats. The sensitivity function $\phi\! \left( q \right) = 1$ is assumed, and the wind field condition is assumed to be:

(16) \begin{align}\left\{ {\begin{array}{*{20}{l}}{{v_w} = 10{\rm{m}}/{\rm{s}}}\\{\varphi = \left( {h - {h_{{\rm{min}}}}} \right)/{h_0} \cdot 2\pi }\\{{h_{{\rm{min}}}} \le h \le {h_{{\rm{max}}}}}\end{array}} \right.\end{align}

where, ${v_w}$ is the size of the wind speed, $\varphi $ is the angle between the wind direction and the north direction, and ${h_0} = 2{\rm{km}}$ , ${h_{{\rm{min}}}} = 18{\rm{km}}$ , ${h_{{\rm{max}}}} = 22{\rm{km}}$ . The simulation results under the ideal wind field condition are shown in Fig. 4. In the region, multiple aerostats can achieve optimal coverage under this algorithm.

Figure 4.

Simulation results in ideal wind.

2.2.3 Motion control of aerostat based on virtual force method

The assumed wind field environment is ideal, and the wind direction can be directed in any direction $\left[ {0,2\pi } \right]$ depending on the flying altitude of the aerostat. The actual wind direction is quite different from the ideal wind. The wind layer in the appropriate direction should be found according to the actual wind field to control the movement of the aerostat toward the centre of mass of the Voronoi.

Figure 5 shows the aerostat motion control strategy based on the virtual force method. A virtual gravitational field centred on the dynamic target point is constructed, and the aerostat is affected by the virtual force ${\vec F_b}$ under the action of the virtual force field. In order to make the aerostat move towards the target point, it is necessary to control the aerostat to reach the desired wind layer by inflating and deflating the air through the altitude adjustment system, and the wind direction of the desired wind layer is consistent with the direction of the virtual force. The altitude ${h_{expect}}$ corresponding to the desired wind layer is defined as:

(17) \begin{align}\begin{array}{*{20}{l}}{{\rm{min}}\,\delta = \delta \!\left( {t,h} \right) = {\rm{arccos}} \left({\frac{{{{\vec F}_b}\left( t \right) \cdot {{\vec k}_w}\left( {t,h} \right)}}{{\left| {{{\vec F}_b}\left( t \right)} \right| \cdot \left| {{{\vec k}_w}\left( {t,h} \right)} \right|}}} \right)}\\{{\rm{subject\;\;to}}\,:\,{h_{{\rm{max}}}} \le h \le {h_{{\rm{max}}}}}\end{array}\end{align}

where, ${\vec F_b}\!\left( t \right)$ is the virtual gravitational force between the aerostat at the time $t$ and the target point, and ${\vec k_w}\!\left( {t,h} \right)$ is the wind direction under the longitude, latitude and altitude $h$ of the aerostat at the time $t$ . This formula indicates that when the wind direction of the expected wind layer is inconsistent with the direction of the virtual force, the wind layer with the wind direction closest to ${\vec F_b}\!\left( t \right)$ is selected. When there are multiple wind layers with the same wind direction, the wind layer with higher wind speed is selected to control the aerostat to move quickly to the target point. In order to move the aerostat to the centre of mass of the Voronoi subregion, the centre of mass of the Voronoi subregion is set as the target point.

Figure 5.

Schematic diagram of aerostat motion control strategy based on virtual force method.

2.2.4 Motion control of aerostat based on ${A^{{*}}}$ algorithm

While the aerostat motion control strategy based on the virtual force method effectively assists the aerostat in identifying the optimal wind direction from its current position to the target location, as well as in navigating back to the restricted area, it fails to account for wind speed at altitude and the duration required for return flight. This section presents a path planning method designed to optimise the return trajectory of an aerostat that has flied outside the restricted area, ensuring it can return in the shortest possible time.

The ${A^{\rm{*}}}$ algorithm is a heuristic graph search algorithm, which searches the minimum cost path in the space in the best priority way, thereby improving the search efficiency. Its core formula is as follows [Reference Yao, Liao, Liu and Yang28–Reference Miyombo, Liu, Mulenga, Siamulonga, Kabanda, Shaba, Xi and Ayodeji30]:

(18) \begin{align}f\!\left( n \right) = g\!\left( n \right) + h\!\left( n \right)\end{align}

where, $f\!\left( n \right)$ is the total cost function, $g\!\left( n \right)$ is the actual consumption from the start node to the current node. $h\!\left( n \right)$ is the estimated value from the current node to the target node, which is called the heuristic function and is the key to finding the optimal path.

${A^{\rm{*}}}$ algorithm is suitable for solving problems where the initial point and target point are known. In this section, the algorithm is designed for the task scenario of aerospace-restricted area coverage, so that it can dynamically search the target point near the restricted area and fly back to the restricted area in the shortest time. The boundary covered by the aerospace-restricted area is a rectangular area. The heuristic function for defining the algorithm is as follow:

(19) \begin{align}h\!\left( n \right) = h\!\left( {{h_i}} \right) = \frac{{{{\left| {\vec f} \right|}^2}}}{{\vec f \cdot \vec v\left( {{h_i}} \right)}}\\[-18pt]\nonumber\end{align}

(20) \begin{align}\vec f = {O_Q} - {p_i}\\[6pt]\nonumber\end{align}

where, ${p_i}$ is the horizontal position coordinate of the aerostat outside the area, ${O_Q}$ is the boundary of the area to be covered, and $\vec f$ is the vertical vector pointing to the boundary of the aerostat location. $\vec v\left( {{h_i}} \right)$ is the wind field at the latitude and longitude of the aerostat that varies with altitude.

In the case of altitude adjustment, the valve and compressor of the aerostat air bag and are controlled to make the aerostat enter the wind layer of different altitudes. It is necessary to limit the search altitude of the aerostat, so that the aerostat can reach the corresponding target air layer after the operation is performed:

(21) \begin{align}{\rm{\Delta }}h = \int {v_h}dt\\[-18pt]\nonumber\end{align}

(22) \begin{align}\left| {{h_j} - {h_i}} \right| \le {\rm{\Delta }}h\\[6pt]\nonumber\end{align}

where, ${\rm{\Delta }}h$ is the altitude difference before and after the altitude adjustment of the aerostat, ${h_1}$ and ${h_2}$ are respectively the altitudes of adjacent wind layers to be explored.