3.1 Parafoil system description

After the parafoil is fully inflated, the aerodynamic force is changed by manipulating the deformation of the canopy, thereby achieving flight along a predetermined path. Since the shape of the parafoil does not change significantly, the system can be approximately treated as a two-body rootless system composed of the canopy and payload. Denote the canopy, parachute ropes, etc. as rigid body $ A$ , the suspension straps, payload, etc. as rigid body $ B$ , and the connecting rope at the intersection point as a spherical joint $ M$ . Then the system has a total of nine degrees of freedom, as shown in Fig. 2.

Figure 2.

Two-body parafoil system.

The following reference systems are required in the specific modelling process, as shown in Fig. 3.

1. (1) Earth-fixed reference frame $ {S}_{e}\left({O}_{e}{X}_{e}{Y}_{e}{Z}_{e}\right)$ .

Figure 3.

Reference system.

The origin $ {O}_{e}$ is at the intended impact point. The $ {O}_{e}{Z}_{e}$ axis is vertically down. The $ {O}_{e}{X}_{e}$ axis lies in the canopy’s longitudinal symmetry plane, perpendicular to $ {O}_{e}{Z}_{e}$ , and points toward the East.

1. (2) Canopy reference frame $ {S}_{p}\left({O}_{p}{X}_{p}{Y}_{p}{Z}_{p}\right)$ .

The canopy reference system is the body-fixed reference system of the canopy. The origin $ {O}_{P}$ is at the canopy’s mass centre. The $ {O}_{P}{X}_{P}$ axis lies in the canopy’s longitudinal symmetry plane, parallel to the canopy’s chord length, and points forward. The $ {O}_{P}{Z}_{P}$ axis is also in the longitudinal symmetry plane of the canopy.

1. (3) Payload reference frame $ {S}_{c}\left({O}_{c}{X}_{c}{Y}_{c}{Z}_{c}\right)$ .

The payload reference system is the body-fixed reference system of the payload. The origin $ {O}_{c}$ is at the payload’s mass centre. The $ {O}_{c}{X}_{c}$ axis is parallel to the payload’s central axis and points forward of the parafoil system. The $ {O}_{c}{Z}_{c}$ axis lies in the payload’s longitudinal symmetry plane.

1. (4) Airstream reference frame $ {S}_{q}\left({O}_{q}{X}_{q}{Y}_{q}{Z}_{q}\right)$ .

The origin $ {O}_{q}$ is at the parafoil system mass centre. The $ {O}_{Q}{X}_{Q}$ axis is opposite the oncoming flow. The $ {O}_{Q}{Z}_{Q}$ axis lies in the parafoil system longitudinal symmetry plane, perpendicular to $ {O}_{Q}{X}_{Q}$ , and points down.

3.2 Dynamics and kinematics equations

To facilitate the design of the parafoil path planning, the Lagrange equations in the form of quasi-coordinates are adopted for dynamic modelling. For the parafoil system shown in Fig. 2, the position components $ {x}_{e}$ , $ {y}_{e}$ and $ {z}_{e}$ of the position vector in the earth-fixed reference frame, the Euler angles $ {\gamma }_{p}$ , $ {\psi }_{p}$ and $ {\theta }_{p}$ of the canopy, as well as the Euler angles $ {\gamma }_{c}$ , $ {\psi }_{c}$ and $ {\theta }_{c}$ of the payload are set as the quasi-coordinates of the parafoil system.

The velocity components $ {v}_{x}$ , $ {v}_{y}$ and $ {v}_{z}$ of the mass centre velocity vector of the system in $ {S}_{e}\left({O}_{e}{X}_{e}{Y}_{e}{Z}_{e}\right)$ , the angular velocity components $ {\omega }_{px}$ , $ {\omega }_{py}$ and $ {\omega }_{pz}$ of the canopy, and the angular velocity components $ {\omega }_{cx}$ , $ {\omega }_{cy}$ and $ {\omega }_{cz}$ of the payload are taken as the quasi-velocities.

With the quasi-coordinates and quasi-velocities defined, the matrix-form representation can now be derived, as shown in Equations (1) and (2).

(1) \begin{align} \boldsymbol{q}={\left[\begin{array}{c@{\quad}c@{\quad}c}\boldsymbol{R} {}& {\boldsymbol{E}}_{c} {}& {\boldsymbol{E}}_{p} \end{array}\right]}^{\mathrm{T}},\boldsymbol{\nu }={\left[\begin{array}{c@{\quad}c@{\quad}c}{\boldsymbol{V}}_{o} & {} {\boldsymbol{\omega }}_{c} {}& {\boldsymbol{\omega }}_{p} {} {} \end{array}\right]}^{\mathrm{T}}, \end{align}

where,

