Nomenclature
- Symbol
-
Definition
- AAR
-
Autonomous aerial refuelling
- LQR, MPC, NMPC
-
Linear quadratic regulator, model predictive control and nonlinear model predictive control, respectively
- UAV
-
Unmanned aerial vehicle
- J
-
NMPC objective function
- m
-
Mass of the receiver (kg)
-
$N_p,N_u$
-
Prediction- and control-horizon lengths, respectively
-
$Q,R,S$
(NMPC)
-
State-error, control-input-increment and terminal-state-error weighting matrices, respectively
- S (aerodynamics)
-
Reference wing area (m2)
- T
-
Sampling period (s)
-
$\vec u,\Delta\vec u$
-
Control-input vector and its increment, respectively
-
$v_x,v_y,v_z$
-
Body-axis velocity components (m/s)
-
$V_a$
-
Airspeed of the receiver (m/s)
-
$x,y,z$
-
Receiver position coordinates in the reference coordinate system (m)
-
$\vec X,\vec X_d,\widehat{\vec X}$
-
Actual, reference and predicted state vectors, respectively
-
$\alpha,\beta$
-
Angle of attack and sideslip angle, respectively (rad)
-
$\delta_a,\delta_e,\delta_r$
-
Aileron, elevator and rudder deflection angles, respectively (rad)
-
$\rho$
-
Air density (kg/m3)
-
$\phi,\theta,\psi$
-
Roll, pitch and yaw Euler angles, respectively (rad)
-
$\omega_x,\omega_y,\omega_z$
-
Body-axis roll, pitch and yaw angular rates, respectively (rad/s)
1.0 Introduction
Autonomous aerial refueling (AAR) technology can expand the combat radius and extend the air holding time of unmanned aerial vehicles (UAVs), and is a multiplier of air power [Reference Rasol, Yuelei and Qing1, Reference Lee and Kim2]. In the process of AAR, it is necessary to study high-precision flight control methods to ensure that the receiver approaches the tanker autonomously and smoothly.
In recent years, the flight control problem of the receiver actively approaching the tanker has attracted the research interest of many scholars at home and abroad, and has achieved fruitful research results [Reference Weijia3]. Burns and Lugo et al. generated the flight path using the Dubins curve method, and combined the nonlinear dynamic inverse control method and vector field theory to achieve the active approach of the receiver to the tanker [Reference Burns, Blue and Zollars4–Reference Owen, Nichols and Colton6]. Some scholars draw on the guidance idea of seeking to guide the receiver to fly autonomously to a safe position behind the tanker by designing guidance laws. For example, Ochi et al. improved the proportional guidance law and established a simple and feasible guidance method by controlling the line-of-sight angle [Reference Ochi and Kominami7]. Smith et al. improved the proportional guidance law based on adaptive terminal guidance, which not only solved the problem of autonomous approach of the receiver and tanker, but also analysed the influence of external disturbances [Reference Smith8]. Yoon and Kung et al. adopted the pure tracking guidance law and its improved method to achieve the autonomous approach of the receiver and tanker, improving the sensitivity and guidance during the approach [Reference Yoon, Kim and Kim9, Reference Kung, Chiang and Chen10]. In addition to drawing on the idea of navigation guidance law, some scholars draw on the idea of guidance control rate to solve the problem of the receiverr’s autonomous approach to the tanker. For example, Jun et al. improved the proportional guidance law based on the backstep control method to achieve terminal heading control and designed the control law of the altitude channel [Reference Jun and Xinming11]. Dogan et al. designed a linear quadratic regulator (LQR) using the multi-input and multi-output state feedback method to achieve the autonomous approach of the receiver and tanker [Reference Dogan and Venkataramanan12]. Schumacher et al. combined the adaptive control method, neural network method and dynamic inverse theory to design the control law of the flight motion of the receiver, achieving autonomous approach to the tanker [Reference Schumacher and Kumar13]. Gong Quanquan and Yuan Suozhong et al. designed the guidance law and control law for the flight motion of the receiver based on the relationship between the position, speed and target line of sight between the tanker and receiver, achieving the autonomous approach of the two [Reference Quanquan, Suozhong and Jin14, Reference Suozhong, Xinhua and Fengying15].
Although the simulation results of the above-mentioned methods demonstrate the effectiveness of the designed guiding control method, it still can not well solve the problems of nonlinearity, strong coupling and multiple constrained variables existing in the dynamical model of the receiver. The model predictive control (MPC) method is a model-based finite-time-domain open-loop optimisation method. It uses the dynamic model of the controlled object to predict the future output of the system and combines rolling optimisation and feedback correction strategies. The performance indicators of the model are minimised within the constraint range to determine the optimal control quantity [Reference Xiaoqi16, Reference Hongbiao, Yu and Xiaoying17]. However, during the rolling optimisation, the model needs to be linearised, and a precise nonlinear model of the system can not be accurately established, which is prone to cause the model mismatch problem [Reference Sixu, Biao and Manjiang18, Reference Hongwei, Chenyu and Lei19]. In view of the deficiencies of the MPC method, nonlinear model predictive control (NMPC) can make up for it well. It not only has good control performance for nonlinear models, but also can explicitly handle system constraints and solve non-convex programming problems [Reference Chen, Liu and Shi20]. Therefore, aiming at the characteristics of nonlinearity and multivariable constraints such as state and input variables in the motion model of the receiver during the receiver autonomous approach to the tanker, an autonomous approach method based on NMPC for tanker/receiver is proposed. The proposed method can quickly achieve the autonomous approach of the receiver and the tanker, and has good control performance.
The organisational structure of this paper is as follows. Section 2 conducts theoretical analysis and mathematical derivation of the in-air motion process of the receiver and establishes a nonlinear dynamic model of the receiver. Section 3 designs a NMPC method. The adaptive gradient descent method is adopted to solve the optimal control quantity online, and the convergence and closed-loop stability of the method are analysed by using the Lyapunov stability theory. Section 4 verified the validity and accuracy of the designed method through numerical simulation; Section 5 summarises and generalises the research work done.
2.0 Nonlinear dynamic model of the unmanned receiver
Based on the background of AAR of UAVs, this paper mainly studies the control methods for the receiver to autonomously approach the tanker. Therefore, establishing an accurate dynamic model of the receiver is the primary task. The receiver system is complex, nonlinear and strongly coupled. Its motion can be divided into two parts: the motion of the centre of mass and the rotation around the centre of mass, namely three linear motions and three angular motions [Reference Weijia3, Reference Bolian21]. In this section, a set of multi-dimensional state variables
$\vec{X}=[x,y,z,v_{x},v_{y},v_{z},\varphi ,\theta ,\psi ,\omega _{x},\omega _{y},\omega _{z}]^{T}$
will be selected. This set of state variables contains the position information of the receiver in the geocentric coordinate system, as well as the velocity, attitude angle and angular velocity information in the body coordinate system, in order to establish the dynamic model of the receiver. Before establishing the dynamic model of the receiver, it is necessary to analyse and model the forces and moments acting on it during its motion process.
2.1 Analysis of force and torque
During the flight of the receiver in the air, The receiver is subjected to different external forces. The common external forces mainly include the gravitational force
$\vec{G}$
, the thrust
$\vec{f}$
of the receiver and the total aerodynamic force
$\vec{A}$
. It is precisely because of the combined effect of these external forces that the receiver can achieve different manoeuvres in the air. The different external forces acting on the receiver have different effects and properties. Next, we will analyse and model these different external forces.
2.1.1 Gravity
The receiver is always affected by the force of gravity during flight, and the direction is always vertically downward. During the AAR process, the flight altitude of the receiver is usually between 4 and 10 km [Reference Bolian21]. The influence of altitude on the gravitational force acting on the receiver can be ignored and assumed as a constant vector. In the geocentric coordinate system, the gravitational force
$\vec{G}$
can be expressed as [Reference Xiaoqi, Changsheng and Wuxing22]:
\begin{equation} \vec{G}=m\left[\begin{array}{l@{\quad}l@{\quad}l} 0 & {} 0 & {} g \end{array}\right]^{T} \end{equation}
where g is the gravitational acceleration, and g = 9.8 m/s 2 is a constant.
According to the transformation relationship between the geocentric coordinate system and the body coordinate system of the receiver [Reference Xiaoqi, Changsheng and Wuxing22], the expression of gravity in the body coordinate system can be obtained:
\begin{align} {\,\vec{\!F}}_{\!G}^{b} & ={C}_{g}^{b}\vec{G}=m\left[\begin{array}{c@{\quad}c@{\quad}c} \cos \theta \cos \psi & {} \cos \theta \sin \psi {}& -\sin \theta \\[2pt] \sin \theta \cos \psi \sin \varphi -\sin \psi \cos \varphi & {} \sin \theta \sin \psi \sin \varphi +\cos \psi \cos \varphi & {} \cos \theta \sin \varphi \\[2pt] \sin \theta \cos \psi \cos \varphi +\sin \psi \sin \varphi {}& \sin \theta \sin \psi \cos \varphi -\cos \psi \sin \varphi {}& \cos \theta \cos \varphi \end{array}\right] \left[\begin{array}{c} 0\\[2pt] 0\\[2pt] g \end{array}\right]\nonumber\\[9pt] {} & =m\left[\begin{array}{c} -g\sin \theta \\[2pt] {g}\cos \theta \sin \varphi \\[2pt] {g}\cos \theta \cos \varphi \end{array}\right] \end{align}
where
${C}_{g}^{b}$
is the rotation matrix from the geocentric coordinate system to the body coordinate system of the receiver; (
$\varphi, \theta, \psi$
) are the roll angle, pitch angle and yaw angle of the receiver, respectively.
The gravitational force acting on the receiver is evenly distributed throughout the body of the receiver. The resultant gravitational force acts on the centre of the body of the receiver and does not generate a gravitational moment. Therefore, the moment of gravity in the body coordinate system of the receiver can be expressed as:
\begin{equation} {\vec{\!M}}_{\!G} {}^{\kern-2.7pt b}=\left[\begin{array}{l@{\quad}l@{\quad}l} 0 {}& 0 & {} 0 \end{array}\right]^{T} \end{equation}
2.1.2 Thrust
The thrust received by the receiver is used to overcome the air resistance and control the flight speed and altitude. Imagine the engine of the receiver is a vector engine. In the body coordinate system of the receiver, the thrust can be expressed as:
\begin{equation} {\vec{\!F}}_{\!f} {}^{\kern-1.8pt b}=\left[f_{x},\ f_{y},\ f_{z}\right]^{T} \end{equation}
When the engine is installed at a certain position away from the centre of the receiver, it will cause the thrust to generate a corresponding thrust torque relative to the origin of the body coordinate system of the receiver. In the body coordinate system of the receiver, suppose the distances from each component of the thrust to the centre of mass of the receiver are l x , l y and l z , respectively. According to the relationship between torque and force, the expression of the torque generated by the thrust is:
\begin{equation} {\vec{\!M}}_{f} {}^{\kern-2.7pt b} =\left[l_{x}\,f_{x},l_{y}\,f_{y},l_{z}\,f_{z}\right]^{T} \end{equation}
2.1.3 Aerodynamics
The receiver will generate aerodynamics under the action of airflow during flight. In the airflow coordinate system, the components of the aerodynamic forces acting on the receiver along the o v x v axis, o v y v axis and o v z v axis are usually represented by drag D, lateral force Y and lift L. The aerodynamic force received by the receiver can be divided into longitudinal aerodynamic force and transverse aerodynamic force.
2.1.3.1 Longitudinal aerodynamics
Longitudinal aerodynamics mainly affect the movement of the receiver in the o
b
x
b
z
b
plane of the body coordinate system, and the involved aerodynamics are drag D and lift L. The main variables affecting longitudinal aerodynamics are the attack angle
$\alpha$
, the pitch angle velocity
$\omega$
y
and the elevator deflection angle
$\delta$
e
of the receiver. In the airflow coordinate system, the expression of the longitudinal aerodynamic force acting on the receiver is [Reference Xiaoqi16]:
\begin{equation} \begin{cases}\displaystyle D={1 \over 2}\rho {V}_{a}^{2}S\left(C_{D}(\alpha)+C_{{D_{{\omega _{y}}}}}{c \over 2V_{a}}\omega _{y}+C_{{D_{{\delta _{e}}}}}\delta _{e}\right)\\[14pt]\displaystyle L={1 \over 2}\rho {V}_{a}^{2}S\left(C_{L}(\alpha)+C_{{L_{{\omega _{y}}}}}{c \over 2V_{a}}\omega _{y}+C_{{L_{{\delta _{e}}}}}\delta _{e}\right) \end{cases} \end{equation}
where
$\rho$
is the density of air, calculated by the formula
$\rho$
= 1.293* (actual pressure/standard physical atmospheric pressure)* (273.15/actual absolute temperature), and under standard conditions,
$\rho$
is usually taken as 1.293 kg/m3;
$V_{a}$
is the airspeed of the receiver; S is the reference area of the receiver wing; c is the average aerodynamic chord length of the receiver wing. For the calculation, formulas of C
D
(
$\alpha$
) and C
L
(
$\alpha$
), can be seen the reference [Reference Xiaoqi16].
According to the transformation relationship between the airflow coordinate system and the body coordinate system of the receiver [Reference Xiaoqi, Changsheng and Wuxing22], the expressions of drag D and lift L in the body coordinate system can be obtained:
\begin{equation} \left[\begin{array}{l} {F}_{A_{x}}^{b}\\[3pt] {F}_{A_{z}}^{b} \end{array}\right]=\left[\begin{array}{c@{\quad}c} \cos \alpha & {} -\sin \alpha \\[3pt] \sin \alpha & {} \cos \alpha \end{array}\right]\left[\begin{array}{l} -D\\[3pt] -L \end{array}\right]={1 \over 2}\rho {V}_{a}^{2}S \left[\begin{array}{l}\displaystyle C_{X}(\alpha)+C_{{X_{{\omega _{y}}}}}(\alpha){c \over 2V_{a}}\omega _{y}+C_{{X_{{\delta _{e}}}}}(\alpha)\delta _{e}\\[14pt] \displaystyle C_{Z}(\alpha)+C_{{Z_{{\omega _{y}}}}}(\alpha){c \over 2V_{a}}\omega _{y}+C_{{Z_{{\delta _{e}}}}}(\alpha)\delta _{e} \end{array}\right] \end{equation}
where
$C_{X}(\alpha )=-C_{D}(\alpha )\cos \alpha +C_{L}(\alpha )\sin \alpha;$
$C_{{X_{{\omega _{y}}}}}(\alpha )=-C_{{D_{{\omega _{y}}}}}\cos \alpha +C_{{L_{{\omega _{y}}}}}\sin \alpha$
;
$C_{{X_{{\delta _{e}}}}}(\alpha )=-C_{{D_{{\delta _{e}}}}}$
$\cos \alpha + C_{{L_{{\delta _{e}}}}}\sin \alpha;$
$C_{Z}(\alpha )=-C_{D}(\alpha )\sin \alpha -C_{L}(\alpha )\cos \alpha$
;
$C_{{Z_{{\omega _{y}}}}}(\alpha )=-C_{{D_{{\omega _{y}}}}}\sin \alpha +C_{{L_{{\omega _{y}}}}}\!\cos \alpha;$
$C_{{Z_{{\delta _{e}}}}}(\alpha )=$
$-C_{{D_{{\delta _{e}}}}}\sin \alpha -C_{{L_{{\delta _{e}}}}}\cos \alpha$
.
2.1.3.2 Lateral aerodynamics
The lateral aerodynamic force mainly affects the movement of the receiver in the o
b
y
b
axis direction of the body coordinate system, and simultaneously influences the rolling and yaw movements of the receiver. The involved aerodynamic force is the lateral force Y. The main variables affecting lateral aerodynamics are the side-slip angle
$\beta$
, the rolling angular velocity
$\omega$
x
, the yaw angular velocity
$\omega$
z
, the aileron rudder deflection angle
$\delta$
a
and the rudder deflection angle
$\delta$
r
of the receiver. In the body coordinate system, the expression of the lateral aerodynamic force acting on the receiver is [Reference Xiaoqi16]:
\begin{equation} {F}_{A_{y}}^{b}={1 \over 2}\rho {V}_{a}^{2}S\left(C_{{Y_{0}}}+C_{{Y_{\beta }}}\beta +C_{{Y_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{Y_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{Y_{{\delta _{a}}}}}\delta _{a}+C_{{Y_{{\delta _{r}}}}}\delta _{r}\right) \end{equation}
where b is the wingspan of the receiver wing;
$C_{{Y_{0}}} , C_{{Y_{\beta }}} , C_{{Y_{{\omega _{x}}}}} , C_{{Y_{{\omega _{z}}}}} , C_{{Y_{{\delta _{a}}}}}$
and
$C_{{Y_{{\delta _{r}}}}}$
are all coefficients.
Therefore, in the body coordinate system, the expression of the total aerodynamic force acting on the receiver is:
\begin{equation} {\overset{\rightharpoonup}{F}}_{\!A} {}^{\!\!b}=\left[\begin{array}{l} {F}_{A_{x}}^{b}\\[4pt] {F}_{A_{y}}^{b}\\[4pt] {F}_{A_{z}}^{b} \end{array}\right]={1 \over 2}\rho {V}_{a}^{2}S\left[\begin{array}{c} \displaystyle C_{X}(\alpha)+C_{{X_{{\omega _{y}}}}}(\alpha){c \over 2V_{a}}\omega _{y}+C_{{X_{{\delta _{e}}}}}(\alpha)\delta _{e}\\[8pt] \displaystyle C_{{Y_{0}}}+C_{{Y_{\beta }}}\beta +C_{{Y_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{Y_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{Y_{{\delta _{a}}}}}\delta _{a}+C_{{Y_{{\delta _{r}}}}}\delta _{r}\\[8pt] \displaystyle C_{Z}(\alpha)+C_{{Z_{{\omega _{y}}}}}(\alpha){c \over 2V_{a}}\omega _{y}+C_{{Z_{{\delta _{e}}}}}(\alpha)\delta _{e} \end{array}\right] \end{equation}
Correspondingly, the expression of the aerodynamic torque acting on the receiver in the body coordinate system can be obtained [Reference Xiaoqi, Changsheng and Wuxing22]:
\begin{equation} {\vec{\!M}}_{A} {}^{\!\!\!b}=\left[\begin{array}{l} {M}_{A_{x}}^{b}\\[3pt] {M}_{A_{y}}^{b}\\[3pt] {M}_{A_{z}}^{b} \end{array}\right]={1 \over 2}\rho {V}_{a}^{2}S\left[\begin{array}{c} \displaystyle b\left(C_{{x_{0}}}+C_{{x_{\beta }}}\beta +C_{{x_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{x_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{x_{{\delta _{a}}}}}\delta _{a}+C_{{x_{{\delta _{r}}}}}\delta _{r}\right)\\[11pt] \displaystyle c\left(C_{{y_{0}}}+C_{{y_{\alpha }}}\alpha +C_{{y_{{\omega _{y}}}}}{b \over 2V_{a}}\omega _{y}+C_{{y_{{\delta _{e}}}}}\delta _{e}\right)\\[11pt] \displaystyle b\left(C_{{z_{0}}}+C_{{z_{\beta }}}\beta +C_{{z_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{z_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{z_{{\delta _{a}}}}}\delta _{a}+C_{{z_{{\delta _{r}}}}}\delta _{r}\right) \end{array}\right] \end{equation}
Based on the above analysis, it can be known that in the body coordinate system, the expressions of the resultant external force
$\,\vec{\!F} {}^{\kern-0.2pt b}$
and the resultant external torque
$\vec{M}^{b}$
acting on the receiver are:
\begin{align} \,\vec{\!F} {}^{\kern0.2pt b} & ={\,\vec{\!F}}_{\!G} {}^{\kern-3pt b}+{\,\vec{\!F}}_{\!f} {}^{\kern-2pt b}+{\,\vec{\!F}}_{\!A} {}^{\kern-3pt b} = mg\left[\begin{array}{c} -\sin \theta \\[2pt] \cos \theta \sin \varphi \\[2pt] \cos \theta \cos \varphi \end{array}\right]\nonumber\\[11pt] & +\left[\begin{array}{c} f_{x}\\[2pt] f_{y}\\[2pt] f_{z} \end{array}\right]+ \displaystyle {1 \over 2}\rho {V}_{a}^{2}S\left[\begin{array}{c} \displaystyle C_{X}(\alpha)+C_{{X_{{\omega _{y}}}}}(\alpha){c \over 2V_{a}}\omega _{y}+C_{{X_{{\delta _{e}}}}}(\alpha)\delta _{e}\\[8pt] \displaystyle C_{{Y_{0}}}+C_{{Y_{\beta }}}\beta +C_{{Y_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{Y_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{Y_{{\delta _{a}}}}}\delta _{a}+C_{{Y_{{\delta _{r}}}}}\delta _{r}\\[8pt] \displaystyle C_{Z}(\alpha)+C_{{Z_{{\omega _{y}}}}}(\alpha){c \over 2V_{a}}\omega _{y}+C_{{Z_{{\delta _{e}}}}}(\alpha)\delta _{e} \end{array}\right] \end{align}
\begin{align} \vec{M}^{b} & = \,{\vec{\!M}}_{\!G} {}^{\kern-3.1pt b}+\, {\vec{\!M}}_{\!f} {}^{\kern-2pt b}+\, {\vec{\!M}}_{\!A} {}^{\kern-2.9pt b}\nonumber\\[5pt] & =\left[\begin{array}{c} l_{x}\,f_{x}\\[2pt] l_{y}\,f_{y}\\[3pt] l_{x}\,f_{z} \end{array}\right]\displaystyle +{1 \over 2}\rho {V}_{a}^{2}S\left[\begin{array}{c}\displaystyle b\left(C_{{x_{0}}}+C_{{x_{\beta }}}\beta +C_{{x_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{x_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{x_{{\delta _{a}}}}}\delta _{a}+C_{{x_{{\delta _{r}}}}}\delta _{r}\right)\\[12pt] \displaystyle c\left(C_{{y_{0}}}+C_{{y_{\alpha }}}\alpha +C_{{y_{{\omega _{y}}}}}{b \over 2V_{a}}\omega _{y}+C_{{y_{{\delta _{e}}}}}\delta _{e}\right)\\[12pt] \displaystyle b\left(C_{{z_{0}}}+C_{{z_{\beta }}}\beta +C_{{z_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{z_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{z_{{\delta _{a}}}}}\delta _{a}+C_{{z_{{\delta _{r}}}}}\delta _{r}\right) \end{array}\right] \end{align}
2.2 Dynamic model of the receiver
Through the analysis of the forces and torques acting on the receiver in Section 2.1, this section will respectively establish the translational dynamics model of the centre of mass and the rotational dynamics model around the centre of mass of the receiver in the body coordinate system based on Newtonr’s laws of motion and the theorem of moments of motion.
2.2.1 Translational dynamics model of the centre of mass
The receiver’s body coordinate system is adopted as the moving reference frame, and its velocity vector relative to the geocentric coordinate system is
$\, \vec{\!V}$
, and the angular velocity vector of rotation is
$\vec{\omega}$
. According to the relationship between the absolute derivative and the relative derivative of a vector expression, in a dynamic coordinate system, the expression for differentiating the velocity vector
$\, \vec{\!V}$
with respect to time is [Reference Weijia3]:
\begin{equation} {d\, \vec{\!V} \over dt}={\delta\, \vec{\!V} \over \delta t}+\vec{\omega}\times \vec{\!V} \end{equation}
where
${d\, \vec{\!V} \over dt}$
represents the absolute derivative of the velocity vector
$\, \vec{\!V}$
in the geocentric coordinate system;
${\delta \, \vec{\!V} \over \delta t}$
represents the relative derivative of the velocity vector
$\, \vec{\!V}$
in the machine coordinate system.
In the body coordinate system,
$\, \vec{\!V}=[v_{x},v_{y},v_{z}]^{T}$
is the velocity vector and
$\vec{\omega}=[\omega _{x},\omega _{y},\omega _{z}]^{T}$
is the angular velocity vector. Substituting them into Equation (13), we can obtain:
\begin{equation} {d\, \vec{\!V} \over dt}=\left[\begin{array}{l} \dot{v}_{x}\\ \dot{v}_{y}\\ \dot{v}_{z} \end{array}\right]+\left[\begin{array}{l} \omega _{x}\\ \omega _{y}\\ \omega _{z} \end{array}\right]\times \left[\begin{array}{l} v_{x}\\ v_{y}\\ v_{z} \end{array}\right]=\left[\begin{array}{l} \dot{v}_{x}+\omega _{y}v_{z}-\omega _{z}v_{y}\\ \dot{v}_{y}+\omega _{z}v_{x}-\omega _{x}v_{z}\\ \dot{v}_{z}+\omega _{x}v_{y}-\omega _{y}v_{x} \end{array}\right] \end{equation}
In the body coordinate system,
$\,\vec{\!F} {}^{\kern-0.2pt b}=[{F}_{\!x}^{b},{F}_{\!y}^{b},{F}_{\!z}^{b}]^{T}$
is the resultant external force acting on the receiver, according to Newtonr’s laws of motion, after calculation and organisation, Equation (14) can be rewritten as:
\begin{equation} \begin{cases}\displaystyle \dot{v}_{x}={{F}_{x}^{b} \over m}+\omega _{z}v_{y}-\omega _{y}v_{z}\\[10pt] \displaystyle \dot{v}_{y}={{F}_{y}^{b} \over m}-\omega _{z}v_{x}+\omega _{x}v_{z}\\[10pt] \displaystyle \dot{v}_{z}={{F}_{z}^{b} \over m}+\omega _{y}v_{x}-\omega _{x}v_{y} \end{cases} \end{equation}
where m represents the mass of the receiver.
Combining Equation (11), the translational dynamics model of the centre of mass of the receiver in the body coordinate system can be obtained:
\begin{align} & \left[\begin{array}{l} \dot{v}_{x}\\ \dot{v}_{y}\\ \dot{v}_{z} \end{array}\right]=\left[\begin{array}{l} \omega _{z}v_{y}-\omega _{y}v_{z}\\ \omega _{x}v_{z}-\omega _{z}v_{x}\\ \omega _{y}v_{x}-\omega _{x}v_{y} \end{array}\right]\nonumber\\[10pt] & +\displaystyle\left[\begin{array}{c} \displaystyle-g\sin \theta +{f_{x} \over m}+{1 \over 2m}\rho {V}_{a}^{2}S\left(C_{X}(\alpha)+C_{{X_{{\omega _{y}}}}}(\alpha){c \over 2V_{a}}\omega _{y}+C_{{X_{{\delta _{e}}}}}(\alpha)\delta _{e}\right)\\[13pt] \displaystyle {g}\cos \theta \sin \varphi +{f_{y} \over m}+{1 \over 2m}\rho {V}_{a}^{2}S\left(C_{{Y_{0}}}+C_{{Y_{\beta }}}\beta +C_{{Y_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{Y_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{Y_{{\delta _{a}}}}}\delta _{a}+C_{{Y_{{\delta _{r}}}}}\delta _{r}\right)\\[13pt] \displaystyle {g}\cos \theta \cos \varphi +{f_{z} \over m}+{1 \over 2m}\rho {V}_{a}^{2}S\left(C_{Z}(\alpha)+C_{{Z_{{\omega _{y}}}}}(\alpha){c \over 2V_{a}}\omega _{y}+C_{{Z_{{\delta _{e}}}}}(\alpha)\delta _{e}\right) \end{array}\right] \end{align}
2.2.2 Rotational dynamics model around the centre of mass
Similarly, taking the body coordinate system of the receiver as the dynamic coordinate system, based on the relationship between the absolute derivative and the relative derivative of the vector expression, the derivative of the angular momentum
$\vec{H}$
with respect to time can be expressed as [Reference Weijia3]:
\begin{equation} {d\vec{H} \over dt}={\delta \vec{H} \over \delta t}+\vec{\omega}\times \vec{H} \end{equation}
where
${d\vec{H} \over dt}$
represents the absolute derivative of angular momentum
$\vec{H}$
in the geocentric coordinate system;
${\delta \vec{H} \over \delta t}$
represents the relative derivative of the angular momentum
$\vec{H}$
in the body coordinate system.
In the body coordinate system, the angular momentum
$\vec{H}$
can be expressed as [Reference Weijia3]:
\begin{equation} \vec{H}=\vec{J}\cdot \vec{\omega}=\left[\begin{array}{c@{\quad}c@{\quad}c} J_{x} & -J_{xy} & -J_{xz}\\[2pt] -J_{xy} & {} J_{y} {}& -J_{yz}\\[2pt] -J_{xz} & {} -J_{yz} {}& J_{z} \end{array}\right]\left[\begin{array}{l} \omega _{x}\\[2pt] \omega _{y}\\[2pt] \omega _{z} \end{array}\right]=\left[\begin{array}{c} J_{x}\omega _{x}-J_{xy}\omega _{y}-J_{xz}\omega _{z}\\[2pt] -J_{xy}\omega _{x}+J_{y}\omega _{y}-J_{yz}\omega _{z}\\[2pt] -J_{xz}\omega _{x}-J_{yz}\omega _{y}+J_{z}\omega _{z} \end{array}\right] \end{equation}
where
$J_{x}$
,
$J_{y}$
and
$J_{z}$
, respectively, represent the moments of inertia of the receiver around the three axes o
b
x
b
, o
b
y
b
and o
b
z
b
of the body coordinate system, as well as the products of inertia
$J_{xy}$
,
$J_{yx}$
,
$J_{xz}$
,
$J_{zx}$
,
$J_{yz}$
and
$J_{zy}$
. And its calculation formulas are as follows [Reference Weijia3].
\begin{equation} \begin{cases} J_{x}=\displaystyle\int (y^{2}+z^{2})dm;\,J_{y}=\int (x^{2}+z^{2})dm;\,J_{z}=\int (y^{2}+x^{2})dm\\[8pt] J_{xy}= \displaystyle\int (xy)dm=J_{yx};\,J_{xz}=\displaystyle\int (xz)dm=J_{zx};\,J_{yz}=\int (yz)dm=J_{zy} \end{cases} \end{equation}
The receiver belongs to the face weight type of UAV, and its inertial products
$J_{xy}$
=
$J_{yx}$
=
$J_{yz}$
=
$J_{zy}$
= 0. Equation (19) can be rewritten as:
\begin{equation} \vec{H}=\vec{J}\cdot \vec{\omega}=\left[\begin{array}{c@{\quad}c@{\quad}c} J_{x} & 0 & -J_{xz}\\[2pt] 0 & {} J_{y} & {} 0\\[2pt] -J_{xz} & {} 0 {} & J_{z} \end{array}\right]\left[\begin{array}{c} \omega _{x}\\ \omega _{y}\\ \omega _{z} \end{array}\right]=\left[\begin{array}{c} J_{x}\omega _{x}-J_{xz}\omega _{z}\\[2pt] J_{y}\omega _{y}\\[2pt] -J_{xz}\omega _{x}+J_{z}\omega _{z} \end{array}\right] \end{equation}
In the body coordinate system,
$\vec{M}^{b}=[{M}_{x}^{b},{M}_{y}^{b},{M}_{z}^{b}]^{T}$
represents the resultant external torque acting on the receiver. According to the theorem of moment of motion, combined with Equations (17) and (20), after calculation and organisation, it can be obtained as follows:
\begin{equation} \left[\begin{array}{l} \dot{\omega }_{x}\\ \dot{\omega }_{y}\\ \dot{\omega }_{z} \end{array}\right]=\left[\begin{array}{c@{\quad}c@{\quad}c} J_{x} {}& 0 & {} -J_{xz}\\[2pt] 0 & {} J_{y} & {} 0\\[2pt] -J_{xz} & {} 0 & {} J_{z} \end{array}\right]^{-1}\left(\left[\begin{array}{c} {M}_{x}^{b}\\[4pt] {M}_{y}^{b}\\[4pt] {M}_{z}^{b} \end{array}\right]-\left[\begin{array}{c} (J_{z}-J_{y})\omega _{y}\omega _{z}-J_{xz}\omega _{x}\omega _{y}\\[4pt] (J_{x}-J_{z})\omega _{x}\omega _{z}+J_{xz}\big({\omega _{x}^{2}}-{\omega _{y}^{2}}\big)\\[4pt] (J_{y}-J_{x})\omega _{x}\omega _{y}+J_{xz}\omega _{y}\omega _{z} \end{array}\right]\right) \end{equation}
Combining Equation (12), the rotational dynamics model of the receiver around the centre of mass in the body coordinate system can be obtained:
\begin{align} \left[\begin{array}{c} \dot{\omega }_{x}\\ \dot{\omega }_{y}\\ \dot{\omega }_{z} \end{array}\right] & =\vec{\Gamma}\cdot \left[\begin{array}{c} J_{xz}\omega _{x}\omega _{y}-(J_{z}-J_{y})\omega _{y}\omega _{z}\\[3pt] J_{xz}\big({\omega _{y}^{2}}-{\omega _{x}^{2}}\big)-(J_{x}-J_{x})\omega _{x}\omega _{z}\\[3pt] -J_{xz}\omega _{y}\omega _{z}-(J_{y}-J_{x})\omega _{x}\omega _{y} \end{array}\right]+\vec{\Gamma}\cdot \left[\begin{array}{c} l_{x}f_{x}\\[2pt] l_{y}f_{y}\\[2pt] l_{x}f_{z} \end{array}\right]\nonumber\\[10pt] & + \displaystyle\vec{\Gamma}\cdot {1 \over 2}\rho {V}_{a}^{2}S\left[\begin{array}{c} \displaystyle b\bigg(C_{{x_{0}}}+C_{{x_{\beta }}}\beta +C_{{x_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{x_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{x_{{\delta _{a}}}}}\delta _{a}+C_{{x_{{\delta _{r}}}}}\delta _{r}\bigg)\\[10pt] \displaystyle c\bigg(C_{{y_{0}}}+C_{{y_{\alpha }}}\alpha +C_{{y_{{\omega _{y}}}}}{b \over 2V_{a}}\omega _{y}+C_{{y_{{\delta _{e}}}}}\delta _{e}\bigg)\\[10pt] \displaystyle b\bigg(C_{{z_{0}}}+C_{{z_{\beta }}}\beta +C_{{z_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{z_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{z_{{\delta _{a}}}}}\delta _{a}+C_{{z_{{\delta _{r}}}}}\delta _{r}\bigg) \end{array}\right] \end{align}
where
$\vec{\Gamma}=\left[\begin{array}{c@{\quad}c@{\quad}c} J_{x} & {} 0 & {} -J_{xz}\\[3pt] 0 & {} J_{y} & {} 0\\[3pt] -J_{xz} & {} 0 & {} J_{z} \end{array}\right]^{-1}=\displaystyle {1 \over J_{x}J_{y}J_{z}-J_{y}{J}_{xz}^{2}} \left[\begin{array}{c@{\quad}c@{\quad}c} J_{y}J_{z} & {} 0 {}& J_{y}J_{xz}\\[3pt] 0 & {} J_{x}J_{z}-{J}_{xz}^{2} & {} 0\\[3pt] J_{y}J_{xz} & {} 0 {}& J_{x}J_{y} \end{array}\right]$
.
To establish the dynamic model of the receiver with complete state variables, it is also necessary to respectively provide the positional motion model of the receiver in the geocentric coordinate system and the angular rotation model in the body coordinate system [Reference Bolian21]:
\begin{equation} \left[\begin{array}{c} \dot{x}\\ \dot{y}\\ \dot{z} \end{array}\right]=\left[\begin{array}{c} \cos \theta \cos \psi v_{x}+\left(\sin \theta \cos \psi \sin \varphi -\sin \psi \cos \varphi \right)v_{y}+\left(\sin \theta \cos \psi \cos \varphi +\sin \psi \sin \varphi \right)v_{z}\\[2pt] \cos \theta \sin \psi v_{x}+\left(\sin \theta \sin \psi \sin \varphi -\cos \psi \cos \varphi \right)v_{y}+\left(\sin \theta \sin \psi \cos \varphi -\cos \psi \sin \varphi \right)v_{z}\\[2pt] \sin \theta v_{x}-\cos \theta \sin \varphi v_{y}-\cos \theta \cos \varphi v_{z} \end{array}\right] \end{equation}
\begin{equation} \begin{cases} \dot{\varphi }=\omega _{x}+\tan \theta (\omega _{y}\sin \varphi +\omega _{z}\cos \varphi )\\ \dot{\theta }=\omega _{y}\cos \varphi -\omega _{z}\sin \varphi \\ \dot{\psi }=\sec \theta (\omega _{y}\sin \varphi +\omega _{z}\cos \varphi ) \end{cases} \end{equation}
Combining Equations (16) and (22) to (24), when
$\vec{X}=[x,y,z,v_{x},v_{y},v_{z},\varphi ,\theta ,\psi ,\omega _{x},\omega _{y},\omega _{z}]^{T}$
is taken as the state variable,
$\vec{u}=[f_{x},f_{y},f_{z},\delta _{a},\delta _{e},\delta _{r}]^{T}$
is taken as the control variable, the dynamic model of the receiver can be expressed as:
\begin{equation} \dot{\vec{X}}(t)=f(\vec{X}(t),\vec{u}(t)) \end{equation}
where
$f(\vec{X}(t),\vec{u}(t))=[A_{3\times 1},B_{3\times 1},C_{3\times 1},D_{3\times 1}]^{T}$
;
\begin{align*}A_{3\times 1}&=\left[\begin{array}{c} \cos \theta \cos \psi v_{x}+(\sin \theta \cos \psi \sin \varphi -\sin \psi \cos \varphi )v_{y}+(\sin \theta \cos \psi \cos \varphi +\sin \psi \sin \varphi )v_{z}\\[3pt] \cos \theta \sin \psi v_{x}+(\sin \theta \sin \psi \sin \varphi -\cos \psi \cos \varphi )v_{y}+(\sin \theta \sin \psi \cos \varphi -\cos \psi \sin \varphi )v_{z}\\[3pt] \sin \theta v_{x}-\cos \theta \sin \varphi v_{y}-\cos \theta \cos \varphi v_{z} \end{array}\right]\!;\\B_{3\times 1}&=\left[\begin{array}{c} \displaystyle\omega _{z}v_{y}-\omega _{y}v_{z}-g\sin \theta +{f_{x} \over m}+{d \over m}\left(C_{X}(\alpha)+C_{{X_{{\omega _{y}}}}}(\alpha){c \over 2V_{a}}\omega _{y}+C_{{X_{{\delta _{e}}}}}(\alpha)\delta _{e}\right)\\[6pt] \displaystyle\omega _{x}v_{z}-\omega _{z}v_{x}+{g}\cos \theta \sin \varphi +{f_{x} \over m}+{d \over m}\left(C_{{Y_{0}}}+C_{{Y_{\beta }}}\beta +C_{{Y_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{Y_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{Y_{{\delta _{a}}}}}\delta _{a}++C_{{Y_{{\delta _{r}}}}}\delta _{r}\right)\\[6pt] \displaystyle \omega _{y}v_{x}-\omega _{x}v_{y}-{g}\cos \theta \cos \varphi +{f_{z} \over m}+{d \over m}\left(C_{Z}(\alpha)+C_{{Z_{{\omega _{y}}}}}(\alpha){c \over 2V_{a}}\omega _{y}+C_{{Z_{{\delta _{e}}}}}(\alpha)\delta _{e}\right) \end{array}\right]\!;\\ C_{3\times 1} & =\left[\begin{array}{c} \omega _{x}+\tan \theta (\omega _{y}\sin \varphi +\omega _{z}\cos \varphi )\\[3pt] \omega _{y}\cos \varphi -\omega _{z}\sin \varphi \\[3pt] \sec \theta (\omega _{y}\sin \varphi +\omega _{z}\cos \varphi ) \end{array}\right];\\[8pt] D_{3\times 1} & =\left[\begin{array}{c} \displaystyle{\big(J_{y}J_{z}J_{xz}-{J}_{y}^{2}J_{xz}+J_{x}J_{y}J_{xz}\big)\omega _{x}\omega _{y}-\big(J_{y}{J}_{z}^{2}-{J}_{y}^{2}J_{z}+J_{y}{J}_{xz}^{2}\big)\omega _{y}\omega _{z} \over J_{x}J_{y}J_{z}-J_{y}{J}_{xz}^{2}}\\[11pt] \displaystyle {\left(J_{x}J_{z}J_{xz}-{J}_{xz}^{3}\right)\left({\omega _{y}^{2}}-{\omega _{x}^{2}}\right)-\left(J_{x}J_{z}-{J}_{xz}^{2}\right)\left(J_{x}-J_{x}\right)\omega _{x}\omega _{z} \over J_{x}J_{y}J_{z}-J_{y}{J}_{xz}^{2}}\\[11pt] \displaystyle {\left(J_{y}{J}_{xz}^{2}-J_{x}{J}_{y}^{2}+{J}_{x}^{2}J_{y}\right)\omega _{x}\omega _{y}-\left(J_{y}J_{z}J_{xz}-{J}_{y}^{2}J_{xz}+J_{x}J_{y}J_{xz}\right)\omega _{y}\omega _{z} \over J_{x}J_{y}J_{z}-J_{y}{J}_{xz}^{2}} \end{array}\right]\\[10pt] & \quad +\vec{\Gamma}\left[\begin{array}{c} \displaystyle l_{x}f_{x}+db\left(C_{{x_{0}}}+C_{{x_{\beta }}}\beta +C_{{x_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{x_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{x_{{\delta _{a}}}}}\delta _{a}+C_{{x_{{\delta _{r}}}}}\delta _{r}\right)\\[11pt] \displaystyle l_{y}f_{y}+dc\left(C_{{y_{0}}}+C_{{x_{\alpha }}}\alpha +C_{{y_{{\omega _{y}}}}}{b \over 2V_{a}}\omega _{y}+C_{{y_{{\delta _{e}}}}}\delta _{e}\right)\\[11pt] \displaystyle l_{z}f_{z}+db\left(C_{{z_{0}}}+C_{{z_{\beta }}}\beta +C_{{z_{{\omega _{x}}}}}{b \over 2V_{a}}\omega _{x}+C_{{z_{{\omega _{z}}}}}{b \over 2V_{a}}\omega _{z}+C_{{z_{{\delta _{a}}}}}\delta _{a}+C_{{z_{{\delta _{r}}}}}\delta _{r}\right) \end{array}\right]. \qquad\qquad\qquad\end{align*}
In addition, during the process of refueling operations, the receiver will be affected by atmospheric turbulence. The Dryden model is a commonly used atmospheric turbulence model in flight control, and its descriptive equation is [Reference Bolian21]:
\begin{equation} \begin{cases} \displaystyle\Phi _{x}={\sigma }_{x}^{2}\cdot {L_{x} \over \pi }\cdot {1 \over 1+(L_{x}\Omega )^{2}}\\[13pt] \displaystyle\Phi _{y}={\sigma }_{y}^{2}\cdot {L_{y} \over \pi }\cdot {1+12(L_{y}\Omega )^{2} \over [1+4(L_{y}\Omega )^{2}]^{2}}\\[13pt] \displaystyle\Phi _{z}={\sigma }_{z}^{2}\cdot {L_{z} \over \pi }\cdot {1+12(L_{z}\Omega )^{2} \over [1+4(L_{z}\Omega )^{2}]^{2}} \end{cases} \end{equation}
where Φ
x
, Φ
y
and Φ
z
are the three-axis spectral functions of atmospheric turbulence in the body coordinate system; Ω is the spatial frequency of atmospheric turbulence;
$\sigma$
x
,
$\sigma$
y
and
$\sigma$
z
are the turbulence intensity at a certain height and
$\sigma_x = \sigma_y = \sigma_z$
; L
x
, L
y
and L
z
are the turbulence scale at a certain height, and L
x
=L
y
=L
z
.
The atmospheric turbulence can be seen as external uncertainty disturbances experienced by the receiver. Therefore, Equation (25) can be rewritten as
\begin{equation} \dot{\vec{X}}(t)=f(\vec{X}(t),\vec{u}(t))+\varepsilon (t) \end{equation}
where
$\varepsilon$
(t) represents external uncertainty disturbances that include atmospheric turbulence.
At this point, through the analysis of the motion process of the receiver, the nonlinear dynamic model of the receiver as shown in Equation (27) is obtained. Our target is to design a controller based on the nonlinear model of Equation (27) to enable the receiver to actively approach the tanker smoothly and rapidly. It is required that the final flight speed and attitude of the receiver when successfully approaching the tanker be the same as those of the tanker, and maintain a certain safety distance and height deviation from the tanker.
3.0 Design of the NMPC method
In the actual process of AAR, we expect the receiver to approach the tanker actively and quickly in a short time and with low fuel consumption. This can not only save the aerial refueling time and improve the success rate of the receiver in performing subsequent tasks, but also reduce the exposure time of the tanker/receiver in the air, ensuring the safety of the two UAVs. Essentially, the approach control problem of the tanker/receiver can be regarded as an optimal control problem, that is, to seek an optimal control quantity
$\vec{u}^{\ast }(t)$
so that the approach control time of the tanker/receiver is short and the fuel consumption is low. This section will adopt the NMPC method to solve the optimal control problem of the active approach of the receiver to the tanker.
3.1 Objective function
NMPC is usually described in two forms: the discrete-time model and the continuous-time model. Compared with the continuous-time model, the discrete-time model is more robust in the case of small sample time intervals [Reference Chen and Tan23]. Therefore, the prediction model in the discrete-time state is used. Assuming the sampling time is T, the dynamic model of the receiver in the discrete-time state in the form of Equation (28) is shown as follows.
\begin{equation} \begin{cases} \vec{X}(k+1)=\vec{X}(k)+T\cdot f(\vec{X}(k),\Delta \vec{u}(k))+\vec{\varepsilon}(k)\\[3pt] \vec{u}(k)=\vec{u}(k-1)+\Delta \vec{u}(k) \end{cases} \end{equation}
where k is the sampling moment;
$\Delta \vec{u}(k)$
is the increment of the control input.
Assuming that
$\vec{\varepsilon}(k)$
remains constant in the predicted time domain, the following form of objective function can be established:
\begin{align} J&={\sum }_{i=1}^{N_{p}-1}\big[\vec{X}_{d}(k+i| k)-\hat{\vec{X}}(k+i| k)\big]^{T}Q\big[\vec{X}_{d}(k+i| k)-\hat{\vec{X}}(k+i| k)\big]\nonumber\\&\quad+{\sum }_{i=0}^{N_{u}-1}\left(\Delta \vec{u}\right)^{T}(k+i| k)R\left(\Delta \vec{u}\right)(k+i| k)\nonumber\\ &\quad+\big[\vec{X}_{d}(k+N_{p}| k)-\hat{\vec{X}}(k+N_{p}| k)\big]^{T}S\big[\vec{X}_{d}(k+N_{p}| k)-\hat{\vec{X}}(k+N_{p}| k)\big] \end{align}
where the first item and the second item are respectively the state cost and control input cost functions, which realise the active approach control from the receiver to the tanker; the third item is the terminal state cost function, which is used to ensure the stability of the NMPC controller;
$[\vec{X}_{d}(k+1| k),\vec{X}_{d}(k+2| k),\cdots ,\vec{X}_{d}(k+N_{p}| k)]^{T}$
is the reference value of the state vector of the receiver at time k, which can be obtained from the output value of the dynamic model of the tanker;
$[\hat{\vec{X}}(k+1| k),\hat{\vec{X}}(k+2| k),\cdots ,\hat{\vec{X}}(k+N_{p}| k)]^{T}$
is the predicted value of the state vector of the receiver at time k;
$[\Delta \vec{u}(k| k),\Delta \vec{u}(k+1| k),\cdots ,\Delta \vec{u}(k+N_{u}-1| k)]^{T}$
is the increment of the optimal control input of the receiver at moment k; N
p
and N
u
are respectively the lengths of the prediction and control time domains and
$N_{u}\leq N_{p};\,Q=diag[q_{1},q_{2},\cdots ,q_{12}]^{T} , R=diag[r_{1},r_{2},\cdots ,r_{6}]^{T}$
and
$S=diag[s_{1},s_{2},\cdots ,s_{12}]^{T}$
are respectively the weight matrices of the corresponding states, control inputs and terminal states.
3.2 Constraint conditions
During the process of the receiver actively approaching the tanker, its movement speed, attitude and the control amount applied to it will all be subject to certain constraints and limitations. At the same time, considering that the two UAVs need to maintain a certain safety distance when approaching and the receiver is located directly behind and below the tanker, the specific constraint conditions are expressed as:
\begin{equation} \vec{u}_{\min }\leq \vec{u}(k+i-1| k)\leq \vec{u}_{\max } \left(i=1,2,\cdots ,N_{u}\right) \end{equation}
where the subscripts min and max are respectively represented as the upper and lower boundaries of the constraints.
\begin{equation} \begin{cases} \, \vec{\!V}_{\min }\leq \hat{\vec{V}}(k+i| k)\leq \, \vec{\!V}_{\max }\\[3pt] \varphi _{\min }\leq \hat{\varphi }(k+i| k)\leq \varphi _{\max }\\[3pt] \theta _{\min }\leq \hat{\theta }(k+i| k)\leq \theta _{\max }\\[3pt] \psi _{\min }\leq \hat{\psi }(k+i| k)\leq \psi _{\max }\\[3pt] \vec{\omega}_{\min }\leq \hat{\vec{\omega}}(k+i| k)\leq \vec{\omega}_{\max } \end{cases} \qquad \left(i=1,2,\cdots ,N_{p}\right) \end{equation}
\begin{equation} \begin{cases} \sqrt{{l_{d}}^{2}-(\Delta y)^{2}-(\Delta z)^{2}}-\xi \leq x_{d}(k+N_{p}| k)-\hat{x}(k+N_{p}| k)\leq \sqrt{{l_{d}}^{2}-(\Delta y)^{2}-(\Delta z)^{2}}+\xi \\[3pt] \xi \leq y_{d}(k+N_{p}| k)-\hat{y}(k+N_{p}| k)\leq \xi \\[3pt] \sqrt{{l_{d}}^{2}-\left(\Delta x\right)^{2}-(\Delta y)^{2}}-\xi \leq z_{d}(k+N_{p}| k)-\hat{z}(k+N_{p}| k)\leq \sqrt{{l_{d}}^{2}-\left(\Delta x\right)^{2}-(\Delta y)^{2}}+\xi \\[3pt] \xi \leq \, \vec{\!V}_{d}(k+N_{p}| k)-\hat{\vec{V}}(k+N_{p}| k)\leq \xi \\[3pt] \xi \leq \varphi _{d}(k+N_{p}| k)-\hat{\varphi }(k+N_{p}| k)\leq \xi \\[3pt] \xi \leq \theta _{d}(k+N_{p}| k)-\hat{\theta }(k+N_{p}| k)\leq \xi \\[3pt] \xi \leq \psi _{d}(k+N_{p}| k)-\hat{\psi }(k+N_{p}| k)\leq \xi \end{cases} \end{equation}
where
$\Delta x=x_{d}(k+N_{p}| k)-\hat{x}(k+N_{p}| k);\,\Delta y=y_{d}(k+N_{p}| k)-\hat{y}(k+N_{p}| k);\,\Delta z=z_{d}(k+N_{p}| k)-\hat{z}(k+N_{p}| k)$
; l
d
is the safe distance between the tanker and receiver;
$\xi$
is a minimum value greater than 0.
Given the initial state
$\hat{\vec{X}}(0| 0)=\vec{X}\left(0\right)$
, combined with the dynamic model and constraint conditions, solving the optimal control variable
$\vec{u}^{*}(t)$
so that the objective function J takes the minimum value.
3.3 Solution of the optimal control quantity
After establishing the prediction model for the motion control of the receiver, in order to obtain the optimal control quantity, the NMPC needs to solve a nonlinear optimisation problem shown in Equation (29) online in each control step. This paper will adopt the adaptive gradient descent method [Reference Lugo-Cárdenas, Flores and Salazar5] to optimise the objective function online in order to obtain the optimal control quantity.
Let the state error
$\vec{e}(k+i| k)=\vec{X}_{d}(k+i| k)-\hat{\vec{X}}(k+i| k),(i=1,2,\cdots ,N_{p})$
, and the objective function (29) can be rewritten as:
\begin{align} J&={\sum }_{i=1}^{N_{p}-1}\vec{e}^{T}(k+i| k)Q\vec{e}(k+i| k)+{\sum }_{i=0}^{N_{u}-1}\left(\Delta \vec{u}\right)^{T}(k+i| k)R\left(\Delta \vec{u}\right)(k+i| k)\nonumber\\ &\quad +\vec{e}^{T}(k+N_{p}| k)S\vec{e}(k+N_{p}| k) \end{align}
Let
$\vec{E}(k)=[\vec{e}(k+1| k),\vec{e}(k+2| k),\cdots ,\vec{e}(k+N_{p}| k)]^{T}$
be the state error vector,
$\Delta \vec{U}(k)=[\Delta \vec{u}(k| k),$
$\Delta \vec{u}(k+1| k),\cdots ,\Delta \vec{u}(k+N_{u}-1| k)]^{T}$
be the control input increment vector,
$M=diag[Q_{1},Q_{2}\cdots ,$
$Q_{{N_{p}}-1},S]^{T}$
is the state weight matrix,
$H=diag[R_{1},R_{2}\cdots ,R_{{N_{u}}}]^{T}$
be the control weight matrix, and then we have
\begin{equation} J=\vec{E}^{T}(k)M\vec{E}(k)+(\Delta \vec{U})^{T}(k)H(\Delta \vec{U})(k) \end{equation}
The adaptive gradient descent method is adopted to optimise the objective function Equation (34), and the control increment vector
$\Delta \vec{U}(k)$
can be obtained. The expression of gradient descent is:
\begin{equation} \vec{U}(k+1)=\vec{U}(k)-\eta {\partial J \over \partial \vec{U}(k)}=\vec{U}(k)+2\eta \vec{P}^{T}(k)M\vec{E}(k)-2\eta H(\Delta \vec{U})(k) \end{equation}
where
$\eta$
is the adaptive optimisation rate of the gradient descent method.
According to Equation (35), the control incremental vector
$\Delta \vec{U}(k)$
can be obtained:
\begin{equation} \Delta \vec{U}(k)=\left(1+2\eta H\right)^{-1}2\eta \vec{P}^{T}(k)M\vec{E}(k) \end{equation}
where the expression of the matrix vector
$\vec{P}(k)$
is
\begin{equation} \vec{P}(k)=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \displaystyle{\partial \hat{\vec{X}}(k+1 | k) \over \partial \vec{u}(k| k)} {}& 0 {}& \cdots {}& 0\\[9pt] \displaystyle{\partial \hat{\vec{X}}(k+2| k) \over \partial \vec{u}(k| k)} & {} \displaystyle{\partial \hat{\vec{X}}(k+2| k) \over \partial \vec{u}(k+1 | k)} & {} \cdots & {} 0\\[6pt] \displaystyle \vdots {}& \vdots {}& \ddots & {} \vdots \\[3pt] \displaystyle {\partial \hat{\vec{X}}(k+N_{p}| k) \over \partial \vec{u}(k| k)} & {} \displaystyle{\partial \hat{\vec{X}}(k+N_{p}| k) \over \partial \vec{u}(k+1 | k)} & {} \cdots {}& \displaystyle{\partial \hat{\vec{X}}(k+N_{p}| k) \over \partial \vec{u}(k+N_{u}-1| k)} \end{array}\right] \end{equation}
It can be deduced from the above that in the process of solving the control increment vector, the inversion and multiplication of matrices are involved, which increases the computational load and complexity of NMPC. Under normal circumstances, we choose the first control quantity
$\vec{u}^{\ast }(k| k)$
in the control quantity
$\vec{U}^{\ast }(k)=[\vec{u}^{\ast }(k| k),\vec{u}^{\ast }(k+1 | k),\cdots ,\vec{u}^{\ast }(k+N_{u}-1| k)]^{T}$
obtained through the optimisation solution as the control input quantity for the next optimisation solution [Reference Sun, Zhang and Yang24]. Facts have proved that this choice can still enable NMPC to achieve satisfactory control effects [Reference Lu25]. Therefore, when the control time domain length N
u
is set to 1, Equation (36) can be rewritten as:
\begin{equation} \Delta \vec{u}(k| k)=\left(1+2\eta R_{1}\right)^{-1}2\eta {\vec{P}_{1}}^{T}(k)Q_{1}\vec{E}(k) \end{equation}
where
$\vec{P}_{1}(k)=\left[\displaystyle{\partial \hat{\vec{X}}(k+1 | k) \over \partial \vec{u}(k| k)},\,{\partial \hat{\vec{X}}(k+2| k) \over \partial \vec{u}(k| k)},\,\cdots ,\,{\partial \hat{\vec{X}}(k+N_{p}| k) \over \partial \vec{u}(k| k)}\right]^{T}\!$
.
To accelerate the convergence speed of the gradient optimisation method,
$\eta$
can adaptively adjust in the real-time control process, namely [Reference Hongbiao, Yu and Xiaoying17]:
\begin{equation} \eta _{c}={\kappa \over \lambda _{\max }\left(\vec{P}_{1}(k){\vec{P}^{T}}_{1}(k)\right)} \end{equation}
where
$0\lt \kappa \lt 2;\,\lambda _{\max }$
is the maximum eigenvalue of
$\vec{P}_{1}(k){\vec{P}^{T}}_{1}(k)$
.
Based on the above analysis, the optimal control increment
$\Delta \vec{u}^{\ast }(k| k)$
can be expressed as:
\begin{equation} \Delta \vec{u}^{\ast }(k| k)=\eta _{c}{\vec{P}_{1}}^{T}(k)\vec{E}(k)=\eta _{c}{\sum }_{i=1}^{N_{p}}\left({\partial \hat{\vec{X}}(k+i| k) \over \partial \vec{u}(k| k)}\vec{e}(k+i| k)\right) \end{equation}
Substituting the optimal control increment
$\Delta \vec{u}^{\ast }(k| k)$
into Equation (41), the optimal control quantity
$\vec{u}^{\ast }(k)$
at time k can be obtained:
\begin{equation} \vec{u}^{\ast }(k)=\vec{u}(k-1)+\Delta \vec{u}^{\ast }(k| k) \end{equation}
The optimal control quantity
$\vec{u}^{\ast }(k)$
is input into the dynamic model of the receiver; the NMPC performs a new round of predictive control based on the motion state of the tanker/receiver at time k + 1 until the control target requirements are met.
3.4 Stability analysis
Stability is a prerequisite for ensuring the normal operation of the control system. There are many methods for analysing the stability of control systems, such as Routhr’s criterion, root locus method and Lyapunovr’s theorem, etc. This paper will adopt the Lyapunov stability theorem to analyse the stability of the designed NMPC.
Selecting the Lyapunov quadratic positive definite function as shown in Equation (42):
\begin{equation} \overset{\frown }{\vec{V}}(k)=\hat{\vec{E}}^{T}(k)\hat{\vec{E}}(k) \end{equation}
When the difference calculation is performed on Equation (42), there is
\begin{align} \Delta \overset{\frown }{\vec{V}}(k) {} &=\overset{\frown }{\vec{V}}(k+1)-\overset{\frown }{\vec{V}}(k)\nonumber\\ {} &=\hat{\vec{E}}^{T}(k+1)\hat{\vec{E}}(k+1)-\hat{\vec{E}}^{T}(k)\hat{\vec{E}}(k)\nonumber\\ {} &=\big[\hat{\vec{E}}^{T}(k)+\Delta \hat{\vec{E}}^{T}(k)\big]\big[\hat{\vec{E}}(k)+\Delta \hat{\vec{E}}(k)\big]-\hat{\vec{E}}^{T}(k)\hat{\vec{E}}(k)\nonumber\\ {} &=\Delta \hat{\vec{E}}^{T}(k)\big[2\hat{\vec{E}}(k)+\Delta \hat{\vec{E}}(k)\big] \end{align}
Since
$\vec{E}(k)=[\vec{e}(k+1 | k),\vec{e}(k+2| k),\cdots ,\vec{e}(k+N_{p}| k)]^{T}$
, then there is
\begin{equation} {\partial \hat{\vec{E}}(k) \over \partial \vec{u}(k| k)}=-\left[{\partial \hat{\vec{X}}(k+1 | k) \over \partial \vec{u}(k| k)},\,{\partial \hat{\vec{X}}(k+2| k) \over \partial \vec{u}(k| k)},\,\cdots ,\,{\partial \hat{\vec{X}}(k+N_{p}| k) \over \partial \vec{u}(k| k)}\right]^{T}=-\vec{P}_{1}(k) \end{equation}
Therefore,
$\Delta \hat{\vec{E}}(k)$
can be expressed as:
\begin{equation} \Delta \hat{\vec{E}}(k)={\partial \hat{\vec{E}}(k) \over \partial \vec{u}(k| k)}\Delta \vec{u}(k| k)=-\eta _{c}\vec{P}_{1}(k){\vec{P}_{1}}^{T}(k)\hat{\vec{E}}(k) \end{equation}
Substituting Equation (45) into Equation (43) and rearranging, we can obtain:
\begin{align} \Delta \overset{\frown }{\vec{V}}(k) {} &=-\eta _{c}\big[\vec{P}_{1}(k){\vec{P}_{1}}^{T}(k)\hat{\vec{E}}(k)\big]^{T}\big[2\hat{\vec{E}}(k)-\eta _{c}\vec{P}_{1}(k){\vec{P}_{1}}^{T}(k)\hat{\vec{E}}(k)\big]\nonumber\\ {} &=-\eta _{c}\hat{\vec{E}}^{T}(k){\vec{P}_{1}}^{T}(k)\vec{P}_{1}(k)\big[2\vec{I}-\eta _{c}\vec{P}_{1}(k){\vec{P}_{1}}^{T}(k)\big]\hat{\vec{E}}(k) \end{align}
From Equation (44), it can be known that
$0\lt \eta _{c}\lambda _{\max }(\vec{P}_{1}(k){\vec{P}^{T}}_{1}(k))=\kappa \lt 2$
, then
$\Delta \overset{\frown }{\vec{V}}(k)\lt 0$
, which can prove that the designed NMPC method is theoretically stable and convergent.
4.0 Numerical simulation and result analysis
To verify the control performance of the NMPC method established in this paper during the active approach of the receiver to the tanker, a certain type of UAV was taken as the controlled object, and numerical simulation and result analysis were carried out in the Matlab/Simulink environment on a high-performance desktop computer. The configuration of the high-performance desktop computer is shown in Table 1 as follows.
Configuration of high-performance computers
Table 1. Long description
A table comparing specifications of high-performance computer components. The table has two rows and four columns. Column headers are Name and Specification model. Row 1: Processor, Intel i7-13700KF, Main frequency, 2.7 GHz. Row 2: Motherboard, Asus ROG Z790-A, Graphics card, Raptor ROG-STRIX-RTX 4084-O16G-W. Row 3: Memory, G.skill (16G*4) DDR5 7200, Hard disk, 2T*2 PCIE solid-state Samsung 990Pro. Row 4: Power supply, ROG STRIX Blowing Snow 1000W, Heat dissipation, ROG Rozen III 360ARGB white water cooler.
Physical parameters of the dynamic model [Reference Slegers, Kyle and Costello26, Reference Xu, Wang and Li27]
Table 2. Long description
A table with two columns and eight rows detailing the physical parameters of a dynamic model. The columns are labeled ‘Variable’ and ‘Value’. Row 1: Quality m, 2400 kg. Row 2: Wing area S, 8 m2. Row 3: Wing span b, 12 m. Row 4: Average aerodynamic chord length of the wing c, 6 m. Row 5: Moment of inertia Jx, 5.0 103 kg m2. Row 6: Moment of inertia Jy, 5.0 103 kg m2. Row 7: Moment of inertia Jz, 8.9 103 kg m2. Row 8: Moment of inertia Jzc, 0.006 103 kg m2.
Aerodynamic coefficient of the dynamic model [Reference Ostowari and Naik28]
Table 3. Long description
A table with ten rows and four columns listing aerodynamic coefficients for a dynamic model. The columns are labeled with variable names and their corresponding values. The table includes the following data: Row 1: C L 0, 0.400; C D 0, -0.100; C Y v, -0.100; C l 0, 2.000. Row 2: C m 0, -0.100; C D e, 0.001; C Y delta a, 0.150; C l delta a, 1.000. Row 3: C Y delta r, -0.100; C m delta a, 0.004; C Y delta b, 0.003; C l delta r, 0.150. Row 4: C l delta b, 1.500; C m delta r, -0.050; C Y 0, 0.050; C l 0, 0.001. Row 5: C m 0, 0.001; C Y beta, 0.002; C l beta, 0.200; C m beta, 1.500. Row 6: C Y p, 0.070; C l p, 0.001; C m p, 0.250; C Y n, 2.000. Row 7: C l p, 0.070; C l n, 0.070; C m n, 0.002; C Y delta t, 0.002. Row 8: C l delta t, 0.003; C m delta t, 1.293 kg/m3.
4.1 Parameter initialization
The tanker is assumed to maintain straight and level flight at a constant speed. In the geocentric coordinate system, the initial positions of the tanker are (4.25, 12, 6.58) km, the flight speed is 25 m/s, and the attitude angles (
$\varphi, \theta, \psi$
) are (0, 0, 0) degrees, respectively. In order to verify the effectiveness of the proposed method, the three initial positions of the receiver are be given as (0, 0, 5) km (condition 1), (1, 2, 3) km (condition 2) and (2, 2.5, 5) km (condition 3), the initial speed is 30 m/s, and the initial attitude angles (
$\varphi$
0
,
$\theta$
0
,
$\psi$
0
) are (0, 0, 66) degrees, respectively. The physical parameters and aerodynamic coefficient initialisations of the dynamic model of the receiver are shown in Tables 2 and 3, respectively.
The conditions that the receiver is subject to during its flight are shown in Table 4 as follows.
Variable constraints [Reference Reinhardt and Johansen29, Reference Mathisen, Gryte and Gros30]
Table 4. Long description
A table with two rows and five columns comparing minimum and maximum values for various flight variables. The columns are labeled Variable, Minimum value, Maximum value, Variable, Minimum value, and Maximum value. The rows list the variables and their respective minimum and maximum values. Row 1: Speed V, 0, 120 m/s, Rolling angle φ, -180 degrees, 180 degrees. Row 2: Pitch angle θ, -180 degrees, 180 degrees, Yaw angle ψ, -180 degrees, 180 degrees. Row 3: Engine thrust f, 0, 124.50 kN, Aileron rudder deflection angle δa, -35 degrees, 35 degrees. Row 4: Elevator deflection angle δe, -35 degrees, 35 degrees, Rudder deflection angle δr, -35 degrees, 35 degrees.
Our goal is that, based on the above parameters and constraints, the receiver will eventually fly at the same speed and in the same direction as the tanker and be located directly below and behind the tanker.
Rendezvous trajectory in condition 1.
4.2 Simulation result analysis
Suppose the length of the prediction domain N
p
= 20 and the length of the control domain N
u
= 1, the receiver eventually flies forward parallel to the ground horizontal line at a constant speed of 25 m/s, and the deviation between the flight position and the position of the tanker is
$\varDelta$
= (−0.05, 0.0, −0.08) km. According to the settings of the parameters in section 4.1, the numerical simulation results are as follows.
Figures 1–6 present the motion trajectory and the three-axis position errors of the tanker/receiver during autonomous approach and rendezvous with three different conditions, respectively. It can be seen from that although the initial positions of the receiver is different, the receiver can quickly achieve autonomous approach and rendezvous with the tanker, and the three-axis errors between the receiver and the tanker are eventually stabilized at −50 m, 0 m and −80 m, respectively, indicating that the receiver is directly below and behind the tanker, meeting the flight position requirements set by the numerical simulation.
The position errors of the three axes of the tanker/receiver.
Rendezvous trajectory in condition 2.
The position errors of the three axes of the tanker/receiver.
Rendezvous trajectory in condition 3.
The position errors of the three axes of the tanker/receiver.
Figures 7–10 present the change curves of flight speed, roll angle, pitch angle and yaw angle during the autonomous approach of the receiver to the tanker in condition 1, respectively. It can be seen that the receiver flies at an initial speed of 30 m/s, after approximately 30 s, it begins to reach the maximum speed of 115 m/s, and then flies towards the tanker at the maximum speed independently. After approximately 100 s, the receiver begins to decelerate and fly towards the tanker. Finally, after around 150 s, it decelerates to 25 m/s and maintains the same flight speed as the tanker. After a period of change, the flight attitude angle of the receiver remains consistent with that of the tanker at about 150 s. That is, the receiver eventually maintains the same flight direction as the tanker, meeting the requirements of the flight speed and direction of the tanker/receiver.
The flight speed of the receiver.
The rolling angle of the receiver.
The pitch angle of the receiver.
The yaw angle of the receiver.
Taking into account the impact of atmospheric turbulence on the tanker and receiver, a simulation is conducted to verify the control performance of the proposed method under turbulence effects. The simulation results are presented in Figures 11 and 12. It can be seen that although atmospheric turbulence can have a certain impact, the receiver is directly below and behind the tanker, meeting the flight position requirements set by the numerical simulation, indicating the control stability of the proposed method.
Rendezvous trajectory under turbulence.
The position errors of the three axes of the tanker/receiver.
To further verify the control performance of the proposed method, the relative distance between the tanker and the receiver is taken as the performance indicator. The simulation comparison results between the proposed method, LQR and MPC are shown in Figure 13. As seen in the figure, the NMPC-based control performance is the best.
The comparison results with different methods.
The simulation results show that the proposed method can meet the requirements for autonomous approach and rendezvous between the receiver and the tanker. The receiver can maintain a position directly behind and below the tanker while flying at the same speed as the tanker. The established method has good control performance, verifying its effectiveness.
5.0 Conclusion
In this paper, an autonomous approach and rendezvous method for unmanned tanker/receiver based on NMPC is established. Through the analysis and modelling of the resultant external force and torque acting on the receiver during its motion process, appropriate state variables and control variables are selected to establish a nonlinear dynamic model of the entire state of the receiver. The objective function with state terminals is designed by applying the NMPC theory. Combined with the constraint conditions, the adaptive gradient descent method is adopted to solve the optimal control quantity online, and the stability analysis is carried out. Numerical simulation verifies the control performance of the established method.
The simulation results show that the established method can quickly achieve the autonomous approach and rendezvous of the tanker and receiver in the air, ensuring that the receiver maintains the same speed and direction as the tanker directly behind and below it. The established method features high control accuracy, stable control process and excellent control performance.
Data and materials availability
The authors confirm that the data supporting the findings of this study are available within the article.
Acknowledgements
This work is supported by the Natural Science Basic Research Program-General Projects-Youth Projects of Shaanxi Provincial Department of Science and Technology under Grant No. 2025JC-YBQN-944, the Natural Science Basic Research Program-General Projects-Surface Projects of Shaanxi Provincial Department of Science and Technology under Grant No. 2025JC-YBMS-746, the Scientific Research Plan Project of the Education Department of Shaanxi Province-Youth Innovation Team Project under Grant No. 23JP071 and the Doctoral Research Foundation of Xihang University under Grant No. 2024KY0210. The authors would like to thank all editors and anonymous reviewers for improving this article.
Author contribution
H.S. was responsible for experimental conceptualisation and design, and was the main contributor to writing the manuscript; H.F., S.S.S. and X.Y.W. reviewed and revised the manuscript, and verified the experimental design; R.H.D. analysed and explained the experimental data. All authors have reviewed the manuscript.
Competing interests
The authors declare no competing interests.
