Airborne wind energy systems are capable of extracting energy from higher wind speeds at higher altitudes. The configuration considered in this paper is based on a tethered kite flown in a pumping orbit. This pumping cycle generates energy by winching out at high tether forces and driving a generator while flying figures-of-eight, or lemniscates, as crosswind pattern. Then, the tether is reeled in while keeping the kite at a neutral position, thus leaving a net amount of generated energy. In order to achieve an economic operation, optimization of pumping cycles is of great interest.

In this paper, first the principles of airborne wind energy will be briefly revisited. The first contribution is a singularity-free model for the tethered kite dynamics in quaternion representation, where the model is derived from first principles. The second contribution is an optimal control formulation and numerical results for complete pumping cycles. Based on the developed model, the setup of the optimal control problem (OCP) is described in detail along with its numerical solution based on the direct multiple shooting method in the CasADi optimization environment. Optimization results for a pumping cycle consisting of six lemniscates show that the approach is capable to find an optimal orbit in a few minutes of computation time. For this optimal orbit, the power output is increased by a factor of two compared to a sophisticated initial guess for the considered test scenario.

Quaternions, optimal control, tethered kite, airborne wind energy, kite power, pumping cycle.

theDOIsuffix \VolumeVV \IssueI \MonthMM \YearYYYY \pagespan1 \ReceiveddateXXXX \ReviseddateXXXX \AccepteddateXXXX \DatepostedXXXX \subjclass[msc2010]49N90

Quaternion-based model for airborne wind energy system]A quaternion-based model for optimal control of the SkySails airborne wind energy system

M. Erhard]Michael Erhard1 G. Horn]Greg Horn M. Diehl]Moritz Diehl

1 Introduction

The idea of airborne wind energy (AWE) is to generate usable power by airborne devices. In contrast to towered wind turbines, airborne wind energy systems are flying, usually connected by a tether to the ground, like kites or tethered balloons. AWE systems exploit the relative velocity between the airmass and the ground and maintain a tension in the tether. They can be connected to a stationary ground station, or to another moving, but non-flying object, like a land or sea vehicle. Power is generated in form of a traction force, e.g. to a moving vehicle, or in form of electricity. The three major reasons why people are interested in airborne wind energy for electricity production are the following:

  • First, like solar, wind power is one of the few renewable energy resources that is in principle large enough to satisfy all of humanity’s energy needs.

  • Second, in contrast to ground-based wind turbines, airborne wind energy devices might be able to reach higher altitudes, tapping into a large and so far unused wind power resource [5]. The winds in higher altitudes are typically stronger and more consistent than those close to the ground, both on- and off-shore.

  • Third, and most important, airborne wind energy systems might need less material investment per unit of usable power than most other renewable energy sources. This high power-to-mass ratio promises to make large scale deployment of the technology possible at comparably low costs.

In order to achieve an economic and competitive operation of airborne wind energy (AWE) setups, optimal control of energy production is the key to success. In this paper, a novel formulation of the dynamics of a specific AWE system is given, and it is shown how it allows one to solve challenging optimal control problems for this system.

The paper is organized as follows: in Section 2, we introduce the main ideas of airborne wind energy, following the lines of [8], and then explain in Section 3 the experimental system developed by the company SkySails, the modeling of which is the main subject of this paper. Starting point is an existing differential equation model of the system from [9], described in Section 4, which suffers from singularities caused by the choice of coordinate system. The main contribution of this paper is presented in Section 5, where a new singularity free model based on quaternions is developed. It is demonstrated by numerical simulation that the old and new model describe identical system dynamics when away from the singularities, while only the new model can be reliably simulated everywhere. It is then shown in Section 6 how the new quaternion-based model can be used to formulate the complex periodic optimal control problem (OCP). In Sect. 7, the discrete formulation of the OCP to be solved in the dynamic optimization environment CasADi [4] is given. Furthermore, resulting complete operational pumping cycles for the SkySails airborne wind energy system are presented. The paper ends with a conclusion in Section 8.

2 Crosswind kite power and pumping cycles

