2.1 Oilfield environment modeling

The oilfield dynamic environment is modeled in Fig. 1. The UAVs need to track the moving target, avoiding the static and dynamic obstacles in the oilfield environment. In Fig. 1, the static obstacles contain oil tanks, water tanks, pumping wells and so on. The dynamic obstacles, such as flying birds, are assumed to move at a constant speed ${v_d}$ in a certain direction. The position at $k$ is defined as ${p_d} {\left( k \right)} = {[{x_d}{\left( k \right)},{y_d}{\left( k \right)},{z_d}{\left( k \right)}]^T}$ . The dynamic obstacles are modeled kinematically in Equation (1). ${\rm{\Delta }}T$ is the sampling time, and ${I_3}$ is the third-order identity matrix.

(1) \begin{equation}\left[ {\begin{array}{*{20}{l}}{{p_d} {\left( {k + 1} \right)}}\\[5pt]{{v_d} {\left( {k + 1} \right)}}\end{array}} \right] = \left[ {\begin{array}{*{20}{l}}{{I_3}}\quad {}{{\rm{\Delta }}T{I_3}}\\[5pt]0 \quad {}{{I_3}}\end{array}} \right]\left[ {\begin{array}{*{20}{l}}{{p_d} {\left( k \right)}}\\[5pt]{{v_d} {\left( k \right)}}\end{array}} \right]\end{equation}

2.2 UAV model

In this work, the ${N_u}$ quad-rotor UAVs track the target moving on the ground corporately. The wind disturbance is not considered here. The UAVs must avoid all the obstacles in the oilfield environment in order to track safely. Figure 2 shows the geometry of the target tracking. The UAV ${_i}$ position vector is ${p_i} = {[{x_i},{y_i},{z_i}]^T}$ ( $i = 1, \ldots ,{N_u}$ ), and the target position vector is ${p_t} = {[{x_t},{y_t},{z_t}]^T}$ . Here ${z_t}$ is zero due to the target moving on the ground.

Figure 1.

The oilfield environment model.

Figure 2.

The geometry of the target tracking.

The relative position vector between UAV ${_i}$ and the target is described in Equation (2).

(2) \begin{equation} {p_{i,t}} = {p_i} - {p_t} = {[ {{x_i} - {x_t},{y_i} - {y_t},{z_i} - {z_t}} ]^T}\end{equation}

The distance between the UAV ${_i}$ and the target is expressed in Equation (3), and its projection on the X-Y plane is written in Equation (4).

(3) \begin{align} {d_{i,t}} = \big| {{p_{i,t}}} \big| = \sqrt {{{( {{x_i} - {x_t}})}^2} + {{( {{y_i} - {y_t}})}^2} + {{( {{z_i} - {z_t}})}^2}} \\[-30pt] \nonumber \end{align}

(4) \begin{equation} {r_{i,t}} = \sqrt {{{( {{x_i} - {x_t}})}^2} + {{( {{y_i} - {y_t}})}^2}} \end{equation}

Assuming that the quad-rotor UAV is a mass point in the tracking process, the UAV dynamic model is shown in Equation (5) [Reference Wang, Li, Han and Yan28] :

