https://doi.org/10.5194/wes-7-2039-2022
https://doi.org/10.5194/wes-7-2039-2022 # Flight trajectory optimization of Fly-Gen airborne wind energy systems through a harmonic balance method

Filippo Trevisi, Iván Castro-Fernández, Gregorio Pasquinelli, Carlo Emanuele Dionigi Riboldi, and Alessandro Croce
Abstract

The optimal control problem for flight trajectories of Fly-Gen airborne wind energy systems (AWESs) is a crucial research topic for the field, as suboptimal paths can lead to a drastic reduction in power production. One of the novelties of the present work is the expression of the optimal control problem in the frequency domain through a harmonic balance formulation. This allows the potential reduction of the problem size by solving only for the main harmonics and allows the implicit imposition of periodicity of the solution. The trajectory is described by the Fourier coefficients of the dynamics (elevation and azimuth angles) and of the control inputs (onboard wind turbine thrust and AWES roll angle). To isolate the effects of each physical phenomenon, optimal trajectories are presented with an increasing level of physical representation from the most idealized case: (i) if the mean thrust power (mechanical power linked to the dynamics) is considered as the objective function, optimal trajectories are characterized by a constant AWES velocity over the loop and a circular shape. This is done by converting all the gravitational potential energy into electrical energy. At low wind speed, onboard wind turbines are then used as propellers in the ascendant part of the loop; (ii) if the mean shaft power (mechanical power after momentum losses) is the objective function, a part of the potential energy is converted into kinetic and the rest into electrical energy. Therefore, the AWES velocity fluctuates over the loop; (iii) if the mean electrical power is considered as the objective function, the onboard wind turbines are never used as propellers because of the power conversion efficiency. Optimal trajectories for case (ii) and (iii) have a circular shape squashed along the vertical direction. The optimal control inputs can be generally modeled with one harmonic for the onboard wind turbine thrust and two for AWES roll angle without a significant loss of power, demonstrating that the absence of high-frequency control is not detrimental to the power generated by Fly-Gen AWESs.

Share
Dates
1 Introduction

Airborne wind energy (AWE) is the branch of wind energy which aims at harvesting energy from the wind using airborne systems. Airborne wind energy systems can be classified according to the flight operations, which are linked to the power generation technique. The flight operations can be divided into crosswind, tether-aligned and rotational, as discussed by . Electrical power can be generated by a fixed or a moving ground station, or, alternatively, it can be directly generated on board and transmitted to the ground through the tether. The wing type, soft or fixed, additionally classifies the AWES. This paper focuses on AWESs based on a fixed-wing with onboard generation, known as Fly-Gen AWESs. However, the methods developed can be applied to other AWE architectures, after an appropriate rework of the dynamic models. These methods are suitable for investigating the optimal trajectories of AWESs and, especially when applied to low-fidelity models, for understanding their physical characteristics. The interpretation of the physical characteristics of optimal trajectories and the analysis of how they are influenced by parameters describing the system and its operation is the main goal of this work. With this aim, solutions are compared with analytical results coming from first-principle models whenever possible.

The first analytical power equation of crosswind AWESs was derived by Loyd (1980), and additional refinements, such as the one proposed by , made an effort to modify analytic equations to include gravitational and centrifugal effects. This kind of analytical models can be used to study how power and other relevant trends approximately scale with design parameters. However, they typically neglect the system dynamics and its effect on power generation. These effects can be studied with dynamical models, ranging from low to high fidelity. Typically, low- to mid-fidelity models are used to investigate optimal trajectories of AWE. Low-fidelity dynamic models are characterized by multiple assumptions, which simplify the models, and by the low computational cost. The quasi-steady model assumes the kite as a point mass in steady state for each point of the loop. This model is validated with experimental data , and it is considered accurate for soft kites, where the inertia is low and the AWES quickly reaches the steady state. A similar approach is considered while deriving the Unicycle model . Also this model, based on a point mass, is developed for soft-wing AWESs and computes the velocity vector via quasi-steady flight equations. The kite orientation is found by a turning law that is derived from lateral force equilibrium and is validated through a number of experiments. The Unifoil model is derived by modification of the Unicycle model in order to be applied to fixed-wing AWESs. Indeed, the quasi-steady assumption is removed, and the turning maneuvers are modeled with a yaw dynamic equation.

Higher-fidelity, but still computationally efficient, dynamic models are developed by as a part of the Lagrangian Kite Flight Simulators (LAKSA) package based on minimal coordinates and by to study the dynamics of multiple AWES configurations. Moreover, thorough Newtonian dynamic models are used to compute reference flight paths and the consequent flight path control for soft-wing AWESs and for fixed-wing AWESs .

The dynamic models just introduced are particularly suitable to be used within optimal control studies for their computational inexpensiveness and for the reduced number of nonlinearities compared to even higher-fidelity codes, such as kiteFAST . The Unicycle and Unifoil models, introduced earlier, are mainly used to compute reference flight paths and for flight path control development . To ease the deployment of optimal control problems for AWE, awebox (awebox2022) is developed and used, for instance by , and , to solve optimal control problems. awebox solves optimal control problems in time, imposing periodicity constraints. A similar optimal control problem is studied by and , where the optimal trajectory is found in time using a discretization by direct collocation and a homotopy strategy based on the relaxation of the dynamic constraints . solve an optimal control problem with an experimentally validated dynamic model of a Ground-Gen AWES. They find that, under some prescribed constraints, circular and figure-of-eight trajectories produce similar mean power and that closed-loop control enhances robustness but decreases power production by about 10 %. Control in all operation phases is studied by and : the present work can be understood as a study of the guidance (or the reference trajectory) used during the power generation phase of their study.

investigates the power losses in a circular trajectory with a dynamical quasi-analytical model. He finds that the causes of power losses are mainly two: the AWES span non-perpendicularity with respect to the incoming wind during the motion and the AWES speed fluctuation over the loop. The Makani team (Tucker2020) studies the flight trajectories of Fly-Gen AWESs with a simplified quasi-analytical approach, aiming at describing their physical characteristics. They run their flight simulator for different trajectories and production strategies to derive analytical expressions, which can describe the consequences of different operational choices. Their production strategy at low wind speed is to convert part of the potential energy into kinetic and part into electrical, when the AWES moves downward. To reduce the potential energy exchange, they suggest to squash the trajectories along the vertical direction. Moreover, they explain that using electrical power to push the AWES upward drastically decreases the overall power production, as power needs to be converted from mechanical to electrical and again from electrical to mechanical, so that the related efficiencies are counted twice. They, in accord with the study for Ground-Gen by , conclude that the electrical conversion losses should be considered when deciding on the production strategy. Following these conclusions, the present work also investigates the influence of the power generation efficiencies on the optimal trajectories.

As the aim of this work is to interpret optimal trajectories in a physical way, a low-fidelity dynamic model, similar to the one proposed by (reformulated for Fly-Gen AWESs), is selected. Instead of solving the dynamics and the optimal control problem in time, the present approach models the problem in the frequency domain, making use of a harmonic balance method, which expands the periodic solution as a Fourier series . Working with the Fourier coefficients and not with the time series themselves allows the potential reduction of the problem size depending on the problem at hand, the search for periodic solutions implicitly and the study of the solution in an intuitive way by looking at the contribution of the different harmonics. To the best of the authors' knowledge, this is the first work on AWESs where an optimal control problem aided by a harmonic balance methodology is formulated.

Even though the frequency-domain formulation can be used for any periodic flight trajectories (i.e., circular and figure of eight), only circular trajectories are analyzed here to limit the paper scope and length. Figure-of-eight trajectories are intended to be analyzed and intensively compared with circular trajectories in a future work.

The paper is organized as follows: in Sect. 2 the flight dynamic model, the harmonic balance and the optimal control statement are introduced. In Sect. 3, the main results from steady-state analytical models are recalled from literature, together with the introduction of some key non-dimensional numbers used later in the analyses. In Sect. 4, the solution obtained with the harmonic balance formulation is validated against the time integration. In Sect. 5, optimal control problems with constant wind inflow and no constraints on the mean elevation angle are analyzed. This extreme idealization allows the understanding of some optimal trajectory characteristics which are also present in more realistic cases. Section 6 focuses on the results of a more realistic optimal control problem where the wind shear and a constraint on the minimum elevation angle are included in the analyses. Finally, in Sect. 7 the results are discussed and the main conclusions summarized.

2 Methodology

## 2.1 Flight dynamic model

Two coordinate systems (Fig. 1a) are defined to derive the equations of motion. The ground coordinate system (denoted by G) is inertial and centered at the ground station: ex points downwind, ez toward the zenith and ey completes the right-handed frame. For convenience, spherical coordinates are used to describe the position of the airborne unit, with Lt the tether length, ϕ the azimuth angle and β the elevation angle. The spherical reference frame (denoted by S) is unequivocally defined at every position with the origin at the AWES center of mass, er pointing outward the sphere in the radial direction, eϕ normal to er and contained on a plane parallel to xy, and ${\mathbit{e}}_{\mathit{\beta }}={\mathbit{e}}_{\mathrm{r}}×{\mathbit{e}}_{\mathit{\varphi }}$. The position p, velocity v and acceleration a vectors projected into the spherical reference frame S are

$\begin{array}{}\text{(1)}& \begin{array}{rl}& \mathbit{p}={L}_{\mathrm{t}}{\mathbit{e}}_{\mathrm{r}},\\ & \mathbit{v}={L}_{\mathrm{t}}\stackrel{\mathrm{˙}}{\mathit{\varphi }}\mathrm{cos}\mathit{\beta }{\mathbit{e}}_{\mathit{\varphi }}+{L}_{\mathrm{t}}\stackrel{\mathrm{˙}}{\mathit{\beta }}{\mathbit{e}}_{\mathit{\beta }},\\ & \mathbit{a}=\left(-{L}_{\mathrm{t}}{\stackrel{\mathrm{˙}}{\mathit{\varphi }}}^{\mathrm{2}}{\mathrm{cos}}^{\mathrm{2}}\mathit{\beta }-{L}_{\mathrm{t}}{\stackrel{\mathrm{˙}}{\mathit{\beta }}}^{\mathrm{2}}\right){\mathbit{e}}_{\mathrm{r}}\\ & +\left({L}_{\mathrm{t}}\stackrel{\mathrm{¨}}{\mathit{\varphi }}\mathrm{cos}\mathit{\beta }-\mathrm{2}{L}_{\mathrm{t}}\stackrel{\mathrm{˙}}{\mathit{\varphi }}\stackrel{\mathrm{˙}}{\mathit{\beta }}\mathrm{sin}\mathit{\beta }\right){\mathbit{e}}_{\mathit{\varphi }}\\ & +\left({L}_{\mathrm{t}}{\stackrel{\mathrm{˙}}{\mathit{\beta }}}^{\mathrm{2}}\mathrm{sin}\mathit{\beta }\mathrm{cos}\mathit{\beta }+{L}_{\mathrm{t}}\stackrel{\mathrm{¨}}{\mathit{\beta }}\right){\mathbit{e}}_{\mathit{\beta }}.\end{array}\end{array}$ Figure 1(a) Ground reference frame G (${\mathbit{e}}_{x}-{\mathbit{e}}_{y}-{\mathbit{e}}_{z}$) and spherical reference frame S (${\mathbit{e}}_{\mathrm{r}}-{\mathbit{e}}_{\mathit{\varphi }}-{\mathbit{e}}_{\mathit{\beta }}$) and (b) sketch for the spanwise unit vector s definition.

