AWEbox : An Optimal Control Framework for Single-and Multi-Aircraft Airborne Wind Energy Systems

: In this paper, we present AWEbox , a Python toolbox for modeling and optimal control of multi-aircraft systems for airborne wind energy (AWE). AWEbox provides an implementation of optimization-friendly multi-aircraft AWE dynamics for a wide range of system architectures and modeling options. It automatically formulates typical AWE optimal control problems based on these models, and ﬁnds a numerical solution in a reliable and efﬁcient fashion. To obtain a high level of reliability and efﬁciency, the toolbox implements different homotopy methods for initial guess reﬁnement. The ﬁrst type of method produces a feasible initial guess from an analytic initial guess based on user-provided parameters. The second type implements a warm-start procedure for parametric sweeps. We investigate the software performance in two different case studies. In the ﬁrst case study, we solve a single-aircraft reference problem for a large number of different initial guesses. The homotopy methods reduce the expected computation time by a factor of 1.7 and the peak computation time by a factor of eight, compared to when no homotopy is applied. Overall, the CPU timings are competitive with the timings reported in the literature. When the user initialization draws on expert a priori knowledge, homotopies do not increase expected performance, but the peak CPU time is still reduced by a factor of 5.5. In the second case study, a power curve for a dual-aircraft lift-mode AWE system is computed using the two different homotopy types for initial guess reﬁnement. On average, the second homotopy type, which is tailored for parametric sweeps, outperforms the ﬁrst type in terms of CPU time by a factor of three. In conclusion, AWEbox provides an open-source implementation of efﬁcient and reliable optimal control methods that both control experts and non-expert AWE developers can beneﬁt from.


Introduction
Airborne wind energy (AWE) is a renewable energy technology that aims at harvesting strong and steady high altitude winds that cannot be reached using conventional wind technology, at a fraction of the material resources [1]. It is based on the principle of one or more tethered autonomous aircraft flying fast crosswind maneuvers. In the majority of AWE concepts, electricity is either produced by onboard turbines on the aircraft and conducted to a ground station through the tether (drag-mode), or in a periodic fashion by reeling out the tether at a high tension to drive a winch at the ground station, and reeling back in at low tension, so as to achieve a net positive energy output over one period (lift-mode). Although there exist many other interesting AWE concepts, e.g., those based on tethered rotorcrafts [2,3], we will limit the scope of this paper to rigid-wing lift-and drag-mode systems. The reader is referred to [4,5] for a recent and comprehensive overview of the different technologies.
The principle of AWE was first investigated in 1980 by Miles Loyd, who derived an upper limit for the power that could be produced by a crosswind AWE system [6]. Since then, in the past two decades, in particular, AWE has gained increasing interest from both academia and industry, leading to significant technological progress and many small-to medium-scale prototypes, the largest of which was based on a 26 m wing span aircraft [7]. While AWE developers are considering a multitude of different designs, most systems are based on a single-aircraft setup. At this moment, AWE technology is still in a pre-commercial stage, with some companies taking the first steps towards market entry [8]. One of the central unresolved challenges for AWE developers is achieving techno-economic performance at utility scale, i.e., designing systems that produce large amounts of electricity at low cost.
Multi-aircraft systems have been proposed and investigated in the literature as a more efficient and cheaper way of producing utility-scale electricity [9][10][11][12]. In a multiaircraft AWE system, two or more tethered aircraft fly tight crosswind maneuvres around a shared main tether, thereby minimizing the latter's crosswind motion and hence, the associated dissipation losses due to aerodynamic drag. These systems can be up to twice as efficient as their single-aircraft counterparts [10], while having superior, modular, upscaling properties [12], intrinsically smooth power output profiles [13], and higher potential power densities in farm configurations [14]. As a consequence of the increased system complexity, this system class has thus far only been investigated in simulation studies.
A crucial condition for the performance of both single-and multi-aircraft systems is finding power-efficient flight paths that satisfy flight envelope constraints and airframe load limits. This is not only necessary for path planning purposes but also for, e.g., offline performance prediction, design optimization, and control strategy design. Optimal control is an evidently suitable path planning technique for AWE, given its natural ability to handle unstable, nonlinear, constrained systems with multiple in-and outputs. In the past decade, it has become an established method in the field, leading to various applications ranging from performance assessment studies [15], over model predictive control [16], and system identification [17], to flight path planning for a real-world soft-kite system [18]. We refer the reader to [19] for a complete overview of applications.
Despite its obvious advantages, optimal control comes with its own challenges: the dependence on an accurate model; the computational burden associated with finding a numerical solution; and a rather complex implementation that heavily relies on expert knowledge. In case an accurate model is unavailable, it is possible to resort to modelfree, adaptive techniques, such as extremum seeking (ES) [20] or iterative learning control (ILC) [21]. However, in [22], a validated reference model was proposed for a lift-mode, rigid-wing single-aircraft system. While identifying the parameters of such a model is a complicated and time-consuming task [17,23], this shows that deriving a physical model that fits the measurements very well is, in principle, possible. In the following, we will focus on the two other challenges mentioned above.
In order to increase computational efficiency, a general model structure based on non-minimal coordinates was proposed in [24], resulting in smooth dynamic equations of low symbolic complexity. Additionally, since the system nonlinearity gives rise to highly non-convex optimization problems, a feasible initial guess is typically needed for fast and reliable convergence of Newton-type optimization solvers. Such an initial guess is, in many cases, not available a priori. Therefore, a homotopy procedure was proposed that produces a close-to-optimal, feasible initial guess based on a generic, naive one [25]. Combined with a direct-collocation-based transcription method, this led to reported computation times of below one minute for a representative power cycle of a lift-mode AWE system with a six-degree-of-freedom aircraft model [26]. Another homotopy variant was investigated in [27] to efficiently compute drag-mode power cycles for large-scale wind data.
In this work, we consider direct collocation approaches with local support, i.e., using several subintervals with each its own, small set of collocation points. However, for smooth trajectories, it is often more efficient to apply a global support approach [28], using Fourier interpolating polynomials for periodic solutions. This approach has been proposed for drag-mode AWE systems in [29], using a frequency-domain formulation. For lift-mode systems, the optimal solution contains sharp transitions between the reel-out and reel-in phases, and it is a priori unclear which method will be more efficient. For future comparison, the global-support approach can be easily implemented using the collocation framework of the proposed toolbox.
While many open-source AWE control and simulation frameworks exist for singleaircraft models [30,31], and even for multi-aircraft models [32], there are only a few available open-source implementations of optimization methods tailored for AWE. The MATLAB library MegAWES [20] provides an implementation of a megawatt-class system model and of a power optimization algorithm based on ES. The Optimal Control Library OpenOCL [33] provides a user-friendly MATLAB interface for formulating and solving OCPs, which can be linked to the optimization-friendly lift-mode model implemented in the framework OpenAWE [34]. However, this framework does not offer homotopy-based initial guess refinement or multi-aircraft models.
In this paper, we present AWEbox, an open-source Python framework for modeling and optimal control of single-and multi-aircraft AWE systems. The contributions of the software package can be summarized as follows: • Usability: the user specifies only high-level modeling and optimization parameters.
AWEbox implements optimization-friendly system dynamics for single-and multiaircraft systems, for various system architectures and combinations of model options. It automatically formulates typical AWE optimization problems and implements and interfaces the algorithms needed to compute a numerical solution efficiently. • Reliability: AWEbox increases reliability by efficiently computing a feasible OCP initial guess via homotopy methods, based on an analytic initial guess defined using a small number of user input parameters. The framework also implements an algorithm to efficiently and reliably perform parameter sweeps. • Extensibility: within the baseline non-minimal-coordinates structure, users can add new or alternative modeling components (e.g., wind model, aerodynamics, etc.) in a straightforward fashion. The homotopy procedure for initial guess refinement can be extended in a modular fashion so that new model components can be introduced without affecting reliability.
AWEbox is freely available and open-source under the GNU LGPLv3, which allows use in proprietary software. The toolbox heavily builds on lower-level open-source software packages such as CasADi [35], a framework for algorithmic differentiation and optimization, and the nonlinear program (NLP) solver IPOPT [36].
The remainder of this paper is structured as follows. Section 2 discusses the multiaircraft modeling procedure, while Section 3 gives an overview of the optimization ingredients used to formulate and numerically solve periodic AWE optimal control problems. Section 4 then outlines the software implementation details. Section 5 presents two case studies that highlight the efficiency and reliability of the implementation, as well as its multi-aircraft capability. Section 6 draws conclusions based on these results and makes suggestions for further research.