Every hobby kite pilot or kite surfer knows this observation: As soon as a kite is flying fast loops in a crosswind direction the tension in the lines increases significantly. The hobby kite pilots have to compensate the tension strongly with their hands while kite surfers make use of the enormous crosswind power to achieve high speeds and perform spectacular stunts. The reason for this observation is that the aerodynamic lift force of an airfoil increases with the square of the flight velocity, or more exactly, with the apparent airspeed at the wing, which we denote by . More specifically,


where is the density of the air, the airfoil area, and the lift coefficient which depends on the geometry of the airfoil.

Thus, if we fly a kite in crosswind direction with a velocity that is ten times faster than the wind speed , the tension in the line will increase by a factor of hundred in comparison to a kite that is kept at a static position in the sky. The key observation is now that the high speed of the kite can be maintained by the ambient wind flow, and that either the high speed itself or the tether tension can be made useful for harvesting a part of the enormous amount of power that the moving wing can potentially extract from the wind field.

The idea of power generation with tethered airfoils flying fast in a crosswind direction was already in the 1970’s and 1980’s investigated in detail by the American engineer Miles Loyd [19]. He was arguably the first to compute the power generation potential of fast flying tethered wings - a principle that he termed crosswind kite power. Loyd investigated (and also patented) the following idea: an airplane, or kite, is flying on circular trajectories in the sky while being connected to the ground with a strong tether. He described two different ways to make this highly concentrated form of wind power useful for human needs, that he termed lift mode and drag mode: while the lift mode uses the tension in the line to pull a load on the ground, the drag mode uses the high apparent airspeed to drive a small wind turbine on the wing.

This paper is concerned with systems operating in lift mode, which has the advantage over the drag mode that that it does not need high voltage electrical power transmission via the tether. The idea is to directly use the strong tether tension to unroll the tether from a drum, and the rotating drum drives an electric generator. As both the drum and generator can be placed on the ground, we call this concept also ground-based generation or traction power generation. For continuous operation, one has to periodically retract the tether. One does so by changing the flight pattern to one that produces much less lifting force. This allows one to reel in the tether with a much lower energy investment than what was gained in the power production phase. The power production phase is also called reel-out phase, and the retraction phase reel-in or return phase. Loyd coined the term lift mode because one uses mainly the lifting force of the wing. But due to the periodic reel-in and reel-out motion of the tether, this way of ground-based power generation is often also called pumping mode; sometimes even the term Yo-Yo mode was used to describe it.

Figure 1: The basic idea of the power generating device is to only build the wing tips of a gigantic windmill, in form of tethered airfoils, i.e. kites. Fast crosswind motion of the kites creates an enormous traction force [graphics courtesy Boris Houska].
Figure 2: SkySails functional prototype setup for kites of sizes ranging from 20 to 40 m (30 m shown here). A tether in the range 150–400 m transfers the forces from the flying system to the ground based main winch, which is attached to a motor/generator. A specific feature of this setup is the single tether and the airborne control pod located under the kite. The steering actuator in the control pod pulls certain lines in order to steer the kite.

It is interesting to compare crosswind kite power systems with a conventional wind turbine, as done in Fig. 2, which shows a conventional wind turbine superposed with an airborne wind energy system. Seen from this perspective, the idea of AWE is to only build the fastest moving part of a large wind turbine, the tips of the rotor blades, and to replace them by automatically controlled fast flying tethered kites. The motivation for this is the fact that the outer 30% of the blades of a conventional wind turbine provide more than half of the total power, while they are much thinner and lighter than the inner parts of the blades. Roughly speaking, the idea of airborne wind energy systems is to replace the inner parts of the blades, as well as the tower, by a tether.

The power that can be generated with a tethered airfoil operated either in drag or in lift mode had under idealized assumptions been estimated by Loyd [19] to be approximately given by


where is the area of the wing, the lift and the drag coefficients, and the wind speed. Note that the lift-to-drag ratio enters the formula quadratically and is thus an important wing property for crosswind AWE systems. For airplanes, this ratio is also referred to as the gliding number; it describes how much faster a glider without propulsion can move horizontally compared to its vertical sink rate.