(2) \begin{align} \left\{\begin{array}{l}\boldsymbol{R}={\left[\begin{array}{c@{\quad}c@{\quad}c}{x}_{e} & {y}_{e} & {z}_{e}\end{array}\right]}^{\mathrm{T}} \\[4pt] {\boldsymbol{V}}_{o}={\left[\begin{array}{c@{\quad}c@{\quad}c}{V}_{ox} & {} {V}_{oy} {}& {V}_{oz}\end{array}\right]}^{\mathrm{T}} \\[4pt] {\boldsymbol{E}}_{p}={\left[\begin{array}{c@{\quad}c@{\quad}c}{\gamma }_{p} {}& {\psi }_{p} & {} {\theta }_{p}\end{array}\right]}^{\mathrm{T}} \\[4pt] {\boldsymbol{\omega }}_{p}={\left[\begin{array}{c@{\quad}c@{\quad}c}{\omega }_{px} {}& {\omega }_{py} & {} {\omega }_{pz}\end{array}\right]}^{\mathrm{T}}. \\[4pt] {\boldsymbol{E}}_{c}={\left[\begin{array}{c@{\quad}c@{\quad}c}{\gamma }_{} {}& {\psi }_{c} & {}{\theta }_{c}\end{array}\right]}^{\mathrm{T}} \\[4pt] {\boldsymbol{\omega }}_{c}={\left[\begin{array}{c@{\quad}c@{\quad}c}{\omega }_{cx} {} & {\omega }_{cy} {}& {\omega }_{cz}\end{array}\right]}^{\mathrm{T}}\end{array}\right. \end{align}

The multi-body dynamics method is adopted to conduct parafoil modelling. The detailed derivation process is available in Zhang et al. [Reference Zhang, Gao, Guo, Feng and Ge14], where relevant steps and explanations are presented for further exploration. By integrating the kinematic and dynamic relationships, the Lagrange equations are established as shown in Equation (3).

(3) \begin{align} \left\{\begin{array}{l}\dot{\boldsymbol{q}}=\boldsymbol{B}\boldsymbol{\nu }\\[4pt] \boldsymbol{M}\left(\boldsymbol{q}\right)\dot{\boldsymbol{\nu }}=\boldsymbol{F} ,\end{array}\right. \end{align}

the state transition matrix $ \boldsymbol{B}$ is expressed as follows,

(4) \begin{align} \boldsymbol{B}=diag\left({\boldsymbol{E}}_{3\times 3}^{\mathrm{T}} \quad {\boldsymbol{D}_{c}^{-1} \quad {\boldsymbol{D}}}_{p}^{-1}\right)\!. \end{align}

(5) \begin{align} {\boldsymbol{D}}_{c}=\left[\begin{array}{c@{\quad}c@{\quad}c}1 {} & -\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_{c} \mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{c} & {} \mathrm{t}\mathrm{a}\mathrm{n}{\theta }_{c} \mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_{c}\\ 0 {} & \mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{c}/\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{c} & {} -\mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_{c}/\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{c}\\ 0 {} & \mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_{c} {} & \mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{c}\end{array}\right]\!, \end{align}

(6) \begin{align} {\boldsymbol{D}}_{p}=\left[\begin{array}{c@{\quad}c@{\quad}c}1 & {} -\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_{p} \mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{p} & {} \mathrm{t}\mathrm{a}\mathrm{n}{\theta }_{p} \mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_{p}\\ 0 {} & \mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{p}/\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{p} {} & -\mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_{p}/\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{p}\\ 0 {} & \mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_{p} & {} \mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_{p}\end{array}\right]\!, \end{align}

the generalised mass matrix $ \boldsymbol{M}$ is expressed as follows,

(7) \begin{align} \boldsymbol{M}=\left[\begin{array}{c@{\quad}c@{\quad}c}({m}_{p}+{m}_{c}){\boldsymbol{E}}_{3\times 3}+{\boldsymbol{M}}_{a} {}& -{\boldsymbol{B}}_{c}^{p}{m}_{c}\tilde{\boldsymbol{R}}_{c} & {} -({m}_{p}{\boldsymbol{E}}_{3\times 3}+{\boldsymbol{M}}_{a})\tilde{\boldsymbol{R}}_{p}\\ {m}_{c}\tilde{\boldsymbol{R}}_{c} {} & {\boldsymbol{I}}_{c} & {} {0}_{3\times 3}\\ -({m}_{p}{\boldsymbol{E}}_{3\times 3}+{\boldsymbol{M}}_{a})\tilde{\boldsymbol{R}}_{p} {}& {0}_{3\times 3} & {} {\boldsymbol{I}}_{p}+{\boldsymbol{I}}_{a}\end{array}\right]\!, \end{align}

the generalised force matrix $ \boldsymbol{F}$ is expressed as follows,

(8) \begin{align} \boldsymbol{F} {}& =\left[\begin{array}{c}-({m}_{c}+{m}_{p}){\boldsymbol{E}}_{3\times 3}\left(\tilde{\boldsymbol{\omega}}_{p}{\boldsymbol{B}}_{d}^{p}{\boldsymbol{V}}_{o}\right)\\ -\tilde{\boldsymbol{\omega}}_{c}\left({\boldsymbol{I}}_{c}{\boldsymbol{\omega }}_{c}\right)-{m}_{c}{\boldsymbol{R}}_{c}\left(\tilde{\boldsymbol{\omega}}_{c}{\boldsymbol{B}}_{d}^{c}{\boldsymbol{V}}_{o}\right)\\ -\tilde{\boldsymbol{\omega}}_{p}\left({\boldsymbol{I}}_{p}{\boldsymbol{\omega }}_{p}\right)-{m}_{p}\tilde{\boldsymbol{R}}_{p}\left(\tilde{\boldsymbol{\omega}}_{p}{\boldsymbol{B}}_{d}^{p}{\boldsymbol{V}}_{o}\right)\end{array}\right]+\left[\begin{array}{c}{m}_{c}\boldsymbol{g}+{m}_{p}\boldsymbol{g}\\ \tilde{\boldsymbol{R}}_{c}\left({m}_{c}\boldsymbol{g}\right)\\ \tilde{\boldsymbol{R}}_{p}\left({m}_{p}\boldsymbol{g}\right)\end{array}\right]\nonumber\\[8pt] {} & +\left[\begin{array}{c}{\boldsymbol{F}}_{p}+{\boldsymbol{B}}_{c}^{p}{\boldsymbol{F}}_{c}+{\boldsymbol{F}}_{AM}\\ {\boldsymbol{M}}_{c}\\ {\boldsymbol{M}}_{p}+{\boldsymbol{M}}_{AM}\end{array}\right]-\left[\begin{array}{c}{m}_{c}{\boldsymbol{B}}_{c}^{p}\big[\tilde{\boldsymbol{R}}_{c}\big(\tilde{\boldsymbol{\omega}}_{c}{\boldsymbol{\omega }}_{c}\big)\big]+{m}_{p}\big[\tilde{\boldsymbol{R}}_{p}\big(\tilde{\boldsymbol{\omega}}_{p}{\boldsymbol{\omega }}_{p}\big)\big]\\ {0}_{3\times 3}\\ {0}_{3\times 3}\end{array}\right]\!,\end{align}

The symbols in this equation are defined in Equations (7) and (8) as follows: Transformation matrices include $ {\boldsymbol{B}}_{e}^{p}$ (from earth-fixed to canopy frame), $ {\boldsymbol{B}}_{e}^{c}$ (from earth-fixed to payload frame), and $ {\boldsymbol{B}}_{c}^{p}$ (from payload to canopy frame). Matrix notations comprise $ {\boldsymbol{E}}_{3\times 3}$ (3 × 3 identity matrix) and $ {0}_{3\times 3}$ (3 × 3 zero matrix). Mass properties are $ {m}_{p}$ (canopy mass) and $ {m}_{c}$ (payload mass), while $ {\boldsymbol{I}}_{p}$ and $ {\boldsymbol{I}}_{c}$ are inertia matrices about the spherical joint. Position vectors $ {\boldsymbol{R}}_{p}$ and $ {\boldsymbol{R}}_{c}$ extend from the canopy and payload mass centres to the spherical joint, respectively. Aerodynamic forces and moments are $ {\boldsymbol{F}}_{p},{\boldsymbol{M}}_{p}$ (at canopy mass centre) and $ {\boldsymbol{F}}_{c},{\boldsymbol{M}}_{c}$ (at payload mass centre), with $ {\boldsymbol{F}}_{AM},{\boldsymbol{M}}_{AM}$ representing added mass effects governed by $ {\boldsymbol{M}}_{a}$ (added mass matrix) and $ {\boldsymbol{I}}_{a}$ (added inertia matrix). $ {\boldsymbol{V}}_{ox},{\boldsymbol{V}}_{oy},{\boldsymbol{V}}_{oz}$ are velocity components of mass centre.

3.3 State equation

From Equations (1) to (8), the complete expression of the dynamic and kinematic equations for the parafoil system has been derived. To conduct the sparse identification of the system more conveniently, the model is further transformed into a nonlinear state equation suitable for the design of the control system, and the system state variables are defined as follows,

(9) \begin{align} \boldsymbol{x}={\left[\begin{array}{c@{\quad}c}\boldsymbol{q} & {} \boldsymbol{v}\end{array}\right]}^{\mathrm{T}}, \end{align}

define $ {\delta }_{l}$ and $ {\delta }_{r}$ as the trailing edges deflection on the left and right of parafoil, which are used as the control inputs and expressed as,

(10) \begin{align} \boldsymbol{u}={\left[\begin{array}{c@{\quad}c}{\delta }_{l} & {} {\delta }_{r}\end{array}\right]}^{T}, \end{align}

subsequently, the state equations are obtained as follows,

(11) \begin{align} \frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{x}=f(\boldsymbol{x},\boldsymbol{u}), \end{align}

For the parafoil dynamic model, Equation (11) can be expressed in matrix form as,

(12) \begin{align} \left[\begin{array}{c}{\dot{x}}_{1}(t)\\ {\dot{x}}_{2}(t)\\ \vdots \\ {\dot{x}}_{18}(t)\end{array}\right]=\left[\begin{array}{c}{f}_{1}(\boldsymbol{x},\boldsymbol{u})\\ {f}_{2}(\boldsymbol{x},\boldsymbol{u})\\ \vdots \\ {f}_{18}(\boldsymbol{x},\boldsymbol{u})\end{array}\right]\!, \end{align}

where $ {x}_{1}\sim {x}_{18}$ correspond to the 18 state parameters at time $ t$ in Equation (2). Given the initial state quantities and control inputs, the motion characteristics of the two-body parafoil system can be obtained via numerical integration methods.

4.1 Data collection

1. 1. Full-state historical data

   Collect the data of the control input, state variables and their derivatives in the time series of the 9 DOF dynamic model of the parafoil system in the simulation and build a historical dataset. These data are arranged according to the state variables to form three data matrices [Reference Slegers and Costello9], $ X$ , $ \dot{X}$ and $ U$ , as shown in Equation (13).

   (13) \begin{align} \left\{\begin{array}{c}X=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\boldsymbol{x}({t}_{1}) & {} \boldsymbol{x}({t}_{2}) {} & {} \cdots {} & {} \boldsymbol{x}({t}_{m})\end{array}\right]\in {\mathbb{R}}^{18\times m}\\[5pt] \dot{X}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \dot{\boldsymbol{x}}({t}_{1}) {} & {} \dot{\boldsymbol{x}}({t}_{2}) & {} {} \cdots & {} {} \dot{\boldsymbol{x}}({t}_{m})\end{array}\right]\in {\mathbb{R}}^{18\times m},\\[5pt] U=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \boldsymbol{u}({t}_{1}) & {} {} \boldsymbol{u}({t}_{2}) {} & {} \cdots {} & {} \boldsymbol{u}({t}_{m})\end{array}\right]\in {\mathbb{R}}^{2\times m}\end{array}\right. \end{align}