AWE Modeling for Optimal Control
In this section, we define the multi-aircraft topologies considered in AWEbox. We outline the optimization-friendly AWE model structure for six-degrees-of-freedom (6DOF) aircraft dynamics as described in [22] and the extensions made for the multi-aircraft case [12,37]. The different additional model components presented in this section are not novel per se but rather collected and integrated from our own previous work as well as from the literature.

Topology
We consider any tree-structured multi-aircraft topology as previously introduced in [12]. Each tree is described by a set of nodes N , where each node n ∈ N represents the end-point of a tether. All tethers in the tree are assumed to be rigid and straight, which is a reasonable assumption if tether tension is high compared to gravity and tether drag. Note that this assumption allows for a very good agreement with experimental flight data of the single-aircraft system in [22]. In case a higher accuracy is desired, the topology can be readily extended by discretizing the tethers into finite elements, adding a node for each segment, as proposed in [10].
Some of the nodes k ∈ K correspond to aircraft nodes, while other nodes l ∈ L are layer nodes, with L := N \K if |N | > 1 and L := N in the single-aircraft case. The parent map P(n) uniquely defines the interlinkage between nodes, and the children map C(n) := {n ∈ N |P(n) = n} returns the set of nodes with parent n. Figure 1 illustrates the proposed notation for some typical examples.

System Dynamics and Variables
The considered topologies require a multi-body modeling approach that should exhibit certain optimization-friendly properties. For one, the dynamics should have a low symbolic complexity to allow for fast repeated numerical evaluation, in particular of its sensitivities. Second, model nonlinearity should be kept low in order to enable fast and reliable use of Newton-type optimization techniques. Third, the model should avoid singularities that might be visited by and crash the optimization algorithm.
In [24], the efficacy of a non-minimal coordinates modeling approach to describe the translational and rotational dynamics of multiple interlinked aircraft is demonstrated. In this approach, each node is considered as a separate rigid body and linked by algebraic constraints. The aircraft orientation is parametrized in a non-singular fashion using the direction cosine matrix (DCM).
The resulting multi-body models are of reasonable complexity and nonlinearity but result in model equations in the form of index-3 DAEs. This representation does not allow for the deployment of classical integration methods within the optimal control problem [38]. Therefore an index-reduction technique is applied, which involves time-differentiation of the algebraic constraints. The resulting model equations for both lift-and drag-mode AWE systems are summarized by the following index-1 DAE: with associated consistency conditions C(x(t)) = 0, ∀t ∈ R. To define the differential state vector x for both lift-and drag-mode systems, we first define the basic multi-aircraft state This state vector firstly contains q andq that are concatenations of the node positions q n ∈ R 3 and velocitiesq n ∈ R 3 respectively, ∀n ∈ N . These are followed by the states specific to aircraft nodes, namely R, ω, δ, which are concatenations of all R k , ω k , δ k , ∀k ∈ K. The DCMs R k := [ê 1,k ,ê 2,k ,ê 3,k ] ∈ R 3×3 contain the chord-wise, span-wise and upwards unit vectors of the aircraft body frames, expressed in the inertial frame {e x , e y , e z }. All DCMs should be orthonormal, i.e., they are constrained to evolve on the 3D manifold defined by c R,k : where the operator P ut is used to select the six upper triangular elements of a matrix. The aircraft angular velocities ω k ∈ R 3 are given in the body frame. The surface deflections δ k = [δ a,k , δ e,k , δ r,k ] ∈ R 3 of aileron, elevator, and rudder, respectively, give control over the aircraft aerodynamics.
The lift-and drag-mode state vector can now be defined as where the tether length l t ∈ R and speedl t ∈ R describe the main tether reel-in and -out evolution. The variable κ is the concatenation of all κ k ∈ R, ∀k ∈ K, which represent the onboard turbine drag coefficients. The controls u lift := (δ,l t ) and u drag := (δ,κ) (5) are given by the concatenation of all aircraft surface deflection ratesδ k ∈ R 3 and by either the tether accelerationl t ∈ R or the concatenation of the turbine drag coefficient derivativeṡ κ k ∈ R. The algebraic variables z := λ describe the concatenation of all Lagrange multipliers λ n ∈ R related to the tether constraints that restrict the position of each node n ∈ N to evolve on a 2D manifold defined by where ∆q(n) := q n − q P(n) , if n / ∈ K, q n + R n r t − q P(n) , if n ∈ K.
In these constraints, r t is the tether attachment point described in the aircraft frame. The variable l n describes the tether length associated with node n and is defined together with the tether diameter d n as with l s and d s the length and diameter of the secondary tethers and l i and d i those of the layer-linking tethers in stacked multi-aircraft configurations.
The ground station is located at the origin of the inertial frame, such that q 0 := 0. The variables θ represent variable system parameters that can be optimized over. In the general stacked multi-aircraft case, they are defined as θ lift := (l s , d) and θ drag := (l t , l s , d) (9) where d := (d t , d s , d i ).
The constant parameters p allow the dynamics to be evaluated for varying model parameters, such as aircraft wing span, wind model parameters, etc. The system parameter values used in the numerical experiments in this paper are listed in Appendix A, Table A1.