The wind velocity is in the positive x-axis direction of G, and projected into the spherical reference frame is

$\begin{array}{}\text{(2)}& \begin{array}{rl}& {\mathbit{v}}_{\mathrm{w}}={v}_{\mathrm{w}}\left(\mathrm{cos}\mathit{\varphi }\mathrm{cos}\mathit{\beta }{\mathbit{e}}_{\mathrm{r}}-\mathrm{sin}\mathit{\varphi }{\mathbit{e}}_{\mathit{\varphi }}-\mathrm{cos}\mathit{\varphi }\mathrm{sin}\mathit{\beta }{\mathbit{e}}_{\mathit{\beta }}\right),\\ & {v}_{\mathrm{w}}\left(h\right)={v}_{\mathrm{w},\mathrm{0}}{\left(\frac{h}{{h}_{\mathrm{0}}}\right)}^{{\mathit{\alpha }}_{\mathrm{s}}}={v}_{\mathrm{w},\mathrm{0}}{\left(\frac{{L}_{\mathrm{t}}}{{h}_{\mathrm{0}}}\mathrm{sin}\mathit{\beta }\right)}^{{\mathit{\alpha }}_{\mathrm{s}}},\end{array}\end{array}$

where the wind speed vw as function of the altitude h is modeled with an exponential law: vw,0 is the reference wind speed at the reference altitude h0 and αs is the wind shear exponent. The relative speed between the AWES and the wind is

$\begin{array}{}\text{(3)}& {\mathbit{v}}_{\mathrm{r}}={\mathbit{v}}_{\mathrm{w}}-\mathbit{v}.\end{array}$

To describe the AWES attitude, a non-sideslip velocity constraint is included in the modeling. Indeed, the wing operates at the highest performance under this condition. To impose this constraint implicitly, the unit vector e1 is defined to point to the direction opposite to the relative wind speed:

$\begin{array}{}\text{(4)}& {\mathbit{e}}_{\mathrm{1}}=-\frac{{\mathbit{v}}_{\mathrm{r}}}{|{\mathbit{v}}_{\mathrm{r}}|}.\end{array}$

The spanwise unit vector s (with origin at the center of mass and pointing in the right-wing span direction) is defined perpendicular to e1 with the procedure illustrated in Fig. 1b. A second vector e3 is defined as a unit vector in a plane parallel to the xz plane with elevation βs (and negative sign):

$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\mathbit{e}}_{\mathrm{3}}=-\left({\mathbit{e}}_{x}\mathrm{cos}\left({\mathit{\beta }}_{\mathrm{s}}\right)+{\mathbit{e}}_{z}\mathrm{sin}\left({\mathit{\beta }}_{\mathrm{s}}\right)\right),\\ & {\mathbit{e}}_{x}=\mathrm{cos}\mathit{\varphi }\mathrm{cos}\mathit{\beta }{\mathbit{e}}_{\mathrm{r}}-\mathrm{sin}\mathit{\varphi }{\mathbit{e}}_{\mathit{\varphi }}-\mathrm{cos}\mathit{\varphi }\mathrm{sin}\mathit{\beta }{\mathbit{e}}_{\mathit{\beta }},\\ & {\mathbit{e}}_{z}=\mathrm{sin}\mathit{\beta }{\mathbit{e}}_{\mathrm{r}}+\mathrm{cos}\mathit{\beta }{\mathbit{e}}_{\mathit{\beta }}.\end{array}\end{array}$

Note that e3 points upwind when βs=0. The unit vector e2 is then defined as

$\begin{array}{}\text{(6)}& {\mathbit{e}}_{\mathrm{2}}=\frac{{\mathbit{e}}_{\mathrm{3}}×{\mathbit{e}}_{\mathrm{1}}}{|{\mathbit{e}}_{\mathrm{3}}×{\mathbit{e}}_{\mathrm{1}}|},\end{array}$

where $|{\mathbit{e}}_{\mathrm{3}}×{\mathbit{e}}_{\mathrm{1}}|$ can take values smaller than one because e3 and e1 are not defined to be perpendicular in general. In this way, e2 is perpendicular to the plane e3e1. Rodrigues' formula is then used to define s through a rotation of ψ around e1, starting from e2:

$\begin{array}{}\text{(7)}& \mathbit{s}={\mathbit{e}}_{\mathrm{2}}\mathrm{cos}\mathit{\psi }+\left({\mathbit{e}}_{\mathrm{1}}×{\mathbit{e}}_{\mathrm{2}}\right)\mathrm{sin}\mathit{\psi }+{\mathbit{e}}_{\mathrm{1}}\left({\mathbit{e}}_{\mathrm{1}}\cdot {\mathbit{e}}_{\mathrm{2}}\right)\left(\mathrm{1}-\mathrm{cos}\mathit{\psi }\right).\end{array}$

With this formulation, s is defined to be always perpendicular to the relative wind, and its components are defined by a unique angle ψ, called hereafter roll angle. When ψ=0, s is perpendicular to e3.

The aerodynamic lift L and the drag D take the standard form

$\begin{array}{}\text{(8)}& \mathbit{L}=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A{C}_{\mathrm{L}}\left|{\mathbit{v}}_{\mathrm{r}}\right|{\mathbit{v}}_{\mathrm{r}}×\mathbit{s},\phantom{\rule{0.25em}{0ex}}\mathbit{D}=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A{C}_{\mathrm{D}}\left|{\mathbit{v}}_{\mathrm{r}}\right|{\mathbit{v}}_{\mathrm{r}},\end{array}$

where ρ is the air density, A is the wing area, and the lift and drag coefficients CL and CD are considered constant. The drag coefficient CD includes the contribution from the tether drag . The gravitational force Fg and the tether force T are

$\begin{array}{}\text{(9)}& {\mathbit{F}}_{\mathrm{g}}=-mg\left(\mathrm{sin}\mathit{\beta }{\mathbit{e}}_{\mathrm{r}}+\mathrm{cos}\mathit{\beta }{\mathbit{e}}_{\mathit{\beta }}\right),\phantom{\rule{0.25em}{0ex}}\mathbit{T}=-T{\mathbit{e}}_{\mathrm{r}},\end{array}$

where m is the AWES mass, g is the gravitational acceleration and T is the norm of the tether force. The thrust produced by the onboard wind turbines Dt is expressed as a linear function of the aerodynamic drag with gain γ:

$\begin{array}{}\text{(10)}& {\mathbit{D}}_{\mathrm{t}}=\mathit{\gamma }\mathbit{D}.\end{array}$

The dynamic equations of motion in compact form read

$\begin{array}{}\text{(11)}& m\mathbit{a}=\mathbit{L}+\mathbit{D}+{\mathbit{D}}_{\mathrm{t}}+{\mathbit{F}}_{\mathrm{g}}+\mathbit{T},\end{array}$

recalling that a is given by Eq. (1).

As the objectives of the optimal control problems are linked to the power production, three different power quantities are defined. The thrust power Pt (i.e., the power linked to the AWES dynamics) is estimated as a dot product of Dt and the relative velocity:

$\begin{array}{}\text{(12)}& {P}_{\mathrm{t}}={\mathbit{D}}_{\mathrm{t}}\cdot {\mathbit{v}}_{\mathrm{r}}.\end{array}$

The shaft power Ps (i.e., the mechanical power that can be converted to electrical power) is modeled using 1D momentum theory (actuator disc) as

$\begin{array}{}\text{(13)}& {P}_{\mathrm{s}}=\left(\mathrm{1}-a\right){P}_{\mathrm{t}}=\left(\frac{\mathrm{1}}{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\mathrm{1}-\mathit{\gamma }{C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}}\right){P}_{\mathrm{t}},\end{array}$

where the induction a is found by setting the thrust given by momentum theory (${\mathbit{T}}_{\mathrm{d}}=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{A}_{\mathrm{t}}\left(\mathrm{4}a\left(\mathrm{1}-a\right)\right)\left|{\mathbit{v}}_{\mathrm{r}}\right|{\mathbit{v}}_{\mathrm{r}}$) equal to Dt, as in , and At is the total turbine area.

Finally, the electrical power exchanged with the grid P takes into account the generator and transmission efficiency ηel:

When power is generated (γ>0), the electrical power distributed to the grid P is lower than the shaft power Ps because of electrical efficiencies. When power from the grid is used, the electrical power requested to the grid P is instead higher in absolute value compared to the shaft power Ps. To model the discontinuity in a continuous optimization framework, the logistic function is used:

$\begin{array}{}\text{(15)}& P={P}_{\mathrm{s}}-\left(\frac{\mathrm{1}-{e}^{-f\mathit{\gamma }}}{\mathrm{1}+{e}^{-f\mathit{\gamma }}}\right)\left(\mathrm{1}-{\mathit{\eta }}_{\mathrm{el}}\right){P}_{\mathrm{s}},\end{array}$

where f is taken equal to 100.

## 2.2 Frequency-domain formulation

Frequency-domain formulations may present clear advantages when solving for periodic solutions of dynamic and control problems. They have the capability of solving for both stable and unstable branches of periodic solutions in an efficient way. Moreover, they potentially use fewer variables to describe the same problems. Since the problem of optimal trajectories for AWESs has a periodic nature, the flight dynamic model just introduced is expressed in the frequency domain. The harmonic balance methodology is then used to transform the differential equations of motion into a set of nonlinear algebraic equations . The equations of motion (Eq. 11) can be written as a set of second-order nonlinear differential equations in the form