Since the first proposal of using tethered wings for energy harvesting by Loyd, academic research and industrial development activities in the field of airborne wind energy (AWE) have been significantly increasing, especially during the last 10 years. The main advantage of the AWE technology is that airborne systems are capable of harvesting energy from higher wind speeds at higher altitudes. In addition, as most setups require less installation efforts compared to conventional wind turbines, e.g. foundations and towers, this technology is a promising candidate to become an additional source of renewable energy in the near future. For an overview on the different geometries which have been developed so far, the reader is referred to [13] and the monograph on AWE [2].

3 The SkySails Power prototype

The geometry under consideration in this paper is based on the SkySails Power prototype, shown in Fig. 2. This setup features an airborne control pod located under the soft kite connected to a ground based winch with motor/generator by a single tether. The functional prototype successfully demonstrated fully automated energy generation by autonomously controlled pumping cycles [12]. {vchfigure} \vchcaptionTrajectory of an experimentally flown pumping cycle [12]. The principle of energy generation is depicted in Fig. 3 and shall be briefly summarized: during the power phase, the kite is flown dynamically in lemniscates (figures of eight). This dynamical crosswind flight leads to high tether forces. By reeling out the tether, electrical energy is produced. At a certain line length, the kite is transferred to the neutral zenith position. At this low-force position, the tether is reeled in (return phase) consuming a certain portion of the previously generated energy. Finally, a reasonable amount of generated net power remains. Hence, the average power output over complete cycles will be the target of the optimization implementation presented in this paper.

The optimal control is based on a dynamical model of the system, and the chosen numerical techniques are similar to the ones used in [16]. In order to allow for optimization of complete pumping cycles, a simple model should be chosen in order to keep the number of optimization variables low. The model developed in this paper is based on a very simple model set up for the SkySails system [9], which describes the basic dynamic properties of the system and has been experimentally validated [11]. The quaternion representation allows for singularity-free equations of motion, and the replacement of trigonometric functions seems to be advantageous for the optimization process.

A further advantage of the considered system is that an experimentally flown trajectory (see e.g. Fig. 3) can be used as initial guess and the optimization result can be directly compared to it in return. It should be noted that the optimization has to be subject to certain geometrical constraints, which will be discussed in this paper. Particularly, topological constraints have to be imposed to preserve the topology of the pumping cycle during the optimization process (the 6 lemniscates given in the initial guess, compare Fig. 3).

4 Setup and reference kite model

In this section, the model from [9], [11] is briefly summarized. The coordinate system is defined in Fig. 3.

Figure 3: Coordinate system definitions of the fixed reference vectors and the body frame vectors , and . The position is determined by the two angles , and the tether length . The wind direction is defined in -direction.

The state vector of the dynamical system consists of the 3-dimensional orientation of the flying system, given by and the tether length .


Note that the tether of the system leads to a relation between orientation and position, which is given in the basis by:


The wind direction is defined in -direction. In order to provide a clear picture of the coordinate system definition, a comparison with the navigation task on earth shall be discussed as illustrated in Fig. 5.

Figure 4: Demonstrative explanation of the used Euler angles by comparing them to angles used for navigation on earth. For that purpose, the earth’s rotation axis has to be tilted by 90 degrees with the north pole towards the wind direction. The angles , correspond to the position on earth and to the bearing. The exact mapping is given in Table 5.

table angle in model ’earth’ navigation quantity earth longitude earth latitude bearing (direction w.r.t. north pole)

Figure 5: Mapping of the Euler angles of the tethered system to earth navigation quantities (compare Fig. 5).

The coordinate system corresponds to an earth coordinate system with symmetry (rotational) axis in wind direction, which could be achieved by rotating the north pole by 90 degrees towards the wind direction. According to the angle mapping given in Table 5, the two angles correspond to longitude and latitude, hence determine the position on earth. The orientation angle corresponds to bearing, i.e. direction w.r.t. north in our earth system. This angle determines the direction of motion (course) and the angle of describes heading directly against the wind. It should be finally noted that due to the ambient wind, the angle does not exactly coincide the kinematic flight direction (course). This effect is fully covered by the developed model and further explained in [12].

The input vector of the model comprises the inputs of the two control actuators: the steering deflection of the airborne control pod, and the winch speed of the tether , respectively.


In the following, the equations of motion shall be summarized for reference purpose only. This model was derived based on the same assumptions as those that will be used in the subsequent section for the derivation of the quaternion representation. As conducted in detail in [9], [11], the equations of motion read:


Here, the airpath speed of the kite, , is given by


The set of equations involves three parameters: lift-to-drag (glide) ratio , steering response proportionality constant and ambient wind speed . Further explanations of parameters including values used later in the paper are summarized in Table 10. {vchtable} \vchcaptionDescription of system parameters and values used for optimization symbol value description 21.0 m projected area of the kite 1.0 aerodynamic force coefficient, given by for aerodynamic lift (drag) coefficient () 0.6 1/s maximal control actuator steering speed 0.7 maximal control actuator steering deflection 5.0 lift-to-drag (glide) ratio. Note, that in contrast to most other systems, the SkySails system is operated at constant glide ratio 0.1 rad/m proportionality constant relating turn rate of the kite to steering actuator deflection 300 m maximal available tether length 0.35 rad minimal elevation angle w.r.t. the horizontal plane 1.2 kg/m air density 5.0 m/s minimal air path speed in order to guarantee flight stability and keep system tethered (this corresponds to free flight velocity) m/s typically chosen as , leads to a windward limit of 10 m/s ambient wind speed Finally, the tether force can be computed based on the air path speed by


with density of air , projected kite area and aerodynamic force coefficient .

5 Quaternion based kite model

In the following, the equations of motion shall be derived in a singularity-free formulation based on quaternions. Quaternions are a well-known tool in aerospace modeling and estimation to avoid the singularities of Euler angles [18]. An alternative singularity-free representation could be based on rotation matrices as natural coordinates directly, as is shown for modeling of AWE systems in [15]. Note that the rotation matrix in the following will be used only as intermediate step in the derivation towards quaternions. The base coordinate system is depicted in Fig. 5a. The time dependence w.r.t. to this reference frame is expressed by rotation matrices , i.e.


The zero-rotation state is defined by , corresponding to , and as illustrated in Fig. 6.

Figure 6: Reference frame axes for the derivation of the quaternion equations of motion. The depicted orientation corresponds to the zero-rotation, i.e.  or .

Note that this seems to be an arbitrary choice a priori. However defining the ambient wind direction, given by , in the kite’s yaw-axis leads to a beneficial simplicity of the equations of motion later.

It should be emphasized that due to the nature of the tethering, position and orientation are related. Hence, the position of the kite is given by


The air path speed vector can be calculated as sum of the ambient wind vector and reversed kinematic velocity vector (apparent wind opposed to motion). Using the time derivative of (15) and introducing the speeds and in the tangent plane, compare Fig. 7, yields


5.1 Model assumptions

For AWE systems, the aerodynamic forces are typically large compared to masses, which in fact is a prerequisite for tethered flight. This is particularly true for the SkySails Power system. Therefore, acceleration effects play a minor role and can be neglected. In addition, the system is assumed to fly always in its aerodynamic equilibrium state. These two simplifications lead to the simple structure of the equations of motion.


Figure 7: The kite is assumed to always fly in an aerodynamic equilibrium state. As aerodynamic model, absence of transversal air flow (top view) and constant lift-to-drag-ratio (side view) is demanded. Figure taken from [11]

The aerodynamics of the system is reduced to two conditions, depicted in Fig. 7 and given in the following.

  1. The kite experiences no transversal wind flow, i.e. no ’side-slip’ angle occurs, compare Fig. 7b.


    Inserting (16) yields


    Thus, the condition results in

  2. The lift-to-drag ratio is constant and given by the parameter . Hence, the kite is assumed to fly always at the same angle of attack. Considering the geometry of Fig. 7c) yields


    Inserting (16) produces


Finally, by considering the motion in the tangent plane, compare Fig. 7, the following relations for the turn rates in kite-body-frame can be derived. {vchfigure} \sidecaption Motion in the tangent plane. Due to the tether, the velocities in the tangent plane and determine and as given in (22) and (23).



Steering of the kite is accomplished by the control actuator in the control pod. Applying an actuator deflection to the steering lines tilts the kite canopy resulting in a curve flight or turn rate around the -axis. The steering behavior can be described by the following equation


This turn-rate law states that the rotation rate is proportional to a constant system parameter , the deflection and the scalar air path speed , which is defined as component of in -direction