Lagrangian Dynamics
The system dynamics (1) can be derived in accordance with the Lagrangian approach proposed in [24]. The system Lagrangian is defined as with c the concatenation of all tether constraints c n for all n ∈ N and with the kinetic energy T and potential energy V defined as Here, m K is the aircraft mass, and the tether mass m t,n := ρ t l n πd 2 n 4 , with ρ t the tether material density, and g is the gravitational acceleration. The kinetic energy related to the aircraft T K,k and to the tethers T t,n [39] are given by T t,n := 1 6 m t,n (q nqn +q P(n)qP(n) +q nqP(n) ) , with J K the aircraft moment of inertia and with the tether velocity at the ground station given byq 0 :=q 1 e t with e t := q 1 l t . (15) Note that for a drag-mode system with constant tether length, this implies thatq 0 = 0.
With the system Lagrangian defined, the translational dynamics read: with F the concatenation of the external forces F n exerted on each of the nodes. The term F corr is a Lagrangian momentum correction term for open systems: This term is non-zero for lift-mode systems since tether mass and energy are entering and leaving the system due to the reeling motion.
The rotational dynamics are projected on a 3D manifold in the aircraft body frame [24] so as to read: with M k the aerodynamic moment exerted on the aircraft and with U the "unskew" operator, i.e., U(skew(a)) = a, as defined in [24]. Next to the dynamic equations, also the holonomic constraints c = 0 need to be enforced. Since these constraints do not explicitly depend on the generalized accelerations q or on the algebraic variables λ, it is not possible to numerically integrate the resulting dynamic equations with standard algorithms. Therefore an index reduction is performed by differentiating c twice with respect to time. Note thatc depends onq.
Because of the index reduction, as well as the over-parametrization of the rotational degrees-of-freedom, the consistency conditions C(x) := (c,ċ, c R ) = 0 must be enforced at an arbitrary time point in the trajectory. These quantities are called invariants, since their value is preserved by the dynamics. System invariants, when not dealt with carefully, can lead to failure of the Linear Independence Constraint Qualification (LICQ) in the context of periodic optimal control. Performing Baumgarte stabilization on the invariants is an effective way to avoid this issue, while simultaneously ensuring that C(x) = 0 is satisfied over the entire time period [40]. Therefore the tether constraint dynamics are augmented with the following Baumgarte stabilization scheme: with κ t a Baumgarte tuning parameter [41].

System Kinematics
The system dynamics also describe the following trivial kinematics. First, as explained in the previous section, the rotational kinematics are augmented with a Baumgarte-type stabilization on the orthogonality conditions [42]: with κ R another tuning parameter. The remaining kinematics together with (16)-(20) then complete the system dynamics summarized by (1). The remaining modeling effort now focuses on the generalized forces F and moments M.

Wind and Atmosphere Model
AWE systems typically operate at altitudes of several hundreds of meters, and the altitude variation within a typical power cycle is of the same order of magnitude [15]. In particular, the multi-aircraft variant, unhindered by the drag losses caused by main tether cross-wind motion, can theoretically operate at arbitrarily high altitudes, wherever the wind power density is highest [10,12]. Therefore, a wind model is needed that accounts for the varying wind power availability with altitude. Within the community, it is common to use one of the following approximations: (a) Logarithmic profile: A logarithmic model [43] is typically used as a very simple wind shear approximation. Assuming steady, laminar flow, the logarithmic model provides us with the following expression for the freestream wind velocity u ∞ (z): which, in this model, is assumed to be aligned with the x-axis in the inertial frame.
Here, u ref is the reference wind speed that is measured at an altitude z ref , whereas z 0 is the surface roughness length, which depends on local terrain characteristics. (b) Power-law profile: Another frequent approximation is given by the power law: where c f is a ground surface friction coefficient. We will use this model in the numerical experiments in this paper, in accordance with case studies in [15,22]. (c) 2D wind data interpolation: The disadvantage of the logarithmic and power-law models is that they are only useful for representing long-term average wind conditions. Realistic wind profiles come in a wide variety of shapes and are subject to strong short-term (hourly, diurnal, seasonal) changes. Furthermore, the approximation accuracy typically breaks down at altitudes relevant to AWE systems [44] while failing to account for the effect of vertical atmospheric stability. Hence, for accurate optimal-control-based power curve and capacity factor estimation, it is often necessary to generate a more detailed but still differentiable wind model based on highly spatially resolved wind speed measurements. To achieve this, we adopt the approach presented in [27,45]. We assume a wind profile that is represented by discrete 2D wind measurements u m,1 , . . . , u m,n lag ∈ R 2 . These measurements correspond to a set of altitudes z 1 , . . . , z n lag . We can then create a smooth wind model to approximate the measured wind profile, by creating an interpolating function based on Lagrange polynomials: with ζ the concatenation of the polynomial coefficients ζ i ∈ R 2 obtained by solving the following optimization problem The cost function is tuned with weight k so that ζ * i ≈ u m,i , ∀i = 1 . . . n lag , while preventing overfitting via the penalization of the second derivative of the interpolating polynomials. The smooth and differentiable wind model is then given by Wind power availability is linear in the air density, and the atmospheric density drop is non-negligible in the altitudes relevant to AWE. Therefore the density variation with altitude ρ(z) is modeled according to the international standard atmosphere model [43]: where R is the universal gas constant. The parameters T 0 and ρ 0 are the temperature and air density at sea level, and T L is the temperature lapse rate. Alternatively, the user can provide a set of density data points at different altitudes. The continuous density profile ρ(z) is then constructed similar to the wind profile with smoothed Lagrange polynomials, cfr. (24) and (25).

Aerodynamic Model
The apparent wind at each aircraft node k ∈ K is defined as We then define the dynamic pressure as q k := 1 2 ρ(q k e z ) u a,k 2 . The aerodynamic forces (in the inertial frame) and moments (in the body frame) on the aircraft wings are then given using with S the aircraft aerodynamic surface and with the aerodynamic coefficients C F,k := (C X,k , C Y,k , C Z,k ) and C M,k := (C l,k , C m,k , C n,k ), which are a function of the angles of attack α k and side-slip angles β k , given by the small-angle approximations α k :=ê 3,k u a,k e 1,k u a,k and β :=ê 2,k u a,k e 1,k u a,k .
The force and moment coefficients C ,k (with ∈ {F,M}) read as The dependence of these coefficients on α k is approximated by second-order polynomials of the form: with the values of the coefficients c ,•,i used in this study given in Table 2 in [22]. The tether drag is modeled as follows. Consider the infinitesimal tether drag force dF td,n (s) := F td,n (s)ds on an infinitesimal segment l n ds, for s ∈ [0, 1], with: with C D,t the tether drag coefficient, and where the segment position and apparent wind speed are given by It is shown in [39,46] that the total drag force can be exactly distributed into contributions on node n and on its parent node P(n), so as to read (1 − s)F td,n (s)ds (36) respectively. In order to be able to numerically evaluate the tether drag, the integrals in (36) are discretized using the midpoint rule. Typically, a number of M td = 5 integration intervals is sufficiently accurate. The generalized forces can now be defined for each node as if n ∈ N \K (37) and the generalized moments are given by the aerodynamic moments, i.e., M k := M A,k , ∀k ∈ K.
In the drag-mode case, also the braking force of the onboard turbines is acting on the aircraft: Note that the tether pulling force and moment exerted on the aircraft are implicitly modeled in the constraint-based dynamics (16) and (18) and should not be considered as part of the generalized forces.