$\begin{array}{}\text{(16)}& \begin{array}{rl}& \mathbf{f}\left(\mathbit{x},\stackrel{\mathrm{˙}}{\mathbit{x}},\stackrel{\mathrm{¨}}{\mathbit{x}},\mathbit{u}\right)=\mathbf{0},\\ & \mathbit{x}={\left[\begin{array}{c}\mathit{\beta }\left(t\right),\phantom{\rule{0.25em}{0ex}}\mathit{\varphi }\left(t\right)\end{array}\right]}^{T},\phantom{\rule{0.25em}{0ex}}\mathbit{u}={\left[\begin{array}{c}\mathit{\psi }\left(t\right),\phantom{\rule{0.25em}{0ex}}\mathit{\gamma }\left(t\right)\end{array}\right]}^{T},\end{array}\end{array}$

where x is the state vector and u is the control vector. By assuming that Eq. (16) accepts periodic solutions, every variable of the state vector is expanded as a Fourier series of order Nx:

$\begin{array}{}\text{(17)}& \begin{array}{rl}& x\left(t\right)\approx \frac{{X}_{\mathrm{0}}}{\mathrm{2}}+\sum _{k=\mathrm{1}}^{{N}_{x}}{X}_{k,\mathrm{s}}\mathrm{sin}\left(k\mathit{\omega }t\right)+{X}_{k,\mathrm{c}}\mathrm{cos}\left(k\mathit{\omega }t\right),\\ & \mathbf{X}={\left[\begin{array}{c}{X}_{\mathrm{0}},\phantom{\rule{0.25em}{0ex}}{X}_{\mathrm{1},\mathrm{s}},\phantom{\rule{0.25em}{0ex}}{X}_{\mathrm{2},\mathrm{s}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots }\phantom{\rule{0.125em}{0ex}}{X}_{\mathrm{1},\mathrm{c}},\phantom{\rule{0.25em}{0ex}}{X}_{\mathrm{2},\mathrm{c}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots }\end{array}\right]}^{T},\end{array}\end{array}$

with $\mathit{\omega }=\frac{\mathrm{2}\mathit{\pi }}{\mathcal{T}}$ being the fundamental frequency of the motion and 𝒯 the period. Alternatively, the state vector can be expressed as

$\begin{array}{}\text{(18)}& \begin{array}{rl}& x\left(t\right)\approx \frac{{A}_{\mathrm{0}}}{\mathrm{2}}+\sum _{k=\mathrm{1}}^{{N}_{x}}{A}_{k}\mathrm{cos}\left(k\mathit{\omega }t-{\mathit{\theta }}_{k}\right),\\ & \mathbf{A}={\left[\begin{array}{c}{A}_{\mathrm{0}},\phantom{\rule{0.25em}{0ex}}{A}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots }\end{array}\right]}^{T},\phantom{\rule{0.25em}{0ex}}\mathbit{\theta }={\left[\begin{array}{c}{\mathit{\theta }}_{\mathrm{1}},\phantom{\rule{0.25em}{0ex}}{\mathit{\theta }}_{\mathrm{2}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots }\end{array}\right]}^{T},\end{array}\end{array}$

where

$\begin{array}{}\text{(19)}& {A}_{\mathrm{0}}={X}_{\mathrm{0}},\phantom{\rule{0.25em}{0ex}}{A}_{k}=\sqrt{{X}_{k,\mathrm{s}}^{\mathrm{2}}+{X}_{k,\mathrm{c}}^{\mathrm{2}}},\phantom{\rule{0.25em}{0ex}}{\mathit{\theta }}_{k}=\mathrm{arctan}\left(\frac{{X}_{k,\mathrm{s}}}{{X}_{k,\mathrm{c}}}\right).\end{array}$

The first and second time derivatives of the state vector can be found analytically:

$\begin{array}{}\text{(20)}& \begin{array}{rl}& \stackrel{\mathrm{˙}}{x}\left(t\right)\approx \sum _{k=\mathrm{1}}^{{N}_{x}}k\mathit{\omega }\left({X}_{k,\mathrm{s}}\mathrm{cos}\left(k\mathit{\omega }t\right)-{X}_{k,\mathrm{c}}\mathrm{sin}\left(k\mathit{\omega }t\right)\right),\\ & \stackrel{\mathrm{¨}}{x}\left(t\right)\approx -\sum _{k=\mathrm{1}}^{{N}_{x}}\left(k\mathit{\omega }{\right)}^{\mathrm{2}}\left({X}_{k,\mathrm{s}}\mathrm{sin}\left(k\mathit{\omega }t\right)+{X}_{k,\mathrm{c}}\mathrm{cos}\left(k\mathit{\omega }t\right)\right).\end{array}\end{array}$

Similarly, the control inputs, assumed to be periodic, can also be expressed as a Fourier series of order Nu:

$\begin{array}{}\text{(21)}& \begin{array}{rl}& u\left(t\right)\approx \frac{{U}_{\mathrm{0}}}{\mathrm{2}}+\sum _{k=\mathrm{1}}^{{N}_{u}}{U}_{k,\mathrm{s}}\mathrm{sin}\left(k\mathit{\omega }t\right)+{U}_{k,\mathrm{c}}\mathrm{cos}\left(k\mathit{\omega }t\right),\\ & \mathbf{U}={\left[\begin{array}{c}{U}_{\mathrm{0}},\phantom{\rule{0.25em}{0ex}}{U}_{\mathrm{1},\mathrm{s}},\phantom{\rule{0.25em}{0ex}}{U}_{\mathrm{2},\mathrm{s}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots }\phantom{\rule{0.125em}{0ex}}{U}_{\mathrm{1},\mathrm{c}},\phantom{\rule{0.25em}{0ex}}{U}_{\mathrm{2},\mathrm{c}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots }\end{array}\right]}^{T},\end{array}\end{array}$

where Nu<Nx because the equations of motion need to be solved at frequencies higher than the control inputs order. By introducing Eqs. (17), (20) and (21) into Eq. (16), the equations of motion can be expanded as a Fourier series of order Nx:

$\begin{array}{}\text{(22)}& \begin{array}{rl}\mathbf{f}\left({\mathbf{X}}_{\mathit{\beta }},{\mathbf{X}}_{\mathit{\varphi }},{\mathbf{U}}_{\mathit{\psi }},{\mathbf{U}}_{\mathit{\gamma }},\mathit{\omega },t\right)& \approx \frac{{\mathbf{F}}_{\mathrm{0}}}{\mathrm{2}}+\sum _{k=\mathrm{1}}^{{N}_{x}}{\mathbf{F}}_{k,\mathrm{s}}\mathrm{sin}\left(k\mathit{\omega }t\right)\\ & +{\mathbf{F}}_{k,\mathrm{c}}\mathrm{cos}\left(k\mathit{\omega }t\right)=\mathbf{0}.\end{array}\end{array}$

The Fourier coefficients of the equations of motion are found numerically by applying the Fourier coefficient definition to the time series, which should have a minimum size of 2Nx+1. A Galerkin methodology is then applied by pre-multiplying Eq. (22) by 1, sin (kωt), and cos (kωt) and subsequently integrating the resulting equation over one period. The result is a set of $\mathrm{2}×\left(\mathrm{2}{N}_{x}+\mathrm{1}\right)$ nonlinear algebraic equations as a consequence of the orthogonality properties of the selected basis of trigonometric functions:

$\begin{array}{}\text{(23)}& \mathbf{F}=\left[\begin{array}{c}{\mathbf{F}}_{\mathrm{0}};\phantom{\rule{0.25em}{0ex}}{\mathbf{F}}_{\mathrm{1},\mathrm{s}};\phantom{\rule{0.25em}{0ex}}{\mathbf{F}}_{\mathrm{2},\mathrm{s}};\phantom{\rule{0.125em}{0ex}}\mathrm{\dots }\phantom{\rule{0.125em}{0ex}}{\mathbf{F}}_{\mathrm{1},\mathrm{c}};\phantom{\rule{0.25em}{0ex}}{\mathbf{F}}_{\mathrm{2},\mathrm{c}};\phantom{\rule{0.125em}{0ex}}\mathrm{\dots }\end{array}\right]=\mathbf{0},\end{array}$

which can be understood as the residuals of the equations of motion expressed in the frequency domain. For given periodic control inputs and a given fundamental frequency, the periodic solution can be found by looking for the Fourier coefficients [Xβ; Xϕ] of the dynamics which solve Eq. (23).

## 2.3 Optimal control problem (OCP)

In this work, the frequency-domain formulation is included within an optimal control problem (OCP). A generic optimization problem can be written as

$\begin{array}{}\text{(24)}& \begin{array}{rl}& {\mathcal{X}}^{*}=\mathrm{arg}\left(\underset{\mathcal{X}}{\mathrm{min}}\phantom{\rule{0.25em}{0ex}}\mathrm{obj}\left(\mathcal{X}\right)\right),\\ & \mathrm{s}.\mathrm{t}.:\phantom{\rule{0.25em}{0ex}}\mathbf{lb}\le \mathcal{X}\le \mathbf{ub}\\ & \mathbf{g}\left(\mathcal{X}\right)\le \mathrm{0}\\ & \mathbf{h}\left(\mathcal{X}\right)=\mathrm{0},\end{array}\end{array}$

where 𝒳 represents the unknown optimization variables, 𝒳* their optimal values, “obj” the objective function, lb and ub the lower and upper bounds of 𝒳, g the inequality, and h the equality constraints. In the present formulation, the optimization variables are the Fourier coefficients of the state variables, of the control inputs and the fundamental frequency:

$\begin{array}{}\text{(25)}& \mathcal{X}=\left[{\mathbf{X}}_{\mathit{\beta }};{\mathbf{X}}_{\mathit{\varphi }};{\mathbf{U}}_{\mathit{\psi }};{\mathbf{U}}_{\mathit{\gamma }};\mathit{\omega }\right].\end{array}$

The negative value of the mean thrust power ${\stackrel{\mathrm{^}}{P}}_{\mathrm{t}}$ (Eq. 12), shaft power ${\stackrel{\mathrm{^}}{P}}_{\mathrm{s}}$ (Eq. 13) or electric power $\stackrel{\mathrm{^}}{P}$ (Eq. 15) over the loop is taken as objective function, where the symbol $\stackrel{\mathrm{^}}{.}$ stands for the mean value over the loop. The equality constraints are the aggregation of the residuals of the equation of motion in the frequency domain F (Eq. 23) and additional physical constraints R in the frequency domain (e.g., certain quantities can be imposed to be constant over the loop):

$\begin{array}{}\text{(26)}& \mathbf{h}\left(\mathcal{X}\right)=\left[\mathbf{F}\left(\mathcal{X}\right);\mathbf{R}\left(\mathcal{X}\right)\right]=\mathbf{0}.\end{array}$