(5) \begin{align} \left\{ {\begin{array}{*{20}{l}}{} {}{{x_i} = {V_i}\ {\rm{cos}}\ {\theta _i}\ {\rm{cos}}\ {\varphi _i}}\\{} {}{{y_i} = {V_i}\ {\rm{cos}}\ {\theta _i}\ {\rm{sin}}\ {\varphi _i}}\\{} {}{{z_i} = {V_i}\ {\rm{sin}}\ {\theta _i}}\\{} {}{{\varphi _i} = u_i^\varphi }\\{} {}{{\theta _i} = u_i^\theta}\\{} {}{{V_i} = u_i^V}\end{array}} \right.\end{align}

The heading angle ${\varphi _i}$ and pitching angle ${\theta _i}$ for the UAV ${_i}$ are written in Equation (6).

(6) \begin{align} \left\{ {\begin{array}{*{20}{l}}{} {}{{\rm{sin}}\ {\varphi _i} = \frac{{{x_{r,i}}}}{{{r_i}}}}\\[5pt]{} {}{{\rm{cos}}\ {\varphi _i} = \frac{{{y_{r,i}}}}{{{r_i}}}}\\[5pt]{} {}{{\rm{sin}}\ {\theta _i} = \frac{{{z_{r,i}}}}{{{d_i}}}}\\[5pt]{} {}{{\rm{cos}}\ {\theta _i} = \frac{{{r_i}}}{{{d_i}}}}\end{array}} \right.\end{align}

where ${V_i}$ is the speed for the UAV ${_i}$ . The minimal speed ${V_{{\rm{min}}}}$ is 0. $u_i^\varphi $ , $u_i^\theta $ , and $u_i^V$ are the control inputs for ${\varphi _i}$ , ${\theta _i}$ and ${V_i}$ . The state vector of the UAV ${_i}$ is denoted as ${{\bf{s}}_i} = {[{x_i},{y_i},{z_i},{\varphi _i},{\theta _i}]^T}$ , and the control input vector of UAV ${_i}$ is ${u_i} = {[u_i^\varphi ,u_i^\theta ,u_i^V]^T}$ . The UAV kinematic constraints are shown as follows:

(7) \begin{align} \left\{ {\begin{array}{*{20}{l}}{} {}{\left| {u_i^\varphi } \right| \le u_{{\rm{max}}}^\varphi }\\[5pt]{\left| {u_i^\theta } \right| \le u_{{\rm{max}}}^\theta }\\[5pt]{\left| {u_i^V} \right| \le u_{{\rm{max}}}^V}\\[5pt] {}{0 \le {V_i} \le {V_{{\rm{max}}}}}\\[5pt] {}{{H_{{\rm{min}}}} \le {z_i} \le {H_{{\rm{max}}}}}\end{array}} \right.\end{align}

where ${H_{{\rm{min}}}}$ and ${H_{{\rm{max}}}}$ represent the minimum and maximum flight altitude of the UAVs, respectively. $u_{{\rm{max}}}^\varphi $ , $u_{{\rm{max}}}^\theta $ and $u_{{\rm{max}}}^V$ are the maximum value for the control input $u_i^\varphi $ , $u_i^\theta $ and $u_i^V$ . ${V_{{\rm{min}}}}$ is the maximum speed.

2.3 Target model

In this work, the target moving on the ground is modeled as a mass point. The next position of the target can be obtained as follows [Reference Song, Cai, Zhou, Zhang, Chen, Li, Deng and Deng10]:

(8) \begin{align} {x_{t,k + 1}} = {v_{t,k}}{\rm{\Delta }}T \, {\rm{sin}} \,{\varphi _{t,k}} + {x_{t,k}} \\[-28pt] \nonumber \end{align}

(9) \begin{equation}{y_{t,k + 1}} = {v_{t,k}}{\rm{\Delta }}T\, {\rm{cos}}\, {\varphi _{t,k}} + {y_{t,k}} \end{equation}

where $\left( {{x_{t,k + 1}},{y_{t,k + 1}}} \right)$ and $\left( {{x_{t,k}},{y_{t,k}}} \right)$ are the position of the target at time $k + 1$ and $k$ , respectively. ${v_{t,k}}$ and ${\varphi _{t,k}}$ are the speed and heading angle of the target at time $k$ , respectively. ${\rm{\Delta }}T$ is the sampling time.

3.1 Description of DMPC

DMPC, as a typical control strategy, is used to predict future states and determine the sequence of the control inputs. The DMPC framework is shown in Fig. 3.

Figure 3.

The DMPC framework.

The DMPC process is described as follows. At time $k$ , ${u_i}{\left[ {k\,:\,k + N - 1 \left] = \right[{u_i}{\left( k \right)}, \ldots ,{u_i}{\left( {k + N - 1} \right)}} \right]}$ is the control-input sequences in future for the UAV ${_i}$ , and $N$ is the predicted steps. The control-input sequence is initialised independently and synchronously. ${{\bf{s}}_i}{\left[ {k\,:\,k + N - 1\left] = \right[{s_i}{\left( k \right)}, \ldots ,{s_i}{\left( {k + N - 1} \right)}} \right]}$ is the future state sequences, and the state information is shared between the UAVs through the communication network. ${J_i}{\left( {s{\left[ {k\,: + N - 1} \right]},{u_i}{\left[ {k\,:\,k + N - 1} \right]}} \right)}$ is the cost function and need to be minimised for obtaining the optimal control-input sequence (Equation (10)). The step $k$ in the optimised control-input sequence is only executed for the accuracy, and the subsequent steps ( ${u_i}{\left[ {k + 1\,:\,k + N - 1} \right]}$ ) are served as the initial control-inputs for the next prediction. The above procedure is repeated until the termination condition is satisfied.

(10) \begin{align} u_i^{\rm{*}}\left[ {k\,:\,k + N - 1} \right] = {\rm{arg \, min}} \, {J_i}\left( {s{\left[ {k\,:\,k + N - 1} \right]},{u_i}{\left[ {k\,:\,k + N - 1} \right]}} \right) \\[-30pt] \nonumber \end{align}

3.2 Cost function

To obtain the optimal control-input sequences, the cost function ${J_i}$ is established in Equation (11) and optimised.

(11) \begin{equation}{J_i}{\left( {s{\left[ {k\,:\,k + N - 1} \right]},{u_i}{\left[ {k\,:\,k + N - 1} \right]}} \right)} = {J_{{\rm{Obs}}}} + {J_{{\rm{Dis}}}} + {J_{{\rm{FOV}}}} + {J_{{\rm{LOS}}}} + {J_{{\rm{safe}}}} + {J_{{\rm{com}}}} + {J_u}\end{equation}

where the cost function includes the obstacle avoidance cost ${J_{{\rm{Obs}}}}$ , the distance cost ${J_{{\rm{Dis}}}}$ , the FOV cost ${J_{{\rm{FOV}}}}$ , the LOS occlusion cost ${J_{{\rm{LOS}}}}$ , the safety cost ${J_{{\rm{safe}}}}$ , the communication cost ${J_{{\rm{com}}}}$ and the control cost ${J_u}$ .

3.2.1 Obstacle avoidance cost ${J_{{\bf{Obs}}}}$

The UAVs must avoid obstacles in the oilfield environment in the tracking for safety. The obstacle avoidance cost is written in Equation (12). ${F_i}{\left(\, j \right)}$ represents the avoidance of the static obstacles, which is expressed in Equation (13) [Reference Tyagi, Kumar and Sujit29]; ${B_i}{\left(\, j \right)}$ represents the avoidance of the dynamic obstacles based on the chance constraint.

(12) \begin{equation}{J_{{\rm{Obs}}}} = {\lambda _{{\rm{Obs}}}}\mathop \sum \limits_{j = k}^{k + N - 1} \mathop \sum \limits_{i = 1}^{{N_u}} [ {{F_i}(\, j) + {B_i}(\, j)}]\end{equation}

where

\begin{equation*}{F_i} = \mathop \sum \limits_{o = 1}^M {P_1}\left( {{h_o} - {z_i}} \right)\end{equation*}

(13) \begin{equation}{P_1} = \left\{\begin{array}{ll}1 {}& (x_{i} - x_{o} )^{2} + (y_{i} - y_{o} )^{2} - r_{o}^{2} \le 0\\[5pt]0 {}& (x_{i} - x_{o} )^{2} + (y_{i} - y_{o} )^{2} - r_{o}^{2} \gt 0\end{array} \right.\end{equation}

where $\left( {{x_o},{y_o}} \right)$ and ${r_o}$ are the position and the radius of the projection centre of the obstacle on the X-Y plane, respectively. $N$ is the predicted steps, and ${N_u}$ is the number of the UAVs. ${\lambda _{{\rm{Obs}}}}$ is the weight of the obstacle avoidance cost. ${h_o}$ is the height of the static obstacles, and $M$ is the number of the static obstacles in the oilfield. ${P_1}$ is a sub-function for determining the position relationship of the UAVs with static obstacles. The UAVs should be located outside the static obstacle for avoiding the collision, i.e. the distance between the UAVs and the static obstacle is larger than the radius of the static obstacle. The UAVs collide with the static obstacles when the distance is less than the radius and the UAV flight height is lower than the height of the static obstacles.

For the dynamic obstacles, the chance constraints have been used in this work including the position chance constraint and the velocity chance constraint, which is reported in Ref. [Reference Wakabayashi, Suzuki and Suzuki27]. In Fig. 4, the position space ${C_{{\rm{pd}}}}$ and the velocity space ${C_{{\rm{vd}}}}$ are defined, respectively. In order to explain conveniently, the point $P$ is introduced as below. A vertical line is drawn from the centre $O$ of the obstacle and intersects with the extension line in the direction of the relative velocity ( ${v_i} - {v_d}$ ). The cross point is defined as $P$ . In order to avoid colliding with the dynamic obstacles, the probability of the UAV entering the position space is limited below a set threshold and suppresses the probability of the point $P$ entering the velocity space. The collision probability for the position chance constraints is described as follows:

(14) \begin{equation}{{\rm{g}}_{{pid}}} = a_{{pid}}^T( {{{\hat p}_i} - {{\hat p}_d}} ) - {b_{{id}}} - er{f^{ - 1}}( {1 - 2{\sigma _{{pid}}}} )\sqrt {2a_{{pid}}^T\left( {{{\rm{\Sigma }}_i} + {{\rm{\Sigma }}_d}} \right){a_{{pid}}}}\end{equation}

(15) \begin{equation}erf{\left( x \right)} = {2 \over {\sqrt \pi }}\mathop\int\nolimits_0^x {e^{ - {t^2}}}dt\end{equation}

(16) \begin{equation}{a_{{pid}}} = \frac{{( {{{\hat p}_i} - {{\hat p}_d}})}}{{{||\, {\hat p}_i} - {{\hat p}_d}||}}\end{equation}

(17) \begin{equation}{B_p} = \frac{{{w_p}}}{{{\rm{g}}_{{pid}}^2}}\end{equation}

where ${r_i}$ and ${r_d}$ are the radius of the UAV ${_i}$ and the dynamic obstacles ( ${b_{{id}}} = {r_i} + {r_d}$ ), respectively. ${\sigma _{{pid}}}$ is the collision probability threshold ( $0 \lt {\sigma _{pid}} \lt 0.5$ ). ${{\rm{\Sigma }}_i}$ and ${{\rm{\Sigma }}_d}$ are the error covariance matrices of the UAV ${_i}$ and the dynamic obstacles, respectively. $erf{\left( . \right)}$ is the standard error function. ${w_p}$ is the weight coefficient.

Figure 4.

Position space ${C_{{\rm{pd}}}}$ and speed space ${C_{{\rm{vd}}}}$ .

The collision probability for the velocity chance constraints is written as below:

(18) \begin{align}{{\rm{g}}_{{vid}}} = a_{{vid}}^T{\hat r_{{vid}}} - {b_{{id}}} - er{f^{ - 1}}( {1 - 2{\sigma _{{\rm{vi}}}}} )\sqrt {2a_{{vid}}^T{\left( {{{\rm{\Sigma }}_i} + {{\rm{\Sigma }}_d}} \right)}{a_{{vid}}}} \\[-25pt] \nonumber \end{align}

(19) \begin{align}{B_v} = \left\{{\begin{array}{l@{\quad}l}{{w_v} \cdot {{\rm{g}}_{{vid}}}} & {}{{\rm{if}}\ \ {{\rm{g}}_{{vid}}} \ge 0}\\0& {}{{\rm{else}}}\end{array}} \right. \\[-28pt] \nonumber \end{align}

(20) \begin{equation}{a_{{vid}}} = \frac{{{{\hat r}_{{vid}}}}}{{| {| {{{\hat r}_{{vid}}}} |} |}}\end{equation}

where ${w_v}$ is the weight coefficient. ${r_{{vid}}}$ is the vector from the intersection $P$ to $O$ . ${\sigma _{{vid}}}$ is the same as ${\sigma _{{pid}}}$ .

The chance-constrained cost function is written in Equation (21).

(21) \begin{equation}B = {B_p} + {B_v}\end{equation}

3.2.2 Distance cost ${J_{{\bf{Dis}}}}$

The distance cost is the distance between the UAV and the target on the horizontal plane. The smaller the distance, the better the tracking performance. The distance cost is expressed as:

(22) \begin{equation}{J_{{\rm{Dis}}}} = {\lambda _{\rm{d}}}\mathop \sum \limits_{j = k}^{k + N - 1} \mathop \sum \limits_{i = 1}^{{N_u}} \sqrt {{{({x_i}{\left(\, j \right)} - {x_t}{\left(\, j \right)})}^2} + {{({y_i}{\left(\, j \right)} - {y_t}{\left(\, j \right)})}^2}}\end{equation}

where $\left( {{x_i},{y_i}} \right)$ and $\left( {{x_t},{y_t}} \right)$ are the coordinates of the UAV and the target on the 2D plane, respectively. ${\lambda _d}$ is the weight of the distance cost.

3.2.3 FOV cost ${J_{{\bf{FOV}}}}$

The FOV cost is used to keep the target in the FOV of the UAV sensor. The target is within the coverage range of the UAV sensors when the detection distance and the field angle constraints are satisfied simultaneously. The FOV cost is defined as [Reference Wu, Wang, Li, Yao, Huang, Su and Yu19]:

(23) \begin{equation}{J_{{\rm{FOV}}}} = \left( {\begin{array}{*{20}{l}}{{\lambda _{m1}}\mathop \sum \nolimits_{j = k}^{k + N - 1} \mathop \sum \nolimits_{i = 1}^{{N_u}} \bigg(1 - {\rm{exp}} \! \bigg[ { - {c_1}\bigg( {\frac{{| {{p_i}(\, j ) - {p_t}(\, j )} |}}{{{R_{{\rm{max}}}}}}{\bigg)^2}} \bigg]} \bigg) ,}\\[7pt]{{{\left| {{p_i}\left(\, j \right) - {p_t}\left(\, j \right)} \right|}^2} \le R_{{\rm{max}}}^2 \cap \left( {{\rm{co}}{{\rm{s}}^{ - 1}}\left( {\frac{{{Z_i}}}{{\left| {{p_i}\left(\, j \right) - {p_t}\left(\, j \right)} \right|}}} \right) \le {\gamma _{{\rm{max}}}}} \right)}\\[7pt]{\mathop \sum \nolimits_{j = k}^{k + N - 1} \mathop \sum \nolimits_{i = 1}^{{N_u}} {\lambda _{m1}} + {\lambda _{m2}}\left| {{p_i}\left(\, j \right) - {p_t}\left(\, j \right)} \right|,}\\[7pt]{{{\left| {{p_i}\left(\, j \right) - {p_t}\left(\, j \right)} \right|}^2} \gt R_{{\rm{max}}}^2 \cup \left( {{\rm{co}}{{\rm{s}}^{ - 1}}\left( {\frac{{{Z_i}}}{{\left| {{p_i}\left(\, j \right) - {p_t}\left(\, j \right)} \right|}}} \right) \gt {\gamma _{{\rm{max}}}}} \right)}\end{array}} \right.\end{equation}

where ${\lambda _{m1}}$ and ${\lambda _{m2}}$ are the weight coefficients. ${c_1}$ is the parameter determining the shape of ${J_{{\rm{FOV}}}}$ . ${R_{{\rm{max}}}}$ and ${\gamma _{{\rm{max}}}}$ represent the maximum detection distance and the field angle of the UAV sensor, respectively. The target is within the coverage range of UAV sensors when the detection distance and field angle constraints are satisfied simultaneously.

3.2.4 LOS occlusion cost ${J_{{\bf{LOS}}}}$

The obstacles in the oilfield environment, including the static and dynamic obstacles, can obscure the UAV’s LOS to the target. Figure 5 shows the LOS obstruction conditions. To determine the relationship between the LOS’s projection on the X-Y plane and the obstacle, a sub-function ${P_1}$ is introduced according to Ref. [Reference Tyagi, Kumar and Sujit29].

(24) \begin{equation}{P_1} = \left\{ {\begin{array}{l@{\quad}l}1 {}& {\dfrac{{\left| {m{x_a} - {y_a} + {y_t} - m{x_t}} \right|}}{{\sqrt {{m^2} + {{( - 1)}^2}} }} - {r_a} \le 0}\\[14pt]0& {}{\dfrac{{\left| {m{x_a} - {y_a} + {y_t} - m{x_t}} \right|}}{{\sqrt {{m^2} + {{( - 1)}^2}} }} - {r_a} \gt 0}\end{array}} \right.\end{equation}

where $\left( {{x_a},{y_a}} \right)$ and ${r_a}$ are the position and the radius of the projection centre of the obstacle on the X-Y plane, respectively. $m$ is the slope of the LOS’s projection on the X-Y plane. The LOS is not obstructed at ${P_1} = 0$ , but LOS is obstructed by the obstacles at ${P_1} = 1$ and ${\theta _{{\rm{Obs}}}} \gt {\theta _{{\rm{LOS}}}}$ .

(25) \begin{equation} \begin{array}{*{20}{l}}{{J_{{\rm{LOS}}}}} {}{ = {\lambda _{{\rm{LOS}}}}\mathop \sum \limits_{j = k}^{k + N - 1} \mathop \sum \limits_{i = 1}^{Nu} {P_1}{\rm{max}}\left( {0,{\rm{tan}}{\left( {{\theta _{{\rm{Obs}}}}} \right)}- {\rm{tan}} {\left( {\theta _{{\rm{LOS}}}^i {\left(\, j \right)}} \right)} } \right)}\\{{\rm{tan}}{\left( {{\theta _{{\rm{Obs}}}}} \right)}} {}{ = \mathop \sum \limits_{a = 1}^A \frac{{{h_a}}}{{\sqrt {{{({x_t} - {x_a})}^2} + {{({y_t} - {y_a})}^2}} - {r_a}}}}\\{{\rm{tan}}{\left( {{\theta _{{\rm{LOS}}}}} \right)}} {}{ = \frac{{{Z_i}}}{{\sqrt {{{({x_t} - {x_i})}^2} + {{({y_t} - {y_i})}^2}} }}}\end{array}\end{equation}

where ${\lambda _{{\rm{LOS}}}}$ is the weight coefficient. $A$ is the number of obstacles in the oilfield, including the static and dynamic obstacles.

Figure 5.

LOS obstruction conditions.

Note: (a) LOS is not obstructed. (b) LOS is obstructed. The LOS’s projection on the X-Y plane intersects with the obstacle and ${\theta _{{\rm{Obs}}}} \gt {\theta _{{\rm{LOS}}}}$ . (c) LOS is not obstructed. The LOS’s projection on the X-Y plane intersects with the obstacle and ${\theta _{{\rm{Obs}}}} \lt {\theta _{{\rm{LOS}}}}$ .

3.2.5 Safety cost ${J_{{\bf{safe}}}}$ and communication cost ${J_{{\bf{com}}}}$

The safety distance cost ${J_{{\rm{safe}}}}$ and the communication distance cost ${J_{{\rm{com}}}}$ are introduced to keep the safety and communication for the UAVs, which is normalised between 0 and 1 by the hyperbolic tangent functions [Reference Yao, Wang and Ji18] :

(26) \begin{equation}{f_{{\rm{Sdi}}}}{\left( {d,{d_{{\rm{min}}}},{\rm{\Delta }}d} \right)} = \frac{1}{2}\left[ {1 - {\rm{tanh}}\left( {\frac{{d - {d_{{\rm{min}}}} - \frac{{{\rm{\Delta }}d}}{2}}}{{\frac{{{\rm{\Delta }}d}}{8}}}} \right)} \right]\end{equation}

(27) \begin{equation}{f_{{\rm{Cdi}}}}{\left( {d,{d_{{\rm{max}}}},{\rm{\Delta }}d} \right)} = \frac{1}{2}\left[ {1 + {\rm{tanh}}\left( {\frac{{d - {d_{{\rm{max}}}} + \frac{{{\rm{\Delta }}d}}{2}}}{{\frac{{{\rm{\Delta }}d}}{8}}}} \right)} \right]\end{equation}

where ${d_{{\rm{min}}}}$ and ${d_{{\rm{max}}}}$ represent the minimum safe distance and maximum communication distance, respectively. $d$ is the distance between the UAVs, and ${\rm{\Delta }}d$ is the buffer distance ( $0 \lt {\rm{\Delta }}d \lt\! \lt {d_{{\rm{min}}}}$ ).

(28) \begin{equation}{J_{{\rm{safe}}}} = {\lambda _{{\rm{safe}}}}\mathop \sum \limits_{j = k}^{k + N - 1} \mathop \sum \limits_{i \ne j}^{{N_u}} {f_{{\rm{Sdi}}}}{( {{d_{{\rm{ij}}}},{d_{{\rm{min}}}},{\rm{\Delta }}d} )}\end{equation}

(29) \begin{equation}{J_{{\rm{com}}}} = {\lambda _{{\rm{com}}}}\mathop \sum \limits_{j = k}^{k + N - 1} \mathop \sum \limits_{i \ne j}^{{N_u}} {f_{{\rm{Cdi}}}}{( {{d_{{\rm{ij}}}},{d_{{\rm{max}}}},{\rm{\Delta }}d} )}\end{equation}

(30) \begin{equation}{d_{{\rm{ij}}}} = \sqrt {{{({x_i} - {x_j})}^2} + {{({y_i} - {y_j})}^2} + {{({z_i} - {z_j})}^2}}\end{equation}

where ${d_{{\rm{ij}}}}$ is the distance between UAV ${_i}$ and UAV ${_j}$ ( $i \ne j$ ), ${\lambda _{{\rm{safe}}}}$ and ${\lambda _{{\rm{com}}}}$ are the weights of the safe and communication costs, respectively.

3.2.6 Control cost ${J_{\bf{u}}}$

UAVs have limited energy in the tracking target. Balancing control inputs helps reduce energy consumption. The balance of the control inputs is introduced as the control cost, which is given by:

(31) \begin{equation}{J_u} = {\lambda _u}\mathop \sum \limits_{j = k}^{k + N - 1} \mathop \sum \limits_{i = 1}^{{N_u}} \left| {\left| {{u_i}\left(\, j \right) - {u_i}\left( {j - 1} \right)} \right|} \right|\end{equation}

where ${\lambda _u}$ is a parameter.

4.1 SSA method

SSA is inspired by the swarming behaviour of the salps with the great performance including superior exploration-exploitation balance, practical simplicity and robust performance [Reference Mirjalili, Gandomi, Mirjalili, Saremi, Faris and Mirjalili30]. The population is divided into the leader and follower groups. The leaders are individuals at the front of the chain, leading the entire population to search for food sources. The other individuals as the followers connect one by one in the chain, following the leaders. The food source is identified, and the information is transferred from one individual to another. The mathematical process is described below.

First, the population is initialised randomly as:

(32) \begin{equation}{x_{W \times D}} = rand\left( {W,D} \right) \times \left( {ub - lb} \right) + lb\end{equation}

where $W$ is the population size and $D$ is the dimension.

In SSA, the position of the leader is updated as:

(33) \begin{equation}x_j^1 = \left\{ {\begin{array}{l@{\quad}l} {{F_j} + {c_1}(( u{b_j} - l{b_j})\ {c_2} + l{b_j}),} & {}{{c_3} \ge 0.5}\\ {{F_j} - {c_1}(( u{b_j} - l{b_j})\ {c_2} + l{b_j}),} {}& {{\rm{else}}}\end{array}} \right.\end{equation}

where $x_j^1$ is the location of the leader in the $j$ th dimension. ${F_j}$ is the location of the food source in the $j$ th dimension. $u{b_j}$ and $l{b_j}$ are the upper and lower bounds in the search space. ${c_2}$ and ${c_3}$ are the random numbers distributed uniformly in the interval $\left[ {0,1} \right]$ . ${c_1}$ decreases adaptively with the iteration and is denoted as Equation (34).

(34) \begin{equation}{c_1} = 2{e^{ - {{(\frac{{4t}}{{{T_{{\rm{max}}}}}})}^2}}}\end{equation}

where $t$ is the current iteration, and ${T_{{\rm{max}}}}$ is the maximum iteration.

In the salp chain, the position of the followers is updated as Equation (35).

(35) \begin{equation}x_j^i = \frac{1}{2}\left( {x_j^i + x_j^{i - 1}} \right)\end{equation}

where $x_j^{i - 1}$ and $x_j^i$ are the positions of the $\left( {i - 1} \right){\rm{th}}$ following individual and the $i{\rm{th}}$ following individual in the $j{\rm{th}}$ dimension, respectively ( $i \ge 2$ ,).

4.2 ISSA method

4.2.1 Sine-tent-cosine chaotic mapping

The initialisation stage is very important to improve the convergence speed and optimisation accuracy for SSA. The population of the traditional SSA is initialised randomly. In Ref. [Reference Haupt and Haupt31], the initial population is distributed uniformly in the search space and is helpful for improving the performance.

In this work, sine-tent-cosine chaotic mapping is used to initialise the salp population. In the sine-tent-cosine chaotic mapping, a continuous chaotic sequence is first generated using the sine function, and then the sequence is segmented with the tent function, and finally, the sequence is further transformed using the cosine function. The mathematical expression is defined as Equation (36).

(36) \begin{equation}x\left( {i + 1} \right) = \left( {\begin{array}{l@{\quad}l}{{\rm{cos}}\left( {\pi {\left( {r\,{\rm{sin}}{\left( {\pi x{\left( i \right)}} \right)} + 2{\left( {1 - r} \right)}x{\left( i \right)} - 0.5} \right)} } \right),} {}&{x{\left( i \right)} \lt 0.5}\\[5pt]{{\rm{cos}}\left( {\pi { \left( {r\,{\rm{sin}}{\left( {\pi x{\left( i \right)}} \right)} + 2{\left( {1 - r} \right)}\!\left( {1 - x{\left( i \right)}} \right) - 0.5} \right)}} \right),} {}& {{\rm{else}}}\end{array}} \right.\end{equation}

where $x{\left( i \right)}$ denotes the input value in the current iteration, $x{\left( {i + 1} \right)}$ is the output value in the next iteration and $r$ is a parameter in the range $\left[ {0,1} \right]$ .

4.2.2 Leader position updating based on differential variance operation

As can be seen from Equation (33), the leaders update their position according to the food source. However, the individuals do not share information in the updating stage. Inspired by the differential variance, the mutation strategy of DE/best/2 is used to improve the convergence speed and computational accuracy. The updating formula is written as Equation (37).

(37) \begin{equation}x_j^1 = {F_j} + {\rm{g}}{\left( {x_j^a - x_j^b} \right)} + {\rm{g}}{\left( {x_j^c - x_j^d} \right)}\end{equation}

where ${F_j}$ is the location of the food source in the $j{\rm{th}}$ dimension, and $g$ is the ratio factor. $x_j^a$ , $x_j^b$ , $x_j^c$ and $x_j^d$ are the different individuals selected from the population.

4.2.3 Follower position updating based on nonlinear adaptive inertia weights

To trade off the exploratory and exploitative capabilities in searching, the nonlinear adaptive inertia weights are introduced to update the position of the followers. The follower position is updated as Equations (38) and (39).

(38) \begin{equation}x_j^i = \frac{1}{2}\left( {x_j^i + \omega {\left( t \right)}{x_j}^{i - 1}} \right)\end{equation}

(39) \begin{equation}\omega {\left( t \right)} = \frac{{0.55}}{{\sqrt {1 + {{\big(\frac{{Kt}}{{{T_{{\rm{max}}}}}}\big)}^{10}}} }} + 0.4\end{equation}

where $K$ is the weight control factor, and $t$ is the current iteration. ${T_{{\rm{max}}}}$ is the maximum iteration. The weights are initialised as the maximum value of 0.95. As the iteration increases, the weights gradually converge to the minimum value of 0.4.

The pseudo-code of ISSA is shown as Algorithm 1. The flowchart of the trajectory optimisation process is shown in Fig. 6.

Figure 6.

The flowchart of the trajectory optimisation process.

Algorithm 1

The pseudo-code of ISSA

5.1 Simulation experiments and results

In this section, the ISSA is proposed as the solver for CCDMPC and compared with the seven typical algorithms. The specific parameter settings are shown in Table 1. The population size is $W = 30$ , and the maximum iteration ${T_{{\rm{max}}}} = 100$ .

Table 1.

Parameter settings for the eight algorithms

Figure 7.

Oilfield environment model.

The target tracking trajectories planned by the eight algorithms are shown in Figs 8–11. All the planned trajectories can avoid static obstacles. However, the tracking trajectory planned by the ISSA method fits the motion trajectory of the moving obstacle more closely, and thus the tracking performance is better, which can be seen in Fig. 11.

Figure 8.

Tracking trajectories planned by the WOA and GWO, respectively, in top-view.

Figure 9.

Tracking trajectories planned by the PSO and HHO, respectively, in top-view.

Figure 10.

Tracking trajectories planned by the CGWO and HJPSO, respectively, in top-view.

Figure 11.

Tracking trajectories planned by the AOA-HHO and ISSA, respectively, in top-view.

To confirm the dynamic obstacle avoidance, the real-time distance between the UAVs and the bird flock is exhibited in Figs 12–15. Assume that the radius of the dynamic obstacles is 200 m and the UAV is regarded as a particle. ${d_{{\rm{US}}}}$ is defined as the Euclidean distance between the UAVs and dynamic obstacles. To prevent collision between the UAVs and dynamic obstacles, it should satisfy ${d_{{\rm{US}}}} \ge $ 200 m. The distance between the UAVs and the bird flock falls below 200 m during 1128–1140s for GWO, and the distance is less than 200 m during 1116–1128 s for CGWO. This indicates that the UAVs collided with the bird flock in this interval. This fact also appears for WOA and PSO. However, the proposed ISSA method has planned the trajectories with the safe distance at the tracking time, as well as HHO, HJPSO and AOA-HHO.

Figure 12.

Real-time distance between the UAVs and the bird flock using CCDMPC+WOA and CCDMPC+GWO.

Figure 13.

Real-time distance between the UAVs and the bird flock using CCDMPC+PSO and CCDMPC+HHO.

Figure 14.

Real-time distance between the UAVs and the bird flock using CCDMPC+CGWO and CCDMPC+HJPSO.

Figure 15.

Real-time distance between the UAVs and the bird flock using CCDMPC+AOA-HHO and CCDMPC+ISSA.

Figures 16–19 show the total tracking time for moving targets. The undetected flag is set to 0 when one UAV can track the target; otherwise, the undetected flag is set to 1. The two UAVs cannot track the moving target at the first 400 and 200 s for WOA and CGWO, respectively, due to the high density of the obstacle distribution. This indicates that the WOA and CGWO methods are difficult in solving the tracking problem with the high density of the obstacles. This fact also appears for the AOA-HHO method and even becomes more serious in the whole tracking process for the PSO method. The GWO, HHO and HJPSO methods show better performance in tracking time. The proposed method, ISSA, can provide the best tracking trajectories for the UAVs in the whole tracking process. The tracking time reaches 1620 s and accounts for 98.2% of the tracking time.

Figure 16.

The total tracking time of the moving target using CCDMPC+WOA and CCDMPC+GWO.

Figure 17.

The total tracking time of the moving target using CCDMPC+PSO and CCDMPC+HHO.

Figure 18.

The total tracking time of the moving target using CCDMPC+CGWO and CCDMPC+HJPSO.

Table 2.

Statistical results of the total tracking time for all the algorithms

ASO: Avoiding static obstacles, ADO: Avoiding dynamic obstacles.

Figure 19.

The total tracking time of the moving target using CCDMPC+AOA-HHO and CCDMPC+ISSA.

Figure 20 is the graph of the distance between the UAVs and the schematic of the flight altitude during target tracking using the CCDMPC+ISSA method, respectively. The distance graph shows UAVs can avoid colliding and keep communicating with each other. In addition, the UAV moves up and down within a limited altitude range, which prevents the UAV from colliding with obstacles in the oilfield environment or obscuring the LOS between the UAV and the target.

Figure 20.

Distance graph between UAVs and flight altitude graph of UAVs.

Table 2 gives statistical results for all the algorithms. The algorithms run independently 20 times. Note that the algorithm with the addition of sine-tent-cosine mapping is referred to as the sine-tent-cosine salp swarm algorithm (TSSA), and the algorithm incorporating both the sine-tent-cosine mapping and differential variance operation is referred to as the sine-tent-cosine differential salp swarm algorithm (TDSSA). The indicators, such as maximum duration, average duration and standard deviation, show the superiority of the algorithms proposed in this paper. This also confirms that the strategies introduced have enhanced the capability of SSA in solving the UAV target tracking problem. In 20 runs, the proposed algorithm guarantees that tracking can avoid dynamic and static obstacles in the environment. The trajectories planned by the HHO and AOA-HHO algorithms are also capable of avoiding obstacles but become weaker than the proposed algorithm in other indicators.