Flight trajectory optimization of Fly-Gen airborne wind energy systems through a harmonic balance method

. The optimal control problem for ﬂight trajectories of Fly-Gen airborne wind energy systems (AWESs) is a crucial research topic for the ﬁeld, 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 coefﬁcients 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 ﬂuctuates 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 efﬁciency. 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 signiﬁcant loss of power, demonstrating that the absence of high-frequency control is not detrimental to the power generated by Fly-Gen AWESs.


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 Vermillion et al. (2021). Electrical power can be generated by a fixed or a moving ground station, or, alter-natively, 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 F. Trevisi et al.: Flight trajectory optimization of Fly-Gen AWESs through a harmonic balance method 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 Trevisi et al. (2020a), 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 midfidelity 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 (van der Vlugt et al., 2019) assumes the kite as a point mass in steady state for each point of the loop. This model is validated with experimental data (Schelbergen and Schmehl, 2020), 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 (Fagiano et al., 2014;Vermillion et al., 2021). 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 (Cobb et al., 2020) 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 Sánchez-Arriaga et al. (2017, 2019; Sánchez-Arriaga and Serrano-Iglesias (2021) 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 (Fechner et al., 2015;Fechner and Schmehl, 2016) and for fixed-wing AWESs (Licitra et al., 2019;Malz et al., 2019;Eijkelhof and Schmehl, 2022).
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 (Jonkman et al., 2018). The Unicycle and Unifoil models, introduced earlier, are mainly used to compute reference flight paths and for flight path control development (Cobb et al., 2020;Fernandes et al., 2021). To ease the deployment of optimal control problems for AWE, awebox (awebox, 2022) is developed and used, for instance by Leuthold et al. (2018), Haas et al. (2019) and De , to solve optimal control problems. awebox solves optimal control problems in time, imposing periodicity constraints. A similar optimal control problem is studied by Horn et al. (2013) and Malz et al. (2020a, b), 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 . Licitra et al. (2019) 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 Rapp et al. (2019) and Todeschini et al. (2021): 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.
Pasquinelli (2021) 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 (Tucker, 2020) studies the flight trajectories of Fly-Gen AWESs with a simplified quasianalytical 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 Stuyts et al. (2015), 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 Fernandes et al. (2021) (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 (Lau et al., 1982;Pierre and Dowell, 1985;Dimitriadis, 2017). 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.

Flight dynamic model
Two coordinate systems (Fig. 1a) are defined to derive the equations of motion. The ground coordinate system (denoted by F G ) is inertial and centered at the ground station: e x points downwind, e z toward the zenith and e y completes the right-handed frame. For convenience, spherical coordinates are used to describe the position of the airborne unit, with L t the tether length, φ the azimuth angle and β the elevation angle. The spherical reference frame (denoted by F S ) is unequivocally defined at every position with the origin at the AWES center of mass, e r pointing outward the sphere in the radial direction, e φ normal to e r and contained on a plane parallel to x-y, and e β = e r ×e φ . The position p, velocity v and acceleration a vectors projected into the spherical reference frame F S are p = L t e r , v = L tφ cos βe φ + L tβ e β , a = −L tφ 2 cos 2 β − L tβ 2 e r + L tφ cos β − 2L tφβ sin β e φ + L tβ 2 sin β cos β + L tβ e β . (1) The wind velocity is in the positive x-axis direction of F G , and projected into the spherical reference frame is where the wind speed v w as function of the altitude h is modeled with an exponential law: v w,0 is the reference wind speed at the reference altitude h 0 and α s is the wind shear exponent. The relative speed between the AWES and the wind is 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 e 1 is defined to point to the direction opposite to the relative wind speed: The spanwise unit vector s (with origin at the center of mass and pointing in the right-wing span direction) is defined perpendicular to e 1 with the procedure illustrated in Fig. 1b. A second vector e 3 is defined as a unit vector in a plane parallel to the x-z plane with elevation β s (and negative sign): e 3 = − (e x cos (β s ) + e z sin (β s )) , e x = cos φ cos βe r − sin φe φ − cos φ sin βe β , e z = sin βe r + cos βe β .
Note that e 3 points upwind when β s = 0. The unit vector e 2 is then defined as where |e 3 × e 1 | can take values smaller than one because e 3 and e 1 are not defined to be perpendicular in general. In this way, e 2 is perpendicular to the plane e 3 -e 1 . Rodrigues' formula is then used to define s through a rotation of ψ around e 1 , starting from e 2 : 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 e 3 .
The aerodynamic lift L and the drag D take the standard form where ρ is the air density, A is the wing area, and the lift and drag coefficients C L and C D are considered constant. The drag coefficient C D includes the contribution from the tether drag (Trevisi et al., 2020a). The gravitational force F g and the tether force T are F g = −mg sin βe r + cos βe β , T = −T e r , 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 D t is expressed as a linear function of the aerodynamic drag with gain γ : The dynamic equations of motion in compact form read 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 P t (i.e., the power linked to the AWES dynamics) is estimated as a dot product of D t and the relative velocity: The shaft power P s (i.e., the mechanical power that can be converted to electrical power) is modeled using 1D momentum theory (actuator disc) as where the induction a is found by setting the thrust given by momentum theory (T d = 1 2 ρA t (4a(1−a)) |v r | v r ) equal to D t , as in Trevisi et al. (2020b), and A t 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 P s 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 P s . To model the discontinuity in a continuous optimization framework, the logistic function is used: where f is taken equal to 100. 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 (Dimitriadis, 2017). The equations of motion (Eq. 11) can be written as a set of secondorder nonlinear differential equations in the form

Frequency-domain formulation
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 N x : with ω = 2π T being the fundamental frequency of the motion and T the period. Alternatively, the state vector can be expressed as where The first and second time derivatives of the state vector can be found analytically: Similarly, the control inputs, assumed to be periodic, can also be expressed as a Fourier series of order N u : where N u < N x 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 N x : 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 2N x + 1. A Galerkin methodology is then applied by premultiplying 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 2 × (2N x + 1) nonlinear algebraic equations as a consequence of the orthogonality properties of the selected basis of trigonometric functions: 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).

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 where X represents the unknown optimization variables, X * their optimal values, "obj" the objective function, lb and ub the lower and upper bounds of X , 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: The negative value of the mean thrust powerP t (Eq. 12), shaft powerP s (Eq. 13) or electric powerP (Eq. 15) over the loop is taken as objective function, where the symbol. 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): 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 (gradientbased) 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-ofeight-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.