Inequality constraints g, expressed in the time domain, can also be included in the problem (e.g., the minimum elevation angle over the loop can be bounded). A graphical representation of the OCP setup is given in Fig. 2. The derivatives of flight dynamic model with respect to the optimization variables can be taken analytically and provided to the solver, allowing for a deep and fast convergence of the solution. The OCP is implemented in a MATLAB® environment and solved with the interior-point algorithm implemented in fmincon. As the chosen optimization algorithm (gradient-based) can only look for local optima, the initial guess may influence the solution. In this work, the initial guesses are taken to be circular trajectories, leading to circular-shaped optimal trajectories. Figure-of-eight trajectories can be implemented as initial guesses, which may lead to figure-of-eight-shaped optimal trajectories. A detailed comparison between these two trajectory types is left for future works. The OCPs are solved on an Intel Core i7-9700 3.0 GHz, 16 GB RAM, system. The computation times of the presented examples require from a few seconds to tens of seconds. For example, OCP A in Sect. 6 takes approximately 25 s, and OCP B takes approximately 12 s.

To compare the results of the optimal control problem with idealized analytical expressions, the main results from a refined version of the Loyd power equation (Loyd1980) are here briefly recalled. The thrust power equation, with the assumption of linear crosswind motion , is

$\begin{array}{}\text{(27)}& \begin{array}{rl}{P}_{\mathrm{t},\mathrm{L}}& =\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A{v}_{\mathrm{w}}^{\mathrm{3}}\mathit{\gamma }\frac{{C}_{\mathrm{L}}}{G}{\left(\mathrm{1}+{\left(\frac{G}{\mathrm{1}+\mathit{\gamma }}\right)}^{\mathrm{2}}\right)}^{\mathrm{3}/\mathrm{2}}\\ & =\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A{v}_{\mathrm{w}}^{\mathrm{3}}\mathit{\gamma }\frac{{C}_{\mathrm{L}}}{G}{\left(\mathrm{1}+{G}_{\mathrm{t}}^{\mathrm{2}}\right)}^{\mathrm{3}/\mathrm{2}},\end{array}\end{array}$

where the system glide ratio (including tether drag) is $G=\frac{{C}_{\mathrm{L}}}{{C}_{\mathrm{D}}}$, and, for readability, a modified glide ratio is defined as ${G}_{\mathrm{t}}=\frac{{C}_{\mathrm{L}}}{{C}_{\mathrm{D}}\left(\mathrm{1}+\mathit{\gamma }\right)}$ by including the drag of the onboard propellers. The shaft power takes into account the onboard wind turbine induction a:

$\begin{array}{}\text{(28)}& {P}_{\mathrm{s},\mathrm{L}}=\left(\mathrm{1}-a\right){P}_{\mathrm{t},\mathrm{L}}=\left(\frac{\mathrm{1}}{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\mathrm{1}-\mathit{\gamma }{C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}}\right){P}_{\mathrm{t},\mathrm{L}}.\end{array}$

Finally, the power generated and sent to the grid takes into account the efficiencies of the electrical conversion:

$\begin{array}{}\text{(29)}& {P}_{\mathrm{L}}={\mathit{\eta }}_{\mathrm{el}}{P}_{\mathrm{s},\mathrm{L}}.\end{array}$

For high G, the power equation simplifies to

$\begin{array}{}\text{(30)}& {P}_{\mathrm{L}}\approx \frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A{v}_{\mathrm{w}}^{\mathrm{3}}{C}_{\mathrm{L}}{G}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{el}}\left(\frac{\mathrm{1}}{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\mathrm{1}-\mathit{\gamma }{C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}}\right)\frac{\mathit{\gamma }}{\left(\mathrm{1}+\mathit{\gamma }{\right)}^{\mathrm{3}}}.\end{array}$

For this expression, the value of γ which maximizes the power is only a function of the non-dimensional quantity ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$. In Fig. 3, the electrical power PL, normalized with the maximum electrical power at ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}=\mathrm{0}$, is plotted as a function of γ and ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$. For increasing values of ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$, the value of γ which maximizes power production decreases. The maximum normalized power as a function of ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$ is shown in Fig. 4, highlighting that the analytical expression predicts a decrease in power production for increasing ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$. Figure 3Normalized power $\frac{{P}_{\mathrm{L}}}{{P}_{\mathrm{L}}\left({C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}=\mathrm{0}\right)}$ as a function of ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$ for high glide ratios G. Figure 4Maximum normalized power $\frac{{P}_{\mathrm{L}}}{{P}_{\mathrm{L}}\left({C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}=\mathrm{0}\right)}$ as a function of ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$ for high glide ratios G.

Table 1Reference values for the examples (values from the Makani MX2 description; ), associated non-dimensional parameters and quantities evaluated with the steady-state model for γ maximizing Eq. (29). The tether force can be evaluated as

$\begin{array}{}\text{(31)}& {T}_{\mathrm{L}}=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A{v}_{\mathrm{w}}^{\mathrm{2}}\frac{{C}_{\mathrm{L}}}{{G}_{\mathrm{t}}}{\left(\mathrm{1}+{G}_{\mathrm{t}}^{\mathrm{2}}\right)}^{\mathrm{3}/\mathrm{2}}.\end{array}$

showed that for high G, neglecting gravity and with constant incoming wind, there exists an opening angle $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ (angle swept by the AWES during the circular trajectory; see Fig. 5) which erases the power losses due to centrifugal forces and that it is only a function of the non-dimensional mass parameter:

$\begin{array}{}\text{(32)}& M=\frac{m}{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A{C}_{\mathrm{L}}{L}_{\mathrm{t}}}.\end{array}$

In this idealized case, the turning radius is $R={L}_{\mathrm{t}}\mathrm{sin}\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ and the revolution period is

$\begin{array}{}\text{(33)}& {\mathcal{T}}_{\mathrm{L}}=\frac{\mathrm{2}\mathit{\pi }R}{{v}_{\mathrm{L}}}=\frac{\mathrm{2}\mathit{\pi }R}{{v}_{\mathrm{w}}{G}_{\mathrm{t}}},\end{array}$

where vL is norm of the AWES velocity.

In addition to the non-dimensional mass parameter, the Froude number, which weights the fluid inertial forces to gravity forces, is used in this work:

$\begin{array}{}\text{(34)}& Fr=\sqrt{\frac{{v}_{\mathrm{w}}^{\mathrm{2}}}{g\cdot {L}_{\mathrm{t}}}},\end{array}$

where the reference velocity is the wind velocity and the reference length is the tether length. By combining the previously introduced non-dimensional parameters, the gravity ratio Gr is defined as

$\begin{array}{}\text{(35)}& {G}_{\mathrm{r}}=\frac{M}{F{r}^{\mathrm{2}}{G}_{\mathrm{t}}^{\mathrm{2}}}=\frac{mg}{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A{C}_{\mathrm{L}}{v}_{\mathrm{w}}^{\mathrm{2}}{G}_{\mathrm{t}}^{\mathrm{2}}},\end{array}$

which represents the ratio between gravitational force and aerodynamic lift, similarly to the one introduced in .

In the following sections, the results will be generalized as a function of the non-dimensional parameters just introduced. Input parameters from the Makani MX2 design (Tucker2020) will be used as reference values to present the results (Table 1).

4 Validation of the frequency-domain formulation against time integration

To make sure the frequency-domain formulation is well implemented and finds solutions which respect the equations of motion, they are compared with the solution coming from a time integration scheme. The model described in Sect. 2.1 is solved with the MATLAB® ode45 integration scheme. After solving the periodic solution with the harmonic balance methodology, the Fourier coefficients of the state and control vector are retrieved. The state vector at t=0 is used as an initial condition for the numerical integration. The control inputs must be computed from their Fourier series at every step of the integration. In Appendix B, a comparison for a circular and a figure-of-eight trajectory is shown. The solution of the dynamics, represented by the azimuth and elevation, for the two cases is equivalent, demonstrating that the frequency-domain formulation is accurate enough to be used in the present optimal control problem framework.

5 Optimal control problems with constant inflow and no elevation angle constraints

As the analysis is limited to circular trajectories, a cylindrical reference frame C, similar to the one employed in and , is used to present the results. A graphical representation of C is given in Fig. 5. The longitudinal axis of C is aligned with the mean elevation angle $\stackrel{\mathrm{^}}{\mathit{\beta }}$. The angle βm denotes the minimum elevation angle and Φ the opening angle. The angular position of the AWES is defined by α, and when α=0 the AWES moves upward (i.e., $\stackrel{\mathrm{˙}}{\mathit{\alpha }}>\mathrm{0}$). Figure 5Cylindrical reference system C used to analyze circular trajectories.

To increase complexity incrementally, the optimal control problems (OCPs) are modified from the most idealized case to a realistic one. For the idealized cases analyzed in this section, uniform incoming wind speed (αs=0) and no minimum elevation angle constraints are considered. In this section, βs (Eq. 5) is set equal to zero, such that e3 points upwind. In this way, when the roll is equal to zero (ψ=0), the span direction is perpendicular to the incoming wind.

## 5.1 Optimizing for the mean electrical power in the absence of gravity

For the most idealized case, the gravity is null g=0, such that Fr→∞ and Gr=0. The objective function is taken as the mean electrical power, given in Eq. (15). By solving the OCP for the example (Table 1), it is found that the solution has constant values over the trajectory and the average power output is equal to the one evaluated with the analytical expression in Eq. (29). Figure 6 shows the evolution of β and ϕ, highlighting that the solution is a circle. Due to the constant values of the solution, quantities such as tether force along the axial symmetry axis, γ, AWES velocity and others can be found with the formulation assuming a crosswind straight motion, as Sect. 3.

For the solution to be optimal, it is found that the AWES span is perpendicular to the wind speed, or, in analytical terms, that ψ=0. Figure 7 shows the optimal opening angle Φ* as a function of a modified non-dimensional mass parameter Mt found by solving a number of OCPs with different G (G∈[10 30]), M (M∈[0.025 0.15]) and ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$ (${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}\in \left[\mathrm{0}\phantom{\rule{0.25em}{0ex}}\mathrm{0.4}\right]$). Figure 7Optimal opening angles Φ* (x) found by solving multiple OCPs and analytical expression (–) as a function of the modified non-dimensional mass parameter Mt.

The values of Φ* can be accurately described by

$\begin{array}{}\text{(36)}& \stackrel{\mathrm{̃}}{\mathrm{\Phi }}=\mathrm{arccos}\left(-\frac{{M}_{\mathrm{t}}}{\mathrm{2}}+\frac{\sqrt{{M}_{\mathrm{t}}^{\mathrm{2}}+\mathrm{4}}}{\mathrm{2}}\right),\phantom{\rule{0.25em}{0ex}}{M}_{\mathrm{t}}=\frac{M}{\mathrm{1}+\frac{\mathrm{1}}{{G}_{\mathrm{t}}^{\mathrm{2}}}},\end{array}$

