2.2.1 Formulation

If only the $m$ -th step (final step) matches the target altitude ( ${h_i}$ = ${h^{\rm{{\star}}}}$ ), ST will be $m$ . If trajectory points of $m$ -th and ( $m$ −1)-th step remain on target, ST becomes ( $m$ −1). Extending this logic: if the trajectory matches the target altitude for ( $m-s$ ) consecutive steps before termination, ST decreases as $m$ increases. The final step $m$ -th must satisfy ${h_i}$ = ${h^{\rm{{\star}}}}$ to minimise ST, making it the highest priority. Subsequent steps, such as ( $m$ −1)-th, ( $m$ −2)-th, are prioritised in reverse order, with previous steps (e.g. 1st step) having the lowest priority. This reasoning ensures that trajectory adjustments prioritise minimising ST by sequentially enforcing alignment from the final step backward.

Since the non-deviation state at $m$ -th step is of utmost importance, adjustments to altitude trajectory points from 1st step to ( $m$ −1)-th step should prioritise maintaining the non-deviation state at step $m$ . Specifically, the total reward loss caused by adjustments of 1st step to ( $m$ −1)-th step must be smaller than the reward gain from adjusting $m$ -th step to non-deviation state. This ensures optimisation tools like deep reinforcement learning can search for trajectories with higher rewards, driving the trajectory toward the non-deviation state at $m$ -th step. Thus, the inequality must hold as:

(2) \begin{align} {\sum_{(n = 1)}^{(n = m - 1)}} |\Delta {r_n}| \lt {(\Delta {r_m})^*} \end{align}

where ${\rm{\Delta }}{r_n}$ denotes the reward change caused by adjusting the altitude of the trajectory point at $n$ -th step. The adjustment can result in total of $A_2^2$ = 4 kinds of possible state transition, where the state can either deviate or remain non-deviating. ${\rm{\Delta }}{r_m}$ represents the reward gain from adjusting $m$ -th step from a deviated to a non-deviated state. The left-hand side of the inequality quantifies the cumulative reward loss from trajectory adjustments at 1st to ( $m$ −1)-th step, while the right-hand side represents the reward gain from trajectory adjustments at $m$ -th step. Given the diverse state transition possibilities at each step and the various distributions of trajectory points in non-deviating states, the values on both sides of the inequality lies within a specific range. To ensure the inequality holds, the maximum value of the left side should be less than the minimum value of the right side. The maximum value of the left side of the inequality (2) can be derived as:

(3) \begin{align} \max\!\left( \sum_{n=1}^{n = m - 1} \left| \Delta r_n \right| \right)&= \sum_{n=1}^{n = m - 1} \max\!\left( \left| \Delta r_n \right| \right) \nonumber \\[3pt] &= \sum_{n=1}^{n = m - 1} \max\!\left( \lambda^{n-1}- \left(1 - \min\!\left(\max\!\left(\left| \frac{h_i - h^*}{h^* - h_0} \right|\!, 0\right)\!, 1\right)\right)\right) \nonumber \\[3pt]&= \sum_{n=1}^{n = m - 1} \lambda^{n-1}.\end{align}

The minimum value of the right side of the inequality (2) can be deduced as:

(4) \begin{align}\begin{array}{ll}{{\rm{min}}\!\left( {{{\left( {{\rm{\Delta }}{r_m}} \right)}^{\ast}}} \right)} & = {\rm{min}}\!\left( {{\lambda ^{m - 1}} - \left( {1 - {\rm{min}}\!\left( {{\rm{max}}\!\left( {\left| {\dfrac{{{h_i} - {h^{\ast}}}}{{{h^{\ast}} - {h_0}}}} \right|\!,0} \right)\!,1} \right)} \right)} \right)\\[11pt]& = {\lambda ^{m - 1}} - 1 \end{array} \end{align}

A schematic illustration can help explain the physical meaning. In Fig. 5, the blue dashed line represents the trajectory where ST does not exist (denoted as $\infty $ ), achieving the maximum reward value. Except for the last trajectory point, which is slightly deviating from the target altitude, the reward value is $0.\dot 9$ that asymptotically approaches 1, all other trajectory points remain non-deviating, with reward values calculated as the exponential value ${\lambda ^{n - 1}}$ for the $n$ -th step.

Figure 5.

Trajectory whose ST does not exist with maximum reward and trajectory whose ST = $m$ with minimum reward.

The red dashed line represents the trajectory where ST = $m$ , achieving the minimum reward value. Except for the last point, which remains non-deviating at the target altitude with a reward value of ${\lambda ^{\left( {m - 1} \right)}}$ , all other trajectory points are deviating and remain the initial altitude, with rewards equal to zero.

