Mario Zanon

Greg Horn

Sebastien Gros

Moritz Diehl

Sébastien Gros

Control of Dual-Airfoil Airborne Wind Energy systems based on nonlinear MPC and MHE

Keywords: airborne wind energy, fast NMPC and MHE, trajectory tracking, dual airfoils

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Over the last years, conventional wind turbines have grown in size and mass up to a scale at which the major challenges are posed by the structural loads [START_REF] Laks | Control of Wind Turbines: Past, Present, and Future[END_REF], [START_REF] Bossanyi | Further Load Reductions with Individual Pitch Control[END_REF]. The main idea behind Airborne Wind Energy (AWE) is to eliminate all the elements of the system which are not essential for power extraction, resulting in a much lighter structure that only involves an airfoil tethered to the ground. In this configuration, higher altitudes can be reached and the swept area is not fixed by the structure of the system, but can be optimized so as to maximise the extracted power. The system is thus free to operate in previously inaccessible portions of the wind field, where higher wind resources can be found.

The potential of this technology has been established in theory [START_REF] Loyd | Crosswind Kite Power[END_REF]. Several architectures have been proposed for AWE systems, ranging from a single tethered airfoil to complex multi-airfoil structures. In [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF], the power extracted by a single tethered airfoil was compared to the one extracted by two airfoils attached to a Y-shaped tether. The study showed the potential of the dual airfoil configuration, which was able to extract more power than the single airfoil for all considered airfoil dimensions.

Tethered airfoils are strongly nonlinear, constrained unstable systems subject to strong perturbations (e.g. wind M. Zanon, G. Horn, S. Gros and M. Diehl are with the Optimization in Engineering Center (OPTEC), K.U. Leuven, Kasteelpark Arenberg 10, B-3001 Leuven-Heverlee, Belgium. mario.zanon@esat.kuleuven.be * This research was supported by Research Council KUL: PFV/10/002 Optimization in Engineering Center OPTEC, GOA/10/09 MaNet and GOA/10/11 Global real-time optimal control of autonomous robots and mechatronic systems. Flemish Government: IOF / KP / SCORES4CHEM, FWO: PhD/postdoc grants and projects: G.0320.08 (convex MPC), G.0377.09 (Mechatronics MPC); IWT: PhD Grants, projects: SBO LeCoPro; Belgian Federal Science Policy Office: IUAP P7 (DYSCO, Dynamical systems, control and optimization, 2012-2017); EU: FP7-EMBOCON (ICT-248940), FP7-SADCO ( MC ITN-264735), ERC ST HIGHWIND (259 166), Eurostars SMART, ACCM. gusts). Nonlinear model predictive control (NMPC) is an advanced optimization-based control technique able to deal with nonlinear systems while satisfying the given constraints (e.g. avoid stall). Moving horizon estimation (MHE) is an optimization-based observer able to take into account the full nonlinear dynamics of the system. For single airfoils, control schemes based on NMPC have been proposed in [START_REF] Diehl | Efficient nmpc of unstable periodic systems using approximate infinite horizon closed loop costing[END_REF], [START_REF] Houska | Optimal Control for Power Generating Kites[END_REF] and based on both NMPC and MHE have been proposed in [START_REF] Gros | Orbit Control for a Power Generating Airfoil Based on Nonlinear MPC[END_REF], [START_REF] Gros | Control of Airborne Wind Energy Systems Based on Nonlinear Model Predictive Control & Moving Horizon Estimation[END_REF], [START_REF] Zanon | Rotational Start-up of Tethered Airplanes Based on Nonlinear MPC and MHE[END_REF].

To the authors' knowledge, no control scheme has been proposed so far for dual airfoils attached to a Y-shaped tether. In this paper, a MHE and NMPC based control scheme for dual airfoils is presented. Each airfoil is modeled as a rigidbody, 6-DOF object interacting with the air mass. The poweroptimal periodic reference trajectory is computed using the method proposed in [START_REF] Gros | A Relaxation Strategy for the Optimization of Airborne Wind Energy Systems[END_REF].

When using general purpose solvers, the computational burden of solving the MHE and NMPC dynamic optimization problem can be excessive for real-time applications, especially when dealing with fast mechanical systems such as tetherd airfoils. In order to address this issue, the realtime iteration (RTI) scheme has been proposed in [START_REF] Diehl | Real-time optimization and Nonlinear Model Predictive Control of Processes governed by differential-algebraic equations[END_REF]. This scheme was successfully implemented in an extension to the open-source software ACADO that exports tailored, efficient C code [START_REF] Houska | An Auto-Generated Real-Time Iteration Algorithm for Nonlinear MPC in the Microsecond Range[END_REF].

This paper is organized as follows. The process model is presented in Section II and the control scheme based on MHE and NMPC is proposed in Section III. Simulation results are presented in Section IV. Future developments and conclusions are proposed in Section V.

II. SYSTEM MODEL

A. Airfoil Kinematics

In a similar way as in [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF], the system is modeled as a tree structure with three nodes coinciding with the joint position (node i = 0) and the airfoil center of mass (nodes i = 1, 2). Each airfoil is considered as a rigid body having 6 degrees of freedom (DOF). An orthonormal right-hand reference frame E is chosen s.t. a) the wind is blowing in the E 1 -direction, b) the vector E 3 is opposed to the gravitational acceleration vector g.

A right-hand orthonormal reference frame e i is attached to each airfoil node (i = 1, 2) s.t. a) the basis vector e x i spans the wing longitudinal axis, pointing in the forward direction and is aligned with the wing chord, b) the basis vector e z i spans the vertical axis, pointing in the upward direction. Each airfoil attitude is given by the rotation matrix

R i = e x i e y i e z i ,
where vectors e x,y,z i are given in E. The position of node i in the reference frame E is given by r i = [x i , y i , z i ] T . Each tether segment i is approximated as a rigid link of length l i that constrains r i to be on the manifold defined by

c 0 = 1 2 ( r 0 2 2 -l 0 2 ) = 0, c i = 1 2 ( r i + R i r T -r 0 2 2 -l i 2 ) = 0, i = 1, 2,
where r T is the position of the tether attachment point in e i .

Because the center of mass of the airfoils does not coincide with the tether attachment point, translational and rotational dynamics are not separable. Using Lagrange mechanics and the methods described in [START_REF] Gros | Nonlinear MPC and MHE for Mechanical Multi-Body Systems with Application to Fast Tethered Airplanes[END_REF], [START_REF] Gros | Airborne Wind Energy, ch. Modeling of Airborne Wind Energy Systems in Natural Coordinates[END_REF] and introducing the algebraic states λ = [ λ 0 , λ 1 , λ 2 ] T , the equations of motion can be described as index-1 differential algebraic equations (DAE):

M G G T 0   r ω λ   = F, r =   r0 r1 r2   , ω = ω1 ω2 , a i = [ a a i a e i a r i κ i ] T , ȧi = u i , i = 1, 2 R i (0) T R i (0) -I 3 = 0, Ṙi = R i (ω i ) × , i = 1, 2 c i (0) = 0, ċi (0) = 0, i = 1, 2, 3 (1)
where I 3 is the 3 × 3 identity matrix, (ω i ) × ∈ SO(3) is the skew symmetric matrix of the rotational velocity ω i and

M =       ∑ 3 i=0 ξ i 1 2 ξ 1 1 2 ξ 2 0 0 1 2 ξ 1 ξ 1 + m 1 I 3 0 0 0 1 2 ξ 2 0 ξ 2 + m 2 I 3 0 0 0 0 0 J 1 0 0 0 0 0 J 2       , G =       r 0 ∇ r 0 c 1 ∇ r 0 c 2 0 ∇ r 1 c 1 0 0 0 ∇ r 2 c 2 0 2P R 1 (∇ R 1 c 1 ) 0 0 0 2P R 2 (∇ R 2 c 2 )       , F =             F 0 -1 2 gµ 0 l 0 I 3 -∑ 2 i=1 1 2 gµ i l i I 3 F 1 -1 2 gµ 1 l 1 I 3 -gm 1 I 3 F 2 -1 2 gµ 2 l 2 I 3 -gm 2 I 3 T 1 -ω 1 × J 1 ω 1 T 2 -ω 2 × J 2 ω 2 ∇ r 0 ċT 0 ṙ0 -∇ r 0 ċT 1 ṙ0 -∇ r 1 ċT 1 ṙ1 -2P R 1 (∇ R 1 ċ1 ) T ω 1 -∇ r 0 ċT 2 ṙ0 -∇ r 2 ċT 2 ṙ2 -2P R 2 (∇ R 2 ċ2 ) T ω 2            
, where ξ i = 1 3 µ i l i I 3 , µ i denotes the density of tether segment i, m i and J i are respectively the mass and rotational inertia of airfoil i. The sum of the forces applied to node i is denoted by F i . Similarly, the sum of the torques applied to airfoil i is denoted by T i . The projcetion opertation P R (•) is defined as P R (A) = U R T A , and U is the unskew operator

U     a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33     = 1 2   a 32 -a 23 a 13 -a 31 a 21 -a 12   , U(a × ) = a.
Note that, with this formulation, the tether tension is readily given by Γ i = λ i l i .

B. Airfoil Aerodynamic Forces

Introducing the relative velocity v i , i.e. the velocity of the airfoil w.r.t. the air mass given by:

v i = ṙi -w i ,
where

w i = [w x i , w y i , w z i ] T ∈ W ⊂ R 3
is the wind velocity field at the location of airfoil i. The relative velocity in the airfoil reference frame is given by V i = R T i v i . The angle of attack α i and sideslip angle

β i α i = -tan V z i V x i , β i = tan V y i V x i .
The aerodynamic forces and torques are given by

F A i = 1 2 ρA v i C L i v i × e y i -C D i v i , T A i = 1 2 ρA v i 2 C R i C P i C Y i T ,
where the aerodynamic coefficients of lift