The validity of the turn-rate law has been experimentally shown for different kites [9], [14], [17]. In addition, a nice derivation based on geometric considerations can be found in [14]. It should be mentioned that correction terms for mass effects [10] and aerodynamic effects [7] might be added, but are neglected in this paper.

5.2 Quaternions

The final formulation of the equations of motion is based on quaternions [18], introduced as


The time evolution as function of turn rates, which are given in the body frame, reads


Mapping of the kite axes to has to be considered for . Regarding (1214) or Fig. 6 yields , and . The quaternion time evolution then results in


The relation to the rotation can be given by


The derivation of the system’s equations of motion is done as follows: The turn rates in kite-body-frame, , and are expressed as functions of and subsequently inserted in (28). In order to achieve this, (29) is used for and inserted into (22), (23) and (24).

5.3 Equations of motion

The state vector for the quaternion formulation is given by


and the input vector by


Performing the derivation steps of the previous section yields the following set of equations of motion:


with air path speed


It should be remarked that if the winch speed is defined relative to the wind speed by substituting , the product of wind speed and time (or cycle time ) becomes an invariant. As a consequence, solutions for different wind speeds are simply obtained by scaling the dynamic’s time axis proportional to .

Although these equations of motion (EQM) preserve the quaternion norm , a stabilization term in order to compensate for numerical inaccuracies should be added. The EQM for the numerical optimization read then


where denotes the equations of motion (3234) and the damping constant is chosen to lead to a slow decay ( 1/s for the parameters given in Table 10).

5.4 Transformation to Euler angles for comparison

For sake of completeness and in order to allow for a simple geometric interpretation of the results in the Cartesian coordinate system (see Fig. 5), transformations from the set of Euler angles to quaternions and back shall be summarized. The sequence of rotations


is considered in quaternion representation


Calculating the quaternion product () [18] of the single rotations yields


The kite position w.r.t.  can be found by inserting (29) into (15) as


The Euler angles can be computed from the quaternions by


where is the four-quadrant arctangent function.

5.5 Simulation example

In order to complete the presentation of the quaternion-based model, a simulation example will be discussed. The simulations have been carried out using the quaternion-based model on the one and the reference model from Sect. 4 on the other hand. Both computations have employed a Runge-Kutta-4 ODE propagation using timesteps of  s, a tether length of  m and the initial conditions and . A constant steering command of has been used to fly next to the singularity of the Euler angle representation. {vchfigure} \vchcaptionSimulation of the simple model ODE, formulated in Euler angles (solid red) and quaternions (dashed blue). The 3d trajectory is shown on the left, the right plot zooms into the vicinity of the singularity, where the Euler angle formulation fails. The points in the right subfigure represent the discrete simulation timesteps of constant duration. The simulation results are depicted in Fig. 5.5. It can be recognized that the propagation of the model (610) fails due to the singularity, while the quaternion formulation (3234) leads to correct numerical simulation results.

6 Optimal control problem formulation

The optimization task is to control power cycles with optimal power output averaged over complete cycles. This quantity is computed as integral over the mechanical power


It should be mentioned at this point that an open loop optimal control problem is set up, which is based on controls as model input, but does not involve any feedback control based on measured values.

6.1 Setup of the optimization problem

As steering speed will also be limited, the steering deflection is introduced as additional state to the model and the pod steering actuator speed, defined as , is used as control, instead. Further, the state with for computing the integral in (43) is added. Thus, the following relations have to be added to complete the set of EQM (3234)


The control vector now reads


and the state


The NLP problem is formulated as minimization of


The second term introduces a penalty in order to achieve a smooth steering behavior. The weighting factor has to be chosen appropriately.

The optimization, i.e. minimization of , has to be carried out by variation of subject to the following constraints. First the optimization is subject to the model equations of motion (3234). In addition, physical, geometrical and topological constraints have to be added, which will be discussed in the subsequent sections.

6.2 Physical constraints

Due to weight and geometrical design considerations for the airborne system, the two limitations of limited steering speed and deflection of the control pod actuator have to be taken into account. Hence, the steering speed limit is given by


and the limited steering deflection is implemented by


respectively. In addition, winching speed is limited