The black dashed arrows indicate the change at each step for these two extreme cases. The goal is to transform ST from non-existence to ST = $m$ , i.e. to reduce ST. The reward changes from 1st to ( $m$ −1)-th step represent a loss, corresponding to the left side of the inequality (2), while the reward change for the last trajectory point represents a gain, corresponding to the right side of the inequality. To achieve this transition, the gain should logically exceed the loss.

In practice, since the WIG is underactuated, the current state affects the state that can be reached in the next step. Therefore, when adjusting one trajectory point, the other points also tend to change accordingly. In most cases, it is difficult to achieve the two extreme scenarios described above. The solid blue line represents a realistic possible trajectory in Fig. 5.

2.2.2 Feasibility

To guaranteed:

(5) \begin{align} \mathop \sum \limits_{n = 1}^{n = m - 1} {\lambda ^{n - 1}} \lt {\lambda ^{m - 1}} - 1 \end{align}

The right-hand side is a geometric sequence with common ratio $\lambda $ , first term ${\lambda ^{m - 2}}$ , and number of terms $\left( {m - 1} \right)$ . Its sum is:

(6) \begin{align}S = \frac{{{\lambda ^{m - 1}} - 1}}{{\lambda - 1}},{\rm{\;\;\;\;}}\lambda \ne 1 \end{align}

So the inequality becomes:

(7) \begin{align}{\lambda ^{m - 1}} - 1 \gt \frac{{{\lambda ^{m - 1}} - 1}}{{\lambda - 1}} \end{align}

Therefore, when $\lambda \gt 2$ , the inequality (2) holds. In practical scenarios, such as the underactuated nature of the WIG discussed in the engineering context, reaching the considered boundaries is often challenging. Therefore, the actual benefit is likely underestimated, while the loss is overestimated. As a result, the condition generally holds even when the exponential base value equals 2. Thus, $\lambda \geqslant 2$ , the inequality (2) is valid.

2.2.3 Extension

After ensuring that the final step, the $m$ -th step, does not deviate from the target altitude, the remaining problem is how to prevent the ( $m$ −1)-th trajectory point from deviation. Building on the preceding analysis, it can be proven likewise that:

(8) \begin{align} \mathop \sum \limits_{n = 1}^{n = m - 2} \!\left| {{\rm{\Delta }}{r_n}} \right| \lt \left( {{\rm{\Delta }}{r_{m - 1}}} \right) \end{align}

which ensures that the reward gain for the ( $m$ −1)-th trajectory point not deviating from the target altitude exceeds the loss from any adjustments in the previous $m$ −2 steps. The physical meaning is illustrated in Fig. 6, depicting extreme cases where the trajectory is ST = $m$ with maximised reward, and ST = $m$ −1 with minimised reward. While both trajectories exhibit no deviation at the final trajectory point, the former only deviates slightly at the ( $m$ −1)-th trajectory point, whereas the latter deviates at all trajectory points except the ( $m$ −1)-th trajectory point. The inequality (8) guarantees that the performance of the first $m$ −2 trajectory points can be sacrificed to facilitate the transition of the ( $m$ −1)-th trajectory point toward not deviating from the target altitude. By extension, this kind of inequality allows sequential adjustment of trajectory points from the end backward, gradually aligning them with the target altitude.

Figure 6.

Trajectory whose ST = $m$ with maximum reward and trajectory whose ST = $m$ –1 with minimum reward.

3.3.1 WIG model

The case is about a small WIG’s altitude change in the longitudinal plane, its model and coordinate frames for motion description, grid for aerodynamic calculation and kinetic analysis are shown in Fig. 8. Table 1 summarises the geometric configuration of the wing components for the WIG. It has three degrees of freedom: horizontal and vertical translations, and pitch. The two control inputs are the elevator deflection and throttle. The state variables include displacement, pitch angle, linear velocity, angular velocity, linear acceleration and angular acceleration. Among these, the flight altitude is extracted. The dynamic and kinematic models for calculating the WIG’s motion response are referenced from the study [Reference Hu, Zhang, Ding, Jiao, Zhang and Zhang4, Reference Hu, Li, Zhang and Zhang5].

Table 1.

Geometric parameters of wing parts in the WIG

Figure 8.

The virtual small WIG and its coordinate frames.

The WIG’s aerodynamics is simulated by vortex lattice method, and the ground effect is simulated by image method. A vortex element and a corresponding control point are assigned to each panel. The total force ${\bf{F}}$ and moment ${\bf{M}}$ with respect to the origin are obtained by summing the individual contributions from all vortex elements, as:

(14) \begin{align} {}&{{{\bf{F}}_i} = \rho {{\rm{\Gamma }}_i}\!\left( {{{\bf{V}}_\infty } + {{\bf{v}}_i}} \right) \times {{\bf{l}}_i}} \nonumber \\[3pt]&{{\bf{F}} \, \,= \mathop \sum \limits_{i = 1}^N {{\bf{F}}_i}} \nonumber \\[3pt] {}&{{\bf{M}} = \mathop \sum \limits_{i = 1}^N {{\bf{F}}_i} \times {{\bf{r}}_i}} \end{align}

where $\rho $ is the air density. ${{\rm{\Gamma }}_i}$ denotes the strength of the vortex, ${{\bf{l}}_i}$ represents the vector of the vortex segment oriented transversely, ${{\bf{V}}_\infty }$ is the freestream velocity vector, and ${{\bf{v}}_i}$ indicates the local velocity disturbance at the midpoint of the vortex segment.

Regarding the aerodynamics of the WIG, when an elevator deflection $\delta $ is executed during flight, the forces acting on the vehicle, as depicted in Fig. 8, can be decomposed into two distinct components. The first component consists of the aerodynamic forces and moment applied to the entire craft: lift ${F_L}$ , drag ${F_D}$ and pitching moment $M$ . These represent the resultant forces and moment at the centre of gravity, derived from the combined contributions of the horizontal stabiliser (lift ${F_{L,h.s}}$ , drag ${F_{D,h.s}}$ , pitching moment ${M_{h.s}}$ ), the vertical tail (lift ${F_{L,v.t}}$ , drag ${F_{D,v.t}}$ , pitching moment ${M_{v.t}}$ ) and the main wing (lift ${F_{L,m.w}}$ , drag ${F_{D,m.w}}$ , pitching moment ${M_{m.w}}$ ). The second component is the thrust $T$ generated by the propulsion system. Therefore, the total resultant forces and moment acting on the WIG craft within this planar force system are expressed as:

(15) \begin{align}\begin{array}{l} {}{{F_x} = {F_D} + T \cdot {\rm{cos}}\!\left( \alpha \right)}\\[6pt] {}{{F_z} = {F_L} + T \cdot {\rm{sin}}\!\left( \alpha \right)}\\[6pt] {}{{M_y} = M + T \cdot {L_2}}\end{array} \end{align}

where $\alpha $ denotes the pitch angle during actual flight operations, and ${L_2}$ represents the vertical distance between the centre of gravity $CG$ and the line of action of the thrust force when the pitch angle is ${0^ \circ }$ . Then, the WIG’s aerodynamics will be fed into the its motion simulation as:

(16) \begin{align}\begin{array}{l} {}{{\dot{\bf u}} = {\bf{u}} \times \omega + {\bf{C}}_{\bf{g}}^{\bf{b}} \cdot \left( {\frac{{\bf{F}}}{m} + {{\bf{g}}^{\bf{n}}}} \right)}\\[10pt] {}{\dot \omega = {{\bf{I}}^{ - 1}}\left( { - \omega \times {\bf{I}}\omega + {\bf{M}}} \right)}\end{array} \end{align}

where ${\bf{F}}$ is [ ${F_x}$ ,0, ${F_z}$ ] and ${\bf{M}}$ is [0, ${M_y}$ ,0]. ${\bf{u}}$ , $\omega $ are line and angle velocity, and ${\dot{\bf u}}$ and $\dot \omega $ are the corresponding accelerations. $m$ represents mass of the WIG, and ${{\bf{g}}^{\bf{n}}}$ signifies the acceleration of gravity. ${\bf{C}}_{\bf{g}}^{\bf{b}}$ is the coordinate frames’ transformation matrix. ${\bf{I}}$ denotes moment of inertia.

For the current validation scenario, the initial altitude of the WIG is ${h_0}{\rm{\;}}$ = 0.3 $c$ , and the target altitude is ${h^{\rm{{\star}}}}$ = 0.6 $c$ , where $c$ is the chord length of the wing, $c$ = 0.4 m. The entire trajectory consists of 24 steps, with each time step lasting 0.02 s. When measuring settling time, the acceptable deviation is $\varepsilon {\rm{\;}}$ = ${10^{ - 3}}\left| {{h^{\rm{{\star}}}} - {h_0}} \right|$ . The acceptable deviation parameter $\varepsilon $ is selected in this study to reduce computational costs associated with high-dimensional problems while maintaining adequate validation of the algorithm’s performance. For practical WIGs’ applications, the allowable deviation should be determined based on specific performance requirements, while also considering sensor errors (including altitude measurement inaccuracies) and sampling rate limitations.