1. 2. Determining key parameters for parafoil path planning

   According to the experience of engineering applications [Reference Yakimenko45], for the path planning of the parafoil system, only a few nonlinear terms are crucial, such as the state variables $ {x}_{d}$ , $ {y}_{d}$ , $ {z}_{d}$ , $ {\psi }_{p}$ and the left and right trailing edge deflections $ {\delta }_{l}$ and $ {\delta }_{r}$ . While roll dynamics significantly affect turning performance, the yaw state serves as an adequate descriptor for autonomous lateral path planning under steady-flight conditions. This simplification aligns with the 4 DOF modeling conventions for large parafoils and prioritises real-time computational efficiency.

   Library terms are pruned based on their physical interpretability in parafoil dynamics. For example, cross-terms like $ {\delta }_{l}\cdot {\psi }_{p}$ have no direct correspondence in force/moment equations, while their effects are implicitly encapsulated by state variables $ {\psi }_{p}$ . A few other terms can be regarded as zero terms. Therefore, the function $ f(\boldsymbol{x},\boldsymbol{u})$ on the right-hand side of Equation (12) can be sparsely represented in a high dimensional nonlinear function space.

1. 3. Truncation of historical data

   The 1st, 2nd, 3rd and 8th rows of the corresponding historical data matrices $ X$ and $ \dot{X}$ are extracted to form new data matrices $ \bar{X}$ and $ \dot{\bar{X}}$ . When all data under the time series $ \{{t}_{1},{t}_{2},\cdots ,{t}_{m}\}$ are expanded, the matrices can be obtained in Equations (14) to (16).

   (14) \begin{align} \bar{X}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{x}_{d}({t}_{1}) {}& {x}_{d}({t}_{2}) {} {}& \cdots {} {}& {x}_{d}({t}_{m})\\[3pt] {y}_{d}({t}_{1}) {} {}& {y}_{d}({t}_{2}) {} {}& \cdots {}& {} {y}_{d}({t}_{m})\\[3pt] {z}_{d}({t}_{1}) {} {}& {z}_{d}({t}_{2}) {} {}& \cdots {}& {} {z}_{d}({t}_{m})\\[3pt] {\psi }_{p}({t}_{1}) {}& {} {\psi }_{p}({t}_{2}) {} {}& \cdots {}& {} {\psi }_{p}({t}_{m})\end{array}\right]\!, \end{align}
   (15) \begin{align} \dot{\bar{X}}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{\dot{x}}_{}({t}_{1}) {} {}& {\dot{x}}_{d}({t}_{2}) {} {}& \cdots {} {}& {\dot{}}_{d}({t}_{m})\\[3pt] {\dot{y}}_{d}({t}_{1}) {} {}& {\dot{y}}_{d}({t}_{2}) {} {}& \cdots {}& {} {\dot{y}}_{d}({t}_{m})\\[3pt] {\dot{z}}_{d}({t}_{1}) {} {}& {\dot{z}}_{d}({t}_{2}) {} {}& \cdots {}& {} {\dot{z}}_{d}({t}_{m})\\[3pt] {\dot{\psi }}_{p}({t}_{1}) {}& {} {\dot{\psi }}_{p}({t}_{2}) {} {}& \cdots {}& {} {\dot{\psi }}_{p}({t}_{m})\end{array}\right]\!, \end{align}
   (16) \begin{align} U=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{\delta }_{l}({t}_{1}) {} {}& {\delta }_{l}({t}_{2}) {} {}& \cdots {} {}& {\delta }_{l}({t}_{m})\\[3pt] {\delta }_{r}({t}_{1}) {} {}& {\delta }_{r}({t}_{2}) {} {}& \cdots {} {}& {\delta }_{r}({t}_{m})\end{array}\right]\!. \end{align}