Note that this bound on the reel-in speed imposes implicitly an constraint on the windward position during the reel-in phase, as can be shown by computing the steady state angle [11],[12]. Further, in order to keep the system tethered, a certain minimum air path speed and hence tether force is required. Recalling (10) yields


Note that this condition is needed to guarantee the validity of the model during the optimization process.

6.3 Geometric constraints

As already introduced, optimization is done for closed pumping orbits, which impose periodic boundary conditions to the optimization problem


We also add one non-periodic boundary condition in order to define the initial value of the “energy state” :


The tether length has to be constrained in order to keep the cycle within a certain range of line length. This is accomplished by setting an upper bound to the state


For a real world system, earth (or water) surfaces are a serious constraint on the equations of motion. Practical safety considerations taking into account flight trajectory deviations due to wind gusts recommend the choice of a certain minimum elevation angle . Geometric considerations yield


Using (42) results in the constraint


It should be finally noted that the operational flight altitude strongly depends on the wind profile and ideally, the optimal flight altitude lies above the given safety limit. However, for the optimal control computations in this paper, the wind field is assumed constant and homogeneous for simplicity and therefore these constraints play a major role.

6.4 Topological constraints

A crucial point in the formulation of the optimization problem is the topology of the trajectory. For operational reasons, basically in order to avoid twisting of the tether, the patterns of choice are lemniscates. Hence, one wishes to preserve the topology of the pattern during the optimization process by imposing certain constraints. If this is not done, the optimizer tries to ’unwind’ certain trajectory features, which in almost all cases leads to a failure of the optimization process.

The topological constraints aim at keeping the configuration of the pattern, e.g. two lemniscates as sketched in Fig. 6.4, or six, as in the actual computations of this paper. {vchfigure} \vchcaptionTopological constraints shown for a pumping cycle consisting of two lemniscates. In order to avoid ’unwinding of the topology’ by the optimizer, the trajectory horizon is divided into stages and additional constraints are imposed on those to preserve the topology. The trajectory can be divided into stages, mainly parts featuring the property flying to the right (left), which corresponds to and , respectively. Formally, an even number of time intervals with duration is introduced. Further, the time points with the property are defined as shown in Fig. 6.4. Note that and that is the overall pumping cycle time.

As is not a state of the system, (7) can be used to formulate an equivalent constraint on . Considering the sign of together with the behavior of the four-quadrant arctangent function in (40), the following conditions can be deduced


The upper (lower) relation corresponds to ().

6.5 Summary of the optimal control problem

The optimal control problem (OCP) can be stated formally as follows:




The set of equations for the ODE model is defined by (2834) and (4445). The boundary conditions are given by the periodicity and initial state constraints (5356), the path constraints by (4952), (57), (59) and the stage wise path constraints by (60).

7 Numerical computation of optimal pumping orbits

In order to numerically solve the optimal control problem, a direct method is used to transfer the continuous time OCP into a nonlinear programming problem (NLP). This is accomplished by using the direct multiple shooting method. The direct multiple shooting method that was originally developed by Bock and Plitt [6] performs first a division of the time horizon into subintervals, that we will call multiple shooting (MS) intervals in the sequel, and it uses piecewise control discretization on these intervals. The continuous time dynamic system is transformed to a discrete time system by using an embedded ordinary differential equation (ODE) solver to solve the ODE on each MS interval individually. Details of the setup of the discrete time control problem will be summarized in the following subsections, and finally numerical results will be presented, that have been obtained by using the optimization environment CasADi [4, 3] and the NLP solver IPOPT [20].

7.1 Multiple shooting and control discretization grid

The choice of time discretization grid size choice is a compromise between keeping the total number of optimization variables low and at the same time reproducing the dynamics of the controls and states in an adequate way. In order to represent controls and states of different timescales, nested grids are introduced as depicted in Fig. 7.1. {vchfigure} \vchcaptionSetup of the discretization grid. The upper level grid is the splitting of the trajectory into stages as introduced in Sect. 6.4. The stage of duration is divided into equidistant multiple shooting (MS) intervals for the states and the winch control. Note that the stages can comprise different numbers of MS intervals in order to account for different phases as power and return phase, respectively. These MS intervals are divided into subintervals for the kite steering actuator speed (control) to allow for a proper consideration of the constraints (49) in the solution. The initial value on teach the control vector on each MS interval consists of components, i.e.