which for high glide ratios coincides with the analytical formulation given in .

Table 2Settings of the two optimal control problems maximizing the mean thrust power considering gravity. V1,s and V1,c are the first Fourier coefficients of the norm of the AWES velocity v. ## 5.2 Optimizing for the mean thrust power considering gravity

Gravity is now included in the modeling, and the objective function is taken as the mean thrust power ${\stackrel{\mathrm{^}}{P}}_{\mathrm{t}}$ (Eq. 12). The results of two slightly different OCPs are shown for the sake of understanding the results, and they are summarized in Table 2. In OCP A, the control inputs are modeled with five harmonics. In OCP B, the time series of the control input ψ is modeled as a constant, and only one harmonic is used for the control input γ. Additionally, the norm of the AWES velocity $v=|\mathbit{v}|$ is imposed to be constant (additional equality constraint). As the control inputs act up to the first harmonic, this constraint is set by imposing the first Fourier coefficients of the AWES velocity to zero ($\mathbf{R}=\left[{V}_{\mathrm{1},\mathrm{s}};{V}_{\mathrm{1},\mathrm{c}}\right]$), while no constraints are imposed on the higher-order harmonics. In Table 2, the mean thrust power (objective function) is also reported and compared with the analytical formulation (Eq. 27). The objective function of the two OCPs is almost the same, showing that the two problems are basically equivalent.

By solving the OCPs, it is found that the optimal solutions have a negative mean elevation of ${\stackrel{\mathrm{^}}{\mathit{\beta }}}_{A}\approx -\mathrm{8.2}$ and ${\stackrel{\mathrm{^}}{\mathit{\beta }}}_{B}\approx -\mathrm{7.8}$. The trajectories, shown in Fig. 8, have a circular shape (a circle with radius $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ is marked as $-\phantom{\rule{0.25em}{0ex}}.$), but it is no longer a perfect circle. Figure 8Optimal trajectory for OCP A (), B ($-\phantom{\rule{0.25em}{0ex}}-$) and a circle with radius $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ ($-\phantom{\rule{0.25em}{0ex}}.$).

Figure 9 shows the trends of the control input ψ as a function of α (see Fig. 5 for definition). For OCP A, it fluctuates with small amplitude about the mean value, which is close to zero. For this reason, it is modeled as a constant in OCP B. The optimal constant value is also close to zero, meaning that the AWES span is perpendicular to the wind speed direction. Since the two OCPs present similar optimal values of power, it is found that the optimal solutions are not sensitive to the fluctuations of the roll angle ψ.

Figure 10 shows the evolution of γ as a function of α. The mean values for the two OCPs are close to the value maximizing Eq. (27), denoted in figure as γL. γ takes values higher than the mean in the descending leg of the loop and negative values in the ascendant leg. This means that in the ascendant leg the onboard wind turbines are operated as propellers. In OCP A γ is modeled with five harmonics, while in OCP B it is modeled with just one. As the objective function of the two optimal solutions is similar, it is found that the trend of γ can be modeled with just one harmonic.

Figure 11 shows the norm of the optimal AWES velocity v for the two OCPs. In OCP B, this value is constrained to be a constant by imposing only its first harmonic null. Since it is not possible to impose constraints at frequencies where the control (γ and ψ) is not acting, higher-order harmonics are not constrained. The mean values of v are similar the one predicted by the steady model. Since the fluctuations for OCP A are small compared to the mean, the influence of the velocity fluctuations on the overall performances is investigated in OCP B, showing that they impact weakly the optimal solution. As the two OCPs are basically equivalent, it is found that optimal trajectories are characterized by a AWES constant velocity. Figure 11Norm of the optimal AWES velocity v as a function of the angular position.

Figure 12 shows the magnitude of the tether force for the two cases. As it scales with the relative wind speed squared, also the tether force has an almost constant trend (the fluctuations are small compared to the mean).

To compare the two OCPs and draw some conclusions, the power output, shown in Fig. 13, is to be analyzed. The mean thrust power output for the two OCPs is slightly higher than Pt,L. This is due to a nonlinear effect induced by the combination of gravity and mean elevation angle different from zero as compared to the idealized case. Indeed, for negative mean elevation, the combination leads to an increase of mean thrust power, while the opposite occurs in case of positive mean elevation. As the effects on power is almost negligible and it does not primarily impact the main outcomes of this paper, a detailed explanation of this phenomenon is avoided here, but the reader can find more details in . The theoretical thrust power, given in Eq. (27), is derived neglecting gravity. However, it approximates well the power output obtained through the OCP, which includes gravity. Figure 13Optimal thrust power production and consumption Pt as a function of the angular position.

As the two analyzed OCPs are almost equivalent, the optimal trajectories are characterized by the perpendicularity of the AWES span with respect to the wind (ψ=0) and a constant AWES velocity. In order to keep the AWES velocity constant over the loop, the onboard wind turbines balance the action of the gravitational force. In the descendent leg, the onboard wind turbines harvest the gravitational potential energy, and in the ascendant leg, that power is given back to the system.

Following these considerations, the power trend, as shown in Fig. 13, can then be approximated as

$\begin{array}{}\text{(37)}& {P}_{\mathrm{t}}\left(\mathit{\alpha }\right)\approx {P}_{\mathrm{t},\mathrm{L}}+{\mathbit{F}}_{\mathrm{g}}\cdot \mathbit{v}\approx {P}_{\mathrm{t},\mathrm{L}}-mg{v}_{\mathrm{w}}{G}_{\mathrm{t}}\mathrm{cos}\mathit{\alpha }.\end{array}$

The onboard wind turbine thrust can be approximated with ${\mathbit{D}}_{\mathrm{t}}\approx \left(\stackrel{\mathrm{^}}{\mathit{\gamma }}+{A}_{\mathit{\gamma },\mathrm{1}}\mathrm{cos}\left(\mathit{\alpha }-{\mathit{\theta }}_{\mathit{\gamma },\mathrm{1}}\right)\right)\mathbit{D}$, where D is constant because the AWES velocity is found to be constant, and from Fig. 10 it is found that ${\mathit{\theta }}_{\mathit{\gamma },\mathrm{1}}\approx \mathrm{180}$. As the thrust power can be written as the product of Dt and the relative wind speed vr (Eq. 12), the amplitude of the first Fourier coefficient of γ, considering Eq. (37), can be approximated by

$\begin{array}{}\text{(38)}& {A}_{\mathit{\gamma },\mathrm{1}}\approx \frac{{\mathbit{F}}_{\mathrm{g}}\cdot \mathbit{v}}{-\mathrm{cos}\mathit{\alpha }\mathbit{D}\cdot {\mathbit{v}}_{\mathrm{r}}}\approx \frac{mg}{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A{C}_{\mathrm{D}}{v}_{\mathrm{w}}^{\mathrm{2}}{G}_{\mathrm{t}}^{\mathrm{2}}}={G}_{\mathrm{r}}G.\end{array}$

Figure 14 shows the comparison of Aγ,1 found numerically by running the OCP (with the settings of OCP B) for a different combination of M (M∈[0.025 0.15]), G (G∈[10 30]), and Fr (Fr∈[0.1 0.2]) and the analytical approximation given in Eq. (38). Figure 14Amplitude of the first Fourier coefficient of γ (x) found by solving multiple OCPs and analytical approximation (–) as a function of GrG.

Figure 15 shows the first Fourier coefficient of elevation β and azimuth ϕ as a function of the non-dimensional parameter Mt, as they represent the width and height of the trajectory. The analytical expression given in Eq. (36) is still a good approximation of the optimal trajectory shape. Figure 15First Fourier coefficient of elevation β and azimuth ϕ (x) found by solving multiple OCPs and analytical expression (–) as a function of the modified non-dimensional mass parameter Mt.

## 5.3 Optimizing for the mean shaft power considering gravity

In this section, the onboard wind turbine induction is included in the power evaluation, and the mean shaft power ${\stackrel{\mathrm{^}}{P}}_{\mathrm{s}}$ is considered as objective function. To present the results, two different OCPs are introduced (Table 3). The mean shaft power for OCP A and B is almost identical, highlighting that one harmonic to model the productive drag and two for the roll are enough. The power for the analytical case is found by maximizing Eq. (28) with respect to γ. The figures in this section refer to OCP B.

Table 3Settings of the two optimal control problems maximizing the mean shaft power considering gravity. Figure 16 shows the trajectory in the βϕ plane. The trajectory deviates from a circular shape, especially along the β axis, and has a mean elevation angle of ${\stackrel{\mathrm{^}}{\mathit{\beta }}}_{B}=-\mathrm{5.5}$, higher than for the case without induction. Figure 16Optimal trajectory (–) and a circle with radius $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$.

Figure 17 shows the roll angle as a function of the angular position in the loop. Even in this case with induction, the fluctuations are relatively small. When the AWES increases the turning radius (approximately between $-\mathrm{90}{}^{\circ }<\mathit{\alpha }<\mathrm{0}{}^{\circ }$ and $\mathrm{90}{}^{\circ }<\mathit{\alpha }<-\mathrm{180}{}^{\circ }$; see Fig. 16), the roll decreases. Figure 17Optimal ψ (blue line) and its mean (red line) as a function of the angular position.

Figure 18 shows γ as a function of the angular position. The mean value is smaller compared to the value maximizing Eq. (28). By comparing the trends with Fig. 10, it is clear that the fluctuations of γ are lower than in the case without induction. The onboard wind turbine induction a has a similar trend to γ, as a and γ are linked through the expression in Eq. (28). When γ takes negative values – in the ascendant leg – the onboard wind turbines are operated as propellers and the induction is negative. In the descendent leg, γ takes values larger than the mean and so does the induction. Higher values of induction result in a lower ratio between shaft power, which is the power the optimizer maximizes, and thrust power, which is the power directly linked to the dynamics. Therefore, high values of γ are not beneficial for the shaft power production.

In Fig. 19, the AWES velocity is shown, highlighting that it fluctuates over the loop. When maximizing the mean thrust power (Sect. 5.2), the optimal AWES velocity over the loop was found to be constant. Here, it is optimal to convert part of the potential energy into electrical energy and part into kinetic energy, letting the velocity fluctuate over the loop. Figure 19Norm of the optimal AWES velocity v as a function of the angular position.

To conclude the analysis of the example, Fig. 20 shows the shaft and the thrust power. As anticipated, when γ takes higher values than the mean, the induction grows and the ratio between shaft and thrust power decreases consequently. Figure 20Optimal shaft power production and consumption Ps and thrust power Pt as a function of the angular position.