C L i , drag C D i , roll C R i , pitch C P i and yaw C Y i are functions of α i , β i
, ω i and of the control surfaces a a i , a e i , a r i [START_REF] Pamadi | Stability, Dynamics, and Control of Airplanes[END_REF]. The drag due to the onboard turbines is given by

F G i = -κ i v i v i ,
where κi = u κ i and u κ i is a control variable. This model assumes that the generated force is opposed to the local velocity.

C. Tether Drag

To compute the tether drag, let's consider tether i connecting nodes i and k. The tether drag model assumes that an infinitesimal portion l i dσ of tether i located at position

r ki (σ ) = σ r i + (1 -σ )r k , σ ∈ [0, 1]
generates an infinitesimal drag force dF ki given by:

dF ki (σ ) = F ki (σ )dσ = - 1 2 ρd k l k C D T v T v T dσ v T = σ ṙi + (1 -σ ) ṙk -[W (z) 0 0] T z = σ z i + (1 -σ ) z k ,
where C D T is the drag coefficient of a cylinder. An infinitesimal displacement of point r ki (σ ) is given by δ r ki (σ ) = σ δ r i + (1σ )δ r k , therefore the contribution of the drag of tether k to the generalized force acting on node i is given by:

F T k,i = 1 0 σ F ki (σ )dσ , (2)
and conversely, the contribution of the drag of tether k to the generalized force acting on node k is given by:

F T k,k = 1 0 (1 -σ )F ki (σ )dσ . (3)
It should be observed that this formulation accounts for the total forces and moments acting on the system due to tether drags. Note that integrals (2) and (3) have to be evaluated numerically. In this paper, a trapezoidal quadrature is used with M = 6 points. The contribution of the tether drags to the generalized forces acting on the generalized coordinates r i is given by:

F T i = F T i,i + ∑ k∈N(i) F T k,i ,
where N(i) associates to each node i the set of nodes connected to it through a tether.

D. Wind Turbulence Model

Similarly to [START_REF] Gros | Control of Airborne Wind Energy Systems Based on Nonlinear Model Predictive Control & Moving Horizon Estimation[END_REF], [START_REF] Zanon | Airborne Wind Energy, ch. Control of Rigid-Airfoil Airborne Wind Energy Systems[END_REF], it is assumed that the wind field at the position of each airfoil i is given by

w i = w t i + w i , i = 1, 2,
i.e. it is the superposition of a turbulent windfield w t and a laminar logarithmic windshear model blowing in the xdirection

w i (z) = w 0 log(z i /z r ) log(z 0 /z r ) 0 0 , i = 1, 2,
where w 0 denotes respectively the wind velocity at height z 0 and z r denotes the terrain roughness.

For control and estimation purposes, it is proposed here to use a rather simple turbulence model, given by ẇt i = -

w t i τ + u t i , i = 1, 2,
where u t i is the forcing term in this first-order differential equation, here modeled as a pseudo-control. While in the estimation problem the optimizer is free to choose the term u t i so as to best fit the measurements, in the control problem, these modes are uncontrollable and, for the prediction horizon, the forcing term is assumed to be 0.

While elaborate stochastic wind models exist in the literature, no turbulence model has been developed specifically for AWE systems.

In the following, the process dynamics and the process initial conditions will be put in the form:

M(X) Ẋ λ = f (X, U) , C (X(0)) = 0,
where f and M and C lump together the process dynamics and consistency conditions given by (1). The control vector is given by U = [ u 1 , u 2 , u t 1 , u t 2 ] ∈ R 14 and the state vector is given by X

= [ r 0 , ṙ0 , X 1 , X 2 ] ∈ R 56 and X i = [ r i , ṙi , R i , ω i , a i , w i ] ∈ R 25 .
C L i Lift coefficient - [-1,1] Γ i Tether tension N [0,∞) Variable Weight r i 8e-8 ṙi 1.5e-6 R i 6e-8 ω i 6e-8 a a i , a e i , a r i 0.1 κ i 9.4e-6 u a i , u e i , u r i 9.4e-6 u κ i 9.4e-8

E. Constraints

The control surfaces deflections a a i , a e i , a r i and their u a i , u e i , u r i rates are constrained by the mechanical design of the airfoil and the chosen actuators.

In addition to the previous constraints, in order to keep the system in the region where the model assumptions are valid, further path constraints need to be added on the lift coefficients C L i to avoid stall and on the tether tension Γ i , to make sure that the tethers are never in compression.

All constraints are summarized in Table II. In the following, all path constraints are lumped together as the inequality constraint function h (X, U) ≤ 0.

F. Available Sensors

The airfoils are assumed to be both equipped with GPS, an inertial measurement unit (IMU), a pitot tube, a variometer, an air probe and control surface encoders. The IMU onboard each airfoil i measures the linear accelerations qi = R T i ri + [ 0 0 g ] T and the rotational velocities ω i , given in frame e i .

All available sensors are listed in Table III with the associated standard deviation σ . Because of actuator noise and inaccuracy, the control inputs computed by the controller, may not be perfectly implemented by the system and are also subject to noise. They are thus added as pseudomeasurements to the cost function. In the following, all measurements will be lumped together in the measurement function y(X, U).

III. CONTROLLER SYNTHESIS

A. Periodic Power Optimal Reference Trajectory

In order to extract the maximum amount of energy from the airmass, the reference trajectory is computed by solving a periodic optimal control problem (OCP) for a given reference windspeed w 0 . The extracted power is given by

P = 2 ∑ i=1 F G i v i .
In order to let the optimizer find the best trajectory, also the period T of the orbit, the tether lengths l i and the tether diameters d i are treated as optimization variables and lumped

together in θ = [ T l 0 l 1 l 2 d 0 d 1 d 2 ].
The resulting periodic OCP is given by minimize

U(•),X(•),λ (•),θ 1 T T 0 P dt, subject to M(X, θ ) Ẋ λ = f (X, U, θ ) , X(0) = X(T ), h(X, U, θ ) ≤ 0, (4)
where the dependence on θ has been made explicit. The system being highly nonlinear, solving the periodic OCP ( 4) is nontrivial and a good initial guess is required in order to be able to solve it. In this paper, this difficulty has been addressed by using the relaxation technique proposed in [START_REF] Gros | A Relaxation Strategy for the Optimization of Airborne Wind Energy Systems[END_REF] and successfully used in [START_REF] Horn | Airborne Wind Energy[END_REF]. The computation of a suitable initial guess involves, itself, the solution of a simpler periodic OCP and a subsequent homotopy procedure [START_REF] Gros | A Relaxation Strategy for the Optimization of Airborne Wind Energy Systems[END_REF].

B. NMPC Formulation

The receding horizon NMPC scheme is formulated using a least squares (LSQ) function penalizing the deviation of the process control inputs and states from the reference trajectories. The NMPC is based on repeatedly (k = 0, 1, ...) solving the dynamic optimization problem:

min U(•),X(•),λ (•) 1 2 t 0 +T P t 0 ( X -X r Q + U -U r R ) dt, s.t. M(X) Ẋ λ = f (X, U) , X(t 0 ) = X(t 0 ), h(X, U) ≤ 0, ( 5
)
where t k = k T s and T s is the NMPC sampling time, T H the NMPC prediction horizon, X r and U r define the state and control reference computed in ( 4) and Q and R are constant positive-definite weighting matrices. Vector X(t k ) is the process state estimated at time instant t k by solving the MHE problem [START_REF] Houska | Optimal Control for Power Generating Kites[END_REF]. Note that the process state estimate must satisfy the consistency condition C X(t k ) = 0. However, the consistency conditions are enforced by the state estimator, and therefore need not appear in the NMPC formulation.

C. MHE Formulation

The MHE scheme is formulated using a least squares (LSQ) function penalizing the deviation of the process control inputs and outputs from the measurements. The MHE estimate of the state X(t k ) is computed by repeatedly (k = 0, 1, ...) solving the following dynamic optimization problem:

min U(•),X(•),λ (•) 1 2 t 0 t 0 -T E y(X, U) -ȳ 2 Q E dt, s.t. M(X) Ẋ λ = f (X, U) , C(X(t 0 )) = 0, h(X, U) ≤ 0, (6)
where y(X, U) is the system measurement function defined in Section II-F, ȳ the corresponding set of measurements and Q E the corresponding covariance matrix.

D. The Direct Multiple Shooting Method

The system dynamics being unstable, problems ( 5) and ( 6) are best solved using simultaneous approaches [START_REF] Biegler | An overview of simultaneous strategies for dynamic optimization[END_REF] such as multiple shooting or collocation. In this paper, a uniform time discretization based on N elements and such that t 0 < t 1 < . . . < t N has been used in the framework of the Direct Multiple Shooting method [START_REF] Bock | A multiple shooting algorithm for direct solution of optimal control problems[END_REF]. The discrete-time formulation is thus obtained by independently integrating the system over each time interval [t k ,t k+1 ] and the path constraints are evaluated on the selected time grid. The basis functions for the control vector parametrization (CVP) have been chosen piecewise constant. The discretization of both ( 5) and ( 6) yields a least-squares NLP, that can be efficiently solved with the generalized Gauss-Newton method. The dimension of the resulting QP is then reduced by condensing [START_REF] Leineweber | An Efficient Multiple Shooting Based Reduced SQP Strategy for Large-Scale Dynamic Process Optimization (Parts I and II)[END_REF].

E. The Real Time Iterations with Shift

The Real Time Iteration scheme is based on solving only a single full Newton-type iteration per sampling instant. The needed computations for the RTI scheme reduce to the computation of the sensitivities of the problem and the solution of a single QP. In this context, a clever initialization of the algorithms is crucial to guarantee contractivity of the scheme [START_REF] Diehl | Real-time optimization and Nonlinear Model Predictive Control of Processes governed by differential-algebraic equations[END_REF]. Based on the solution at the previous time instant, the initial guess for (5) and ( 6) is obtained by shifting the state and control vectors X and U in time.