The idea of multiple shooting is to solve the ODE on small intervals with duration starting with initial values The used scheme is sketched in Fig. 63. {vchfigure} \vchcaptionPrinciple of the implemented shooting scheme for . The ODE model is solved on intervals by integration steps, starting at and ending at . The trajectory is ’connected’ by imposing constraints, e.g. . Note that in contrast to the usual implementation of multiple shooting, several piecewise constant control steps are taken per MS interval for one of the controls. The integration is performed from an initial state to the final state , using the controls . Here, we are utilizing integration steps of the classical Runge-Kutta integrator of order 4 (RK4), each with a step size . In order to give a formal representation, the single RK4 integration step shall be defined as


The result of the RK4 integration steps on the interval with duration is then given by


It should be emphasized that separate piecewise constant control steps on the subgrid are provided for the integration steps of the MS interval.

In summary, the optimization vector can be subdivided into subvectors with all variables belonging to stage , which are given by


7.2 Summary of the discretized optimal control problem

The complete discretized OCP can be stated as follows




where next denotes the subsequent index pair in time. Note that the second term penalizes changes in winch speed (i.e. acceleration) and has been added to the discretized OCP only. The discretized path constraints evaluate the constraints (4952), (57), (59) on the multiple shooting grid, and the discretized stage wise path constraints evaluate (60) on the grid.

7.3 Numerical results

In the following, numerical results for optimal complete pumping cycle will be presented. The parameters are given in Table 10 and are related to the SkySails functional prototype of Fig. 2. As initial guess, an experimentally flown trajectory is used. The topology of 6 lemniscates leads to and the (initial) cycle time of  s is divided up into intervals. The number of RK4 integration steps is chosen as on each MS integration interval. Thus, the optimization vector contains a total of 2762 optimization variables. The 3d trajectories at different iteration steps are shown in Fig. 8.

Figure 8: Evolution of pumping cycles during the optimization process. Starting with the experimentally flown trajectory (top left), the 3d trajectories for different iteration steps (number given in figure) to the optimal solution (bottom right) are shown.

It can observed that optimization starts far from optimum and runs through some crude path configurations towards the optimal solution. This behavior could be attributed to the interior point method of the used IPOPT solver.

The optimization figures of merit are given by the Loyd factor , which is defined by


This factor corresponds to the Loyd limit, which corresponds to the generated power for a kite flying ’maximal crosswind’ (at ) continuously, i.e. no retraction phase, compare (2) and [19]. For the derived model, this quantity is given by


The according time series for the controls and power, given as , are presented in Fig. 9.

Figure 9: Time series for controls and and power, computed as tether force times reeling speed . The upper row shows the experimentally flown pumping cycle, which is used initial guess. The lower row gives the optimal solution. Iteration numbers for the intermediate steps are annotated in the respective subplots.

In summary, it can be stated that the optimizer increased the power generation of the pumping cycle by more than a factor of two from to . This is partly due to the fact that the optimized figure lies deeper in the crosswind region. Further, transfer and return phases are carried out more efficiently. The lemniscates are deformed to minimal curves fulfilling the imposed constraints. For sake of a fair assessment, it should be noted that the experimental trajectory was taken from a test flight aiming at flight controller development rather than optimal power output. Nevertheless, valuable suggestions for improvement of even these ’classical’ control setups can be extracted from such optimization results, see e.g. the winch controller implementation for the transfer phase in [12].

8 Conclusion