Power Output
For lift-mode systems, the generated power is the product of the main tether force with the tether speed. The pulling force using a tether n experienced at node n is given by the expression F t,n := −λ n ∇ q n c n . Note that a positive multiplier corresponds to a positive pulling force. The power transferred through tether n is then given by P t,n := F t,nq n . For the main tether, this expression can be simplified to P t,1 := −λ 1 l tlt . The mechanical power that arrives at the ground station is given by P lift := −P t,1 .
In drag-mode systems, electrical power is generated by the onboard turbines and transferred to the ground station through the tethers. Each aircraft k ∈ K generates an amount of electrical power P turb,k := η turb κ k u a,k 3 , with η turb the onboard turbine efficiency. Note that for the case of power consumption, i.e., κ < 0, the efficiency needs to be inverted. This can be implemented using the logistic function, as proposed in [29]. The total power output generated by the drag-mode system is then given by P drag := ∑ k∈K P turb,k .

Optimization Ingredients
In this section, we discuss all the necessary ingredients to formulate, discretize, and reliably solve power optimization problems for the system model described in the previous section. We state the periodic optimal control problem formulation in continuous time, and we discuss common system constraints. We explain the transcription method to convert the problem into an NLP, and we summarize the interior-point solution strategy used by IPOPT to solve it. Then we describe how the initial guess is constructed, and how it can be refined using two different homotopy methods that are tailored for interior-point NLP solvers. Finally, we discuss a third homotopy method that is tailored to perform parameter sweeps with interior-point NLP solvers.

Problem Formulation for Periodic Orbits
The main goal of the toolbox is to facilitate automated computation of dynamically feasible, power-optimal periodic orbits for both lift-and drag-mode systems while satisfying a set of relevant system constraints. In order to achieve this, we formulate a periodic optimal control problem of a free time period T, which has the distinctive property that the system state at the initial and final time of the OCP time horizon can be chosen freely by the solver but must be equal. Given that some key system parameters θ, such as the tether diameters and lengths, have a huge impact on the system power output and the optimal flight trajectories, they are included as optimization variables as well.
Let the optimization variables be defined as w := (x(·), u(·), z(·), θ, T). Then we can compute a power-optimal state and control trajectories and a corresponding system design θ for given parameters p by solving the following continuous-time optimization problem: The Lagrange cost term is given by the sum of the negative power output and a penalty on the controls in order to mitigate actuator fatigue, as well as on the side slip angle and the angular accelerations in order to avoid aerodynamic side forces and aggressive maneuvers: withŵ(t) := (u(t), β(t),ω(t)) and W a constant diagonal weighting matrix. The variables β andω are the vertical concatenations of the side slip angles β k and angular accelerationsω k , ∀k ∈ K. Proper tuning of the weighting matrix W is necessary to achieve fast convergence of the optimization algorithm as well as to obtain a locally unique solution. We refer the reader to the open-source code for the weighting factors used in the numerical experiments in this study. The function ψ : R n x → R is used to impose a technical constraint that removes the phase invariance inherent to periodic OCPs. For lift-and drag-mode systems, this function is different and reads as either The inequality constraints h are discussed in the following section. Note that the consistency conditions C(x(t)) = 0 are not enforced at any given time within the time horizon of the OCP. In combination with the periodicity constraint (39d), this would lead to LICQ deficiency for all feasible trajectories. There exist several technical solutions for this issue [40]. In the dynamic correction approach chosen here, Baumgarte stabilization is applied to the consistency conditions in the system dynamics, as previously mentioned in Section 2.3. Therefore the dynamics of C are exponentially stable, and since with the value of periodicity, it holds that C(x(0)) = C(x(T)), the only feasible periodic state trajectories are those where C(x(t)) = 0, ∀t ∈ [0, T].

System Constraints
A particular feature of OCP (39) is that it has an economic cost function, which is not lower bounded, as opposed to tracking cost functions [47]. OCPs with an economic cost function tend to have extreme solutions in the absence of constraints. In the context of AWE power optimization, it is, therefore, crucial to impose constraints that avoid a violation of the flight envelope, and that preserve the structural integrity of the airframe and the tether.
The flight envelope consists of upper and lower bounds on the angle-of-attack α (to avoid stall) and the side-slip angle β (to avoid additional drag and preserve model validity) for all aircraft in the system. Additionally, the stress in the tethers should not exceed the yield strength with a certain safety factor f s : Here, the tether force magnitude can be simplified to F t,n = λ n l n , following the definition in Section 2.7. The aircraft orientation is also constrained in order to avoid collision of the airframe with the tether, which might occur during sharp turns in transition maneuvers: where γ max is the maximum angle between the tether vector and the upwards unit vector of the aircraft body frame, which should be set lower than at most π/2. In the multi-aircraft case, the following anti-collision constraints might be included: where f b is safety factor in multiples of the wing span b.
Along these nonlinear constraints, variable bounds are typically imposed on variables such as flight altitude, tether length, speed and acceleration, aircraft angular velocity, control surface deflections and their rates, etc. One pair of variable bounds that is crucial in the context of periodic optimal control are the bounds on the time period T. Since the OCP will be discretized in a discrete number of numerical integration intervals, the integration accuracy is variable along with T. Therefore, T should be bounded from above to guarantee an acceptable simulation accuracy. Also, by translating a priori knowledge on the optimal value of T into variable bounds, we narrow the search space and exclude many possible local solutions, which typically increases reliability and speeds up the convergence of the NLP solver.

Problem Transcription
The continuous-time OCP (39) has an infinite number of variables and constraints. Hence, we apply direct optimal control to transcribe the OCP to an NLP. We choose transcription by direct collocation, which is a fully simultaneous approach, where the numerical simulation variables are treated as variables in the optimization problem [26]. We choose this approach for the following reasons.
First, fully simultaneous optimal control is characterized by faster contraction rates of the Newton-type iterations compared to simultaneous and sequential optimal control, in particular for highly nonlinear and unstable systems [48]. Second, in the fully simultaneous case, the simulation problem is solved directly by the NLP solver, which is typically more robust than the rootfinder used in standard available numerical integrators. Finally, since OCP (39) is highly non-convex, the NLP solver benefits from computing the Newton step using exact Hessian information. The NLP Hessian becomes considerably cheaper to evaluate in the fully simultaneous approach.
Although the resulting direct collocation NLP is comparably large, it is also sparse. In combination with a sparsity-exploiting NLP solver, direct collocation is a highly efficient transcription method for the models presented in this paper.
In direct collocation, the time horizon is divided into N (usually equidistant) intervals described by [t i , t i+1 ], where 0 < t 0 < t 1 < · · · < t N = T. The control trajectory is parameterized as a piecewise constant functionũ(t) := u i if t ∈ [t i , t i+1 ). The state trajectory is parametrized by piecewise polynomials of order M + 1, i.e.,x(t) : Note that it holds thatx i (t i + ∆t i τ j ) = x i,j . The state derivative is given by the derivative of the polynomials, i.e.,x The algebraic variables are also discretized in each i'th time interval as z i,j , and allocated to the collocation points τ 1 , . . . , τ M .
Let us now define x i := x i,0 , X i := [x i,1 , . . . , x i,M ] and Z i := [z i,1 , . . . , z i,M ]. Then, for given state vector x i at the start of each interval, the collocation variables X i and Z i are uniquely determined by enforcing the system dynamics (1) at the grid points τ 1 , . . . , τ M : The state transition from one interval node to the next is given by the equation The system of Equation (48) corresponds to that of an implicit Runge-Kutta integration scheme, where the choice of collocation grid points τ uniquely defines the Butcher-Tableau of the specific integration method. Here, we choose as collocation grid points the roots of Gauss-Radau polynomials, more specifically those corresponding to the Radau IIa integration scheme because of its high order accuracy and its excellent stability properties (A-and L-stability), which is particularly relevant for DAE systems [49].
Further, the inequality constraints are imposed on the interval nodes, and the Lagrange term in the cost function can be computed via a quadrature rule [26]: where the quadrature weights are given by The NLP resulting from discretizing the OCP (39) using direct collocation is then formulated as min w with the decision variables summarized by w := (x 0 , X 0 , Z 0 , u 0 , x 1 , . . . , u N−1 , x N , θ, T).
For the remainder of this text, we will write NLP (52) in more compact form as the parametric NLP

