Optimal Control Algorithms with Adaptive Time-Mesh Refinement for Kite Power Systems

This article addresses the problem of optimizing electrical power generation using kite power systems (KPSs). KPSs are airborne wind energy systems that aim to harvest the power of strong and steady high-altitude winds. With the aim of maximizing the total energy produced in a given time interval, we numerically solve an optimal control problem and thereby obtain trajectories and controls for kites. Efficiently solving these optimal control problems is crucial when the results are used in real-time control schemes, such as model predictive control. For this highly nonlinear problem, we derive continuous-time models—in 2D and 3D—and implement an adaptive time-mesh refinement algorithm. By solving the optimal control problem with such an adaptive refinement strategy, we generate a block-structured adapted mesh which gives results as accurate as those computed using fine mesh, yet with much less computing effort and high savings in memory and computing time.


Introduction
In the last few years we have seen a fast increase in installed wind power capacity. The installed cumulative power grew by 12.6% in 2016, reaching a total of 486.8 GW [1]. Despite the significant and increasing number of wind farm installations (at low heights), most of the existing wind energy is observed at great heights and is still unexplored. Far from the Earth's surface, the wind is more powerful and, more importantly, much more consistent. Nevertheless, the conventional wind turbines mounted on towers cannot harvest such winds. To exploit high-altitude winds, a growing number of research groups in Europe and in the U.S. have recently been proposing and developing airborne wind energy systems (AWES), mainly using kites or rigid wing systems. The huge potential of the research in this topic has already been recognized and the economic viability of the concept is starting to be explored by a few start-up companies, including AmpyxPower (the Netherlands, [2]), KiteGen (Italy, [3]), KitePower (the Netherlands), and Makani Power (now part of Google X, USA, [4]), to name a few. In Portugal, Omnidea [5] has been developing an AWES based on the Magnus effect. Overviews of AWES research and technologies are given in [6][7][8].
Most of the systems being developed are based on exploiting the powerful crosswind kite power described by Loyd in 1980 [9]. These systems use the fact that aerodynamic lift is proportional to the square of the (apparent) wind velocity F lift = 1 /2c L (α)Av 2 a .
Therefore, maximum power extraction occurs when the kite is flying at high speeds in a plane orthogonal to the wind direction (the kite behaves much like the tip of the blade in a turbine). In kite power systems with a generator on the ground, the generation of electrical power is done using a controlled kite (tethered wing) by unwinding a cable coiled around a drum connected to a generator. Energy is generated when the cable is unwinding. However, since the length of the cable is finite, eventually the cable has to be coiled back. Power is consumed when coiling back the cable, but the whole cycle is done in such a way that the balance-the total energy produced-is positive.
We consider the optimal control problem with a continuous-time model of the kite power system, further developing the work reported in [10,11]. The solution to this problem provides the trajectory of the kite and the corresponding control laws that maximize the total energy produced in a given time interval. The efficient solution of optimal control problems is necessary so that the results may be used in a real-time control schemes such as model predictive control (see e.g., [12][13][14][15]). We solve this problem in a time-mesh that is adaptively refined to achieve a desired level of accuracy. The adaptive mesh refinement (AMR) strategy used considers, in the optimal control problem (OCP) context, several refinement criteria in a multi-level scheme [16]. A major feature of the strategy is to consider the local error of the dual variables as a refinement criterion. This error can be computed efficiently by comparing the solution of a linear differential equation system-the adjoint equation of the maximum principle-with the numerically obtained multipliers.
Details of this technique and its use in other nonlinear systems are reported in [16,17]. The refinement strategy resulted in higher accuracy levels with a lower overall computational time when compared with the use of traditional meshes with equidistant spacing in the continuous-time OCP, and with a priori discretized versions of the OCP. In particular, the time to obtain a solution with a similar level of accuracy represented 4% to 30% of the time needed using traditional equidistant meshes. The benefits of using an adaptive time mesh are particularly more evident in highly nonlinear systems, as is the case of the controlled kite and other nonholonomic systems. The technique led to the use of finer meshes in time intervals where sharp turns of the kite occur. A simplified problem in two dimensions was solved with the AMR 29 times faster than with an equidistant mesh leading to the same level of accuracy. The complete 3D problem was successfully solved using the AMR strategy, while when using an equidistant mesh leading to the same level of accuracy, the maximum number of iterations of the nonlinear optimization solver was exceeded.