The elevator’s deflection $\delta $ is constrained by both the deflection angle limits of $ \pm {10^ \circ }$ and the angle-of-attack limits of $ \pm {10^ \circ }$ . Assuming the current pitch angle is ${\alpha _c}$ , with horizontal and vertical velocities denoted as ${V_x}$ and ${V_y}$ respectively, the corresponding quasi-static angle-of-attack is ${\alpha _q}$ = ${\rm{arctan}}\!\left( {{V_y}/{V_x}} \right)$ . Thus, the feasible range for the elevator deflection falls within the interval [max( $ - {10^ \circ } - {\alpha _c} - {\alpha _q}$ , $ - {10^ \circ }$ ), min( ${10^ \circ } - {\alpha _c} - {\alpha _q}$ , ${10^ \circ }$ )] to prevent flow separation and stall. The other control input is engine thrust, and its maximum is $T$ = 7.06 N.

In practical engineering applications, factors such as sensor noise (which affects state acquisition accuracy), actuator delays (impacting timely control output), model uncertainties (causing deviations between expected and actual actions), and safety constraints (further limiting the manoeuvrability of the underactuated WIGs) do indeed introduce significant challenges. While these external constraints are not explicit requirements of the core algorithm proposed in this study, the algorithm’s performance is undoubtedly influenced by these practical factors. Typically, when the effects of such practical factors are modeled as Gaussian distributions, the proposed method would yield a single trajectory that maximises the expected reward in a statistical sense under such distributions. Provided these practical factors are appropriately modeled without ambiguity, the proposed method guarantees a unique solution.

3.3.2 Result

Conventional multi-objective function approaches typically involve testing multiple weight values through trial and error, to ensure the trajectory appear the smallest ST and the smallest IAE during the state transition. In this study, the weight $k$ is respectively set to multiple values between 0 and ${10^2}$ for constructing different reward functions.

To demonstrate the statistical significance of the deep reinforcement learning results, the 25th–75th percentile range is now displayed for each trajectory point based on eight training runs, indicating that 50% of the learning outcomes fall within this range. The resulting trajectories obtained via deep reinforcement learning are shown in Figs 9 and 10. Table 2 summarises their ST performance, including the 25th−75th percentile range of ST data from eight optimisation trials. Figure 9 presents a composite view of all trajectories, while Fig. 10 displays a magnified view of the trajectories corresponding to increasing weight value $k$ in ascending order, to highlight specific details. It can be found that as the weight of ST term increases, ST decreases gradually until stabilising at 13 steps.

Table 2.

Settling time (ST) of trajectories by multi-objective reward functions with varying $k$

Figure 9.

Overall trajectories by multi-objective reward functions ${R_{mo}}$ with varying $k$ .

Figure 10.

Magnified views of trajectories by multi-objective reward functions ${R_{mo}}$ with varying $k$ .

In the proposed reward function, $\lambda $ is required to be at least 2. $\lambda $ is sequentially set as 2, 3, 4 in ${r_{sc}}$ to construct the reward functions and put into the deep reinforcement learning for the corresponding trajectories, as shown in Fig. 11. The first two subfigures illustrate all trajectories overlaid in a single plot, whereas the following subfigures provide zoomed-in views of trajectory segments corresponding to progressively increasing base value $\lambda $ . Table 3 compares the settling time of trajectories achieved by the proposed reward function across different parameter configurations. It can be found that when $\lambda \geqslant 2$ , the corresponding trajectories stabilise at the 13-th step.

Figure 11.

Overall and magnified views of trajectories by the proposed reward function ${R_{sc}}$ with varying  $\lambda $ .

Unlike the conventional multi-objective function approach which requires multiple adjustments to the weight $k$ , the proposed reward function can directly reach the extreme without trial-and-error iterations. The three trajectories of $\lambda = 2$ , $\lambda = 3$ and $\lambda = 4$ should theoretically be identical, because under the same condition that ST has been minimised, the trajectory points before ST are guided solely by AE function. However, the resulting trajectory points do not perfectly coincide. The reason is that the WIG’s state space is very complex. Due to limited computational resources, deep reinforcement learning performs a limited search, making it challenging to achieve the globally maximum reward. Nevertheless, the trend of overall ST minimisation remains consistent. Therefore, a simpler example with a smaller state space is used in subsequent verifications to minimise IAE.