In Fig. 4 the dependence of the analytical expression of the shaft power on ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$ is shown. In Fig. 21, the dependence of the optimal mean shaft power is analyzed as a function of the same non-dimensional coefficient for three different Froude numbers (i.e., three wind speeds). In the current example, Fr=0.1 corresponds to vw=5.4 m s−1, Fr=0.15 corresponds to vw=8.1 m s−1 and Fr=0.2 corresponds to vw=10.8 m s−1. For increasing Froude number, the solution gets closer to the analytical formulation because the power fluctuations gradually lose impact on the mean power production. Indeed, the aerodynamic forces become dominant with respect to the gravitational force. Figure 21Optimal shaft power production normalized with the analytical expression of thrust power as a function of ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$. Figure 22Optimal values of $\stackrel{\mathrm{^}}{\mathit{\gamma }}$ and Aγ,1, normalized with GrG as a function of ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$ for different Froude numbers. Figure 23Optimal values of Aβ,1 and Aϕ,1, normalized with the analytical expression of the opening angle $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$, as a function of ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$ for different Froude numbers.

In Sect. 5.2, it was found that the amplitude of the first Fourier coefficient of γ for ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}=\mathrm{0}$ (i.e., optimizing for the thrust power) can be approximated by ${A}_{\mathit{\gamma },\mathrm{1}}\approx {G}_{\mathrm{r}}G$. In Fig. 22, the trends of $\frac{{A}_{\mathit{\gamma },\mathrm{1}}}{{G}_{\mathrm{r}}G}$, with γ being modeled with a single harmonic, for the three investigated Froude numbers, are shown as a function of ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$. For ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}\to \mathrm{0}$, trends are close to 1. For increasing Fr, the curves collapse to a unique curve. In particular, at ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}=\mathrm{0.23}$, which is the value for the example, the ratio $\frac{{A}_{\mathit{\gamma },\mathrm{1}}}{{G}_{\mathrm{r}}G}\to \mathrm{0.54}$ for increasing Fr. For increasing values of ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$, the ratio $\frac{{A}_{\mathit{\gamma },\mathrm{1}}}{{G}_{\mathrm{r}}G}$ tends to zero, highlighting the fact that a less fluctuating value of γ over the loop is beneficial. The plot shows also the value of $\stackrel{\mathrm{^}}{\mathit{\gamma }}$ as a function of ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$, highlighting that for increasing Fr the trends collapse to the value maximizing Eq. (28), indicated as γL.

Finally, Fig. 23 shows the ratio of the first Fourier coefficient of the elevation angle β and the azimuth angle ϕ with the opening angle $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ evaluated with Eq. (36). For ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}\to \mathrm{0}$, values are close to 1, as noted in Fig. 15. As ${C}_{\mathrm{D}}\frac{A}{{A}_{\mathrm{t}}}$ increases, the values of Aβ,1 decrease more than Aϕ,1, showing that the optimal trajectory no longer has a circular shape and the height decreases more than the width. This effect is visible also in the example in Fig. 16. At low Froude numbers (i.e., low wind speeds) this effect is more evident.

## 5.4 Optimizing for the mean electrical power considering gravity

In this section, the electrical efficiency is included into the optimal control problem, and the mean electrical power is considered as objective function. Two OCPs, whose characteristics are given in Table 4, are introduced to present results. The power for OCP A and B is almost identical, highlighting that one harmonic to model the productive drag and two for the roll are enough. The power for the analytical case is found by maximizing Eq. (29) with respect to γ.

Table 4Settings of the two optimal control problems maximizing the mean electrical power considering gravity. Figure 24 shows the trajectory in the βϕ plane. The trajectories deviate from the circular trajectory with opening angle $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ (Eq. 36), especially along the β axis, and have a mean elevation angle of ${\stackrel{\mathrm{^}}{\mathit{\beta }}}_{A}=-\mathrm{5.1}$ and ${\stackrel{\mathrm{^}}{\mathit{\beta }}}_{B}=-\mathrm{4.9}$. Figure 24Optimal trajectory for OCP A (), B ($-\phantom{\rule{0.25em}{0ex}}-$) and a circle with radius $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ ($-\phantom{\rule{0.25em}{0ex}}.$).

Figure 25 shows the roll angle as a function of the angular position in the loop. As in the case maximizing mean shaft power, the fluctuations are relatively small. Figure 27Norm of the optimal AWES velocity v as a function of the angular position.

Figure 26 shows γ as a function of the angular position. In OCP A, the time evolution of γ is modeled with five harmonics. In the ascendant leg, γ takes null values, meaning that the power is neither spent nor consumed. Indeed, spending power drastically reduces the overall power production because of the conversion efficiency from electrical to thrust power. This is also highlighted by . In OCP B, the time evolution of γ is modeled with just one harmonic. The trend is, however, similar to OCP A, with the minimum value being slightly negative. This means that Aγ,1 is similar to $\stackrel{\mathrm{^}}{\mathit{\gamma }}$.

Figure 27 shows the norm of the AWES velocity over the loop, showing that the trend is similar for the two OCPs, and, as noted in Sect. 5.3, it is optimal to convert part of the potential energy into electrical energy and part into kinetic energy.

The electrical power as a function of the angular position is shown in Fig. 28. As expected when analyzing the trend of γ, the electrical power is null in the ascendant leg and larger than the mean in the descending part. In OCP B, the mean power is slightly lower than in OCP A but the power fluctuations are lower, which could be beneficial from a grid and power smoothing perspective. Due to the electrical efficiency, the wind turbines are not used as propellers anymore. This results in an even more squashed trajectory (Fig. 24) with respect to the previous Sect. 5.3 (Fig. 16) to partially limit the potential energy exchange into kinetic energy and thus the AWES speed fluctuation. Figure 28Optimal electrical power production as a function of the angular position.

One could try to investigate how the optimal values evolve for an increasing wind speed. Figure 29 shows the mean value of γ and its first Fourier coefficient Aγ,1 as a function of GrG (see Eq. 35 for definition) for OCP B. As noted when analyzing Fig. 26, at vw=6 m s−1 Aγ,1 is slightly larger than $\stackrel{\mathrm{^}}{\mathit{\gamma }}$. As wind speed increases up to approximately 8.5 m s−1, Aγ,1 keeps being similar to $\stackrel{\mathrm{^}}{\mathit{\gamma }}$, meaning that the minimum value of γ is close to zero and so is power (Pmin≈0). If the wind speed increases again, Aγ,1 becomes lower than $\stackrel{\mathrm{^}}{\mathit{\gamma }}$, meaning that power is always generated over the loop (Pmin>0). The main effect of the electrical efficiency on the OCP is to prevent the onboard wind turbines from being operated as propellers. Therefore, when the value of Aγ,1 which maximizes the shaft power Ps is larger than $\stackrel{\mathrm{^}}{\mathit{\gamma }}$, results are expected to be modified with respect to Sect. 5.3. Instead, when Aγ,1 is lower than $\stackrel{\mathrm{^}}{\mathit{\gamma }}$, trends are expected to be equal to the analyses in Sect. 5.3. Indeed, when analyzing Fig. 22, it is found that $\frac{{A}_{\mathit{\gamma },\mathrm{1}}}{{G}_{\mathrm{r}}G}\to \mathrm{0.54}$ for high Fr. This means that for low GrG, the first Fourier coefficient of γ can be approximated with ${A}_{\mathit{\gamma },\mathrm{1}}\approx \mathrm{0.54}{G}_{\mathrm{r}}G$, as shown in Fig. 29. Figure 29Optimal values of $\stackrel{\mathrm{^}}{\mathit{\gamma }}$ and Aγ,1 as a function of GrG.

In Fig. 30, the mean power normalized with the power evaluated with Eq. 29 is shown as a function of GrG for a case with ηel=1, which is equivalent to the case in Sect. 5.3, and for a case with ηel=0.8, as in this section. The two curves for wind speed lower than 8.5 m s−1 diverge. A low electrical efficiency ηel not only decreases the power output as in Eq. (29), but also decreases the efficiency with respect to the analytical approximation due to its effect on the dynamics. Figure 30Normalized electrical power for a case with ηel=0.8 () and ηel=1 ($-\phantom{\rule{0.25em}{0ex}}-$) as a function of GrG.

To conclude, Fig. 31 shows the evolution of the first Fourier coefficient of the elevation angle Aβ,1 and of the azimuth Aϕ,1 as a function of GrG. For high wind speed (i.e., low GrG), their value is similar to the approximation given in Eq. (36). As GrG increases, for ηel=1, Aϕ,1 stays almost constant, while Aβ,1 decreases. Smaller Aβ,1 means smaller vertical height, which results in lower potential energy converted into electrical and kinetic energy over the loop. For ηel=0.8 after vw=8.5 m s−1, both Fourier coefficients decrease rapidly, meaning that smaller loops are performed. Smaller loops are therefore beneficial at low wind speed as they decrease the energy fluctuations. Figure 31Optimal values of Aβ,1 and Aϕ,1 normalized with the analytical expression of the opening angle $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ for a case with ηel=0.8 (–) and ηel=1 (– –) as a function of GrG.

6 Optimal control problem considering gravity, wind shear and elevation constraint

In this section, the wind shear is included in the problem. The reference altitude is h0=100 m, the wind shear exponent is αs=0.2 and the reference wind speed is ${v}_{\mathrm{w},\mathrm{0}}=\mathrm{6}$ m s−1. To make the problem more realistic, a constraint on the minimum elevation angle of βm=10 (which is equivalent to a constraint on the minimum flight altitude) is included. The value of βs needed to compute e3 and then the spanwise unit vector s (Eq. 5) is taken as the mean elevation angle ${\mathit{\beta }}_{\mathrm{s}}=\stackrel{\mathrm{^}}{\mathit{\beta }}$. With this definition, the case of no roll (ψ=0) is obtained when the wing span is in the plane perpendicular to the mean elevation angle, as in .

Two OCPs are solved and they are summarized in Table 5. OCP A features five harmonics to model the control inputs, while OCP B has one harmonic to model the onboard wind turbine thrust and two for the roll, as in the previous sections. The two optimizations have similar mean electric power outputs, meaning that they are almost equivalent.

Table 5Settings of the two optimal control problems maximizing the mean power considering gravity and wind shear. Figure 32 shows the trajectory for OCP A and B and compares it with a circle of radius $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ centered at an elevation of $\stackrel{\mathrm{^}}{\mathit{\beta }}=\mathrm{arctan}\sqrt{{\mathit{\alpha }}_{\mathrm{s}}}$, as this formulation identifies the elevation of the center of the wind power window . As for the cases analyzed in Sect. 5.3 and 5.4, the trajectory is squashed along the vertical direction. The constraint on the minimum elevation angle is not used, as the trajectory of both cases is always strictly higher than βm=10. Figure 32Optimal trajectory for case A (), case B ($-\phantom{\rule{0.25em}{0ex}}-$) and a circle with radius Φ ($-\phantom{\rule{0.25em}{0ex}}.$).