Kite Power System Model
The AWES considered in the paper is a kite generator system (KGS) involving a controlled kite module (CKM) connected to a ground generator module (GGM) by a single tether (see Figure 1). Among the recent AWES proposals, the KGS type is one of the most simple, yet appears to be one of the most promising. One of its main advantages is having the generator on the ground, reducing not only safety risks but also the building costs since it does not require heavy infrastructures. The CKM comprises the kite and the airborne control unit (ACU) while the GGM comprises the generator unit and the cable drum unit.
There is a single tether going from the GGM to the ACU and the line tension in the described models below refers to the tension in this single tether. From the ACU to the kite, three or more lines are attached. The main line, which takes most of the load, is attached to the leading edge of the kite in a bridle configuration. The other two lines attached to each corner of the trailing edge are used to control the roll angle (ψ) and angle of attack (α). By pulling both of these lines we increase the angle of attack and a differential pull in these lines controls the roll of the kite.

2D Kite Model
We start by addressing a simplified 2D model of a kite power system moving in a vertical plane containing the wind vector. We consider two coordinate systems: Global G: A Cartesian coordinate system (x, z) where x is aligned according to the wind direction v w = (v w , 0) − basis ( e x , e z ), Local L: A polar coordinate system (r, β) centred at the kite position-basis ( e r , e β ). and, therefore, the rotation matrix from L coordinate system to G is R LG = e r e β = cos(β) −sin(β) sin(β) cos(β) and the rotation matrix from G coordinate system to L is R GL = R −1 LG = R LG .

Acting Forces
Considering the mass of the kite (m) and the total force acting on it (F), Newton's law reads In (1) F can be decomposed into tether tension at the ground ( F th ), gravity ( F grav ), and aerodynamic force ( F aer (α)), i.e., The tether force is completely aligned with the first coordinate of the local coordinate system while the kite gravitational force is completely aligned with the last coordinate of the global coordinate system (see Figure 2). We consider the case in which the drag force is aligned with the apparent wind velocity and the lift force has its upward normal direction (see Figure 2).
where R 90 is the 90 • anticlockwise rotation matrix. Thus In a polar coordinate systemṗ Considering the aerodynamic force directions defined by the apparent wind, we can write and therefore, Consideringṙ = v t ∈ R, we haver = 0 and therefore we can define the state x = r, β,β , the control u = (v t , α), and the dynamic equatioṅ

3D Model
We model the forces acting on the kite in a spherical coordinate system positioned at the center of the kite, similarly to [18,19]. We consider three coordinate systems: Global G: An inertial Cartesian coordinate system (x, y, z) where the origin is on the ground at the point of attachment of the tether and x is aligned according to the wind direction v w = (v w , 0, 0), on the basis of ( e x , e y , e z ). We consider that the kite is positioned in a point p with coordinates (x, y, z). Local L: A non-inertial spherical coordinate system (r, φ, β) on the basis of ( e r , e φ , e β ) (Figure 3a). Body B: A non-inertial Cartesian coordinate system attached to the kite body on the basis of ( e 1 , e 2 , e 3 ). e 1 coincides with the kite's longitudinal axis pointing forward, e 2 in the kite transversal axis points to the left wing tip, and e 3 in the kite vertical axis points upwards (Figure 3b).
Considering the position the rotation matrix from L coordinate system to G is and the rotation matrix from G coordinate system to L is  We consider the apparent wind velocity v a = v w −ṗ. Assume that its radial component v a,r is always strictly positive and that the kite body is at all times positioned in such a way that its longitudinal axis is aligned with the apparent wind velocity, that is e 1 = −v a / v a . Let ψ be the roll angle measuring rotation around the e 1 axis. We consider thatẽ 2 = e 2 is initially(for ψ = 0) in the plane τ, tangent to a sphere centred at the origin (containing the axis e φ and e β ). We have that e 2 ⊥ e r , andẽ 2 ⊥ e 1 . We can then defineẽ 2 = e r × e 1 e r × e 1 . Finally, we consider the kite body has an anti-clockwise rotation of ψ around the e 1 axis: the roll angle. We assume that the roll angle ψ can be controlled directly. In fact, since the kite has some mass, the roll angle cannot be selected arbitrarily at each instant. We would have to control the angular acceleration and thereby the angular velocity and roll angle. However, since the rotation of the kite is much faster than its translational movement in the defined operational range, we can consider, as a simplification, that the roll angle is a directly actuated control variable. For example, in a two-line kite where d is the distance between attachment points and ∆r is the relative difference between the lengths of each line, we have sin ψ = ∆r/d (cf. [19]). Using Rodrigues' formula to rotateẽ 2 by ψ around e 1 , we obtain e 2 =ẽ 2 cos ψ + ( e 1 ×ẽ 2 ) sin ψ + e 1 ( e 1 ·ẽ 2 )(1 − cos ψ) (6) and finally, we define e 3 forming a right-handed coordinate system e 3 = e 1 × e 2 .
As before, the total force acting on it can be decomposed as where In the local coordinate system (8) where the second term is − 1 /m F inert with F inert representing the inertial forces (centrifugal and Coriolis) in the local coordinate system. Now, we can write We assume that the tether accelerationr can be controlled directly by a t . Denoting by T the tension on the tether at the base, we have T = F r − ma t . Defining the state x = r, φ, β,ṙ,φ,β and the control u = (a t , α, ψ), the dynamic equation iṡ