Solution Strategy
There are two common solution approaches for inequality-constrained nonlinear programs such as (53): sequential quadratic programming (SQP) methods and interiorpoint (IP) methods [50]. SQP methods are based on iteratively solving a series of convex quadratic programs (QP) that are local approximations of the NLP. IP methods, on the other hand, perform iterations directly on a relaxed version of the Karush-Kuhn-Tucker (KKT) system corresponding to NLP (53), which reads as where 1 denotes a vector of ones. Together with the conditions µ ≥ 0, s ≥ 0, the KKT system (54) for barrier parameter τ = 0 gives the first-order necessary conditions of optimality. However, in this case, the KKT system is non-smooth due to the complementarity condition diag(s)µ = 0, and therefore difficult to solve with Newton-type methods. Therefore, in IP methods, the iterations generally start on a smooth KKT system related to a barrier parameter τ 0 > 0, which is then gradually reduced to a smaller value τ f > 0, so that the final solution approximates the exact solution of (54) up to sufficient accuracy. It holds that w * τ − w * 2 = O(τ), where w * τ and w * are the solutions to the KKT system for τ > 0 and for τ = 0, respectively.
The advantage of IP methods is that the iterations are computationally cheaper compared to those of SQP methods: per iteration, only one linear system has to be solved, as opposed to one QP of equal size. Additionally, because IP methods start iterating on a problem with relaxed inequality constraints, and only gradually tighten these constraints, they are particularly robust in case little or no a priori knowledge on the active set of the optimal solution is available, as is typically the case for AWE systems.
In this work, we use the interior-point NLP solver IPOPT [36] in combination with the linear solver MA57 [51]. IPOPT implements a particularly reliable algorithm that implements a filter line search method for globalization [52]. The algorithm also exploits the sparsity of the direct collocation NLP which makes it particularly efficient for this application.

Circular Initial Guess Construction
In order to efficiently converge to a solution of a highly nonlinear, non-convex NLP, even a robust NLP solver such as IPOPT typically requires a good initial guess. Therefore we propose here a circular flight trajectory initialization based on a limited number of user-defined parameters π 0 : whereq 0 is the aircraft flight speed, N 0 l the number of loops, l 0 t the initial tether length, θ 0 e the (average) elevation angle of the main tether, and θ 0 c the trajectory cone angle with respect to the average main tether vector. The angle ϕ 0 denotes the phase angle with which the periodic initial guess can be shifted in time. The parameter θ 0 is a direct guess for the system parameters θ.
Building on the parameters π 0 , we then define a stationary tether frame as e 1 := cos(θ 0 e ) · e x + sin(θ 0 e ) · e z , e 2 := e 1 × e x e 1 × e x , e 3 := e 1 × e 2 , after which we can define for each aircraft k a frame that is rotating about the main tether: with the rotation angle ϕ k (t) for each aircraft defined as ϕ k (t) := ϕ 0 + ω 0 t + 2π(k − P(k) − 1)/|C(P(k))\L| , (58) and with the rotation radius and speed, and the time period of one loop defined as respectively.
In the general multi-aircraft case, the node positions and (angular) velocities are then initialized at each time point on the collocation grid t i by In the single-aircraft case (K = L = {1}), the aircraft position is initialized using (61). The layer node velocities are set to zero.
The aircraft DCMs are initialized so that the initial guess meets the flight envelope constraints rather than exactly satisfying the kinematic relation (20). The apparent wind speed for each aircraft at time t i is given by with u ∞ (·) as the user-defined wind profile. The DCM is then initialized to have a zero angle of attack and a zero side-slip angle: e 0 3,k,i ←ê 1,k,i ×ê 2,k,i , ∀k ∈ K .
The tether multipliers are trivially initialized as λ 0 n ← 1 Nm −1 , ∀n ∈ N , to ensure a strictly positive tether force. All remaining states and controls are initialized as zero. Finally, the initial overall cycle period is set to T 0 ← T 0 l N 0 l . The initial guess is summarized with the vectorw 0 . In the following, we will refer to the method which usesw 0 as an initial guess for solving P f as "NH" (no homotopy).

Homotopy-Based Initial Guess Refinement
Even the educated initial guess defined in the previous section often leads to very slow convergence or even solver failure when solving P f (p). In order to increase computation speed and improve reliability, we propose a refinement procedure based on homotopy methods that reliably produces a close-to-optimal, feasible initial guess based on the analytic user-defined initialization.
The basic idea is to first solve a trivial version of the intended NLP, and then to repeatedly compute the solution while updating the NLP in a controlled and smooth way to the full nonlinear final problem. Homotopy methods (also known as continuation methods) are widely used in the field of non-convex optimization when little or no a priori knowledge of the location of the optimal solution is available [53,54]. Homotopy methods were originally introduced in the field of AWE optimization in [25]. In this paper, we generalize this approach for multiple homotopy stages and discuss particularities when using interior-point methods.
First we construct a homotopy problem H c (p, φ), with homotopy parameters φ ∈ R n φ and φ i ∈ [0, 1], ∀i ∈ {1, . . . , n φ }. Note that φ can be multidimensional to allow for a stepwise introduction of distinct model nonlinearities or couplings. The homotopy problem is defined as with the NLP functions Φ H , G H and H H defined such that H c (p, 1) = P 0 (p) and H c (p, 0) = P f (p). Here, P 0 (p) is a simplified problem that is trivial to optimize for a large set of initial guesses, and P f (p) is the target optimization problem defined in (53). It can be shown that, if H c (p, φ) satisfies the LICQ and second-order sufficient conditions (SOSC) for all p and φ, there exists a unique and piecewise smooth homotopy path w * (p, φ) between the optimal solutions w * (p, 0) and w * (p, 1) [54]. Algorithm 1 (CIPH) describes a classic procedure to follow the homotopy path w * (p, φ). First, we provide an initial guessw 0 , which is the approximate solution of the initial problem P 0 (p). Then, for each step i in the multi-step homotopy, we reduce the homotopy parameter φ i from one to zero with an increment 1 γ in a total of γ iterations. At every iteration, the homotopy problem H c (p, φ) is solved up to a certain (low) accuracy level, while the NLP solver is warm-started with the solution of the previous iteration. To improve performance, the maximum number of NLP iterations can be limited in this stage.

Algorithm 1 Classic Interior-Point-based Homotopy (CIPH)
Require:w 0 , p, γ > 0 Output: The output of the homotopy, then, is an approximate solutionw f to the intermediate problem H c (p, 0), which can be used as an initial guess for solving P f (p) up to high accuracy. If the LICQ and SOSC conditions are fulfilled, there exists a high enough value of γ to guarantee convergence of this algorithm [54] (Theorem 5.2).