2. 4. Library of candidate nonlinear terms

   Combined with the sparsity analysis, by introducing the sparse coefficient $ {\Xi }$ and library of candidate nonlinear terms $ {\Theta }$ , the dynamic equation of the parafoil can be expressed as Equations 17,

   (17) \begin{align} \dot{\bar{X}}={\Xi }{{\Theta }}^{T}(\bar{X},U), \end{align}
   where $ {\Theta }$ is composed of simple combinations of certain items from the state variable data matrix $ \bar{X}$ and the control input data matrix $ U$ .

   In this paper, only closely related combination terms are set for the parafoil system to obtain the library of candidate nonlinear terms. When the terms in the library are expanded according to the time series $ \{{t}_{1},{t}_{2},\cdots ,{t}_{m}\}$ , it can be expressed as Equation (18):

   (18) \begin{align} {{\Theta }}^{T}(\bar{X},U)=\left[\begin{array}{c@{\quad}c@{\quad}c}1 {} {}& \cdots {}& {} 1\\[3pt] {x}_{d}({t}_{1}) {} {}& \cdots {} {}& {x}_{d}({t}_{m})\\[3pt] {y}_{d}({t}_{1}) {} {}& \cdots {}& {} {y}_{d}({t}_{m})\\[3pt] {z}_{d}({t}_{1}) {}& {} \cdots {}& {} {z}_{d}({t}_{m})\\[3pt] {\psi }_{p}({t}_{1}) {} {}& \cdots {}& {} {\psi }_{p}({t}_{m})\\[3pt] \mathrm{c}\mathrm{o}\mathrm{s}({\psi }_{p})({t}_{1}) {}& {} \cdots {}& {} \mathrm{c}\mathrm{o}\mathrm{s}({\psi }_{p})({t}_{m})\\[3pt] \mathrm{s}\mathrm{i}\mathrm{n}({\psi }_{p})({t}_{1}) {} {}& \cdots {}& {} \mathrm{s}\mathrm{i}\mathrm{n}({\psi }_{p})({t}_{m})\\[3pt] {\delta }_{l}({t}_{1}) {} {}& \cdots {} {}& {\delta }_{l}({t}_{m})\\[3pt] {\delta }_{r}({t}_{1}) {}& {} \cdots {} {}& {\delta }_{r}({t}_{m})\end{array}\right]\!, \end{align}

   Based on the above conditions, the corresponding sparse coefficient matrix is a 4 × 9 matrix, and its representation form is as Equation (19):

   (19) \begin{align} {\Xi }=[{\xi }_{1} \, {\xi }_{2} \, {\xi }_{3} \, {\xi }_{4}{]}^{\mathrm{T}}\in {\mathbb{R}}^{4\times 9}. \end{align}