Optimal Control Problem
We consider the problem of optimizing power production in a given time period [t 0 , t f ] (see e.g., [20][21][22] for a reference on optimal control and on the corresponding numerical methods). The instant power production is given by P(t) =ṙT and the energy in the interval is We consider two problems: (P 1 ): The production cycle, when the terminal state is free; (P 2 ): The reel-out and reel-in cycle, by imposing that the terminal state should be near the initial one: Considering t ∈ [t 0 , t f ], the problem (P 1 ) can be stated as: T dt (12) subject to dynamic constraintṡ the left end-point constraint and bounded-state constraints In (P 2 ), in order to concentrate the production in the beginning of the time interval, we change the objective function to e −tṙ T , and thus, this problem can be stated as follows: Optimization-based control has been used in AWES with two main purposes. On the one hand, we can address the problem of planning the best path or trajectories that should be aimed to maximize energy production, which can later be implemented using some other control method. This is an open-loop (OL) optimal control problem. On the other hand, we can use optimal control within a model predictive control (MPC) framework for real-time control of the kite (i.e., devising the actuator values in real-time taking into account current measures of the state of the system). This has been shown to be a successful method to address this challenging problem, e.g., [18]. This is ultimately a closed-loop (CL) control problem.
The two problems are intimately related. A sequence of OL optimal control problems might be used to construct the CL feedback control law, which forms the concept of MPC. In addition, the CL control problem would be to track the optimal trajectory or path obtained by the OL planning problem.
We are directly interested in solving the OL planning problem. The performance as well as the robustness of the optimizer is certainly important when solving such OL problems. However, such performance and robustness become crucial when the OL optimal control problems have to be solved within an MPC framework to generate adequate feedback law to command the actuators in real time. The optimizer has to guarantee that a feasible solution is provided (ideally an optimal one) and that such solution is obtained within the real-time window frame constraint.
Efficient methods to obtain a solution to the OL optimal control problems are discussed in the next section.