The initial value embedding consists in keeping the initial state in the optimization variables. This makes it possible to simulate the system, compute the sensitivities and thus run most of the computations before the current estimated state X(t 0 ) becomes available. Once the estimated state is available, the computation of the new controls can be done in a very short time by solving the QP prepared in the previous phase. See [START_REF] Diehl | Real-time optimization and Nonlinear Model Predictive Control of Processes governed by differential-algebraic equations[END_REF], [START_REF] Houska | An Auto-Generated Real-Time Iteration Algorithm for Nonlinear MPC in the Microsecond Range[END_REF], [START_REF] Diehl | A real-time iteration scheme for nonlinear optimization in optimal feedback control[END_REF] for a detailed description of the RTI scheme.

F. ACADO Code Generation

In order to meet the real-time requirements, the code generation tool of ACADO [START_REF] Houska | An Auto-Generated Real-Time Iteration Algorithm for Nonlinear MPC in the Microsecond Range[END_REF], [START_REF] Ferreau | High-speed moving horizon estimation based on automatic code generation[END_REF] has been used. This tool exports an efficient algorithm based on Direct Multiple Shooting and the RTI scheme. The resulting C code exploits the structure of the specific problem and avoids all irrelevant or redundant computations. The effectiveness of code generation in terms of reduced computational time, has been shown in [START_REF] Houska | An Auto-Generated Real-Time Iteration Algorithm for Nonlinear MPC in the Microsecond Range[END_REF], [START_REF] Ferreau | Model Predictive Control Algorithms for Applications with Millisecond Timescales[END_REF].

IV. SIMULATIONS A. NMPC and MHE Tuning

Both NMPC and MHE have been discretized with a sampling time T s = 125 ms. The NMPC prediction horizon was set to N NMPC = 20 elements and the MHE estimation horizon was set tp N MHE = 10 elements.

The NMPC weighting matrices Q and R were chosen diagonal, with each entry specified in Table II, accordingly with the state it corresponds to. Note that the units of the weights are defined consistently with the variables they correspond to, so as to yield a dimensionless cost.

The MHE matrix Q E was chosen so as to match the variance of the measurement noise, given in Table I. Based on the same variance, Gaussian noise has been added to all measurements in the proposed simulations.

B. Simulation Results

In this Section, the simulation results obtained for the model proposed in Section II and the control algorithm proposed in Section III are presented.

The investigated scenario considers tracking a power optimal periodic reference trajectory in a turbulent wind. The reference trajectory and the simulated trajectory are displayed in Figure 1. It can be seen that the airfoil is well stabilized around the reference even in the presence of strong turbulences.

The wind turbulences estimated by the MHE are displayed in Figure 2, where it can be seen that the estimates closely match the real turbulences. Note that the wind profile used in simulations has not been generated using the simple model described in Section II, but with a Von Karman model. Since is assumes a level flight at constant velocity, the Von Karman model is not the most appropriate turbulence model for AWE systems, and should not be used for certification purposes. Given the lack of wind models specific to AWE systems, the Von Karman wind model is still useful to test the control algorithm. In particular, it allows one to show that the component of a complex turbulent wind which is relevant to the system dynamics can still be well estimated with a simple model that does not rely on specific assumptions on the stochastic properties of the wind.

In particular, as shown in Figure 3, the power spectral density (PSD) of the estimation error is much lower than the one of the original signal at low frequencies. The signal-tonoise ratio (SNR) is thus high at low frequencies. At higher frequencies, where the SNR decreases, also the PSD of the turbulences is low.

For the purpose of illustrating the ability of NMPC to deal with constraints, the bounds on the control surfaces and their rates have been artificially tightened so as to activate them often in the simulation scenario. The resulting trajectories are displayed in Figure 4, together with the bounds.

The simulations were run on a 2.8 GHz CPU and the computational times are consistently below T s = 125 ms, allowing for a real-time implementation.

V. CONCLUSION & FURTHER DEVELOPMENTS

This paper has presented a control scheme for a dual airfoil AWE system in drag mode. The system model is a highly nonlinear index 1 DAE which consists of 56 differential states, 3 algebraic states and 14 controls.

The computation of the reference trajectory involves the solution of an involved periodic optimal control problem. In order to be able to solve this problem, particular care has been taken in the computation of a suitable initial guess.

The control technique, based on Nonlinear Model Predictive Control (NMPC) and Moving Horizon Estimation (MHE), was tested in turbulent wind simulations. The MHE scheme is able to estimate the state based on the fusion of the data coming from all sensors while also estimating the wind turbulences. Based on these estimates, the NMPC scheme stabilizes the system on the periodic reference trajectory while satisfying the imposed constraints. The efficient implementation of the optimization routines results in a scheme which is fast enough for a real-time implementation.

The choice of a tracking cost function guarantees that the system will be stabilized, but is suboptimal in the sense of maximizing the energy extracted from the airmass. Future research will aim at further investigating the stability properties of NMPC schemes based on cost functions that directly maximize the extracted energy.

Dual airfoil AWE systems can also operate in pumping mode. Future work will aim at developing a control scheme also for these systems.

Fig. 1 .

 1 Fig. 1. Dual airfoil trajectory. The airfoil trajectories are displayed in red line and the joint position in thick black line. The reference trajectory is displayed in dashed blue line and the black arrow shows the reference wind direction.

Fig. 2 .

 2 Fig.2. Real wind turbulence in dashed line and turbulence wind speed estimated by the MHE in thick line for the two airfoils. As the estimation error is small, the dashed line is covered by the continuous line.

Fig. 3 .

 3 Fig. 3. Power spectral density of the turbulences and of the estimation error (top graph). Signal-to-noise ratio for the estimated signal (bottom graph).

Fig. 4 .

 4 Fig. 4. Control surfaces and their time rates.

TABLE II SYSTEM

 II CONSTRAINTS AND WEIGHTS.


Mario Zanon
email: mario.zanon@esat.kuleuven.be

Sébastien Gros
email: grosse@chalmers.se

Johan Meyers
email: johan.meyers@mech.kuleuven.be

Moritz Diehl
email: moritz.diehl@esat.kuleuven.be

Airborne Wind Energy: Airfoil-Airmass Interaction

) however, do not model the interaction between the airfoils and the airmass. In this paper, the impact of the airfoil-airmass interaction on the extracted power is studied. As this phenomenon is complex to model, a blade element-momentum approach is proposed and the problem is solved by means of optimal control techniques.

INTRODUCTION

To overcome the major difficulties posed by the growing size and mass of conventional wind turbine generators [START_REF] Laks | Control of Wind Turbines: Past, Present, and Future[END_REF]; [START_REF] Bossanyi | Further Load Reductions with Individual Pitch Control[END_REF], the Airborne Wind Energy (AWE) paradigm proposes to eliminate the structural elements not directly involved in power generation. An emerging consensus recognizes crosswind flight as the most efficient approach to Airborne Wind Energy [START_REF] Loyd | Crosswind Kite Power[END_REF]. Crosswind flight extracts power from the airflow by flying an airfoil tethered to the ground at a high velocity across the wind direction. Power can be generated by (a) performing a cyclical variation of the tether length, together with cyclical variation of the tether tension or (b) by using on-board turbines, transmitting the power to the ground via the tether. In this paper, option (b) is considered, as investigated by e.g. [START_REF] Power | Makani Power Homepage[END_REF].

Because it involves a much lighter structure, a major advantage of power generation based on crosswind flight over conventional wind turbines is that higher altitude can be reached and a larger swept area be achieved, thereby reaching wind resources that cannot be tapped into by conventional wind turbines [START_REF] Diehl | Windenergienutzung mit schnell fliegenden Flugdrachen: eine Herausforderung für die Optimierung und Regelung -Wind Power via Fast Flying Kites: a Challenge for Optimization and Control[END_REF].

This research was supported by Research Council KUL: PFV/10/002 Optimization in Engineering Center OPTEC, GOA/10/09 MaNet and GOA/10/11 Global real-time optimal control of autonomous robots and mechatronic systems. Flemish Government: IOF/KP/SCORES4CHEM, FWO: PhD/postdoc grants and projects: G.0320.08 (convex MPC), G.0377.09 (Mechatronics MPC); IWT: PhD Grants, projects: SBO LeCoPro; Belgian Federal Science Policy Office: IUAP P7 (DYSCO, Dynamical systems, control and optimization, 2012-2017); EU: FP7-EMBOCON (ICT-248940), FP7-SADCO ( MC ITN-264735), ERC ST HIGHWIND (259 166), Eurostars SMART, ACCM.

Main tether

Secondary tether

Secondary tether

Fig. 1. Schematic of a dual-airfoil AWE system (from [START_REF] Payne | Self-Erecting Windmill[END_REF], Fig. 3).

Unfortunately, the drag due to the motion of the tether during crosswind flight has a significant impact on the system performance. To tackle this issue, the dual-airfoil design was first introduced in [START_REF] Payne | Self-Erecting Windmill[END_REF] and later investigated in e.g. [START_REF] Houska | Optimal Control for Power Generating Kites[END_REF]; [START_REF] Podgaets | Flight control and stability of a multiple kites tethered system[END_REF]; [START_REF] Williams | Modeling of Optimal Power Generation using Multiple Kites[END_REF]. The key idea of the dual-airfoil desing is to fly two airfoils connected on a single main tether (see Figure 1) in a balanced manner. As a result, only the shorter secondary tethers move at a high velocity and generate drag, while the motion of the main tether is negligible.

In order to fully exploit the gains of reducing the tether drag, the system and the planned trajectory must be carefully designed.

In [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF], a trajectory design method based on optimal control has been proposed to maximize the extracted power. This results in trajectories that a) balance the forces on the main tether, so as to minimize its motion, b) the generator drag forces maximize the system efficiency, c) the tether lengths are chosen so as to achieve the best trade-off between reaching higher altitude and adding airborne mass, and d) the tether diameters must be selected so as to achieve the best trade-off between reducing the drag and withstanding the forces in the system.