Penalty-Based Homotopy
The fixed-step continuation approach described in the previous section is simple to implement and works well in practice [27,45,55]. Nevertheless, it has two drawbacks. First, the choice of a fixed homotopy parameter step renders the algorithm less robust than if an adaptive-step strategy were used. Second, in terms of computational efficiency, γ · n φ NLPs need to be solved by default even when larger steps would be feasible.
Of course, adaptive step size strategies for the homotopy path following exist and are well-established [53,54]. However, they increase the complexity of the algorithm as well as the number of hyperparameters to tune. Therefore we propose a simple but effective variation of Algorithm 1, which can be used in particular when the chosen NLP solver is a well-globalized solver. The idea is to use the underlying globalization routines (e.g., line-search) of the NLP solver to choose a suitable homotopy parameter step size.
The resulting homotopy strategy is penalty-based and builds on the reformulation H p of NLP (68) to read: In this formulation, the parameters φ are treated as decision variables with a high linear penalty S ∈ R n φ + . The homotopy path is now parametrized by the bounds on φ, i.e., Algorithm 2 (PIPH) describes the alternative homotopy procedure. The lower bounds φ are successively set to zero for each homotopy stage, allowing the NLP solver to find a path for the homotopy parameter φ i in stage i, while simultaneously applying correction steps to the decision variables w. Afterward, the problem is solved again with φ i = 0 to ensure completion of the homotopy stage.

Algorithm 2 Penalty-based Interior-Point-based Homotopy (PIPH)
Because of the high linear penalty on φ, the NLP solver will take the largest possible parameter step that is acceptable to the line-search filter, hence providing both robustness and speed. Additionally, only 2 · n φ NLPs need to be solved instead of the γ · n φ NLPs in the classic continuation homotopy. This can allow for a significant speed-up even if the number of iterations per NLP solve is naturally higher.
Note that the convergence of Algorithm 2 is only guaranteed for small enough updates of the parameterφ. In practice, however, convergence is almost always achieved for jumps from 1 to 0.

Interior-Point-Based Homotopy
The homotopy methods presented above are based on the idea of solving a sequence of closely related problems, where the solution of each problem is used to warm-start the next. However, because an interior-point NLP solver by default starts iterating on the relaxed KKT problem (54) (with a high barrier parameter τ), it is unable to exploit the (active set) information contained in the initial guess, if it is the solution to the non-smooth KKT problem. To circumvent this issue, we apply the following barrier strategy [27,56,57] throughout the homotopy:

1.
The initial problem P 0 is solved from an initial barrier parameter τ 0 to an intermediate τ i < τ 0 , so that the KKT system remains sufficiently smooth.

2.
The homotopy problem H p is repeatedly solved for constant barrier parameter τ i . 3.
The final problem P f is solved from τ i to a final value τ f < τ i . Using this strategy, the Newton iterations quickly converge from one intermediate problem to the next during the homotopy stage.

Homotopy Design
In this paper, we propose two homotopy stages (n φ = 2). The initial problem P 0 (p) thus comprises two alterations with respect to the final problem P f (p). Firstly, the aerodynamic forces and moments in the model are replaced with the direct force controls F f,k ∈ R 3 and moment controls M f,k ∈ R 3 for all k ∈ K, which are then added to the control vector u. This step relaxes the nonlinearities and couplings related to aerodynamics [25]. Secondly, the initial problem does not optimize the average power output but rather the tracking error with respect to the user-generated initial guess.
The homotopy problems H c (p, φ) and H p (p,φ) are then constructed by replacing F A,k and M A,k in (37) with as well as by changing the stage cost function to withx 0 (t) the initial state trajectory guess. Additionally, in order to reduce the initial degrees of freedom, the system parameters are fixed to their initial values until the second homotopy step. The system parameters are thus only free optimization variables when the cost function transitions from tracking error to power output: Substituting Equations (70)-(72) into the model, cost function and constraints, we obtain after repeated discretization with direct collocation the functions Φ H (w, p, φ), G H (w, p, φ) and H H (w, p, φ).

Parametric Sweep Warmstarting
Once a solution for NLP (53) has been found, it is often interesting to investigate the sensitivity of the optimal solution with respect to one or more of the model parameters p. A typical example is when we compute the NLP solution for different values of u ref (in the case of a logarithmic or power-law wind profile) to compute a power curve for a particular AWE system. One approach is to apply Algorithms 1 or 2 to compute a solution for all parameter values based on the same initial guess. However, in case the distance between the different parameter values is small, it is more efficient and more reliable to compute an initial guess for one problem from the solution of the previous one.
Algorithm 3 (SIPH) describes how an initial guess for each problem in the set of NLPs P f (p i ), for i = 1, . . . , p, can be generated efficiently. It starts based on the solutionw f,0 of the homotopy problem H c (p 0 , 0) for an initial set of parameters p 0 . This initial solution can be computed using CIPH or PIPH. We assume that the sequence of parameter vectors p 0 , . . . , p p is ordered so as to minimize the distance from one parameter set to the next, as proposed in [27]. Then, we can compute the initial guessw f,i for problem P f (p i ) from the guessw f,i−1 for the previous problem P f (p i−1 ), by updating the parameter vector p from one value to the next via linear interpolation in γ p steps and by recursively solving the problem H c (p, 0). We employ the same barrier strategy as in Section 3.8 and keep the barrier parameter at a constant value τ i while solving H c , to guarantee a smooth transition from one problem to the next.

The AWEbox Software Package
The goal of the AWEbox software package is to provide a user-friendly interface that facilitates the automatic construction of the optimization-friendly dynamics (1). It formulates the power optimization problem (39) and reliably finds a numerical solution. The toolbox is written in Python 3 and relies heavily on the following software packages: CasADi, an opensource symbolic framework for algorithmic differentiation and nonlinear optimization [35]; the interior-point NLP solver IPOPT [36]; and (optionally) the linear solver MA57 [51]. The six main classes and the basic structure of the package are shown in Figure 2, including the dependencies on the external packages. Starting at the lowest level, the Model-class takes the user-provided modeling options and assembles the according state, control, and algebraic variable vectors. Then, the dynamics (1), relevant constraints, and intermediate model outputs are constructed as CasADi Function objects. Table 1 gives an overview of the main modeling options implemented in AWEbox. Central here is the use of CasADi to compute the partial derivatives of the system Lagrangian in (16). Finally, the Model class can also be used in standalone mode, e.g., in case the user is interested in obtaining the dynamics for simulation purposes only.
The NLP class the dynamics and constraints receives from a Model instance and constructs the NLP functions Φ H , G H and H H as CasADi Function objects, using the direct collocation approach presented in Section 3.3.
From a practical viewpoint, it is essential for the convergence of the NLP solver that all variables, equations, and cost terms are properly scaled. Therefore, AWEbox implements a heuristics-based scaling procedure based on the system parameters and the user-defined initialization. We refer the reader to the open-source implementation for the scaling factors obtained in the numerical experiments in this study.