A Multi-Level Adaptive Mesh Refinement Algorithm
The solution of nonlinear optimal control problems such as the ones above can only be obtained numerically. Numerical methods to solve these problems include dynamic programming-based algorithms (solving the Hamilton-Jacobi-Bellman equation), indirect methods (applying the maximum principle to convert the problem into a boundary value problem) and direct methods (discretizing and solving the resulting static nonlinear programming problem).
Over the past decades, direct methods have become increasingly useful when computing the numerical solution of nonlinear optimal control problems (OCPs) [21]. These methods directly optimize the discretized OCP without using the maximum principle and they are known to provide a robust approach for a wide variety of optimal control problems, namely when regularity of the minimizing trajectories can be established [23]. When applying direct collocation methods, the control and the state are discretized in an appropriately chosen mesh of the time interval. Then, the continuous-time OCP is transcribed into a finite-dimensional nonlinear programming problem (NLP) which can be solved using widely available software [24]. Most frequently, in the discretization procedure, regular time meshes having equidistant spacing are used. However, in some cases, these meshes are not the most adequate to deal with nonlinear behaviours. One way to improve the accuracy of the results while maintaining reasonable computational time and memory requirements is to generate a time-mesh with different time steps. The best location for the smaller steps sizes is, in general, not known a priori, so the mesh is refined iteratively. In a mesh-refinement procedure the problem is solved, typically, in an initial coarse uniform mesh in order to capture the basic structure of the solution and of the error. Then, this initial mesh is repeatedly refined according to a chosen strategy until some stopping criterion is attained. Several mesh-refinement methods employing direct collocation methods have been described in recent years [21,[25][26][27].
As stated in [16], the multi-level adaptive mesh refinement process starts with discretizing the time interval [t 0 , t f ] in a coarse mesh used to solve the NLP problem associated to the OCP in order to catch the main structure of the solution. The subintervals of time where the desired accuracy is not achieved are refined taking into account different levels of refinement in a single iteration, i.e., they are divided into smaller subintervals according to user-defined levels of refinement. This procedure adds more node points to the subintervals in higher levels of refinement, corresponding to higher errors, and it adds less node points to those in lower refinement levels ( Figure 4).
In order to proceed with the mesh refinement strategy, we have to define a refinement criterion and a stopping criterion. As refinement criterion we consider the estimate of the relative error of the adjoint multipliers (dual variables). At each refinement step, the local error of the multipliers ε q is evaluated for each time instant t and the procedure selects which time intervals should be further refined: where q MP are the multipliers which are solution of the differential equation system given by the maximum principle [20,28], and q KKT are the multipliers obtained by applying the Kuhn-Tucker conditions to nonlinear optimization problem which results from the transcription of the optimal control. This criterion is chosen because these multipliers give sensitivity information. Furthermore, q MP are solutions to a linear differential equation system, which can be easily solved in a faster way and with higher accuracy.
For the stopping criterion, we consider a threshold for the relative error of the trajectory. At each refinement iteration, the the L ∞ norm of the relative error of the primal variables (ε x ) is computed and this information is taken into account when deciding if the refinement procedure should continue. We do not need to compute ε x (t) for all t-as in the case of the refinement criteria-but just an estimate of ε which is much faster to obtain.
Since the proposed procedure increases the number of nodes, a greater computational time would be expected. To decrease the CPU time, when going from a coarse mesh to a refined one progressively, the previous solution is used as a warm start for the next iteration. This action proved to be vital in the decreasing of the overall computational time.
This algorithm was developed and tested in [16] for solving a problem with a nonholonomic system and the reported results show that the proposed strategy can solve the problems, spending only 4% of the time required by the standard methods, and therefore can be used in real-time optimization-based control schemes.
The results obtained by applying the AMR procedure to kite power systems are described in the next section.

Numerical Results
The proposed algorithm was implemented in MATLAB R2016b combined with the Imperial College London Optimal Control Software-ICLOCS-version 0.1b [29].
We consider the simulation parameters of the kite system defined in Table 1 with the aerodynamic coefficients (adapted from [18]).

2D Problem Results
We solve the simpler problem restricted to two dimensions (φ = 0) and we impose that the terminal state should be near the initial one. We consider the state x = r, β,β , the control u = (v t , α), the dynamic equationẋ and The problem (P 2 ) in 2D can be stated as: This problem (P 2 ) is solved with an adaptive mesh refinement strategy [16] with two meshes: π ML : The mesh generated by the adaptive refinement strategy π F : The equidistant spacing mesh considering the lowest ∆t of π ML The results can be seen in Figures 5 and 6a,b, and the efficiency of the optimal control algorithms is reported in Table 2 where we can see that the use of the adaptive mesh refinement strategy reduces the computational time by about 29 times.

3D Problem Results
Considering t ∈ [0, 20], the problem (P 1 ) in 3D can be stated as follows: 20] x(0) = x 0 = (10, 0, 80 • , 3, 0, 0) 20] This problem is also solved with the adaptive mesh refinement strategy (AMR). The results can be seen in Figures 7-10. We tried to solve the same problem using an equidistant-spaced time-mesh in order to compare the results. However, in such a mesh we could not make the optimizer converge within a reasonable limit for the maximum number of iterations (3000).

Conclusions
We address an optimal control problem of generating electricity through kite power systems. The obtained solution comprises the trajectories and controls for the kite that maximize the total energy produced in a given interval. We report results that show that with an adaptive mesh refinement (AMR) strategy the problem can be solved with high level of accuracy and (in simplified versions) much faster. A simplified problem in two dimensions was solved with the AMR 29 times faster than with an equidistant mesh leading to the same level of accuracy.
The complete 3D problem was successfully solved using the AMR strategy, while when using the equidistant mesh leading to the same level of accuracy, the maximum number of iterations of the nonlinear optimization solver was exceeded.
The use of KPSs to generate electrical power at an utility level will require wind farms with multiple kites sharing the same space. A collision avoidance mechanism would then have to be devised. The optimal control formulation is an adequate framework to deal with the trajectory constraints resulting from the collision avoidance restrictions. See [30] for further details.
Author Contributions: All authors contributed in all topics of this article.

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