The study proposed in [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF] provides a comparison between single and dual airfoil systems for different airborne areas under the assumption of no interaction between the airfoils and the airmass, i.e. the presence of the airfoils does not reduce the velocity of the airmass. It is well known [START_REF] Manwell | Wind Energy Explained: Theory, Design and Application[END_REF]; Bianchi, Fernando D. and Battista, Hernn de and Mantz, Ricardo J. ( 2007) that this hypothesis does not hold true for conventional wind turbines and the Betz law provides an upper bound to the energy that can be extracted from the airmass. Such analysis does not directly apply to AWE systems, and, to the authors' knowledge, no study is available in the literature, which considers interaction between the airfoil and the airmass. In this paper, the effect of such interaction is taken into account in the optimal control problem by introducing an interaction model which adapts the blade element-momentum theory to AWE systems.

This paper is organized as follows. First a multiple-airfoil model is proposed in Section 2. Section 3 describes the powergeneration optimization problem, the solution approach used to compute power-generating trajectories, and the software used to perform the optimization. Section 4 proposes a comparison between optimal power generation based on single and dual airfoils for different approximations of the airfoil-airmass interaction. Finally, Section 5 concludes the paper and outlines further developments.

SYSTEM MODEL

The system is modeled following the procedure proposed in [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF]. Each airfoil has been modeled as a point mass and the tethers have been modeled with a finite element model. An orthonormal right-handed reference frame e = {e x , e y , e z } attached to the ground is chosen to generate the Cartesian coordinate system defining the positions of the airfoils. The frame e is chosen s.t. a) the wind is blowing in the e x -direction, b) the vector e z is opposed to the gravitational acceleration vector g. The origin of the coordinate system coincides with the attachment point of the main tether to the ground.

System Architecture

The system is described as a set of N nodes n ∈ {0, ..., N} with associated coordinate vectors X n ∈ R 3 . The fixed node X 0 = [0, 0, 0] T stands for the attachment point of the AWE system to the ground. The subset A ⊂ {1, ..., N} of the set of nodes describes the nodes associated to the airfoils. Assuming a tree structure, to each node n ∈ {1, ..., N} a single tether k = n is associated, and the parent node i to which the tether is attached is defined by the map i = T (n). See Figures 2 and3 for an illustration. The system architecture is then defined by the number of nodes N, the set A, and the map T . The proposed formulation allows for tree-like system architectures only.

In the following, the notation X n = [x n , y n , z n ] T of the node coordinate vectors X n is used. The position of the node n is given by P n = x n e x + y n e y + z n e z . Each tether k = 1, ..., N has associated length l k and diameter d k .

Airfoil model

For any node n ∈ {0, ..., N}, we define the velocity relative to the airmass: v n = ( ẋn -W ) e x + ẏn e y + żn e z , where W ∈ R is the local wind velocity in the e x direction. A generalization of this formulation to a 3D wind field is straightforward. If n ∈ A, lift and drag forces act on the node. In the following formulas, C n L and C n D are the lift and drag coefficients of the airfoil, ρ is the air density and S the airfoil surface.

Introducing the roll angle ψ i describing the tilting of the lift force around the axis v n , the lift force can be defined by [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF]:

F n L = 1 2 ρSC n L cos(ψ i ) f n L v n -sin(ψ i )e n T v n 2 , e n T = v n × e n r v n × e n r , f n L = e n T × v n .
By definition F n L is always orthogonal to v n , and lies in the plane defined by {e n r , v n } if ψ i = 0. The airfoil drag force is opposed to the relative velocity, and is readily given by:

F n D = - 1 2 ρSC n D v n v n .
The drag generated by the onboard turbines can be modeled as:

F n G = -κ n v n v n , where κn (t) = u κ n (t) , u κ n (t)
is a control variable and we assume that the generated force is opposed to the relative velocity. The resulting aerodynamic power is:

P n = v T n F n G = -κ n v n 3 .
19th IFAC World Congress Cape Town, South Africa. August 24-29, 2014

The resulting airfoil aerodynamic force acting on airfoil n is given by

F n A = F n L + F n D + F n G .
In this model, it is assumed that the time-derivative of the lift coefficient is directly controlled, and the drag coefficient C n D is approximated by [START_REF] Pamadi | Performance, Stability, Dynamics, and Control of Airplanes[END_REF]; [START_REF] Cook | Flight Dynamics Principles[END_REF]

: C n D = C 0 D + C I D (C n L ) 2
, where C 0 D and C I D are the airfoil drag and induceddrag coefficients respectively.

The kinetic and potential energy functions associated with the airfoil dynamics are:

T n A = 1 2 M A Ẋn 2 , V n A = M A gz n ,
where M A is the airfoil mass, and the Lagrange function for the airfoils reads:

L A = T A -V A , T A = ∑ n∈A T n A , V A = ∑ n∈A V n A .

Wind and atmosphere model

Assuming a laminar wind flow with a logarithmic wind shear model blowing uniformly in the e x -direction, the free-flow windspeed W ∞ at altitude z is given by Manwell, J. F., Mc-Gowan, J. G. and Rogers, A. L. ( 2009):

W ∞ (z) = W 0 log (z/z r ) log (z 0 /z r ) , (1)
where W 0 ∈ R is the wind velocity at altitude z 0 and z r is the ground roughness.

To account for the drop of density with altitude the following atmospheric model is introduced Andrews (2010):

T (z) = T 0 -T L z, P(z) = P 0 1 - T L z T 0 gMa RT L , ρ(z) = PM a RT ,
where T 0 is the sea level standard temperature, T L is the temperature lapse rate, P 0 is the pressure at sea level, M a is the molar air density and R is the universal gas constant.

Tether model

In the proposed formulation, the tethers are modeled with a lumped mass Finite Element Model. For a rigid tether k ∈ {1, ..., N} of length l k , density ρ c , diameter d k , we define N k elements linked by massless rigid links, where link k, j connects elements k, j and k, j + 1. Note that with this notation, the position of the endpoint X k,N k of each tether k coincides with the position X k of node k. In the proposed model, all links have the same length l k, j = l k /N k and each element k, j with 1

< j < N k has mass m k, j = m k /(N k +1), while m k,0 = mk, N k = m k /(2(N k + 1))
. The tether kinetic and potential energy functions read:

T k T = N k ∑ j=1 1 2 m k, j Ẋk, j 2 , V k T = N k ∑ j=1 m k, j gz k, j ,
where m k, j is the mass associated with each element and Ẋk, j and z k, j are respectively its velocity and height.

The tether drag on each tether section k, j is given by:

F k, j S = - 1 2 ρ k, j d k l k, j C T v k, j v k, j ,
where C T is the drag coefficient of a cylinder, l k, j is the length of link k, j and v k, j is the velocity of its midpoint, computed as

v k, j = W z k, j + z k, j+1 2 + Ẋk, j + Ẋk, j+1 2 ,
where W is the windspeed at the midpoint's altitude. The lift generated by the tethers is not considered in this formulation.

The contribution of the tether drags to the generalized forces acting on the generalized coordinates X k, j is given by:

F k, j T = F k, j S + F k, j+1
S 2 .

Generalized forces

The vector of generalized forces

F = F T 1,1 , ..., F T N,N N T
, where F k, j ∈ R 3 is the vector of generalized forces acting on the vector of generalized coordinates X k, j , is resulting from the sum of the various contributions coming from tether drags and airfoil aerodynamic forces. Though this summation can be performed very intuitively, it can be formulated as the following systematic construction. For any k ∈ {1, ..., N}, j ∈ {1, ..., N k } , F k, j is given by:

F k, j =      F k, j T if j ∈ 1, ..., N k -1 F k, j T + ∑ F k c ,0 T if T (k c ) = k and j = N k F k, j T + F k A if k ∈ A and j = N k 2.6 System model
In the following, the generalized coordinate vector X = X T 1,1 , ..., X T N,N N T of the system is used. The system is considered as a set of independent tethers and airfoils, with associated Lagrange functions. The tethers introduce a set of constraints in the system configuration, given by:

G k, j = 1 2 (X k, j+1 -X k, j ) T (X k, j+1 -X k, j ) -l 2 k, j = 0,
(2) for k = 1, ..., N, j = 1, ..., N k . The system Lagrange function reads:

L = T -V -λ T G, T = T A + T T , V = V A + V T ,
where λ ∈ R K is the vector of Lagrange multipliers associated to the constraints G. Using the Lagrange equation de Jalón, J. G. and Bayo, E. (1994

) d dt ∂ L ∂ Ẋ -∂ L ∂ X = F,
it can be verified that the system dynamics are given by the following index-3 DAEs:

T Ẋ Ẋ Ẍ + G T X λ + V X = F, G = 0, ( 3
)
where λ is the DAE algebraic state,

G X = ∂ G ∂ X , T Ẋ Ẋ = ∂ 2 T ∂ Ẋ2 and V X = ∂ V ∂ X .
For any t 0 ∈ R, equation (3) can be reformulated as an index-1 DAE by performing index reduction, which yields G(t) = 0, Ġ(t 0 ) = 0, G(t 0 ) = 0. The resulting equations read (together with the consistency conditions):

T Ẋ Ẋ G T X G X 0 Ẍ λ = F -V X - ∂ ∂ X G X Ẋ Ẋ , (4)
G(t 0 ) = 0, Ġ(t 0 ) = G X Ẋ t=t 0 = 0. ( 5
)
19th IFAC World Congress Cape Town, South Africa. August 24-29, 2014

P + P - w A w w 8 8 -
Fig. 4. Schematic of the airflow through an actuator disk (grey area) of area A. W ∞ is the free flow velocity, W the flow velocity at the actuator disk, W -∞ the velocity after the actuator disk. P + and P -stand for the pressure just before and after the actuator disk. It can be verified that the tension in tether k is readily given by:

Γ k = λ k l k .

Blade-Element-Momentum Theory for AWE Systems

The blade-element-momentum (BEM) theory for conventional wind turbines considers the interaction between a blade element and the airflow. The main assumption is that, given the annulus swept by the considered blade element, the change of momentum of the air which passes through it is due only to the force exerced by the blade element on the airmass Burton, T., Sharpe, D., Jenkins, N. and Bossanyi, E. ( 2001).

In the following, the BEM theory is formulated for AWE systems. Note that, in this case, no analytic solution can be computed, as the annulus swept by the AWE system depends on the trajectory flown by the airfoil. The derivation will thus be included in the optimization problem as a constraint which bounds the maximum power that the system can extract from the air mass by adapting the wind velocity at the swept annulus.

As shown in Figure 4, the wind field velocity W (z) at the swept annulus is less than the free flow velocity W ∞ (z), and depends on the average thrust F D exerted on the air mass by the actuator disk. This effect is conventionally represented by the axial induction factor a, defined through W ≡ (1a)W ∞ , 0 < a < 1.

The change of air mass momentum in the annular control volume (cf. Figure 4) results from the thrust force F D acting on the annulus, such that

F D = ρAW (W ∞ -W -∞ ) , (6)
where W -∞ is the steady-state flow velocity in the far wake of the annulus (see Figure 4). Locally around the annular actuator disk, we can further express (since the velocity is continuous over the disk)

F D = A P + -P -,
(7) such that the difference of pressure on the annulus is given by

P + -P -= ρW (W ∞ -W -∞ ) .
(8) From Bernoulli's principle, the conservation of the air mass energy along streamlines before and after the actuator disk requires that

1 2 ρW 2 + P + = 1 2 ρW 2 ∞ + P 0 , 1 2 ρW 2 + P -= 1 2 ρW 2 -∞ + P 0 ,
where P 0 is the pressure in the free flow, resulting in

P + -P -= 1 2 ρ W 2 ∞ -W 2 -∞ . (9)
By combining ( 8) and ( 9), it can be verified that

W -∞ = 2W -W ∞ = (1 -2a)W ∞ . ( 10
)
Using ( 6) and ( 10), the power extracted from the flow at is then given by:

P D (a) = F D W = 2ρAa(1 -a) 2 W 3 ∞ . ( 11
)
This result is essentially the same as the classical Betz result for wind turbines (albeit there for a circular area instead of an annulus). The Betz limit occurs for a = 1/3, leading to a power coefficient C p ≡ P/(ρW 3 ∞ /2) = 16/27, but this presumes an ideal case with frictionless power extraction. In the context of an AWE system, drag forces on the wing, an the possibility to fly the wing at different radii R and speeds, leads to an optimization problem in which a = 1/3 will not be recovered as an optimum. For an AWE system, P D is also given by summing up all the average power dissipated by the system P E , i.e. the total average power resulting from the drag yielded by the tether, airfoil and turbines:

P D = P E = 1 T p T p 0 ∑ i∈A P i + v T i F i D + N ∑ n=1 v T n F n T dt. ( 12
)
Variable P E depends on the wind field velocity W (z), taken as W (z) = (1a)W ∞ (z) in this paper. As a result, within Betz' framework the implicit equality:

P D (a) = P E ((1 -a)W ∞ ) , (13)
must hold for some value of a. In the framework of optimal control a is treated as a decision variable and the implicit relationship (13) as an equality constraint.

Because of the wind shear, the wind velocity on the annulus swept by the AWE system is not uniform, hence it is proposed here to evaluate (11) using the integral:

P D (a) = 2ρa(1 -a) 2 A W 3 ∞ dA. (14)
It can be verified that for a single airfoil, P D can be computed using:

P D (a) = 2ρa(1 -a) 2 T p 0 z tip ẏtip W ∞ (z tip ) 3 dt, (15)
where the coordinate vector X tip = x tip y tip z tip T is the airfoil outer-tip position in frame e. For airfoil n ∈ A, X n tip is given by:

X n tip = A R 2S sin(ψ i ) f n L v n -1 + cos(ψ i )e n T ,
where A R is the airfoil aspect ratio, and w s = A R /S is the airfoil wingspan. For dual-airfoil systems, both airfoils fly the same trajectories, and the area A can be computed based on the trajectory of any of the two airfoils.

Model assumptions & discussion

The proposed model is based on the following assumptions:

(1) the tethers are modelled with a lumped-mass finite element model (2) the lift forces are orthogonal to the airfoil transversal axis (3) the airfoils have a perfect yaw control, resulting in no sideslip (4) the time-derivatives of the lift coefficient and roll angle are controlled and actuation is instantaneous (5) the airfoil-airmass interaction is modeled with a blade element-momentum approach

The proposed model construction can straightforwardly accommodate a 6-DOF airfoil description (including the rotational dynamics) and more elaborate aerodynamic models. Yet in this paper a simple model was preferred, so as to reduce the complexity of the presentation.

The proposed framework for taking into account the interaction between the airfoils and the airmass involves a rather simple formulation and can fail to accurately capture all aerodynamic phenomena, yet, it allows to introduce a model of the interaction in the optimization problem. The improvement of the airfoil-airmass interaction model is the subject of ongoing research.

OPTIMIZATION PROBLEM

The airfoil trajectories as well as the tether lengths and sections are manipulated so as to maximize the system average power generation over an orbit of free duration T p . The periodicity of the system is guaranteed by satisfying the boundary conditions:

Z T (X(0) -X(T p )) = 0, (16) where matrix Z is introduced to satisfy the Linear Independence Constraint Qualification (LICQ) [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF].

In order to ensure that the tethers are always under tension but that their resistance is never exceeded, the constraints:

γ f s π 4 d 2 k ≥ Γ k (t) = λ k l k ≥ 0, ∀t, k = 1, ..., N, (17)
are imposed, where γ is the tether yield strength and f s the safety factor. Moreover, the following bounds are proposed:

0 ≤ C i L ≤ 1, -5 s -1 ≤ Ċi L ≤ 5 s -1 , -80 • ≤ ψ i L ≤ 80 • , -5 s -1 ≤ ψi L ≤ 5 s -1 , -1000 kg/(ms) ≤ κi ≤ 1000 kg/(ms), ∀t, i ∈ A.
(18) The periodic power optimization problem reads:

P = max U,X,θ ,a,T p 1 T p T p 0 ∑ i∈A P i dt, (19)
s.t. (4) -(5), ( 13), ( 16) -(18) where

U i = { Ċi L , ψi L , κ i }, i ∈ A, θ k = {l k , d k }, k = 1, .
.., N. Note that T p is an optimization variable, thus the duration of the orbit will be adapted by the optimizer so as to maximize the average power. In order to be able to treat this problem, a time transformation is introduced (Diehl, 2001, p. 27).

For the dual airfoils the additional constraint X 2 -X 3 2 ≥ w 2 s ensures that no collision occurs between the airfoils. Such a constraint is though never active for an optimal trajectory.

Solution approach

The Optimal Control Problem ( 19) is large-scale and highly non-convex and therefore requires a good initial guess to be tackled by derivative-based optimization. However, no such guess is readily available. To address this issue, a complex procedure is needed to compute an initial guess for problem (19).

For the dual-airfoil system, solving (19) on a full orbit yields quasi-identical trajectories for the two airfoils, hence ( 19) was solved on a half orbit instead, using the periodicity conditions X 2 (0) = X 3 ( 1 2 T p ), X 3 (0) = X 2 ( 1 2 T p ) so as to match the terminal state of one airfoil with the initial state of the other. For both the single and dual-airfoil problems, the control input profiles were discretized using a piecewise-constant parametrization having 20 intervals per full orbit. One collocation element has been used per control interval.

Methods & Software

Because dynamics (4) are unstable, a simultaneous optimal control technique is required to optimize the system model. In this paper, the discretization of the model dynamics (4) was based on a direct collocation approach [START_REF] Biegler | Nonlinear Programming[END_REF], where the model simulation, constraints and optimization are handled simultaneously in a large-scale sparse Nonlinear Program (NLP). Collocation approaches provide a straightforward way to deal with implicit index-1 DAE systems [START_REF] Biegler | Nonlinear Programming[END_REF].

The problem transcription was performed using the opensource optimization framework CasADi [START_REF] Andersson | CasADi -A symbolic package for automatic differentiation and optimal control[END_REF]. The resulting NLP was solved using the interior-point solver IPOPT 3.10.1 [START_REF] Wächter | On the Implementation of a Primal-Dual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming[END_REF] using WSMP as a linear solver.

NUMERICAL RESULTS

The parametric study aims at assessing what is the impact of the airfoil-airmass interaction on the extracted power by means of optimization. In particular, it is of most interest to assess if including a model of this interaction with the airmass (see Section 2.7) makes the single airfoil more advantageous than the dual.

The study was based on airfoils having a maximum gliding ratio L/D = 25. The tethers are assumed to be made of Dyneema R , which has a very high stiffness and yield strength for a low density. The fixed model parameters can be found in [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF]. For the single airfoil the tether has been discretized using 5 segments. For the dual airfoils the main tether has been discretized using 20 segments, while 5 segments were used for both secondary tethers. This results in 41 states for the single airfoil and 207 states for the dual airfoils.

Using the method proposed in Section 3.1, a solution to (19) for the single-airfoil system using parameter values S = 500 m 2 and M A /S = 20 kg/m 2 , and a solution for the dual-airfoil system using parameter values S = 250 m 2 and M A /S = 20 kg/m 2 are computed. The choice of using large-scale systems stems from the results computed in [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF]. The biggest studied system is the one for which dual airfoils are less advantageous, thus the one for which the single airfoil is most competitive.