4.2 Solution of the sparse coefficient matrix

4.2.1 Sparse coefficient estimation via LASSO-regularised regression

Based on the historical data and the library of candidate nonlinear terms, the sparse coefficient matrix $ {\Xi }$ can be obtained via sparse regression. Some coefficients in the matrix $ {\Xi }$ are shrunk to zero to reduce the model complexity, thereby achieving the purpose of feature selection.

The basic idea of Least Absolute Shrinkage and Selection Operator (LASSO) regression is adopted to approximate and solve the sparse coefficients, as shown in Equation (20). LASSO is a linear regression algorithm used for feature selection and model regularisation. It constrains the coefficients in the model by applying the LASSO penalty and shrinks some coefficients to zero. The LASSO regression equation is used to identify the sparse matrix with the best model performance in the library of candidate nonlinear terms,

(20) \begin{align} {\xi }_{k}=\underset{{\xi }_{k}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\left\|{\dot{\bar{X}}}_{k}-{\xi }_{k}{{\Theta }}^{T}(\bar{X},U)\right\|}_{2}^{2}+\lambda {\left\|{\xi }_{k}\right\|}_{1}, \end{align}

where $ {\dot{\bar{X}}}_{k}$ represents the $ k$ -th row of $ \dot{\bar{X}}$ ( $ k= 1,2,3,4$ ), $ {\xi }_{k}$ represents the $ k$ -th row of $ {\Xi }$ and $ \lambda $ is the weight term for controlling the sparsity adjustment.

Combining the summaries in Section 4.1, the basic principle of the sparse identification of the parafoil system can be obtained, as shown in Fig. 4. Here, the solid dots represent the coefficients corresponding to the key terms selected by the sparse identification, and the hollow dots represent the zero terms.

Figure 4.

Principles of sparse identification for parafoil systems.

4.2.2. Expanding historical data quantity

In this section, a method for solving the sparse matrix of the parafoil with data expansion is proposed to improve the accuracy of the sparse identification model and expand the amount of historical data.

Randomly set the initial values of the state variables and control inputs for $ P$ groups (training sample sizes). Perform simulations using the high-fidelity 9 DOF dynamic model established in Section 3. After data organisation, $ P$ groups of data matrices $ \bar{X}$ , $ \dot{\bar{X}}$ , and $ U$ are obtained.

Then the generalised model shown in Equation (17) can be expanded and expressed as Equation (21). The expression of the block matrix can be systematically derived by drawing on the definitions of the library of candidate nonlinear terms and the sparse coefficient matrix, as shown in Equation (22).