Options Variants
Topology single-aircraft, i.e., (1, 1) multi-aircraft, e.g., (1, 2) stacked multi-aircraft, e.g., ( The NLP functions are then passed on to the Optimization class, where their first-and second-order derivatives are constructed using CasADi, which also provides the interface to IPOPT. The Optimization class then constructs the initial guess from Section 3.5 and implements both Algorithms 1 and 2 to prepare the homotopy-based initial guess for solving Problem (53). It is also possible to warm-start of the solver with a user-provided initial guess. Finally, Problem (53) is solved up to high accuracy. The default linear solver for computing the Newton step within IPOPT is MUMPS, but in general, higher performance in terms of speed and reliability is reached using the solver MA57, which has to be installed separately.
On a higher level, the central class with which the user interacts is the Trial class, which knits together the functionality of the lower-level classes. To start with, the user can specify modeling options, physical parameters, discretization options, initialization parameters, etc., as in the (non-exhaustive) example given in Listing 1.
Based on the specified options, the user can create a Trial object, and build the system dynamics, constraints, and NLP functions, including derivatives, as shown in Listing 2. In this example, the power optimization is then solved using the penalty-based homotopy. The Trial class then performs some quality checks on the numerical accuracy of the solution, e.g., by checking consistency condition satisfaction. The class also contains some basic plotting functionality for visualizing the optimal solution. The high-level class Sweep, which builds on the Trial class, can be useful for parametric sweeps, as illustrated by Listing 3. This class builds the parametric NLP functions and their derivatives only once, and implements Algorithm 3 for warm-starting the neighboring NLP problems. The MPC class uses the Trial class and the lower level classes to construct the tracking MPC problem as defined in [58]. The class takes as an input the optimal solution of Problem (39) to construct a periodic reference on the MPC time grid. It also takes care of correct initialization, an initial guess, and periodic reference shifting. The MPC problem can then be recursively solved using IPOPT with the warm-starting strategy from [57]. The main goal of this class is not to provide highly efficient numerical solvers aimed at embedded optimization, such as those implemented in the software packages acados [59] or PolyMPC [60]. Rather, this class provides a reliable controller that conveniently allows for offline closed-loop simulations. Listing 4 provides an example of how such a controller can be constructed and evaluated based on (solved) Trial-object. Although the focus here is reliability and not computational efficiency, the user can also code-generate and compile the MPC solver functions using CasADi for use in an external codebase or for embedded application.

Numerical Results
This section discusses two numerical case studies that highlight the contributions of the AWEbox software package. In the first case study, we discuss and compare the computational performance and robustness of the homotopy algorithms CIPH and PIPH, while solving a single-aircraft lift-mode reference problem. In the second case study, we compute a power curve for a dual-aircraft lift-mode system and compare the performance of the algorithms PIPH and SIPH.
In this paper, we focus mainly on algorithmic developments and performance. However, for validation purposes, we also provide an example in the toolbox repository which re-implements the case study performed in [15] with a very good agreement of results. Small remaining discrepancies can be explained by tuning differences since the exact regularization weights and scaling factors used in this case study are not available. For brevity, we do not discuss these results further in this text.

Single-Aircraft Case Study
The first reference problem aims at finding an optimal power cycle for a lift-mode single-aircraft system, with K = {1}, L = {1}, and P(1) = 0. The aircraft parameters are taken from the Ampyx AP2 reference model presented in [22]. We adopt the same wind profile and atmosphere model as presented in [15]. We assume a "reinforced" version of the AP2 airframe since the real-world airframe load limits lead to an overly pessimistic average power output estimate. Therefore, compared to the OCP in [15], the airspeed limits and tether force limits are omitted and replaced only by a tether stress constraint, while the tether diameter d t is no longer fixed and is treated as an optimization variable. Table A1 in Appendix A summarizes the model parameter values of this reference problem, while Table 2 lists all variable bounds and path constraints. We construct the NLPs (68) and (69) using N = 100 intervals with Radau collocation polynomials of the order M = 4, and the controls are discretized using a piecewise constant parameterization. The resulting NLPs have 15,334 variables, 14,323 equality constraints, and 600 inequality constraints. We solve the problem on an Intel Core i7 2.5 Ghz, 16 GB RAM.
The homotopy meta-parameters are experimentally tuned to minimize the associated CPU time. The intermediate homotopy barrier parameter is chosen as τ i = 10 −2 . For CIPH, the number of parameter update steps per stage are γ 1 = 10 and γ 2 = 1. For PIPH, the homotopy parameter penalties are S 1 = 10 2 and S 2 = 1.
In the following, we wish to investigate the performance and robustness of CIPH and PIPH, compared to the case where the user-provided circular initial guess is applied without refinement ("no homotopy" -NH). For this purpose, the reference problem described above is solved for each method for a set of 100 uniformly sampled initialization parameters π 0 from the set defined by π 0 lb ≤ π 0 ≤ π 0 ub . In the NH-case, performance heavily depends on the a priori knowledge of the user. To account for this fact, we introduce two different users. "User A" is an AWE developer with little a priori knowledge of the location of the optimal solution. Therefore, this user has samples from a wide initialization parameter set. "User B," on the other hand, is a control engineer who is familiar with the system and its optimal behavior for the given conditions. Therefore, User B samples from a parameter set that is defined by a range that is a factor 3 smaller than that of User A, centered around the average parameters as evaluated at the solution of interest. Table 3 summarizes the sampling range for all initialization parameters, for both User A and B. The initial number of loops is chosen to be N 0 l = 1.

Description Variable Min (A) Max (A) Min (B) Max (B) Units
Flight  Figure 3 shows the CPU timing results resulting from the initialization sampling by User A and User B. For User A, NH leads to highly variable CPU timings, ranging from a peak timing of up to 12 min down to a minimum of 15 s. In two cases, NH does not converge as it exceeds the maximum number of iterations of the NLP solver. The minimum NH-timing is 50% lower than the best timings of the homotopy methods. Hence, it is possible for User A to "get lucky" and converge to a solution very fast without initialization refinement. However, the peak NH timing is 8 times higher than the worst PIPH timing and almost six times higher than the worst CIPH timing. The average NH timing is a factor of 1.7 times higher than in the PIPH case and a factor of 1.3 higher than in the CIPH case. Therefore, User A benefits significantly from CIPH/PIPH in terms of expected computational performance and, in particular, in terms of timing consistency. PIPH is, on average, 13% faster than CIPH, while the peak timing is 30% lower. For User B, with much better a priori knowledge, the computation times of NH significantly improve compared to user A: average timings are reduced by a factor of 2.4, to a value slightly lower compared to CIPH/PIPH for User B. The peak NH timing is reduced by a factor of 1.5, which is still a factor 5.5 larger than compared to CIPH/PIPH. Thus, while User B has a slightly better-expected performance in the NH case, he or she can still profit from the improved timing consistency provided by CIPH/PIPH. The difference in timings for the CIPH and PIPH methods is almost negligible. The average timings of these methods do not change much compared to the timings obtained for User A. This highlights the property that by pre-structuring the optimization path, the homotopy methods are not able to exploit a priori user knowledge to achieve better average performance.
Overall the PIPH/CIPH CPU timings range between 30 and 100 s. This is comparable to the CPU timing range reported in [26] for similar model complexity and identical collocation grid (but excluding homotopy timings).
The NLP (53) has multiple local solutions, and the choice of optimization algorithm influences the frequency with which certain solutions are found by the optimizer. In the experiments for User A, a total of nine different local solutions were found. Figure 4 shows the dominant, circular, optimal solution, while Figure 5 shows as an example the third most frequent optimal solution, which is characterized by the well-known lemniscate flight pattern. Table 4 summarizes, for each method, the frequency of local solutions.   The homotopy methods almost always converge to the main solution of interest: out of 100 trials, 100 for PIPH, and 98 for CPIH. In the NH case, on the other hand, this is only the case for 71 trials, while failing to converge in two cases. Hence, the homotopy methods not only improve performance and reliability for User A, but they are also more stable in terms of optimization outcome. For User B, all methods always converge to the main solution.
When comparing the different local solutions, we notice that average power output increases up to 22% with respect to the main solution for solutions with longer optimal time periods T * . The solutions with a longer time period typically consist of more than one loop, which leads to a better ratio of reel-out vs. reel-in time, and, thus, a higher "pumping efficiency." This is in line with the results reported in [26].
Note that for increasing time period T * , consistency condition satisfaction decreases. This is because the consistency condition trajectory is the periodic solution to the stable uncontrolled dynamics of the invariants. Hence, as simulation accuracy decreases, consistency conditions are moving away from the theoretically optimal solution of a constant zero value. For this reason, AWEbox automatically computes the consistency conditions for each solution and gives out a user warning once a threshold is reached. The user can then increase the number of collocation intervals and the integration order, or lower the upper bound on T if applicable.
To give more structure to the problem, we propose the following modification to the OCP. We divide the time horizon into two separate intervals with associated time variables T 1 and T 2 , and we define the total time period as T := T 1 + T 2 . We then impose that the first interval is a single reel-out phase, and the second one a single reel-in phase: In the discrete time grid, 70 time intervals are allotted to the reel-out phase, and 30 intervals to the reel-in phase. The resulting NLP has 33,464 variables, 31,550 equality constraints, and 1402 inequality constraints.
The intermediate barrier parameter is tuned manually to be τ i = 10 −4 for both PIPH and SIPH. The PIPH tuning is the same as in the single-aircraft case. SIPH performs a homotopy with γ p = 10 steps for every new parameter value. Additionally, the maximum number of NLP solver iterations is limited to 100 for both methods.
We search for solutions with three loops, i.e., N 0 l = 3. The reason for this is that the resulting trajectories fit well inside the time period bounds defined in Table 2, for all considered wind speeds. The remaining initialization parameters are set toq 0 = 50 m s , The secondary tether diameter is initialized under the assumption that the secondary tether force equals the main tether force divided by two.
The power curve for the proposed dual-aircraft system is obtained in the following manner. First, the optimal trajectory and design are computed with PIPH for a reference wind speed of u ref = 10 m s . The resulting optimal design is given by l * s = 142.9 m, d * t = 4.3 mm and d * s = 3.2 mm. The average power output isP * = 42.0 kW. Note that this is more than a factor of four higher than the single-aircraft solutions in the first case study, while the number of aircraft has only doubled. The power per wing surface area is thus more than doubled as a result of the reduced main tether drag and higher flight altitude. This is in line with the results reported in [10,12].
The optimal design parameters are then fixed, and the NLP is re-solved for u ref ranging from 0 m s to 20 m s . This is done once with PIPH, every time starting from the identical user-defined initial guess. Then, it is done once using SIPH in two separate sweeps: once downwards and once upwards starting from the solution for u ref = 10 m s . Figure 6 shows the power curve obtained with SIPH, and, additionally, for each wind speed, the optimal time period. Similar to the power curve computed in [37], we identify three operational regions. In the first region of zero wind speed up to the cut-in wind speed u ref = 3 m s , power is consumed to keep the system airborne. The upper bound on T is active here, as the aircraft glide downwards about an almost vertical rotation axis during the reel-out phase. In the reel-in phase, potential energy is injected back into the system as the aircraft fly slow upward trajectories. In the second operating region, power grows cubically until the design wind speed is reached. In the third region, power output still grows with the wind speed, but cubic growth is curtailed in order to satisfy the tether stress constraints. The main strategy to limit power output here is to increase the tether reel-out speed so as the decrease the available wind. The optimal time period increases with respect to the design wind speed, as the reel-out speed increases, while the reel-in speed is constrained and cannot grow proportionally. Figures 7 and 8 illustrate the reel-out and reel-in trajectories for u ref = 18 m s .    SIPH outperforms PIPH at every single wind speed (where PIPH converges), but particularly at low wind speeds, when the optimal solution diverges significantly from the user-defined initial guess. Up until the wind speed of 15 m s , the average CPU time is 5 min and 23 s for SIPH and 15 min and 20 s for PIPH.

Discussion
In this work, we presented AWEbox, an open-source Python toolbox for the modeling and optimal control of single-and multi-aircraft AWE systems. We discussed the general multi-aircraft modeling structure, optimization ingredients, and implementation details needed to efficiently compute power-optimal orbits for a wide range of system architectures and modeling options. In particular, we proposed and implemented two interior-pointbased homotopy method variants, in order to increase the performance and reliability of the optimization algorithms. These methods produce a feasible initial guess for the underlying NLP solver, based on an analytic initial guess shaped by the software user. In a numerical experiment, a reference single-aircraft problem was solved for a large set of different initial guesses.
The penalty-based homotopy method reduced the average and peak CPU timing with a factor of 1.7 and 8, respectively, compared to the case when no homotopy method was applied by a user with little a priori knowledge. With good a priori knowledge available, the homotopy methods did not improve performance, but still, the peak CPU timing was reduced by a factor of 5.5. Overall, computation times were in the range of 30-100 s, which is competitive with those reported in the literature. Additionally, the penalty-based homotopy method consistently led to the same local solution, whereas the no-homotopy case resulted in different local solutions in 29 out of 100 cases.
In a second case study, we computed a power curve of a dual-aircraft AWE system and compared the performance of the penalty-based homotopy method of the previous case study with that of a third homotopy method tailored for parametric sweeps with interiorpoint NLP solvers. The penalty-based method was not able to converge to a solution for all wind speeds, while the sweep method succeeded in doing so while outperforming the penalty-based method, on average, by a factor of three in terms of CPU timings. The average CPU timing per NLP solution was about 5 min.
Future work might entail model accuracy improvements, in particular concerning tether and induction modeling The tether model can be readily improved to include sagging by discretizing the tethers into finite elements, linked by algebraic constraints as proposed in [10]. On the problem formulation side, it can be investigated in which cases globalsupport direct collocation using Fourier approximations can improve efficiency in terms of CPU time and memory usage. Efficient problem formulations and implementations that include stability and robustness considerations would be a useful contribution, in particular for multi-aircraft systems. Finally, efficient algorithms that enable simultaneous trajectory and design optimization with expensive models (e.g., aero-elastic models) could lead to faster and more accurate system design loops.  useful user feedback. The authors would like to thank Reinhart Paelinck for his enthusiasm for and inspiring discussions on multi-aircraft AWE systems. We acknowledge support by the Open Access Publication Fund of the University of Freiburg.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:  Table A1. System parameters.