The computed trajectories for the dual and single airfoil are displayed in Figure 5. It can be seen that the dual airfoils fly at a much higher altitude and at much longer tether lengths than the single airfoil. The trajectories obtained in this study are comparable to the ones obtained in [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF]. The trajectory of each airfoil is smooth and almost circular.

The most noticeable difference from the study proposed in [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF] can be seen in the power extracted from the airmass. In Figure 6 the average power extracted by the single and dual airfoil systems is compared. Two cases are considered: a) the airfoil-airmass interaction is modeled and b) the interaction is not modeled. In both cases, the dual airfoils are more advantageous, but it can be seen that the airfoil-airmass interaction has a stronger impact on the dual airfoils. As it can be seen from Figure 6, the dual airfoils extract 39% less power when the airfoil-airmass interaction is modeled, while the power drop for the single airfoil is 27%.

CONCLUSION & FURTHER DEVELOPMENTS

This paper has proposed a refinement of the results obtained in [START_REF] Zanon | Airborne Wind Energy Based on Dual Airfoils[END_REF]. An airfoil-airmass interaction model has been proposed, in order to investigate if the dual airfoils are indeed more advantageous then the single airfoil.

The results show that the power drop due to the interaction with the airmass is higher for the dual airfoil system. Nevertheless, the dual airfoils still extract more power than the single airfoil.

More elaborate models of the interaction between the airfoil and the airmass are the subject of ongoing research.

Fig. 2 .

 2 Fig. 2. Schematic of the single-airfoil architecture, with N = 1, A = {1}, T (1) = 0.

Fig. 3 .

 3 Fig. 3. Schematic of the dual-airfoil architecture, with N = 3, A = {2, 3}, T (1) = 0, T (2) = 1, T (3) = 1.

Fig. 5 .

 5 Fig. 5. Comparison between the optimal trajectories of the single and dual airfoils for a total wing surface S tot = 500 m 2 . The trajectories are shown as thick lines. The available wind power P w = ρW 3 ∞ /2 is plotted as a dashed line.

19th IFAC World Congress Cape Town, South Africa. August 24-29, 2014


Pascal Belan

Stéphane Gauthier

come

Optimal indirect taxation with a restricted number of tax rates


Zifan Liu
email: zifan.liu@prism.uvsq.fr

Nahid Emad
email: nahid.emad@prism.uvsq.fr

Ben Soufian
email: soufian.benamor@prism.uvsq.fr

Michel Amor
email: michel.lamure@univ-lyon1.fr

Lamure

S B Amor

M Lamure

Towards modeling of epidemic spread: eigenvalue computation

Keywords: Epidemic, pagerank, Complex networks, Non hermitian eigenvalue problem, Sparse matrix

published or not. The documents may come

Introduction

Dynamic complex systems appear in many areas such as physics, biology, and computer networks etc. In the domain of health research, quick response and effective control of widely spreading health crises stay a big challenge for public health officials as well as scientists. In order to simulate the epidemic spreading, such as H1N1 outbreak in France, traditional models need hundreds of experiments and compute the expected outcome by averaging. In addition, these experiments should be adjusted on a daily basis during the initial outbreak. For example, Network Dynamics and Simulation Science Laboratory (NDSSL) has proposed a parallel simulation model Simdemics [START_REF] Bisset | Urgent computing for interaction based socio-technical simulations[END_REF][START_REF] Bisset | Epifast: A fast algorithm for large scale realistic epidemic simulations on distributed memory systems[END_REF], designed to scale to the entire United states (300 Million people). This solver generates roughly 300 GB of data and is expected to increase as more details is added. One run takes 3000 cpu hours on a 1.5TF machine of 448 cores and one experiment takes 100 to 300 runs. As a result, They could only expect 1 to 4 experiments per year. That's not fast enough to answer urgent requests during the beginning phase of outbreak.

There are two deterministic processes being used in the analysis of epidemic spreading, the susceptible-infective-recovered (SIR) process and the susceptible-infective-susceptible (SIS) process, where the difference is whether or not the disease confers lifelong immunity. Most of the existing work in this area has focused on considering the effect of underlying network structure on epidemic dynamics by using tools from statistical physics [START_REF] Barabási | Statistical mechanics of complex networks[END_REF]. The new trend of complex network based models recognize the individual-level randomness and network topology as significant factors on the dynamic of epidemics, which introduces stochastic aspects in the modeling [START_REF] Bowman | Immunization of networks with community structure[END_REF][START_REF] Newman | The spread of epidemic disease on networks[END_REF]. Adaptive network models can be seen as an offspring of this tendency, the idea of which is that individuals usually attempt to reduce their risk of infection by eliminating contact with contagious individuals [START_REF] Bornholdt | Topological evolution of dynamical networks: Global criticality from local dynamics[END_REF][START_REF] Holme | Dynamics of networking agents competing for high centrality and low degree[END_REF][START_REF] Gross1 | Adaptive coevolutionary networks: a review[END_REF]. In other words, the particularity of this approach is to model infection spreading in a population with evolving contacts. However, the above models all need heavy simulations when facing with a large population because more details about the population should be used. Consequently, it is difficult for them to establish effective vaccination campaigns during the beginning phase of an urgent outbreak. We argue that these difficulties could be well handled within an eigenvalue model and discuss its computation using technology of high performance computing.

Before introducing the proposed model, a brief review of epidemiological models is given in the next section. The proposed pagerank-like model is presented in section 3. Section 4 gives the reason why we choose IRAM as computation method to solve the underlying eigenproblem. Parallel resolution of this eigenproblem is discussed in section 5. Section 6 is devoted to numerical simulation.

Review of epidemiological models

One common used measure of infectivity is the epidemic threshold λ c , which is the minimum infectiousness that a disease has to have in order to invade a network. Each susceptible node is infected with rate ν if it is connected to one or more infected nodes. At the same time, infected nodes are cured with rate δ, defining an effective spreading rate λ = ν/δ. If the value of λ is above the threshold, λ ≥ λ c , the infection spreads and becomes persistent. Below it, λ < λ c , the infection dies out exponentially fast [START_REF] Romualdo | Epidemic spreading in scale-free networks[END_REF][START_REF]Modeling community containment for pandemic influenza: A letter report[END_REF].

The homogeneous models (SIR or SIS) assume that the population mixed at random, so that each individual has an equal chance of coming into contact with any other individual [START_REF] Erdös | On the evolution of random graphs[END_REF]. However in real world, it is not rare to find the different mixing rates between the population subgroups [START_REF] Barabási | Emergence of scaling in random networks[END_REF]. So a direct improvement is to avoid the random-mixing assumption. Models that include underlying network structure achieve this goal by assigning to each individual a finite set of contacts [START_REF] Newman | The structure and function of complex networks[END_REF][START_REF] Newman | Assortative mixing in networks[END_REF].

The difference between various network models depends on how individuals are distributed in space and how connections are formed. The Paul Erdös and Alfréd Rényi model (E-R model) consists of n nodes, joined by edges which are placed between pairs of nodes chosen with equal probability p. By using ideas drawn from percolation theory [START_REF] Newman | Spread of epidemic disease on networks[END_REF], it is found that there is a non-zero epidemic threshold,

λ c = 1 < k > (1)
where < k > is the average connectivity. The degree distribution of this model is Poisson and the epidemic growth rate in such a network is reduced in comparison with the random mixing model [START_REF] Keeling | Networks and epidemic models[END_REF]. Despite being one of the oldest and best studied models of a network, E-R model differs from real networks in two crucial ways: it lacks network clustering and it has an unrealistic Poissonian degree distribution [START_REF] Newman | Random graphs as models of networks[END_REF]. By the following, we introduce two other network models to remedy these two inconvenients. Models based on Small-world networks got into our sight by the work of Watts and Strogatz [START_REF] Watts | Collective dynamics of small-world networks[END_REF][START_REF] Watts | Small worlds: the dynamics of networks between order and randomness[END_REF]. Small worlds can be formed by adding a small number of random connections to a lattice1 . And Newman's work [START_REF] Newman | Renormalization group analysis of the small-world network model[END_REF] gives another effective way of constructing this kind of networks. Random networks display low clustering but short path lengths since there are many long-range links, whereas small-world networks have high clustering and short path lengths. These characteristics have important implications in the context of epidemics: the high level of clustering means that most infection occurs locally, but short path lengths means that epidemic spreading through the network is rapid [START_REF] Watts | Collective dynamics of small-world networks[END_REF]. By applying the percolation theory to small-world networks to calculate the threshold, it is found that [START_REF] Moore | Epidemics and percolation in small-world networks[END_REF]:

λ c = 1 + 12φ + 4φ 2 -1 -2φ 4φ = 1 -4φ + O(φ 2 ) (2)
where, it is assumed that each individual is linked to its two nearest neighbors and on average to φ randomly chosen other individuals.

In 1998, a project to map the World Wide Web has revealed a surprising fact that a few highly connected pages are essentially holding the World Wide Web together. Counting how many Web pages have exactly k links showed that the degree distribution followed a power-law. Following researches observed many real world networks that display this phenomenon, while small worlds, random networks have a power-law degree distribution. Scale-free networks can be constructed dynamically by adding new individuals to a network one by one with preferential mechanism. The major contribution of this model is the heterogeneity in numbers of contacts and the existence of hubs (the most highly connected nodes in the network) [START_REF] Cohen | Efficient immunization strategies for computer networks and populations[END_REF]. Hubs in a network play a pivotal role in the spread and maintenance of infection. Research suggests that the simultaneous elimination of as few as 5 to 15 percent of all hubs can crash a system. Despite some practical difficulties, immunizations targeting hubs could be interesting [START_REF] Madar | Immunization and epidemic dynamics in complex networks[END_REF][START_REF] Albert | Error and attack tolerance of complex networks[END_REF][START_REF] Cohen | Resilience of the internet to random breakdowns[END_REF]. Further research indicated the absence of epidemic threshold in scale-free networks [START_REF] Romualdo | Epidemic spreading in scale-free networks[END_REF]. That is, even weakly contagious viruses will spread and persist in the system.