The roll angle ψ is shown in Fig. 33. In LT-GliDe the flight stability of AWESs is studied by linearizing the equations of motion with respect to a fictitious steady-state condition, where the AWES moves in a circular trajectory with a constant velocity. This steady state is characterized by the AWES span being perpendicular to the mean elevation angle direction. In this section, this condition is identified by ψ=0. The roll fluctuations, shown in Fig. 33, are limited in amplitude and might be considered within the linear bounds of the linearization validity of LT-GliDe. More analyses to prove this will be carried out in future works.

Figure 34 shows the evolution of the onboard wind turbine thrust coefficient as a function of the angular position for the two cases. The trends are similar to the analyses in Sect. 5.4. It is optimal to use the onboard wind turbines only to generate power and not as propellers. Even if the trends of γ for OCP A and B are quite different, the overall power production is similar, meaning that power production is not sensitive to harmonics of γ higher than one.

Figure 35 shows the wind speed that the AWES encounters over the loop due to the wind shear. Clearly, at the top of the loop (α=90), the wind speed is the highest, sweeping approximately 1.5 m s−1 over the trajectory. In this section, the mean wind speed over the loop is used to evaluate the Froude number Fr (Eq. 34) and consequently the gravity ratio Gr (Eq. 35). These numbers will be used later in this section to generalize results.

Figure 36 shows the AWES velocity as a function of the angular position. As discussed in the previous sections, in the descending leg the AWES converts the potential energy into electrical energy, producing power, and kinetic energy, accelerating. In the climbing leg instead, electrical power is not spent and kinetic energy is transformed into potential energy. Figure 36Norm of the optimal AWES velocity v as function of the angular position.

Figure 37 shows the power production as a function of the angular position. When looking at power and tether force (not shown here as it follows the AWES velocity trend squared) to characterize the operations of a real system, the maximum power and the maximum and minimum tether force would be constrained not to exceed some given values. To properly include these constraints, additional control inputs, useful to model the de-powering of the AWES (e.g., the lift coefficient), shall be considered in the analysis.

As carried out in the previous section, trends are studied as a function of the Froude number for the optimal control problem B. Figure 38 shows the dependence of $\stackrel{\mathrm{^}}{\mathit{\gamma }}$ and Aγ,1 as a function of GrG. $\stackrel{\mathrm{^}}{\mathit{\gamma }}$ decreases when GrG increases (i.e., the wind speed decreases). At low wind speed, $\stackrel{\mathrm{^}}{\mathit{\gamma }}$ takes low values so that the AWES speed over the loop is higher, which is beneficial to stay airborne. Aγ,1 for low GrG has a linear trend, as noted in Sect. 5.4. When Aγ,1 is equal to $\stackrel{\mathrm{^}}{\mathit{\gamma }}$, the minimum power production over the loop is null. For lower wind speeds (i.e., higher GrG), it is no longer optimal to increase Aγ,1 because the onboard wind turbines would be used as propellers with a high penalty on the mean power production. Figure 38Optimal values of $\stackrel{\mathrm{^}}{\mathit{\gamma }}$ and Aγ,1 as a function of GrG. Figure 39Optimal values of Aβ,1 and Aϕ,1 normalized with the analytical expression of the opening angle $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ for a case with ηel=0.8 () and ηel=1 ($-\phantom{\rule{0.25em}{0ex}}-$) as a function of GrG.

To conclude, Fig. 39 shows the ratio of the first Fourier coefficient of β and ϕ with respect to the analytical expression of the opening angle $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ for a case maximizing electrical power (ηel=0.8) and shaft power (ηel=1). After the cusp in Fig. 38, the two trends diverge, and for ηel=0.8 smaller loops are optimal so that the exchange of potential energy over the loop is reduced.

7 Conclusions and discussion

In this work, a novel methodology to study optimal trajectories for Fly-Gen AWESs is introduced. The chosen low-fidelity dynamic model is characterized by 2 degrees of freedom (the AWES is modeled as a point mass with constant tether length) and two control inputs. The degrees of freedom are the elevation and the azimuth angle. The control inputs are the roll angle, defined as the rotation around the relative velocity direction, and the onboard wind turbine thrust coefficient. An optimal control problem is formulated in the frequency domain through a harmonic balance method. Working with the Fourier coefficients of the time series, instead of the time series themselves, allows the potential reduction of the problem size, the implicit imposition of periodicity and the acquisition of an intuitive understanding of the results by analyzing the harmonic contributions. Moreover, the analytical gradient of the objective function and the constraints with respect to the optimization variables can be provided to the solver, allowing for a deep and fast convergence of the optimal solutions.

The MX2 design from is taken as a reference AWES to introduce the results. To isolate the effects of each physical phenomenon, results are presented with an increasing level of complexity from the most idealized case, and they are compared with analytical solutions from literature, whenever possible. A set of idealized case studies with no constraint on the minimum elevation angle and uniform wind inflow are initially studied. If gravity is neglected, the solution is steady, and it can be described by analytical expressions. If gravity is considered, three different optimal control problems, characterized by three different objective functions, are solved.

• i.

If the mean thrust power (mechanical power neglecting onboard wind turbine induction) is the objective function, the optimal trajectories are circular, have a constant AWES velocity and the wing span is perpendicular to the incoming wind. To obtain this condition, all the potential energy is converted into electrical energy by the onboard wind turbines. At low wind speed, onboard wind turbines are then used as propellers in the ascendant part of the loop. The optimal power, the trajectory shape and the production strategy can be accurately approximated with analytical expressions.

• ii.

If the mean shaft power (mechanical power considering onboard wind turbine induction) is the objective function, the potential energy, in the descending leg, is partially converted into electrical energy and partially into kinetic energy. This is because the power conversion penalizes solutions with high onboard wind turbine induction. Therefore, the velocity fluctuates over the loop, and the trajectories are squashed along the vertical direction to decrease the potential energy exchange.

• iii.

If the mean power electrical provided to the grid is the objective function (i.e., the electrical efficiency is included), the onboard wind turbines never operate as propellers. If operated as propellers, power would be converted from mechanical into electrical while descending and from electrical into mechanical while ascending, leading to large power losses due to the electrical efficiency. This effect is found only at low wind speed, when propelling the AWES in the climbing leg maximizes the mean shaft power. Past a given wind speed, using the onboard wind turbines as propellers does not maximize the mean shaft power, and the influence of the electrical efficiency on the production strategy vanishes.

When the wind shear and a constraint on the minimum elevation angle are included in the optimal control problem for maximizing the electrical power, trends are similar to the case with uniform inflow. Therefore, the power production strategy does not heavily depend on the wind shear. For the analyzed example, the constraint on the minimum elevation angle is not active.

For all the analyzed cases, additional analytical approximations characterizing the solution are introduced. These approximations are found by modeling the control inputs with the lowest number of harmonics. The onboard wind turbine thrust can be modeled with just one harmonic and the roll with two harmonics without loss of generality of the results.

The results of this work align with the discussions in and have a strong mathematical foundation, as the trajectory and the control inputs are found by solving optimal control problems. These methods are planned to be applied, with appropriate modifications, to other AWE architectures and to other trajectory types. A comparison between circular and figure-of-eight trajectories is foreseen. Finally, the physical understanding and methods proposed here are envisaged to be incorporated into the design, analysis and optimization framework T-GliDe , with the aim of improving the power estimation and including an optimal control module.

Appendix A: Nomenclature
 Latin symbols Lt Tether length A Wing area At Total onboard wind turbine area Aγ,1 Amplitude of the first Fourier harmonic of γ Aβ,1 Amplitude of the first Fourier harmonic of β Aϕ,1 Amplitude of the first Fourier harmonic of ϕ a Onboard wind turbine induction CL Lift coefficient CD Drag coefficient T Norm of the tether force vw Norm of the wind velocity v Norm of the AWES velocity v AWES velocity vr Relative wind speed M Non-dimensional mass parameter Mt Modified non-dimensional mass parameter m AWES mass G System glide ratio
 Gt Modified system glide ratio Gr Gravity ratio g Inequality constraints g Gravitational acceleration h Equality constraints R Additional equality constraints in the frequency domain Nx Order of the Fourier series of the state variables Nγ Order of the Fourier series of the control input γ Nψ Order of the Fourier series of the control input ψ Fg Gravitational force Pt Thrust power Ps Shaft power P Electrical power Fr Froude number Greek symbols β Elevation angle ϕ Azimuth angle ψ Roll angle ρ Air density γ Onboard wind turbine factor 𝒯 Revolution period ω Revolution frequency α Angular position in the loop ηel Electrical conversion efficiency Φ Opening angle of the trajectory 𝒳 Optimization variables Symbols ⋅L Quantity evaluated with the steady state (Loyd) model $\stackrel{\mathrm{^}}{\cdot }$ Mean value
Appendix B: Figures of comparison between frequency-domain formulation and time integration Figure B1Azimuth and elevation of the trajectory found with the harmonic balance method and the time integration scheme for a circular-shaped trajectory. Figure B2Time series of the control inputs provided to the harmonic balance method and the time integration scheme for a circular-shaped trajectory. Figure B3Azimuth and elevation of the trajectory found with the harmonic balance method and the time integration scheme for a figure-of-eight-shaped trajectory. Figure B4Time series of the control inputs provided to the harmonic balance method and the time integration scheme for a figure-of-eight-shaped trajectory.

Code availability

The code was developed for this publication. For inquiries about it, please reach out to the corresponding author.

Data availability

Author contributions

FT, ICF and GP conceptualized the study and the research methods. FT developed the research methods. FT and ICF developed the code. FT produced the results and wrote the draft version of the paper. CEDR and AC supervised the research. ICF, GP, CEDR and AC reviewed the draft version.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Financial support

The work by PoliMI had no external funding and was therefore self-funded by the research team. The work by ICF was carried out under the framework of the GreenKite-2 project (PID2019-110146RB-I00) funded by MCIN/AEI/10.13039/501100011033.

Review statement

This paper was edited by Roland Schmehl and reviewed by Nikolaus Vertovec and two anonymous referees.

References

Argatov, I., Rautakorpi, P., and Silvennoinen, R.: Apparent wind load effects on the tether of a kite power generator, J. Wind Eng. Indust. Aerodynam., 99, 1079–1088, https://doi.org/10.1016/j.jweia.2011.07.010, 2011. a