In conclusion, a singularity free model for tethered kite dynamics based on quaternions has been derived, which can be broadly applied for simulation and optimization purposes. It has been demonstrated that optimization runs of complete pumping cycles within a few minutes are feasible with this model. This opens up the way to further extended and systematic studies of AWE efficiencies. It should be finally remarked that these might most likely involve further extensions to the model to take into account e.g. mass effects, tether drag and a wind profile. {acknowledgement} We gratefully thank Mario Zanon for valuable discussions on the implementation of the OCP. This research was supported by KUL PFV/10/002 OPTEC; Eurostars SMART; IUAP P7 (DYSCO); EU: FP7-TEMPO (MCITN-607957), H2020-ITN AWESCO (642682) and by the ERC Starting Grant HIGHWIND (259166).


  1. Corresponding author, e-mail: michael.erhard@skysails.de, Phone +49 40 70299 203, Fax +49 40 70299 222


  1. U. Ahrens, M. Diehl, and R. Schmehl (eds.), Airborne Wind Energy, Green Energy and Technology (Springer Berlin Heidelberg, 2013).
  2. J. Andersson, A General-Purpose Software Framework for Dynamic Optimization, PhD thesis, Arenberg Doctoral School, KU Leuven, Department of Electrical Engineering (ESAT/SCD) and Optimization in Engineering Center, Kasteelpark Arenberg 10, 3001-Heverlee, Belgium, October 2013.
  3. J. Andersson, J. Ãƒâ€¦kesson, and M. Diehl, Casadi: A symbolic package for automatic differentiation and optimal control, in: Recent Advances in Algorithmic Differentiation, edited by S. Forth, P. Hovland, E. Phipps, J. Utke, and A. Walther, Lecture Notes in Computational Science and Engineering Vol. 87 (Springer Berlin Heidelberg, 2012), pp. 297–307.
  4. C. Archer and K. Caldeira, Global assessment of high-altitude wind power, Energies 2, 307–319 (2009).
  5. H. Bock and K. Plitt, A multiple shooting algorithm for direct solution of optimal control problems, in: Proceedings 9th IFAC World Congress Budapest, (Pergamon Press, 1984), pp. 242–247.
  6. S. Costello, G. François, and D. Bonvin, Real-Time Optimization for Kites, in: Proceedings of the 5th IFAC-PSYCO Workshop, (2013), pp. 64–69.
  7. M. Diehl, Airborne Wind Energy: Basic Concepts and Physical Foundations, in: Airborne Wind Energy, edited by U. Ahrens, M. Diehl, and R. Schmehl (Springer, 2013).
  8. M. Erhard and H. Strauch, Control of towing kites for seagoing vessels, Control Systems Technology, IEEE Transactions on 21(5), 1629–1640 (2013).
  9. M. Erhard and H. Strauch, Sensors and navigation algorithms for flight control of tethered kites, in: Control Conference (ECC), 2013 European, (July 2013), pp. 998–1003.
  10. M. Erhard and H. Strauch, Theory and experimental validation of a simple comprehensible model of tethered kite dynamics used for controller design, in: Airborne Wind Energy, edited by U. Ahrens, M. Diehl, and R. SchmehlGreen Energy and Technology (Springer Berlin Heidelberg, 2013), pp. 141–165.
  11. M. Erhard and H. Strauch, Flight control of tethered kites in autonomous pumping cycles for airborne wind energy, Control Engineering Practice 40(0), 13 – 26 (2015).
  12. L. Fagiano and M. Milanese, Airborne wind energy: An overview, in: American Control Conference (ACC), 2012, (June 2012), pp. 3132–3143.
  13. L. Fagiano, A. Zgraggen, M. Morari, and M. Khammash, Automatic crosswind flight of tethered wings for airborne wind energy: Modeling, control design, and experimental results, Control Systems Technology, IEEE Transactions on PP(99), 1–1 (2013).
  14. S. Gros and M. Diehl, Modeling of airborne wind energy systems in natural coordinates, in: Airborne Wind Energy, edited by U. Ahrens, M. Diehl, and R. SchmehlGreen Energy and Technology (Springer Berlin Heidelberg, 2013), pp. 181–203.
  15. B. Houska and M. Diehl, Optimal Control for Power Generating Kites, in: Proc. 9th European Control Conference, (Kos, Greece,, 2007), pp. 3560–3567, (CD-ROM).
  16. C. Jehle, Automatic flight control of tethered kites for power generation, Master thesis, Technical University of Munich, Germany, 2012.
  17. J. B. Kuipers, Quaternions and rotation Sequences: a Primer with Applications to Orbits, Aerospace, and Virtual Reality (Princeton University Press, 1999).
  18. M. L. Loyd, Crosswind kite power, Journal of energy 4(3), 106–111 (1980).
  19. A. Wächter and L. Biegler, On the Implementation of a Primal-Dual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming, Mathematical Programming 106(1), 25–57 (2006).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description