(21) \begin{align} \left\{\begin{array}{c}{\dot{\bar{X}}}^{(1)}={\Xi }{{\Theta }}^{T}({\bar{X}}^{(1)},{U}^{(1)})\\ {\dot{\bar{X}}}^{(2)}={\Xi }{{\Theta }}^{T}({\bar{X}}^{(2)},{U}^{(2)}),\\ \cdots \\ {\dot{\bar{X}}}^{(P)}={\Xi }{{\Theta }}^{T}({\bar{X}}^{(P)},{U}^{(P)})\end{array}\right. \end{align}

(22) \begin{align} \left[\begin{array}{c}{\dot{\bar{X}}}^{(1)}\\ {\dot{\bar{X}}}^{(2)}\\ \cdots \\ {\dot{\bar{X}}}^{(P)}\end{array}\right]=\boldsymbol{\Xi }\left[\begin{array}{c}{{\Theta }}^{T}({\bar{X}}^{(1)},{U}^{(1)})\\ {{\Theta }}^{T}({\bar{X}}^{(2)},{U}^{(2)})\\ \cdots \\ {{\Theta }}^{T}({\bar{X}}^{(P)},{U}^{(P)})\end{array}\right]\!, \end{align}

The left-hand side of the equation is arranged into the matrix $ {\dot{\bar{X}}}_{E}$ , and the right-hand side is arranged into the product of the sparse coefficient matrix and the matrix $ ({{\Theta }}_{E}{)}^{\mathrm{T}}$ , thereby obtaining the expression of the sparse coefficient matrix of the parafoil with data expansion, which corresponds to Equation (23).

(23) \begin{align} {\Xi }={\dot{\bar{X}}}_{E}{\left(({{\Theta }}_{E}{)}^{T}\right)}^{+}, \end{align}

the symbol “+” denotes the pseudo-inverse of a matrix.

By integrating the method elaborated in Section 4.2.1, the LASSO regression is employed to derive the sparse coefficient matrix $ {\Xi }$ . Then, an algorithm for solving the sparse coefficient matrix of the parafoil with data expansion is formulated. The specific details are presented in Algorithm 1.

Algorithm 1:

Parafoil sparse coefficient matrix solution algorithm

This paper employs 10-fold cross-validation to search for the optimal regularisation parameter $ \lambda $ in logarithmic space: 100 logarithmically spaced candidate values are generated within $ [1{0}^{-6},10]$ . The dataset is randomly partitioned into 10 subsets. For each candidate $ \lambda $ , the LASSO model is trained on nine subsets, and the prediction Mean Squared Error (MSE) is computed on the remaining subset. The $ \lambda $ that minimises the average MSE is selected. This process is optimised independently per state dimension to adapt to varying sparsity requirements of the dynamics. A post-processing step truncates sparse regression coefficients with absolute values below $ 1{0}^{-6}$ to zero, eliminating numerical noise-induced perturbations and enhancing model stability and interpretability.

Table 1.

The main parameters of parafoil system

Table 2.

Initial simulation condition of parafoil system

5.1 Influence of training sample sizes

The mechanism by which the amount of data affects the accuracy of the sparse identification model is determined by taking the parafoil system shown in Table 1 as an example and conducting an error analysis. In line with the process depicted in Fig. 4 and Algorithm 1, the number of groups $ P$ of state variables and control inputs is set to 1, 10, 100 and 1,000 respectively. As a result, four approximate sparse identification models with different data volumes are obtained. Under the same simulation conditions, the differences between the high-fidelity model and the four sparse identification models are compared. The initial simulation conditions are set according to the experiments, as detailed in Table 2, and the control variables are shown in Fig. 5.

Figure 5.

Trailing edge deflection on the left/right of parafoil.

It is worth noting that the parameters of the parafoil system used for training the model and validation remain consistent. However, to ensure generalisation ability, we implemented a strict data separation strategy: for training data, randomly generated initial states and control inputs are used; for testing, completely independent initial states (Table 2) are adopted, the control inputs are smooth (Fig. 5), and the samples that do not completely overlap with the training set are included.

Eventually, the comparison curves for the horizontal and vertical displacements, as well as the yaw angle, are acquired and shown in Figs 6 and 7. Analysis shows that:

1. (1) As the simulation time progresses, the error in the horizontal displacement gradually accumulates, manifesting a tendency of fluctuating growth. With the continuous increase in the groups $ P$ of the initial random state variables, the corresponding RMSE of the horizontal displacement steadily diminishes. When $ P$ reaches 1,000, the RMSE attains its minimum value.

2. (2) The error in the vertical displacement demonstrates a linearly ascending trend as the simulation time lengthens. Moreover, subsequent to an increase in $ P$ , the RMSE gradually drops.

3. (3) The error in the yaw angle decreases concomitantly with the increase in the groups $ P$ , and the RMSE also gradually decreases.

Figure 6.

Error analysis of displacement.

Figure 7.

Error analysis of yaw angle.

In summary, a preliminary conclusion can be derived: the initial data volume exerts a significant influence on the accuracy of the sparse identification model. Specifically, the greater the data volume, the higher the accuracy. Nevertheless, it is evident that an augmentation in the data volume will inevitably lead to an increase in computational time consumption. Consequently, in practical applications, it is necessary to carefully consider and strike a balance between the model accuracy and the available computational resources.

In the subsequent analysis, the number of groups $ P$ of random state variables is set to 1,000, and the sparse coefficient matrix of the key parameters of the parafoil system is obtained in Equation (24),

(24) \begin{align} {{\Xi }}^{T}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{\xi }_{1} {}& {\xi }_{2} {}& {\xi }_{3} {}& {\xi }_{4}\end{array}\right]=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}-0.141 & {} -0.060 {} & 10.858 {} & 0\\[3pt] 0 & {} 0 {} & 0 {} & 0\\[3pt] 0 {} & 0 {} & 0 & {} 0\\[3pt] 0 {} & 0 {} & 0 {} & 0\\[3pt] 0 & {} 0 & {} 0 & {} 0\\[3pt] 28.05 {} & -0.149 {} & -0.014 & {} 0\\[3pt] 0.0560 & {} 28.05 & {} 0.015 & {} 0\\[3pt] 0 & {} 0 & {} 0 {} & 0.14\\[3pt] 0.081 & {} -0.068 & {} -0.014 & {} -0.14\end{array}\right]\!. \end{align}