To summarize this section, we focus on the assumptions used for each model and on the existence of epidemic threshold. Recently, an interesting study has proved the close relationship between the epidemic threshold of a network and the largest eigenvalue of network's adjacency matrix, which can subsume many previous known threshold for special case graphs (E-R, BA lished between adjacent individuals in a 2-D grid.

power-law, homogeneous) [START_REF] Wang | Epidemic spreading in real networks: An eigenvalue viewpoint[END_REF].

The proposed model

We propose to make use of Google's pagerank model [START_REF] Page | The pagerank citation ranking: Bringing order to the web[END_REF] by analogy. An individual in a social graph is analogous to a webpage in a web graph. The common concept between pagerank model and our epidemic model is the random walk. In pagerank model, the surfer (or walker) starts from a random page, and then selects one of the outlinks from the page in a random fashion. Each page has two states as being visited by surfer or not. The pagerank (importance) of a specific page represents the probability that the surfer is present at this page. In our epidemic model, the virus could be viewed as a walker and its propagation could be viewed as a path that consists of a succession of random steps. Each individual has two states as being infected or not. The pagerank (importance) of a specific individual represents the probability that the virus reaches this individual during the course of epidemic. To use mathematical formalism, let G = (V, E) be a directed graph with individuals set V and outlinks set E. The graph might be directed. That means, if there is a link i → j in graph G where i, j ∈ V , j → i is not necessarily true. For example, blood disease could only happen from the donators to the acceptors. Suppose that the graph G has n individuals with degree d = (d 1 , d 2 , ..., d n ), where d j is the number of links individual j has to other individuals. At each time step t, the virus has a state s t ∈ V , indicating which individual it is on at time t. If s t = i, then at time step t + 1, the virus moves to one of his neighbors j chosen uniformly at random from all of i's neighbors.

P j,i = P [s t+1 = j|s t = i] = 1 d i if i -→ j 0 otherwise (3)
The probability P (s t = i) for i ∈ V depends only on s t-1 and not on s t-2 , s t-3 , ..., so that {s t } is a Markov chain. {s t } is characterized by its initial state and a transition matrix P , given by P j,i = P (s t = j|s t-1 = i) with P j,i ∈ [0, 1] for all i, j ∈ V and j∈V P j,i = 1. According to Frobenius theorem, λ = 1 is one of the eigenvalue of the transition matrix P and is the biggest eigenvalue. Thus, there is a stationary distribution for the final state of epidemic spread. Let x i be the probability that individual i is infected during epidemic and we write the stationary distribution as x = (x 1 , x 2 , ..., x n ) for the whole population. This infection vector x is independent of the starting distribution and has the relationship: P x = 1 * x = x. To sum up, the infection vector x's implication in social graph is similar to that of pagerank vector in web graph. The problem consists, as a result, to find the dominant eigenvector with 1 as eigenvalue for the transition matrix P of the social graph. An example of 5 individuals has been given in Figure 1. The set

= {0 → 1, 1 → 0, 1 → 2, 2 → 0, 2 → 1, 2 → 4, 3 → 1, 3 → 2, 3 → 4, 4 → 0}.
In matrix formulation, these links give the transition matrix P as:

P =       0 1/2 1/3 0 1 1 0 1/3 1/3 0 0 1/2 0 1/3 0 0 0 0 0 0 0 0 1/3 1/3 0       (4)
The first difficulty with this model is the existence of dangling individuals [START_REF] Bryan | The $25,000,000,000 eigenvector: The linear algebra behind google[END_REF], which containing no outlinks. These individuals will result in one or more columns of zeros in transition matrix P . For our example in Figure 1, if we delete the link 0 → 1, then the first column of matrix P will contain only zeros. Several ideas have been proposed to deal with this problem [START_REF] Ipsen | Pagerank computation, with special attention to dangling nodes[END_REF][START_REF] Eiron | Ranking the web frontier[END_REF]. Research by the initial pagerank paper [START_REF] Page | The pagerank citation ranking: Bringing order to the web[END_REF] indicates that the pagerank could be calculated by removing the links to dangling pages from the web graph. However, theoretically this process might generate new dangling pages and iteratively remove all pages from the graph. In the context of epidemic spread, dangling individuals could be considered as deadends for virus' random walk process. So a nature thinking is to add a loop with probability 1 to these persons themselves. By this way, diagonal elements corresponding to dangling individuals in matrix P are filled with 1. We adopt this simple solution.

The second difficulty with this model is the problem of non-unique rankings. The phenomenon of "small-world" reveals the clustering effect in social networks. Since very few links exist between clusters, some isolated clusters will break the strong connectivity of graph. It can be shown that the transition matrix P will not yield a unique ranking vector x with such isolated clusters [START_REF] Bryan | The $25,000,000,000 eigenvector: The linear algebra behind google[END_REF]. The common solution is to add a jumping vector to the random walk process:

A = αP + (1 -α)ve T ( 5
)
where A is disease transition matrix, v is the teleportation vector, e is the vector [1, ..., 1] T and α, the damping factor, is a positive parameter smaller than 1. The virus has a small probability (1 -α) to jump from any individual to any other individual in a social graph. This would happen, for example, when an infected person (virus carrier) meets someone outside his normal contacts. Considering the preferential attachment of social scale-free networks [START_REF] Barabási | Emergence of scaling in random networks[END_REF], we choose v to be proportional to individuals' degree and normalizes it by "1-norm". In short, there is a small probability that an individual establishes a new temporary link with someone who already has many links.

v =    v 1 . . . v n    , v i = d i i=n i=1 d i (6)
All column sums of the new transition matrix A are 1, so A is still a stochastic matrix with dominant eigenvalue equal to 1: Ax = x. In the example of Figure 1, d = (1, 2, 3, 3, 1), i=4 i=0 d i = 10. By taking α = 0.9, we have:

A = αP + (1 -α)ve T = 0.9 *       0 1/2 1/3 0 1 1 0 1/3 1/3 0 0 1/2 0 1/3 0 0 0 0 0 0 0 0 1/3 1/3 0       + 0.1 *       1/10 2/10 3/10 3/10 1/10             1 1 1 1 1       T =       0.
      (7)
It must be noticed that our proposed model is based on SIS epidemiological process where, at a given time, each individual can be susceptible (S) or infected (I). So we suppose that individuals recover with no immunity to the disease, that is, individuals are immediately susceptible once they have recovered. Furthermore, the characteristic of individuals as well as that of virus are not taken into account. The model depends only on the status of social graphs at a given time. Infection vector x could help health officials to decide the relative importance of different agents in a population facing an epidemic. This is especially useful when the resource of vaccination are limited during the beginning phase of an urgent outbreak. Preorities should be given to those individuals with bigger ranking in x. In addition, a fast computation of this vector could be expected thanks to the efficient implementation of eigenvalue algorithm.

Computation methods

An efficient solution to a large sparse eigenvalue problem strongly depends on the proper choice of iterative methods. Our first objective is to choose the best method to calculate the dominant eigenvector.

A lot of research found that the damping factor α strongly affects the convergence of iterative methods [START_REF] Gleich | Fast parallel pagerank: a linear system approach[END_REF][START_REF] Haveliwala | The second eigenvalue of the google matrix[END_REF]. So another special attention has been paid to investigate how the convergence of the proposed algorithm is influenced by this degree of teleportation.

Algorithms based on matrix-vector product (MVP) might be advantageous. Suppose x is a vector of p-norm 1, Ax = αP x + (1 -α)v(e T x) where e T x is a scalar. So the MVP of A is expressed as MVP of a sparse matrix P plus a vector. Otherwise, direct computation of Ax for a dense matrix A (see example in [START_REF] Holme | Dynamics of networking agents competing for high centrality and low degree[END_REF]) is bottlenecked by memory on large graphs. The simple power method can be used to solve this non-hermitian eigenvalue problem. The idea is to write the initial vector x 0 as a linear combination of n j=1 α j v j , where v j are eigenvectors of A. Without loss of generality, suppose λ 1 is the dominant eigenvalue, we have:

x k = Ax k-1 = A 2 x k-2 = ... = A k x 0 = A k n j=1 α j v j = n j=1 α j A k v j = n j=1 λ k j α j v j = λ k 1 (α 1 v 1 + n j=2 (λ j /λ 1 ) k α j v j ) (8)
For j > 1, |λ j /λ 1 | < 1, so that (λ j /λ 1 ) k → 0, leaving only the pagerank eigenvector v 1 . Employing the known eigenvalue 1 as a shift offers potential for accelerating convergence. Besides, based on the fact that the smallest eigenvalue of A is the reciprocal of the largest eigenvalue of A -1 , shifted inversed version of power method is described in Algorithm 1. The convergence rate of this algorithm is |(λ j -µ)/(λ k -µ)|, where λ j and λ k are eigenvalues of A such that |λ j -µ| -1 is the largest and |λ k -µ| -1 is the second largest eigenvalue of (A -µI) -1 . On one hand, the closer the damping factor is to 1, the closer the matrix A is to the original transition matrix P . If the largest eigenvalue is not well separated from other eigenvalues, then this rate will be close to 1 and the convergence of the method will be very slow.

On the other hand, inverse iteration requires to solve a linear system for the large matrix (A -µI), adding both computational and storage complexity.

To remedy the previous difficulties, Krylov subspace method can be an alternative. Krylov subspace method allows approximation of an eigenpair (λ, x) of A by a Ritz-elements pair (µ ∈ C, x (m) ∈ K m ) where the subspace K m is defined by

K m = Span{x 0 , Ax 0 , ..., A m-1 x 0 } (9)
As shown in power method, the vectors A j x point more and more in the direction of the dominant eigenvector as j increases. This means that for Input: A(n × n): the disease transition matrix with each column sum as 1, x 0 : the initial guess of the infection eigenvector, µ: the shift nearest to the desired eigenvalue (here µ = 1). Output:

x 1 x = x 0 / x 0 ; 2 for j = 1, 2, 3, ... until λ converge do 3 Solve (A -µI)w = x; 4 λ ← µ + x * w w * w ; 5 i ← i max (w); 6
x ← x/(e T i w); 7 end 8 return x;

Algorithm 1: Inverse Iteration Method large m, most of vectors in Krylov subspace will point in about the same direction. The basis is thus ill-conditioned. Generally, it is desired to build a well-conditioned basis {w 1 , w 2 , ..., w m } for K m . One popular choice is an orthogonal basis which leads to Arnoldi procedure. This method approximates k eigenpairs of A by those of a matrix of order m, where k ≤ m ≪ n. This matrix, designated by H m is obtained by orthogonal projection of A onto an m-dimensional subspace K m . Let W m be the matrix whose column w 1 , w 2 , ..., w m constitute an orthogonal basis of K m . The problem is to find µ ∈ C and y ∈ C m such that (H m -µI)y = 0 [START_REF]Modeling community containment for pandemic influenza: A letter report[END_REF] where the matrix H m of dimension m × m, is defined by

H m = W * m AW m . Note that W *
m is the transpose conjugate of W m and x = W m y. Therefore, some eigenvalues of A can be approximated by the eigenvalues of the matrix H m . Solving the problem ( 10) is relatively easy thanks to the Hessenberg structure of matrix. The basic Arnoldi procedure with the orthogonalization refinement [START_REF] Daniel | Reorthogonalization and stable algorithms for updating the gram-schmidt qr factorization[END_REF] is described in Algorithm 2. The sequence of statements in if-clause assures that the new direction f j+1 is numerically orthogonal to the previously computed directions, i.e., to the columns of W j+1 . A larger η sets more strict condition for orthogonality [START_REF] Demmel | Templates for the solution of algebraic eigenvalue problems: a practical guide[END_REF]. Yet a disturbing aspect of Arnoldi method is that the number of iterations (the size of the subspace Input: A(n × n): the disease transition matrix with each column sum as 1, w 0 : the starting vector, m: the size of subspace. Output: w 1 , w 2 , ..., w m : an orthogonal basis of Krylov subspace, H(m × m): the m × m Hessenberg matrix.

1 w 1 = w 0 / w 0 ; 2 v = Aw 1 ; 3 h = w * 1 v; f 1 = v -w 1 h; 4 W 1 ← [w 1 ]; H 1 ← [h]; 5 for j = 1, 2, 3, ..., m -1 do 6 β j = f j ; 7 w j+1 = f j /β j ; 8 W j+1 ← [W j , w j+1 ]; H j ← H j β j e * j ; 9 v = Aw j+1 ; 10 h = W * j+1 v; 11 f j+1 = v -W j+1 h; 12 if ( f j+1 < η h ) then 13 s = W * j+1 f j+1 ; 14 f j+1 = f j+1 -W j+1 s; 15 h = h + s; 16 end 17 H j+1 ← [ H j , h];
18 end Algorithm 2: m-Step Arnoldi Procedure m) needed to compute the desired eigenpair is unknown, except with a very fortunate choice of starting vector w 0 . The basic algorithm increases m until the dominant eigenvalue of A are found. For storage, in addition to A, the method keeps m vectors of length n and an m × m Hessenberg matrix, which gives O(nm + m 2 /2). For computation complexity, matrix-vector product costs about (2m * nnz) operations, where nnz is the number of nonzero elements in A. The modified Gram-Schmidt procedure costs O(m 2 n) operations. Since the size n of a social network may attain millions or even billions of individuals, increasing m causes both storage and computational overhead.

One way to avoid this difficulty is by restarting techniques. The idea is to restart the iteration with a vector that has been preconditioned so that it is more nearly in a k-dimensional invariant subspace which contains the dominant eigenvector. As stated before, the most consuming part in Arnoldi procedure is the MVP due to the very large size of the matrix. Implicitly restarted Arnoldi method (IRAM) [START_REF] Sorensen | Implicitly restarted arnoldi/lanczos methods for large scale eigenvalue calculations[END_REF][START_REF] Sorensen | Implicit application of polynomial filters in a k-step arnoldi method[END_REF], which uses truncated implicitly shifted QR iterations [START_REF] Watkins | The qr algorithm revisited[END_REF], provides a means to reduce the number of MVP needed by Arnoldi procedure from m to r = m-k, with (m-k) the number of shifts used in QR iterations. The sequential algorithm of IRAM is described in Algorithm 3. Each of these shift cycles, in step 6 -10 of Algorithm 3, results in the implicit application of a polynomial in A of degree r to the starting vector [START_REF] Sorensen | Implicitly restarted arnoldi/lanczos methods for large scale eigenvalue calculations[END_REF].

w 0 ← ψ(A)w 0 , where ψ(λ) = r j=1 (λ -µ j ) (11)
The roots of this polynomial are the shifts used in the QR process and these should be selected to filter unwanted information from the starting vector and hence from the Arnoldi factorization. In our case, to find the dominant eigenpair, we could choose as shifts the r eigenvalues with smallest moduli from the spectrum of H m . Concerning the stopping criteria of Algorithm 3, define the vector x = W m y to be a Ritz vector associated with Ritz value µ.

AW m y -W m H m y = Ax -µx = f m e * m y (12)
By using the backward error associated with IRAM [START_REF] Bennani | Stopping criteria for eigensolvers[END_REF], convergence test in step 3 of Algorithm 3 is: f m e * m y < A F ǫ where e m is m th vector

Input:

A(n × n): the disease transition matrix with each column sum as 1, w 0 : the starting vector, m: the size of subspace, r: the number of shifts and m = r + k.

Output:

x: the dominant eigenvector associated with eigenvalue 1.

j R j = H m -µ j I; 8 H m = Q * j H m Q j ; 9 Q = QQ j ; 10 end 11 β k = H m (k + 1, k); σ k = Q(m, k); 12 f k = w k+1 β k + f m σ k ; 13 W k = W m Q(:, 1 : k); H k = H m (1 : k, 1 : k); 14 begin with the k step Arnoldi factorization AW k = W k H k + f k e * k , apply r additional steps of the Arnoldi procedure to obtain a new m-step Arnoldi factorization AW m = W m H m + f m e * m 15 end
Algorithm 3: IRAM of the canonical basis of C m and ǫ is the tolerance. Since the dominant eigenvalue 1 is known, when a Ritz value has converged to 1, IRAM iterations could be stopped.

Analysis of parallelism

We look for individuals or group of people most likely to spread the disease. We have seen that these individuals or groups are determined according to the network status of individuals at a given time. To take into account the networks of realistic size and to find quick solutions to stop the spread of disease, we must perform effective parallel implementations of the methods.

IRAM consists of four main tasks. First, the projection phase: step 1, 2 in Algorithm 3, manipulates the n-sized data sets for sparse MVP. The second phase which includes implicitly shifted QR iterations, step 4-12 in Algorithm 3, acts on m-sized data sets. The third phase constructing the r additional steps of Arnoldi factorization, step 13 -14 in Algorithm 3, manipulates on n-sized data sets for sparse MVP as well. At last the convergence test, step 3 in Algorithm 3, deals with n-sized data sets to calculate f m . Because phase one and three constitute the most expensive part of the algorithm, we propose to distribute them between processors and to run phase two and four redundantly on all processors.

Stochastic simulation using the infection vector

In this experiment, we use a small graph ba to simulate the real world epidemic spread with distribution of vaccination. We consider people who receive vaccination being permanently immunized against viruses. For larger graph, parallelization will be needed due to the memory requirement but the implementation of such parallel simulator is not the objective of the test.

We assume a universal infection rate ν, a jumping rate α (damping factor) and a curing rate δ for every individual. Before each simulation, we randomly choose a set of infected individuals. Then the propagation of virus proceeds by time step. During each time step, an infected individual infects each of its neighbors with probability ν. And this infected individual also passes the disease to another random chosen non-neighbor by probability α. Additionally, every infected individual is cured with probability δ. The result is the average over 10 runs and it is presented in Figure 2. Here, we compare three cases. First of all, without distribution of vaccination, we try to give the worst case for time evolution of infection. Secondly, with random distribution of vaccination, we begin the simulation by distributing vaccination to a random chosen group of individuals. Then, we simulate time evolution of infection. Thirdly, with distribution of vaccination using our pagerank-like model, we calculate the infection vector for the underlying social graph and then distribute vaccination to individuals with big ranking in the vector.

The figure verifies the absence of epidemic threshold in scale-free networks. Without interventions, the epidemic will always enter an endemic state. The second curve, in top-down order from the figure, shows that random distribution of vaccination could not prevent the virus from entering the endemic state. However, distributing vaccination to individuals with big ranking in our pagerank-like model makes the epidemic die out quickly. This simple experiment confirms the important implication of our infection vector for the control of epidemic spread.

Figure 1 :

 1 Figure 1: Small social graph of 5 individuals

1 w 1 5 Q

 15 = w 0 / w 0 ; 2 compute the m-step Arnoldi factorization: AW m = W m H m + f m e * m ; 3 while not converge do 4 compute the spectrum of H m (σ(H m )) and select r shifts µ 1 , µ 2 , ..., µ r ; = I m ; 6 for j = 1, 2, ...r do 7 QR factorization: Q

  Random distribution of vaccination Distribution of vaccination using our model

Figure 2 :

 2 Figure 2: Time series of infection in an 7010-node power-law social graph ba, with α = 0.85, ν = 0.2 and δ = 0.24

Lattices display high clustering but long path lengths because connections are estab-

Acknowledgements

We would like thank Fabrcio Benevenuto from Federal University of Ouro Preto for the twitter graph, Kim Capps from Yahoo! Labs for his help to get access to Alta Vista web graph.