Steady-state model
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 (Loyd, 1980) are here briefly recalled. The thrust power equation, with the assumption of linear crosswind motion (Trevisi et al., 2020b), is where the system glide ratio (including tether drag) is G = C L C D , and, for readability, a modified glide ratio is defined as G t = C L C D (1+γ ) by including the drag of the onboard propellers. The shaft power takes into account the onboard wind turbine induction a: Finally, the power generated and sent to the grid takes into account the efficiencies of the electrical conversion: For high G, the power equation simplifies to For this expression, the value of γ which maximizes the power is only a function of the non-dimensional quantity C D A A t . In Fig. 3, the electrical power P L , normalized with the maximum electrical power at C D the value of γ which maximizes power production decreases. The maximum normalized power as a function of C D A A t is shown in Fig. 4, highlighting that the analytical expression predicts a decrease in power production for increasing C D A A t . The tether force can be evaluated as Trevisi et al. (2020a) showed that for high G, neglecting gravity and with constant incoming wind, there exists an opening angle˜ (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: In this idealized case, the turning radius is R = L t sin˜ and the revolution period is where v L 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: 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 G r is defined as which represents the ratio between gravitational force and aerodynamic lift, similarly to the one introduced in Pasquinelli (2021).
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 (Tucker, 2020) will be used as reference values to present the results (Table 1).

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.

Optimal control problems with constant inflow and no elevation angle constraints
As the analysis is limited to circular trajectories, a cylindrical reference frame F C , similar to the one employed in Trevisi et al. (2020aTrevisi et al. ( , 2021 and Pasquinelli (2021), is used to present the results. A graphical representation of F C is given in Fig. 5. The longitudinal axis of F C is aligned with the mean elevation angleβ. 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.,α > 0). 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 e 3 points up-  wind. In this way, when the roll is equal to zero (ψ = 0 • ), the span direction is perpendicular to the incoming wind.

Optimizing for the mean electrical power in the absence of gravity
For the most idealized case, the gravity is null g = 0, such that F r → ∞ and G r = 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 M t found by solving a number of OCPs with The values of * can be accurately described bỹ which for high glide ratios coincides with the analytical formulation given in Trevisi et al. (2020a).

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 powerP 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 = |v| is imposed to be constant (additional equality constraint). As the control inputs act up to the first har-   monic, this constraint is set by imposing the first Fourier coefficients of the AWES velocity to zero (R = [V 1,s ; V 1,c ]), 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β A ≈ −8.2 • and β B ≈ −7.8 • . The trajectories, shown in Fig. 8, have a circular shape (a circle with radius˜ is marked as − .), but it is no longer a perfect circle. 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 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 P t,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 Pasquinelli (2021). 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.
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 The onboard wind turbine thrust can be approximated with D t ≈ (γ +A γ ,1 cos(α−θ γ ,1 ))D, where D is constant because the AWES velocity is found to be constant, and from Fig. 10 it is found that θ γ ,1 ≈ 180 • . As the thrust power can be written as the product of D t and the relative wind speed v r (Eq. 12), the amplitude of the first Fourier coefficient of γ , considering Eq. (37), can be approximated by 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 F r (F r ∈ [0.1 0.2]) and the analytical approximation given in Eq. (38). Figure 15 shows the first Fourier coefficient of elevation β and azimuth φ as a function of the non-dimensional parameter M t , 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.

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 P 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. 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β B = −5.5 • , higher than for the case without induction. 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 −90 • < α < 0 • and 90 • < α < −180 • ; see Fig. 16), the roll decreases. 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.  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.
In Fig. 4 the dependence of the analytical expression of the shaft power on C D A A 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, F r = 0.1 corresponds to v w = 5.4 m s −1 , F r = 0.15 corresponds to v w = 8.1 m s −1 and F r = 0.2 corresponds to v w = 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.
In Sect. 5.2, it was found that the amplitude of the first Fourier coefficient of γ for C D A A t = 0 (i.e., optimizing for the thrust power) can be approximated by A γ ,1 ≈ G r G. In Fig. 22, the trends of A γ ,1 G r G , with γ being modeled with a single harmonic, for the three investigated Froude numbers, are shown as a function of C D A A t . For C D A A t → 0, trends are close to 1. For increasing F r, the curves collapse to a unique curve. In particular, at C D A A t = 0.23, which is the value for the example, the ratio A γ ,1 G r G → 0.54 for increasing F r. For increas-  ing values of C D A A t , the ratio A γ ,1 G 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γ as a function of C D A A t , highlighting that for increasing F r 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˜ evaluated with Eq. (36). For C D A A t → 0, values are close to 1, as noted in Fig. 15. As C D A A t increases, the values of A β,1 decrease more than A φ,1 , showing that Table 4. Settings of the two optimal control problems maximizing the mean electrical power considering gravity.
OCP N x N γ N ψ Size X Size hP (kW) T (s)  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.

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 γ . Figure 24 shows the trajectory in the β-φ plane. The trajectories deviate from the circular trajectory with opening an-gle˜ (Eq. 36), especially along the β axis, and have a mean elevation angle ofβ A = −5.1 • andβ B = −4.9 • . 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 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 Tucker (2020). 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γ . 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.
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 G r G (see Eq. 35 for definition) for OCP B. As noted when analyzing Fig. 26, at v w = 6 m s −1 A γ ,1 is slightly larger than γ . As wind speed increases up to approximately 8.5 m s −1 , A γ ,1 keeps being similar toγ , meaning that the minimum value of γ is close to zero and so is power (P min ≈ 0). If the wind speed increases again, A γ ,1 becomes lower than γ , meaning that power is always generated over the loop (P min > 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 P s is larger thanγ , results are expected to be modified with respect to Sect. 5.3. Instead, when A γ ,1 is lower thanγ , trends are expected to be equal to the analyses in Sect. 5.3. Indeed, when analyzing Fig. 22, it is found that A γ ,1 G r G → 0.54 for high F r. This means that for low G r G, the first Fourier coefficient of γ can be approximated with A γ ,1 ≈ 0.54G r G, as shown in Fig. 29.  In Fig. 30, the mean power normalized with the power evaluated with Eq. 29 is shown as a function of G r G 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.
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 G r G. For high wind speed (i.e., low G r G), their value is similar to the approximation given in Eq. (36). As G r G 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 v w = 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.

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 h 0 = 100 m, the wind shear exponent is α s = 0.2 and the reference wind speed is v w,0 = 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 e 3 and then the spanwise unit vector s (Eq. 5) is taken as the mean elevation angle β s =β.
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 Trevisi et al. (2021). 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. Figure 32 shows the trajectory for OCP A and B and compares it with a circle of radius˜ centered at an elevation of β = arctan √ α s , as this formulation identifies the elevation of the center of the wind power window (Argatov et al., 2011). 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 • . The roll angle ψ is shown in Fig. 33. In LT-GliDe (Trevisi et al., 2021) 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 circu-  lar 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 F r (Eq. 34) and consequently the gravity ratio G r (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 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γ and A γ ,1 as a function of G r G.γ decreases when G r G increases (i.e., the wind speed decreases). At low wind speed,γ takes low values so that the AWES speed over the loop is higher, which is beneficial to stay airborne. A γ ,1 for low G r G has a linear trend, as noted in Sect. 5.4. When A γ ,1 is equal toγ , the minimum power production over the loop is null. For lower wind speeds (i.e., higher G r G), 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.
To conclude, Fig. 39 shows the ratio of the first Fourier coefficient of β and φ with respect to the analytical expression of the opening angle˜ 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.

Conclusions and discussion
In this work, a novel methodology to study optimal trajectories for Fly-Gen AWESs is introduced. The chosen lowfidelity 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 Tucker (2020) 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 Tucker (2020) 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 (Trevisi et al., 2022), with the aim of improving the power estimation and including an optimal control module.