awebox: Modelling and optimal control of single- and multiple-kite systems for airborne wind energy, https://github.com/awebox/awebox, last access: 27 May 2022. a

Cobb, M. K., Barton, K., Fathy, H., and Vermillion, C.: Iterative Learning-Based Path Optimization for Repetitive Path Planning, With Application to 3-D Crosswind Flight of Airborne Wind Energy Systems, IEEE T. Control Syst. Technol., 28, 1447–1459, https://doi.org/10.1109/TCST.2019.2912345, 2020. a, b

De Schutter, J., Leuthold, R., Bronnenmeyer, T., Paelinck, R., and Diehl, M.: Optimal control of stacked multi-kite systems for utility-scale airborne wind energy, in: 2019 IEEE 58th Conference on Decision and Control (CDC), 11–13 December 2019, Nice, France, 4865–4870, https://doi.org/10.1109/CDC40024.2019.9030026, 2019. a

Dimitriadis, G.: Time Integration, in: Introduction to Nonlinear Aeroelasticity, in: chap. 3, John Wiley & Sons, Ltd, 63–111, https://doi.org/10.1002/9781118756478.ch3, 2017. a, b

Eijkelhof, D. and Schmehl, R.: Six-degrees-of-freedom simulation model for future multi-megawatt airborne wind energy systems, Renew. Energy, 196, 137–150, https://doi.org/10.1016/j.renene.2022.06.094, 2022. a

Fagiano, L., Zgraggen, A. U., Morari, M., and Khammash, M.: Automatic Crosswind Flight of Tethered Wings for Airborne Wind Energy: Modeling, Control Design, and Experimental Results, IEEE T. Control Syst. Technol., 22, 1433–1447, https://doi.org/10.1109/TCST.2013.2279592, 2014. a

Fechner, U. and Schmehl, R.: Flight path control of kite power systems in a turbulent wind environment, in: 2016 American Control Conference (ACC), 6–8 July 2016, Boston, MA, USA, 4083–4088, https://doi.org/10.1109/ACC.2016.7525563, 2016. a

Fechner, U., van der Vlugt, R., Schreuder, E., and Schmehl, R.: Dynamic model of a pumping kite power system, Renew. Energy, 83, 705–716, https://doi.org/10.1016/j.renene.2015.04.028, 2015. a

Fernandes, M. C. R. M., Paiva, L. T., and Fontes, F. A. C. C.: Optimal Path and Path-Following Control in Airborne Wind Energy Systems, in: Advances in Evolutionary and Deterministic Methods for Design, Optimization and Control in Engineering and Sciences, edited by: Gaspar-Cunha, A., Periaux, J., Giannakoglou, K. C., Gauger, N. R., Quagliarella, D., and Greiner, D., Springer International Publishing, Cham, 409–421, https://doi.org/10.1007/978-3-030-57422-2_26, 2021. a, b

Gros, S. and Diehl, M.: Modeling of Airborne Wind Energy Systems in Natural Coordinates, in: Airborne Wind Energy, Green Energy and Technology, edited by: Ahrens, U., Diehl, M., and Schmehl, R., Springer, Singapore, 181–203, https://doi.org/10.1007/978-3-642-39965-7_10, 2013. a

Gros, S., Zanon, M., and Diehl, M.: A relaxation strategy for the optimization of Airborne Wind Energy systems, in: 2013 European Control Conference (ECC), 17–19 July 2013, Zurich, Switzerland, 1011–1016, https://doi.org/10.23919/ECC.2013.6669670, 2013. a

Haas, T., De Schutter, J., Diehl, M., and Meyers, J.: Wake characteristics of pumping mode airborne wind energy systems, J. Phys.: Conf. Ser., 1256, 012016, https://doi.org/10.1088/1742-6596/1256/1/012016, 2019. a

Horn, G., Gros, S., and Diehl, M.: Numerical trajectory optimization for airborne wind energy systems described by high fidelity aircraft models, in: Airborne Wind Energy, edited by: Ahrens, U., Diehl, M., and Schmehl, R., Springer, 205–218, https://doi.org/10.1007/978-3-642-39965-7_11, 2013. a

Jonkman, J., Hayman, G., Mudafort, R., Damiani, R., Wendt, F., Google Incorporated, and USDOE: KiteFAST, https://www.osti.gov//servlets/purl/1786962 (last access: 27 May 2022), 2018. a

Lau, S. L., Cheung, Y. K., and Wu, S. Y.: A variable parameter incrementation method for dynamic instability of linear and nonlinear elastic systems, J. Appl. Mech. Trans. ASME, 49, 849–853, https://doi.org/10.1115/1.3162626, 1982. a

Leuthold, R., De Schutter, J., Malz, E. C., Licitra, G., Gros, S., and Diehl, M.: Operational Regions of a Multi-Kite AWE System, in: 2018 European Control Conference (ECC), 12–15 June 2018, Limassol, Cyprus, 52–57, https://doi.org/10.23919/ECC.2018.8550199, 2018. a

Licitra, G., Koenemann, J., Bürger, A., Williams, P., Ruiterkamp, R., and Diehl, M.: Performance assessment of a rigid wing Airborne Wind Energy pumping system, Energy, 173, 569–585, https://doi.org/10.1016/j.energy.2019.02.064, 2019. a, b

Loyd, M.: Crossswind Kite Power, J. Energy, 4, 106–111, 1980. a, b

Malz, E., Koenemann, J., Sieberling, S., and Gros, S.: A reference model for airborne wind energy systems for optimization and control, Renew. Energy, 140, 1004–1011, https://doi.org/10.1016/j.renene.2019.03.111, 2019. a

Malz, E., Hedenus, F., Göransson, L., Verendel, V., and Gros, S.: Drag-mode airborne wind energy vs. wind turbines: An analysis of power production, variability and geography, Energy, 193, 116765, https://doi.org/10.1016/j.energy.2019.116765, 2020a. a

Malz, E., Verendel, V., and Gros, S.: Computing the power profiles for an Airborne Wind Energy system based on large-scale wind data, Renew. Energy, 162, 766–778, https://doi.org/10.1016/j.renene.2020.06.056, 2020b. a

Pasquinelli, G.: An Engineering Model for Power Generation Estimation of Crosswind Airborne Wind Energy Systems, MS thesis, Politecnico di Milano, Milano, https://doi.org/10.13140/RG.2.2.31882.39363, 2021. a, b, c, d

Pierre, C. and Dowell, E. H.: A Study of Dynamic Instability of Plates by an Extended Incremental Harmonic Balance Method, J. Appl. Mech., 52, 693–697, 1985. a

Rapp, S., Schmehl, R., Oland, E., and Haas, T.: Cascaded Pumping Cycle Control for Rigid Wing Airborne Wind Energy Systems, J. Guidance Control Dynam., 42, 2456–2473, https://doi.org/10.2514/1.G004246, 2019. a

Sánchez-Arriaga, G. and Serrano-Iglesias, J.: Modeling and Natural Mode Analysis of Tethered Multi-Aircraft Systems, J. Guidance Control Dynam., 44, 1199–1210, https://doi.org/10.2514/1.G005075, 2021. a

Sánchez-Arriaga, G., García-Villalba, M., and Schmehl, R.: Modeling and dynamics of a two-line kite, Appl. Math. Model., 47, 473–486, https://doi.org/10.1016/j.apm.2017.03.030, 2017. a

Sánchez-Arriaga, G., Pastor-Rodríguez, A., Sanjurjo-Rivo, M., and Schmehl, R.: A lagrangian flight simulator for airborne wind energy systems, Appl. Math. Model., 69, 665–684, https://doi.org/10.1016/j.apm.2018.12.016, 2019. a

Schelbergen, M. and Schmehl, R.: Validation of the quasi-steady performance model for pumping airborne wind energy systems, J. Phys.: Conf. Ser., 1618, 032003, https://doi.org/10.1088/1742-6596/1618/3/032003, 2020.  a

Stuyts, J., Horn, G., Vandermeulen, W., Driesen, J., and Diehl, M.: Effect of the Electrical Energy Conversion on Optimal Cycles for Pumping Airborne Wind Energy, IEEE T. Sustain. Energ., 6, 2–10, https://doi.org/10.1109/TSTE.2014.2349071, 2015. a

Todeschini, D., Fagiano, L., Micheli, C., and Cattano, A.: Control of a rigid wing pumping Airborne Wind Energy system in all operational phases, Control Eng. Pract., 111, 104794, https://doi.org/10.1016/j.conengprac.2021.104794, 2021. a

Trevisi, F., Gaunaa, M., and Mcwilliam, M.: The Influence of Tether Sag on Airborne Wind Energy Generation, J. Phys.: Conf. Ser., 1618, 032006, https://doi.org/10.1088/1742-6596/1618/3/032006, 2020a. a, b, c, d, e

Trevisi, F., Gaunaa, M., and Mcwilliam, M.: Unified engineering models for the performance and cost of Ground-Gen and Fly-Gen crosswind Airborne Wind Energy Systems, Renew. Energy, 162, 893–907, https://doi.org/10.1016/j.renene.2020.07.129, 2020b. a, b

Trevisi, F., Croce, A., and Riboldi, C. E. D.: Flight Stability of Rigid Wing Airborne Wind Energy Systems, Energies, 14, 7704, https://doi.org/10.3390/en14227704, 2021. a, b, c

Trevisi, F., Riboldi, C. E. D., and Croce, A.: Sensitivity analysis of a Ground-Gen Airborne Wind Energy System design, J. Phys.: Conf. Ser., 2265, 042067, https://doi.org/10.1088/1742-6596/2265/4/042067, 2022. a

Tucker, N.: Airborne Wind Turbine Performance: Key Lessons From More Than a Decade of Flying Kites, in: The Energy Kite Part I., edited by: Echeverri, P., Fricke, T., Homsy, G., and Tucker, N., Makani Technologies LLC, 93–224, https://x.company/projects/makani/# (last access: 27 May 2022), 2020. a, b, c, d, e, f

van der Vlugt, R., Bley, A., Noom, M., and Schmehl, R.: Quasi-steady model of a pumping kite power system, Renew. Energy, 131, 83–99, https://doi.org/10.1016/j.renene.2018.07.023, 2019. a

Vermillion, C., Cobb, M., Fagiano, L., Leuthold, R., Diehl, M., Smith, R. S., Wood, T. A., Rapp, S., Schmehl, R., Olinger, D., and Demetriou, M.: Electricity in the air: Insights from two decades of advanced control research and experimental flight testing of airborne wind energy, Annu. Rev. Control, 52, 330–357, https://doi.org/10.1016/j.arcontrol.2021.03.002, 2021. a, b