Table 3.

Settling time (ST) of trajectories by the proposed reward functions with varying  $\lambda $

3.3.3 Notes

For the conventional approach, achieving the target trajectory requires a trial-and-error process of adjusting the weighting coefficient across values such as $k = 0$ , $k = {10^{ - 4}}$ , $k = {10^{ - 3}}$ , $k = {10^{ - 2}}$ , $k = {10^{ - 1}}$ , $k = {10^0}$ and $k = {10^1}$ and $k = {10^2}$ . As $k$ increases, the trajectory gradually evolves and stabilises at the target. This process exemplifies the time-consuming tuning that conventional methods require. In stark contrast, the proposed method directly generates this target trajectory without the need for any iterative parameter adjustments, thus eliminating this costly trial-and-error phase.

Table 4.

Optimisation iteration counts for trajectory generation by multi-objective ${R_{mo}}$ and the proposed ${R_{sc}}$ reward function

Figure 12 comparatively plot altitude-distance trajectories by multi-objective ${R_{mo}}$ and the proposed ${R_{sc}}$ reward function in the spatial domain. The core objective of the proposed method is not primarily to achieve the fastest possible obstacle avoidance, but rather to prioritise rapid stabilisation at the target altitude (minimising settling time) first, followed by minimisation of the overall deviation. An overly aggressive climb would typically cause significant overshoot, resulting in increased settling time and larger deviation, which contradicts our design goals. The proposed method is applicable to both climbing and descending manoeuvres, achieving fast stabilisation with minimal deviation in scenarios which involves climbing pre-obstacle, descending post-obstacle or adjusting altitude based on mission efficiently. To prevent obstacle collisions, sufficient safety distance must be maintained before initiating avoidance manoeuvres. Similarly, adequate safety altitude must be preserved during descent to prevent water impact.

Figure 12.

Altitude-distance trajectories by multi-objective ${R_{mo}}$ and the proposed ${R_{sc}}$ reward function in the spatial domain.

Table 4 presents the parameter tuning cost of the conventional methods along with lower cost of the proposed method. The conventional multi-objective reward function requires eight parameter-tuning attempts and approximately 67,552 iterations to identify the target trajectory. In contrast, the proposed reward function achieves the same objective in only 7,522 iterations without requiring any parameter tuning. This provides substantial advantage in computational efficiency when comparing single optimisation against multiple tuning attempts.

The variations in elevator deflection and engine thrust for trajectories obtained by both the conventional method and the proposed method are illustrated in Fig. 13. While the differences are minor during the initial phase, significant discrepancies emerge during the stabilisation phase.

Figure 13.

Control input by multi-objective ${R_{mo}}$ and the proposed ${R_{sc}}$ reward function in the time domain.

3.4.1 Model

A simple second-order response model for single-degree-of-freedom altitude adjustment is used for validation. The model’s motion response is defined as:

(17) \begin{align}\begin{array}{l} {}{a = \frac{{F - mg}}{m}}\\[5pt] {}{v = v + a{\rm{\Delta }}t}\\[5pt]{h = h + v{\rm{\Delta }}t}\end{array} \end{align}

where $m$ = 1 kg, $F \in \left[ {0,20} \right]$ N, $g\;$ = 10 ${\rm{m}}/{{\rm{s}}^2}$ and ${\rm{\Delta }}t\;$ = 0.1 s. $a$ is the acceleration, and $v$ is the velocity. The initial altitude is ${h_0}$ = 0 m and ${v_0}{\rm{\;}}$ = 0 m/s. The target altitude is ${h^{\rm{{\star}}}}$ = 0.5 m, and the acceptable deviation ratio is $\varepsilon \;$ = ${10^{ - 2}}\left| {{h^{\rm{{\star}}}} - {h_0}} \right|$ .

3.4.2 Result

The control group is set up using a single-objective function ${r_s}$ that only considers the minimisation of ST. Since ${r_s}$ ensures the fastest attainment of the target altitude without constraining the performance of the transition, it results in multiple trajectories with varying IAE values. Through deep reinforcement learning, total of 200 trajectories are obtained for ${r_s}$ to reach the extremum, while ${r_{sc}}$ yields a single trajectory. Their comparison is shown in Fig. 14. It can be observed that 200 trajectories produced by ${r_s}$ are dispersed. The proposed reward, under the condition of achieving the minimum ST, directly obtains the trajectory with the minimum IAE.

Figure 14.

Trajectories obtained by the single-objective function ${r_s}$ and the proposed ${r_{sc}}$ .