Combined with Equation (24), the sparse identification model of the parafoil system is constructed. This model is completely data-driven, and its parameters are directly obtained through model reduction. Its form is as Equation (25):

(25) \begin{align} \left\{\begin{array}{c@{\quad}c}{\dot{x}}_{d} {} ={\Theta }(\bar{X},U){\xi }_{1}\\[3pt] {\dot{y}}_{d} {} ={\Theta }(\bar{X},U){\xi }_{2}\\[3pt] {\dot{z}}_{d} {} ={\Theta }(\bar{X},U){\xi }_{3}\\[3pt] {\dot{\psi }}_{p} {} ={\Theta }(\bar{X},U){\xi }_{4}\end{array}\right.\!\!\!\!\!\!. \end{align}

5.2 Model performance comparison

In the real implementation of path planning, the model selection process is multifaceted. It requires a comprehensive consideration of the model’s prediction accuracy for the parafoil’s path and its computational efficiency. Additionally, the running time of each model is a crucial determinant. This holistic approach empowers decision makers to make more informed choices when selecting the most suitable model for specific path planning. By doing so, the goal is to ensure that the selected model attains optimal performance, striking an effective balance between high-precision prediction accuracy and efficient computational resource utilisation. The primary objective of this section is to evaluate the single simulation running time and position prediction capabilities of three distinct models: the high-fidelity model, the sparse identification model and the 4 DOF dynamic model (which is prevalently employed in path planning for the parafoil system).

5.2.1 Computational efficiency evaluation

Firstly, 1,000 computer simulations were conducted to compare the simulation times for the three models, as presented in Fig. 8. The analysis reveals that the running time of the high-fidelity model is approximately 0.053 s, while those of the sparse identification model and the 4 DOF dynamic model are around 0.0015 s and about 0.0014 s. The performance disparities are remarkable. Evidently, the running time of the high-fidelity model is substantially longer than that of the other two models. In contrast, the computational times of the sparse identification model and the 4 DOF dynamic model are comparable and relatively short. It implies that, due to its high complexity and large computational load, the high-fidelity model is ill-suited for the rapid design demands of parafoil path planning. Conversely, the other two models exhibit distinct advantages in computational efficiency and are more appropriate for path planning scenarios with high real-time requirements, which is consistent with the findings of previous research.

Figure 8.

The average running time of 1,000 simulations for different models.

5.2.2 Computational accuracy evaluation

Based on the flight test data of a certain type of parafoil system, a detailed comparison was made between the simulated values and the experimental values. A total of four airdrop tests were carried out for this type of parafoil system, and the third typical flight data was selected for model validation. The reason is that the wind field uncertainty of this group of tests is relatively low, which is suitable for verifying the model accuracy. In the test, all sensors were calibrated by the factory and time – synchronised through a 4Hz hardware trigger synchronisation mechanism. After using the Kalman filter to remove noise and outliers, the filtered data was resampled at 10Hz to calculate the RMSE.

The computational accuracy of different models was further evaluated. Figure 9(a) shows the contrast curve of the horizontal position trajectory of the parafoil system, Fig. 9(b) shows the contrast curve of the system height varying with time and Fig. 9(c) shows the contrast curve of the parafoil yaw angle varying with time. Overall, the simulation results of the models are generally in line with the trend of the test values, demonstrating a high degree of coincidence.

The models were further analysed quantitatively. Specifically, the statistical results of the mean error (ME) and RMSE were derived for the sparse recognition model and the 4 DOF dynamic model, which is most prevalently utilised in path planning, in comparison with the test values. The results are presented in Table 3.

The ME provides an intuitive indication of the overall deviation, both in terms of direction and magnitude, between the model’s predicted values and the actual values. It serves as a preliminary assessment of the model’s accuracy. As can be seen from the results, regardless of whether it is the error in horizontal position, altitude or yaw angle, the sparse identification model outperforms the 4 DOF dynamic model.

The RMSE comprehensively considers the magnitude and distribution of errors, reflecting the model’s precision and stability, and is an effective indicator for comparing model performance. In terms of horizontal displacement, height and yaw angle, the RMSE sparse recognition model is significantly better than the 4 DOF model, which indicates that the accuracy of measurement or calculation is improved.

Based on the analysis of running time, it can be concluded that the high-fidelity model is not suitable for rapid trajectory planning design. The proposed sparse identification model significantly outperforms the conventional 4 DOF dynamic model in predictive accuracy, with remarkable RMSE reductions of 12.96% (horizontal position), 54.44% (height) and 37.96% (yaw angle) in trajectory prediction. Combining enhanced accuracy and model simplicity, it is particularly suitable for path planning applications requiring both computational efficiency and high precision.

Table 3.

Error comparison and reduction rates

Figure 9.

Dynamic model validation: Experimental Data vs. Model Prediction.

5.3 Parafoil path planning using sparse identification model

To validate the practical application of the sparse identification model in the path planning task, an optimal control algorithm for the parafoil system is designed. Subsequently, fixed-point path planning of the parafoil is executed under diverse initial conditions.

5.3.1 Path planning description

The path planning problem of the parafoil system essentially aims to find a path from the deployment point to the target point, subject to specific constraint conditions and achieving specific performance indicators. (The path planning problem for the parafoil system essentially involves determining a route from the deployment point to the target point, while adhering to specific constraints and achieving desired performance metrics.) It can be formulated as the following optimisation problem:

Given the initial and terminal moments $ {t}_{0}$ and $ {t}_{f}$ , the objective is to find the manipulation control function $ u(t)$ that minimises the following performance function:

(26) \begin{align} J=F({\boldsymbol{x}}_{0},{t}_{0};\,{\boldsymbol{x}}_{f},{t}_{f}). \end{align}

Meanwhile, the parafoil dynamic equations and the terminal constraints for landing must be satisfied.

(27) \begin{align} \dot{\boldsymbol{x}}(t)=f (\boldsymbol{x}(t),\boldsymbol{u}(t)),\mathrm{}t\in [{t}_{0},{t}_{f}], \\[-30pt] \nonumber \end{align}

(28) \begin{align} \boldsymbol{G}({\boldsymbol{x}}_{0},{t}_{0};\,{\boldsymbol{x}}_{f},{t}_{f})=0. \\[-15pt] \nonumber \end{align}

For programming and calculation convenience, the time interval is normalised. Let

(29) \begin{align} \tau =(t-{t}_{0})/({t}_{f}-{t}_{0}). \end{align}

Here, $ \tau \in \left[\mathrm{0,\ 1}\right]$ .

After the transformation, Equation (26–29) can be converted into

(30) \begin{align} J=F({\boldsymbol{x}}_{0},0;\,{\boldsymbol{x}}_{f},1), \\[-30pt] \nonumber \end{align}

(31) \begin{align} \boldsymbol{x}(\tau)={\Delta }tf (\boldsymbol{x}(\tau ),\boldsymbol{u}(\tau )),\tau \in [\mathrm{0,1}], \\[-30pt] \nonumber \end{align}

(32) \begin{align} \boldsymbol{G}({\boldsymbol{x}}_{0},0;\,{\boldsymbol{x}}_{f},1)=0, \end{align}

the discrete time step is defined as $ {\Delta }t={t}_{f}-{t}_{0}$ .

In the path planning algorithm presented in this paper, the control variables are used to construct the performance index function, that is

(33) \begin{align} J={\int }_{0}^{{t}_{f}}\,\left({\left|{\delta }_{s}\right|}^{2}\right)dt, \end{align}

where $ {t}_{f}$ represents the landing time of the parafoil, $ {\delta }_{s}$ is the symmetric trailing edge deflections of the parafoil, $ {\delta }_{s}=0.5({\delta }_{r}+{\delta }_{l})$ .

In this paper, a discretisation method based on Lagrangian polynomial interpolation is adopted to solve the path planning problem, and the classical fourth-order Runge-Kutta algorithm is applied for the integral operation [Reference Conway46, Reference Hager47]. For the described nonlinear programming problem, the Sequential Quadratic Programming (SQP) algorithm can be used for the solution. Since this method is quite commonly used, its detailed description will not be repeated here.

5.3.2 Path planning result

A path planning simulation is conducted for the sparse identification model. The parameters of the parafoil are shown in Table 1. The initial conditions for the simulation are set as follows:

1. (1) Target points: Assume that its position is the zero point in the northeast-up coordinate system, that is, (0, 0, 0) m.

2. (2) Deployment points: The airdrops are carried out at altitudes of 2,000 meters, 3,000 meters and 4,000 meters, respectively, and the position coordinates are taken as (1,200, 1,000, 2,000) m, (2,200, 1,500, 3,000) m and (3,200, 2,000, 4,000) m, respectively.

Figure 10 shows the three-dimensional trajectories obtained with the optimal control method under different airdrop release conditions; Fig. 11 shows the optimal yaw angles.

Figure 10.

Comparative analysis of 3D flight path under varying deployment points.

Figure 11.

Optimal yaw angle mapping for different planned path.

The simulation results show that at the three altitudes, the landing time of the system is about 166.4 s, 249.7 s and 332.9 s, respectively. When the deployment altitude is low, the system does not have enough time to use up the remaining altitude, resulting in a relatively large yaw angle. This behaviour aligns with the actual operation of the parafoil system, confirming that the sparse identification model can accurately describe the movement characteristics of the parafoil. In addition, even though the differences in airdrop release conditions, the approximate model can still be effectively used to plan the flight path to the desired target point according to the requirements. This indicates that the approximate model is well-suited for parafoil path planning and achieves the expected results.