Federico Bribiesca Argomedo
email: federico.bribiesca-argomedo@gipsa-lab.fr

Emmanuel Witrant

Christophe Prieur

Sylvain Brémond

Rémy Nouailletas

Jean-François Artaud

Lyapunov-Based Distributed Control of the Safety Factor Profile in a Tokamak Plasma

A real-time model-based controller is developed for the tracking of the distributed safety-factor profile in a tokamak plasma. Using relevant physical models and simplifying assumptions, theoretical stability and robustness guarantees were obtained using a Lyapunov function. This approach considers the couplings between the poloidal flux diffusion equation, the time-varying temperature profiles and an independent total plasma current control. The actuator chosen for the safety factor profile tracking is the Lower Hybrid Current Drive, although the results presented can be easily extended to any non-inductive current source. The performance and robustness of the proposed control law is evaluated with a physics-oriented simulation code on Tore Supra experimental test cases.

Introduction

Real time control of internal profiles within a tokamak plasma is a key issue to achieve (and maintain) in a safe manner high-performance operation. Given the high uncertainty in online profile reconstruction and measurements, as well as in the modeling of transport phenomena inside the plasma, controlling these internal profiles is a very challenging task and necessitates robust feedback approaches. In particular, an adequate control of the safety factor profile or q-profile, which is determined by the relationship between the toroidal and the poloidal components of the magnetic field, is very important for the plasma discharge MHD stability and possible enhanced energy confinement. An overview of emerging and existing challenges of tokamak plasma control is given in [START_REF] Pironti | Fusion, tokamaks, and plasma control: an introduction and tutorial[END_REF][START_REF] Walker | Emerging applications in tokamak plasma control[END_REF]. For a discussion on advanced tokamak scenarios, refer for instance to [START_REF] Taylor | Physics of advanced tokamaks[END_REF][START_REF] Gormezano | High performance tokamak operation regimes[END_REF] and [START_REF] Wolf | Internal transport barriers in tokamak plasmas[END_REF].

In this article, the problem of controlling the safety factor profile (or q-profile) by controlling the gradient of the poloidal magnetic flux profile in a tokamak plasma is considered. The control design is based on the distributed control-oriented model proposed in [START_REF] Witrant | A controloriented model of the current control profile in tokamak plasma[END_REF], whereas it is assesed in simulation using the physics-oriented code METIS, a module of the CRONOS suite of codes, suitable for closed-loop control simulations [START_REF] Artaud | The CRONOS suite of codes for integrated tokamak modelling[END_REF]. The poloidal magnetic flux, denoted ψ(R, Z), is defined as the flux per radian of the magnetic field B(R, Z) through a disc centered on the toroidal axis at height Z, having a radius R and surface S, as depicted in Fig. 1. As the safety factor is basically scaling as the ratio of the normalized radius to the poloidal magnetic flux gradient, controlling this latter variable allows controlling the safety factor profile, which is our main objective. This is a challenging problem for several reasons:

• the evolution of the q-profile is governed by the resistive diffusion of the magnetic flux, which is a parabolic equation with spatially-distributed rapidly time-varying coefficients that depend on the solution of another partial differential equation related to heat transport; • the control action is distributed in the spatial domain but nonlinear constraints are imposed on its shape (with only a few engineering parameters being available for control, strong restrictions on the admissible shape are imposed); • important uncertainties exist in most measurements, estimations and models;

• nonlinear source terms appear in the evolution equation (in particular the bootstrap current).

The problem of poloidal magnetic flux profile control is closely related, via the Maxwell equations, to the control of plasma current profile. Some previous works show the possibility of controlling profile shape parameters, for instance on Tore Supra: [START_REF] Wijnands | Feedback control of the current profile on Tore Supra[END_REF], where the current profile shape is characterized by the internal inductance and the central safety factor value and experimental results are presented; and [START_REF] Barana | Feedback control of the lower hybrid power deposition profile on Tore Supra[END_REF], where the control of the width of the lower hybrid power deposition profile is shown and experimentally validated. Also, [START_REF] Imbeaux | Real-time control of the safety factor profile diagnosed by magneto-hydrodynamic activity on the Tore Supra tokamak[END_REF] proposes a discrete real-time control of the steadystate q-profile, considering some possible operating modes. Other works consider the distributed nature of the system and use discretized linear models identified around experimental operating point. An example of such work can be found in [START_REF] Laborde | A model-based technique for integrated real-time profile control in the JET tokamak[END_REF], where a model based on a Galerkin projection was used to control multiple profiles in JET; in [START_REF] Moreau | Plasma control in Tore Supra[END_REF], where a reduced order model is used to control some points in the q-profile and in [START_REF] Moreau | Plasma models for real-time control of advanced tokamak scenarios[END_REF], where the applicability of these identification and integrated control methods to various tokamaks is presented.

Specific contributions from the automatic control research community have also started to appear, dealing with simplified control-oriented models that retain the distributed nature of the system. An example is [START_REF] Ou | Optimal tracking control of current profile in tokamaks[END_REF] and related works, where an infinite-dimensional model, described by partial differential equations (PDEs), is used to construct an extremum-seeking controller for the current profile, considering fixed shape profiles for the current deposited by the RF antennas and for the diffusivity coefficients. Other PDE-control approaches, related to Tore Supra, can also be mentioned: [START_REF] Gahlawat | Control and verification of the safety-factor profile in tokamaks using sum-of-squares polynomials[END_REF], where sum-of-square polynomials are used to construct a Lyapunov function considering constant diffusivity coefficients; [START_REF] Gaye | Sliding mode stabilization of the current profile in tokamak plasmas[END_REF], where a slidingmode controller was designed for the infinite-dimensional system, considering timeinvariant diffusivity coefficients; [START_REF] Bribiesca Argomedo | Polytopic control of the magnetic flux profile in a tokamak plasma[END_REF], where a polytopic linear parameter-varying approach is used to build a common Lyapunov function guaranteeing stability of the discretized system with time-varying diffusivity profiles and finally [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF], where an infinite dimensional Lyapunov function is constructed to guarantee the stability and robustness of the controlled system, considering distributed time-varying diffusivity coefficients as well as non-linear shape constraints in the actuation profiles due to the use of two engineering parameters (the power and the parallel refraction index of the lower-hybrid antennas) to control the safety factor profile.

Some references dealing with the control of parabolic partial differential equations (in particular, diffusion equations) can be mentioned as: [START_REF] Zuazua | Approximate controllability for semilinear heat equations with globally lipschitz nonlinearities[END_REF] and, more recently [START_REF] Coron | Global steady-state controllability of one-dimensional semilinear heat equations[END_REF], studying controllability aspects of semilinear heat equations; [START_REF] Krstic | Adaptive boundary control for unstable parabolic PDEs-part I: Lyapunov design[END_REF] and subsequent papers where a Lyapunov function is used on a heat equation with unknown destabilizing parameters; [START_REF] Smyshlyaev | Closed-form boundary state feedback for a class of 1-D partial integro-differential equations[END_REF], where the control of diffusion-convection equations with constant diffusion and distributed convection coefficients is treated; [START_REF] Smyshlyaev | On control design for PDEs with space-dependent diffusivity or time-dependent reactivity[END_REF], where the case of distributed (but not time-varying) diffusion coefficients is studied; and [START_REF] Vazquez | Control of Turbulent and Magnetohydrodynamic Channel Flows Boundary Stabilization and State Estimation[END_REF], where constant diffusion coefficients and distributed time-varying convection coefficients are considered.

The main contributions of this paper are:

• the use of a physically relevant control-oriented infinite-dimensional model (see [START_REF] Witrant | A controloriented model of the current control profile in tokamak plasma[END_REF]) for the design of a distributed control law to track the gradient of the magnetic flux profile (and consequently the q-profile) by means of LH current drive, with particular care given to time-varying effects and possible extension to arbitrary non-inductive current sources;

• building upon [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF] and [START_REF] Bribiesca Argomedo | D1-input-to-state stability of a time-varying nonhomogeneous diffusive equation subject to boundary disturbances[END_REF], the consideration of time varying diffusivity coefficient profiles in the control design, guaranteeing the stability of the system and its robustness with respect to several common sources of errors and unmodeled dynamics;

• the inclusion of previously neglected couplings between the total plasma current control and the magnetic flux profile control; • real-time optimization that includes the nonlinear constraints imposed by the current deposit profiles while preserving the theoretical stability and robustness guarantees;

• the validation of the proposed control approach using the METIS code [START_REF] Artaud | METIS user's guide[END_REF].

The main advantages of the chosen Lyapunov approach for this problem are:

• no assumption of fixed shape of the control inputs is required (this allows to consider, for example, the shift of the location of the peak current density deposited by the LH antenna when changing the engineering parameters);

• no assumption of fixed or slowly time-varying shape of the diffusivity profiles is made (this allows to consider the fast time-variant nature of the temperature profile and the changes in its shape when heating actuators are used);

• since it is not based on an identified linear model, the validity domain for the analysis is larger (as illustrated in the simulation results);

• the extensive robustness analysis presented in [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF] provides a clear idea of the impact of the tuning parameters and the profile estimation errors on the behaviour and stability of the system; • the resulting control law is easy to implement and requires little online computing ressources (ideal for real-time implementation); • the resulting control law can be easily extended to include other non-inductive actuators (e.g. Electron Cyclotron Current Drive antennas).

All the numerical simulations were performed using Tore Supra parameters, as this work aims to assess the capabilities of this novel safety factor control approach before experimental implementation. Tore supra was chosen for its capability to produce long lasting plasma discharges, offering a unique opportunity to develop and test plasma safety factor profile control schemes on relevant timescales. The extension of these results to other tokamaks will be discussed.

In Section 2, the simplified, infinite-dimensional, control-oriented model (based on [START_REF] Witrant | A controloriented model of the current control profile in tokamak plasma[END_REF]) used for the control design is presented, along with the main related physical hypotheses. The couplings between the distributed system and the finite-dimensional subsystem regulating the total plasma current are formulated and an appropriate change of variables suitable for the stability analysis is introduced. In Section 3 the Lyapunov-based control approach is shown and discussed. In Section 4, the proposed control law is implemented and evaluated for Tore Supra with METIS.

Reference Model

In this section, we couple a simplified 1D model of poloidal magnetic flux diffussion with a 0D model of the time-evolution of the total plasma current via the external boundary condition of the 1D model.

Poloidal (1D) Magnetic Flux Transport Model

In this article, the q-profile regulation is done by controlling the gradient of the poloidal magnetic flux profile.

For the development of a suitable control law, a simplified model for the magnetic flux profile ψ in its one-dimensional representation is considered as in [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics[END_REF]:

∂ψ ∂t = η C 2 µ 0 C 3 ∂ 2 ψ ∂ρ 2 + η ρ µ 0 C 2 3 ∂ ∂ρ C 2 C 3 ρ ∂ψ ∂ρ + η V ρ B φ0 F C 3 j ni (1)
where ρ . = φ πB φ 0 is an equivalent radius indexing the magnetic surfaces, φ is the toroidal magnetic flux, B φ0 is the toroidal magnetic field at the geometric center of the plasma, η is the parallel resistivity of the plasma, j ni represents the current density profile generated by non-inductive current sources (source term), µ 0 is the permeability of free space, F is the diamagnetic function, C 2 and C 3 are geometric coefficients and V ρ is the derivative of the plasma volume with respect to the spatial variable ρ. Some other variable definitions are given in Table 1.

Neglecting the diamagnetic effect caused by poloidal currents and using a cylindrical approximation of the plasma geometry (ρ << R 0 , where R 0 is the major plasma radius) the coefficients in (1) simplify as follows:

F ≈ R 0 B φ0 , C 2 = C 3 = 4π 2 ρ R 0 , V ρ = 4π 2 ρR 0 (2)
Variables Description Units ψ

Poloidal magnetic flux profile T m 2 φ Toroidal magnetic flux profile T m 2 q Safety factor profile q . = dφ/dψ R 0

Location of the magnetic center m B φ0

Toroidal magnetic field at the center T ρ Equivalent radius of the magnetic surfaces m a

Location of the last closed magnetic surface m (minor radius of the torus)

r Normalized spatial variable r . = ρ/a t Time s V Plasma Volume m 3 F Diamagnetic Function T m C 2 , C 3 Geometric coefficients η Parallel resistivity Ωm η Normalized diffusivity coefficient η /µ 0 a 2 s -1 µ 0 Permeability of free space: 4π × 10 -7 Hm -1 n Electron average density m -3 j ni Non-inductive effective current density Am -2 j Normalized non-inductive effective current density T m 2 µ 0 a 2 R 0 j ni I p Total plasma current A V loop Toroidal loop voltage V η lh LH current drive efficiency Am -2 W -1 P lh
Lower Hybrid antenna power W N Hybrid wave parallel refractive index

I Ω Ohmic current A V Ω Ohmic voltage V Table 1. Variable definition
Defining a normalized spatial variable r = ρ a , where a (assumed constant) is the equivalent (minor) radius of the last closed magnetic surface, the simplified model is obtained, as [START_REF] Witrant | A controloriented model of the current control profile in tokamak plasma[END_REF][START_REF] Artaud | The CRONOS suite of codes for integrated tokamak modelling[END_REF]:

∂ψ ∂t (r, t) = η (r, t) µ 0 a 2 ∂ 2 ψ ∂r 2 + 1 r ∂ψ ∂r + η (r, t)R 0 j ni (r, t) (3)
with the boundary conditions ∂ψ/∂r(0, t) = 0 and ∂ψ/∂r

(1, t) = -R 0 µ 0 I p (t)/(2π) or ∂ψ/∂t(1, t) = -V loop (t)/(2π)
where I p is the total plasma current and V loop is the toroidal loop voltage, with initial condition:

ψ(r, t 0 ) = a 2 B φ0 1 r s q(s, t 0 ) ds + ψ(1, t 0 ) (4)
For the purposes of this article, j ni is considered as having two main components:

• the auto-induced bootstrap current j bs (produced by trapped particles in the "banana" regime); • the LHCD (Lower Hybrid Current Drive) current deposit j lh .

Using the previous hypotheses and approximating the toroidal magnetic flux as φ(ρ, t) ≈ -B φ0 a 2 ρ 2 /2 (see [START_REF] Witrant | A controloriented model of the current control profile in tokamak plasma[END_REF]), we can consider the safety factor profile to be related to the poloidal magnetic flux as q(ρ, t) = -B φ0 a 2 ρ/ψ ρ .

Normalizing the constants in (3) by defining η . = η /µ 0 a 2 and j . = µ 0 a 2 R 0 j ni to simplify notations, an equilibrium ψ, j, I p , η is defined as the stationary solution of:

0 = η r rψ r r r + ηj r ∀r ∈ (0, 1) (5)
with boundary conditions ψ r (0) = 0 and ψ r (1) = -R 0 µ 0 I p /(2π); where, for any function ξ depending on the independent variables r and t, ξ r and ξ t are used to denote ∂ ∂r ξ and ∂ ∂t ξ respectively. These equations imply that the time-derivative of ψ is constant over the spatial interval [0, 1] (and, in general, different from zero when

V loop = 0).
Around this equilibrium, neglecting the nonlinear dependence of the bootstrap current on the state, the dynamics of the system is given by:

ψt = η r r ψr r + η j, ∀(r, t) ∈ (0, 1) × [0, T ) (6)
with boundary conditions ψr (0, t) = 0 ψr (1, t) = -R 0 µ 0 Ĩp (t)/(2π) and initial condition ψ(r, 0) = ψ0 (r); and where the dependence of ψ, j and η on (r, t) is implicit; Ĩp . = I p -I p and 0 < T ≤ +∞ is the time horizon. For all variables, a tilde represents the difference between the actual value and the equilibrium ( ξ . = ξ -ξ). Our focus is on the evolution of z . = ∂ ψ/∂r, with input u . = j, as defined by:

z t (r, t) = η(r, t) r [rz(r, t)] r r + [η(r, t)u(r, t)] r , ∀(r, t) ∈ (0, 1) × [0, T ) (7)
with Dirichlet boundary conditions:

z(0, t) = 0 z(1, t) = - R 0 µ 0 Ĩp (t) 2π (8)
and initial condition:

z(r, 0) = z 0 (r) (9)
where z 0 . = ∂ ψ0 /∂r.

0D total plasma current dynamic model

Assuming I lh = η lh P lh /(R 0 n) (where η lh is the efficiency of the LH current drive, P lh is the power delivered by the LH antennas and n is the electron line average density) and using a transformer model as in [START_REF] Kazarian-Vibert | Full steady-state operation in Tore Supra[END_REF], the evolution of Ĩp . = I p -I p around an equilibrium (I p , P lh , N , V Ω , I Ω ) can be considered, neglecting the variations of bootstrap current, as given by:

L p M M L Ω ˙Ĩ p ˙Ĩ Ω = -R p 0 0 -R Ω Ĩp ĨΩ + η lh Rp nR0 0 0 1 Plh ṼΩ (10)
with initial condition Ĩp (0) = ĨΩ (0) = 0.

Let us define the matrices A, B and D as follows: 10) can be rewritten, around the equilibrium and adding an integrator to reject constant disturbances, as:

A =    - LΩRp LpLΩ-M 2 MRΩ LpLΩ-M 2 0 MRp LpLΩ-M 2 - LpRΩ LpLΩ-M 2 0 -1 0 0    , B =    - M LpLΩ-M 2 Lp LpLΩ-M 2 0    , D =    LΩη lh Rp (LpLΩ-M 2 )nR0 - Mη lh Rp (LpLΩ-M 2 )nR0 0    Equation (
ζ = Aζ + B ṼΩ + D Plh (11)
where ζ . = Ĩp ĨΩ E T , and E is the integral of the tracking error of I p .

For simplicity in the calculations and proofs, this article considers constant matrices A, B and C. A more advanced approach, using a linear parameter-varying (LPV) formulation (such as [START_REF] Bribiesca Argomedo | Polytopic control of the magnetic flux profile in a tokamak plasma[END_REF] for the total plasma current) can avoid this assumption to extend our results to the time-varying case.

Since L p L Ω -M 2 > 0, matrix A has two eigenvalues with negative real parts (corresponding to the physical system) and one zero eigenvalue (corresponding to the integrator).

The second boundary condition in ( 8) can be written as:

z(1, t) = Cζ (12)
where C . = -R0µ0 2π 0 0 . The relation between the 0D and 1D models is concisely expressed in equations [START_REF] Imbeaux | Real-time control of the safety factor profile diagnosed by magneto-hydrodynamic activity on the Tore Supra tokamak[END_REF] and [START_REF] Moreau | Plasma control in Tore Supra[END_REF].

Feedback Control Approach

Control Problem

For a given reference q-profile, our aim is to ensure a proper tracking using both the LH antenna and poloidal coils. The main control problem can thus be stated as:

• To guarantee the exponential stability of the origin, in the topology of the L 2 norm, defined in the usual manner (see for instance [START_REF] Curtain | An Introduction to Infinite-Dimensional Linear Systems Theory[END_REF]), of solutions of system ( 7)-( 9), with boundary condition given by ( 10)- [START_REF] Moreau | Plasma control in Tore Supra[END_REF], both in open-loop (with u = 0) and by closing the loop with a controlled input u(•, t) corresponding to a change in non-inductive current actuation (in particular LH current);

• to be able to adjust (in particular, to accelerate) the rate of convergence of the system using the controlled input; • to guarantee the stability of the system in presence of a large class of errors. In particular actuation errors, estimation/measurement errors and state disturbances have been considered in [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF].

This problem is illustrated in Figure 2.

The chosen approach to solve this problem is the construction of a strict Lyapunov function. Using this Lyapunov function, a simple control law that guarantees the stability of the system and allows attenuating possible disturbances is presented, based on the results presented in [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF].

Proposed Lyapunov Function

Let us consider the following candidate Lyapunov function: (i) there exists a twice continuously differentiable positive weight f :

V = V 1 + V 2 (13) V 1 = 1 2 1 0 f (r)z 2 (r, t)dr V 2 = 1 2 ζ T P ζ with f : [0, 1] → [ε, ∞) a
[0, 1] → [ε, ∞)
as defined in [START_REF] Moreau | Plasma models for real-time control of advanced tokamak scenarios[END_REF] such that the function V 1 (z) is a strict Lyapunov function for system (7)-( 9) with homogeneous boundary conditions, verifying for some positive constant α 1 , V1 ≤ -α 1 V 1 ; (ii) an independent control loop regulates the total plasma current while ensuring, for some symmetric positive definite matrix P ∈ R 3×3 , some matrix K ∈ R 1×3 and some positive constant α 2 :

P [A + BK + α 2 I 3 ] - R 2 0 µ 2 0 8π 2 (f (1) + f ′ (1)) η(1, t)I 3 ≺ 0 ( 14
)
where • ≺ 0 denotes the negative definiteness of a square matrix, I 3 is the 3 × 3 identity matrix, f ′ represents d dr f ; then the function V (z, ζ) is a strict Lyapunov function for the interconnected system, satisfying for some positive constant γ:

V ≤ -min{α 1 , α 2 }V (z, ζ) + γ sup 0≤τ ≤t | Plh (τ )|, ∀t ∈ [0, T ) (15)
The proof of this result is given in Appendix A.

Remark 3.3 It can be shown, based on the previous theorem, that for constant values of Plh the exponential convergence of the interconnected system can be achieved with the same rate as that of the infinite-dimensional system if the rate of exponential convergence of the plasma current control is faster by at least -R 2 0 µ 2 0 /(8π 2 ) (f (1) + f ′ (1)) η(1, t) than that of the magnetic flux profile. This term accounts for the impact of the coupling between the 0D (finite-dimensional) plasma circuit equation and the 1D (infinite-dimensional) magnetic-flux diffusion equation by means of the total plasma current (boundary condition). Even if this condition is not verified, the error in the infinite-dimensional state will always remain bounded as long as Plh remains bounded and, furthermore, this error will be inversely proportional to the rate of convergence of this subsystem (α 2 ), which is why a fast convergence of the total plasma control is desirable.

For the rest of this article, condition ( 14) is assumed to be verified for some α 2 >> α 1 , and therefore α = α 1 and we can focus only on the evolution of the infinite dimensional subsystem. This assumption is physically justified by the fact that Ohmic current can be generated at a much faster time scale than that of the poloidal magnetic flux diffusion.

Control Law

Based on the constructed Lyapunov function and particularly on equation (A.8) (in the Appendix), a feedback controller was designed, as detailed in [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF], based on a constrained optimization problem. This controller attempts to maximize the rate convergence of the Lyapunov function to zero while respecting the actuator limitations (the cost function is the control-dependent term appearing in the upper bound of the time-derivative of the Lyapunov function given by equation (A.8)).

At each time step, a couple of antenna parameters (P * lh , N * ) is chosen as follows:

(P * lh , N * ) = arg min

(P lh ,N )∈U 1 0 f (r) ηu(P lh , N ) r zdr (16)
thus maximizing the convergence rate with an admissible control action, and subject to the constraint:

1 0 f (r) ηu(P * lh , N * ) r zdr ≤ 0 (17)
to guarantee preserving the stability of the system. Here U . = [P lh,min , P lh,max ] × [N ,min , N ,max ] is the set of admissible values for the engineering parameters of the LH antennas (P lh is the LH power and N the Hybrid wave parallel refractive index). It is important to note that a feasible starting point for the optimization problem always exists (it corresponds to setting u to zero, and thus keeping the open-loop response of the system). For more details, see [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF].

The constraints imposed by the achievable shapes of the LH current deposit are taken into account in the shape of the u function in the problem. Other constraints

Type of error

Type of expected upper bound for the L 2 norm of z Unmodeled dynamics (w)

ce -αt z 0 L 2 + c t 0 e -α(t-τ ) w(•, τ ) L 2 dτ Actuation error (ε u ) ce -αt z 0 L 2 + c 2 t 0 e -α(t-τ ) ε u (•, τ ) H 1 dτ Profile estimation error (ε z ) ce -αt z 0 L 2 + c 3 t 0 e -α(t-τ ) ε z (•, τ ) L 2 dτ Resistivity estimation error (ε η ) ce -α2t z 0 L 2 Table 2.
Expected ISS gain bounds in a simple case, from [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF] (such as a variation rate on the engineering parameters) can also be taken into consideration while solving this problem.

Remark 3.4

The robustness of the proposed scheme with respect to a variety of disturbances and errors was studied in detail in [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF] and extended to the case where there is a non-vanishing disturbance in the boundary condition in [START_REF] Bribiesca Argomedo | D1-input-to-state stability of a time-varying nonhomogeneous diffusive equation subject to boundary disturbances[END_REF].

In particular, Input-to-State Stability (ISS) results were obtained with respect to the considered disturbances. Input-to-State Stability means to guarantee a bounded gain (in some appropriate sense) between the inputs of the system (in this case the errors in estimations of the poloidal magnetic flux profile and resistivity and the disturbances produced by the non-homogeneous boundary conditions and actuation errors) and the controlled output. For a detailed account of ISS results in the finitedimensional case, see [START_REF] Sontag | Input to state stability: Basic concepts and results[END_REF].

A schematic representation of the ISS concept is shown in Figure 3. An example of typical trajectories in ISS systems are shown for the Open-Loop (O.L.) and Closed-Loop (C.L.) cases. It can be seen that, in the presence of disturbances, the system does not converge to zero, but the error remains bounded. The feedback control usually seeks to attenuate the effect of these disturbances. The expected robustness results are provided in Table 2, which summarizes the results presented in [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF]. To interpret this table adequately, the errors should be considered to affect the evolution of the system according to the following equation:

z t = η r [rz] r r + [η(u(η, ẑ) -ε u )] r + w (18)
with the control action u being a function of the estimations of z and η (ẑ . = z -ε z and η . = η -ε η , respectively) and w being a state disturbance that can represent unmodeled dynamics.

Remark 3.5

The constant c appearing in Table 2 arises from a change in norms in the L 2 space; the constant c 2 is determined by the bounds on the resistivity and its spatial derivative; the constant c 3 depends on the "gain" of the controller (once the constraints are satisfied); the constant α is the unactuated rate of convergence of the coupled system; and finally the constant α 2 is a function of the gain of the controller and the bounds on the error of estimation of η and its spatial derivative (for small enough estimation errors, the system remains exponentially stable). For more details, see [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF].

Application

Simulation Scheme

Although rapid simulations using the simplified model described in [START_REF] Witrant | A controloriented model of the current control profile in tokamak plasma[END_REF] were used in the early stages of controller development for tuning purposes, the proposed scheme is validated on a more complex simulation scheme. The Matlab/Simulink interface of the METIS code developed by the french CEA [START_REF] Artaud | METIS user's guide[END_REF] was included in a flexible platform, developed to easily simulate the Tokamak evolution under different assumptions (e.g.: prescribed total plasma current vs. independent control loop using V Ω ) and different actuator models (such as the LH current deposit profile). This platform includes a controller subsystem, a local plasma current control loop and the METIS interface.

The METIS code solves the full version of the resistive diffusion equation, equivalent to (1) (without the simplifying hypotheses made for the control formulation), as presented in [START_REF] Artaud | The CRONOS suite of codes for integrated tokamak modelling[END_REF], but on a 21-point grid. It computes the MHD equilibrium (used to index the isoflux surfaces) using a fast solver based on moments of the Grad-Shafranov equation, which is not taken into account in the controloriented model. The shapes of the current sources are based on simplified analytical formulations (unlike the functions used for the optimization, which are based on scaling laws as in [START_REF] Witrant | A controloriented model of the current control profile in tokamak plasma[END_REF]). The temperature profile (required to compute the resistivity profile) is computed solving time-independent transport equations and using a 0D solver for computing the energy content of the plasma (this is the main difference with respect to CRONOS, which also solves the temperature transport equation). For more details, see [START_REF] Artaud | METIS user's guide[END_REF].

Closed-loop

Tracking with Approximated Equilibrium (ψ, j, I p , η)

General Description

The presented simulation results are based on general parameters of Tore Supra shot TS-31463, but with a much larger variation range for P lh and N (in order to better illustrate the robustness of the controller when the equilibrium is poorly known or when large variations around the known equilibrium are needed). The shot was simulated with METIS.

The chosen inputs to the controller were the spatial derivative of the poloidal magnetic flux profile, as well as a reference profile (which could be generated from a reference safety factor profile), and an estimate of the parallel resistivity profile. Although the poloidal magnetic flux profile (considered available in Tore Supra from Equinox reconstructions in real time [START_REF] Blum | Reconstruction of the equilibrium of the plasma in a Tokamak and identification of the current density profile in real time[END_REF]) and parallel resistivity may not be perfectly known, the robustness of Lyapunov-based controllers with respect to estimation errors in both profiles has already been studied in [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF] and summarized in Table 2. Plasma parameter estimation could also be done using other methods, such as those found in [START_REF] Felici | Real-time physics-model-based simulation of the current density profile in tokamak plasmas[END_REF]. The controller outputs were chosen to be N and P lh .

Since, in practical applications, the full set of available or desired equilibria might not be easily known (in particular due to nonlinearities such as the bootstrap current and couplings with the temperature equation), a single approximated equilibrium point was estimated, simulating in METIS the Open-Loop behaviour of the system for fixed values of the N and P lh inputs. The ability of the controller to reach or approach other desired profiles and stabilize the system around the corresponding equilibrium was tested.

The ramp-up phase of the simulated shot was done in open-loop (as far as N and P lh are concerned, using only an independent control loop for I p ) based on TS-31463. The controller was activated at t = 9 s. Since Real-Time implementation is desired, an offline computation was done to construct a table with profiles of LH current deposit as a function of the input parameters to the controller using scaling laws as in [START_REF] Witrant | A controloriented model of the current control profile in tokamak plasma[END_REF]. The real-time optimization algorithm can then perform a constrained gradient-descent, using the estimated values of the resistivity profile and the state, and find a suitable control action in less than two hundred microseconds (using MATLAB functions and an Intel processor running at 2.43 GHz).

The global parameters during the flat-top phase are: constant total plasma current of 580 KA; constant line-average electron density of 14.5 × 10 18 m -2 ; noninductive heating and current drive using LHCD with P lh ∈ [0.17, 3.5] MW and N ∈ [1.70, 2.30]; and constant toroidal magnetic field at the center of the plasma of 3.69 T. 4.2.2. First case: Independent I p control, large variations of P lh , temperature profile disturbed by ICRH heating. This test case was chosen to illustrate the robustness of the controller under non-ideal circumstances:

• the total plasma current I p is independently controlled using the poloidal field coils. Since it is no longer considered to be perfectly followed, it becomes a source of disturbances in the plasma edge (boundary condition of the partial differential equation); • during the shot, reference profiles are chosen far from the equilibrium used for the controller design (corresponding to values of P lh and N of 2.76 MW and 2, respectively). This highlights the contribution of non-linear terms in the equations that were neglected near the equilibrium;

• a disturbance in the temperature profile is introduced in the form of 1.5 MW of power from ion-cyclotron radio heating (ICRH) antennas, which cannot be compensated with the LH actuator, thus rendering the target q-profile inaccessible to the controller;

• the model given to the controller for the LH current deposit is based on a gaussian profile approximation with parameters determined by scaling laws, as described in [START_REF] Witrant | A controloriented model of the current control profile in tokamak plasma[END_REF]. This does not correspond to the internal METIS model (where, even though the LH efficiency is calculated using scaling laws depending on plasma and wave parameters [START_REF] Goniche | Lower hybrid current drive efficiency on Tore Supra and JET[END_REF], the shape is based on Landau absorption, accessibility and caustics [START_REF] Artaud | METIS user's guide[END_REF]);

• the general parameters were taken from experimental measurements of shot TS-31463 and therefore introduce measurement noise to all the variables used to compute the control action.

The tracking of the desired q-profiles and the engineering parameters prescribed by the controller for the LH antenna are presented in Figure 4. First, Figure 4 (a) depicts the evolution of the safety factor at six points (corresponding to r = 0, 0.1, 0.2, 0.3, 0.5 and 0.8) during the control window (9 s ≤ t ≤ 23 s). The solid lines represent the simulated evolution of the closed-loop system while the dashed lines represent the reference values given to the controller (calculated from an open-loop simulation without disturbances). It can be seen that the tracking is satisfactory regardless of the radial position. For t between 13.4 s and 16.6 s (and again after 22.6 s) the tracking error does not converge to zero, which corresponds to the introduction of a disturbance caused by ICRH. Nevertheless, the system remains stable and the error small, even though the reference profile is no longer achievable. As soon as the disturbance is removed, the tracking errors are promptly reduced. No undershoot or overshoot is noticeable in the tracking of the reference profile (which is a desirable property to avoid triggering unwanted magnetohydrodynamic modes). There is a lag between the reference profile and the response of the system, which is to be expected since: (a) an error between both signals has to appear before the feedback controller can act, and (b) the rate at which the gradient-descent optimization algorithm is allowed to modify the controlled inputs is limited to filter out noise and to prevent sudden changes to the engineering parameters of the LH antennas. Figure 4 (b) presents the evolution of the engineering parameters used to track the safety factor profile. In solid lines, the parameters N and P lh set by the controller are shown. The parameters used to construct the reference profile are depicted with dashed lines. These figures show that the controller is able to properly reconstruct the engineering parameters required to obtain the desired safety factor (with the saturation on the rate of change of the parameters revealed by the constant slope during big changes in reference). A mark at the beginning of the control action shows the approximated equilibrium around which the controller was designed. The dashed-dotted line shows the ICRH power injected into the system (only applied to the closed-loop system and not taken into account for the reference generation, as would be the case with unexpected disturbances encountered during actual tokamak operation). These results show the robustness of the proposed Lyapunov-based approach to changes in operating conditions. Although the value of N is modified during closed-loop operation, the changes are very small. Finally, even though there is some overshoot in the control parameters chosen by the controller, as mentioned before, these do not cause overshoots in the tracked safety factor profiles.

The tracking of the total plasma current and the resulting loop voltage are depicted in Figure 5. Figure 5 (a) shows the plasma current tracking efficiency of a well-tuned proportional-integral (PI) controller despite constant disturbances (in this case, induced by changes in LH current). During the shot, the current driven by the LH antenna varies considerably (since the LH power is driven between 3 and 1 MW). Nevertheless, the tracking error in the plasma current remains small and goes to zero once the safety factor profile stabilizes. This error acts as a small disturbance in the boundary condition of the partial differential equation that eventually vanishes. The resulting loop voltage can be seen in Figure 5 (b). This figure shows the capability of the controller to handle both the non-inductive current drives and the addition of inductive current. The general shot conditions are the same as those chosen for the previous example. Nevertheless, the safety factor reference profile is generated by changing the N parameter between 1.75 and 2.25 while P lh is held constant at 2.7 MW. The results of the simulation can be seen in Figures 6 and7. Although the safety factor reference is different from that used in the previous example (for instance, the variations of q at r = 0.5 are much more important in the second simulation, while the central safety factor varies less), the tracking remains satisfactory, as seen in Figure 6 (a). As in the previous example, the tracking error does not converge to zero when the ICRH disturbance is introduced (between 13.4 s and 16.6 s, and after 22.6 s), yet the system remains stable and the error small. Figure 6 (b) shows a good reconstruction of the original engineering parameters used to generate the reference, except when the ICRH disturbance is present. The effect of the disturbance is attenuated by an offset in the equilibrium N values. As in the previous example, the overshoots present at some points in the control profile do not induce overshoots for the safety-factor tracking.

Figure 7 shows that the total plasma current tracking (boundary condition) is better than in the previous case, which is to be expected since the variations of LH power are much smaller. This has a direct impact on the tracking of the safety factor profile near the edge (r = 0.8 in Figure 6 (a)). Finally, Figure 7 (b) shows that the non-inductive control effort to maintain the prescribed total plasma current is smaller than in the previous case. 4.2.4. Third case: Independent I p control, reference step, P lh disturbance, temperature profile disturbed by ICRH heating. This last case was chosen to illustrate the interest of closing the loop for the system. The general conditions are the same as those in the first case. The objective is to stabilize the system and switch between two wellknown equilibrium values (i.e. such that, without disturbances, the open-loop system will converge to the reference). In the first two cases, one of the main objectives was to illustrate the possibility of tracking references with imperfect knowledge of the actual equilibrium values. In this case, the objective is to compare the performance and disturbance rejection of the system in closed-loop and open-loop. Two different disturbances are applied to the system:

• a disturbance in the temperature profile is introduced as an additional 1.5 MW of power from ion-cyclotron radio heating (ICRH) antennas, which cannot be compensated with the LH actuator, thus rendering the target q-profile inaccessible to the controller; • a disturbance in the LH power applied to the system (vs. the one prescribed by the controller) of 1 MW (this can be compensated by increasing the LH actuation).

The results of the simulation can be seen in Figures 8 (a) and (b). It can be seen that, in the absence of disturbances, the closed-loop system behaves similarly to the open-loop one (albeit with small differences between 10 and 13 s). The ICRH disturbance is attenuated, particularly at values of r between 0.1 and 0.3 (where the control action is more present). Finally, the LH power disturbance is rejected by the controller, while the open-loop system converges to a safety factor profile far from the reference.

This simulation the interest of using a feedback control with respect to using an open-loop strategy only.

Conclusions and Perspectives

In this article, a control law was designed for the safety factor profile tracking in a tokamak, via the gradient of the poloidal magnetic flux. This control law is set using a physically-relevant simplified model. Based on specific simplifying assumptions, theoretical stability was guaranteed despite large classes of disturbances and actuation/estimation errors. The couplings between the 0D total plasma current control and the 1D magnetic flux profile evolution (via the boundary condition), as well as actuator saturations, were included in the theoretical developments.

The chosen non-inductive current actuator, based on Tore Supra capabilities, was the LH current drive. Nevertheless, the theoretical basis is independent from the actual form of the current deposit and the proposed controller design strategy can easily be applied using other non-inductive current sources (such as ECCD and neutral-beam injection). The controller was validated using the METIS code [START_REF] Artaud | METIS user's guide[END_REF]. The tracking of the full safety-factor profile shows promising results as safety factors far from the estimated equilibrium could be adequately reached, with central safety factors varying both above and below the critical q = 1 value. The robustness of the proposed scheme with respect to unmodeled disturbances in the temperature profile was tested by heating the plasma with ICRH power.

The chosen controller solves a reduced online optimization problem and relies on some off-line calculations to run in real-time. The average time required by the control algorithm to compute the control values at each time step was under 160 µs (using using a Matlab R function running on a processor at 2.54 GHz), well below the simulation sampling time. Future experimental validation of the control scheme in Tore Supra is expected.

The proposed method is based on a first-principle driven approach that is, as such, applicable to any tokamak. To use these results on other tokamaks, a new weight for the presented Lyapunov function may be required. The heuristic procedure presented in [START_REF] Bribiesca Argomedo | D1-input-to-state stability of a time-varying nonhomogeneous diffusive equation subject to boundary disturbances[END_REF] allows to compute such weights in a simple manner for relatively general resistivity profiles. If needed, the procedure used to prove the main result in [START_REF] Bribiesca Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF] can be applied to the resistive diffusion equation without using the cylindrical hypothesis (the terms C 2 , C 3 , F and 1/ρ in equation (1) could be computed from a 2D MHD equilibrium).

Figure 1 .

 1 Figure 1. Coordinates (R, Z) and surface S used to define the poloidal magnetic flux ψ(R, Z).

Figure 2 .

 2 Figure 2. Diagram representing the coupling between the finite-dimensional and infinite-dimensional subsystems and the main control problem.

  twice continuously differentiable positive function and P = P T ∈ R 2×2 a symmetric positive definite matrix. Remark 3.1 The weighting function f (r) is computed to guarantee the strict decrease of the Lyapunov function. A constant f (r) is not, in general, enough to guarantee this decrease for all admissible values of the resistivity profile. Theorem 3.2 If the following two conditions are satisfied:

Figure 3 .

 3 Figure 3. Schematic depiction of the Input to State Stability property and closed-loop improvements.

  Tracking of the safety factor profile. Dashed line: reference q value; full line: obtained q value.

  Controller action and ICRH disturbance.

Figure 4 .

 4 Figure 4. Safety factor profile tracking and radio-frequency antenna control parameter evolution.

  Tracking of the total plasma current.

  Resulting V loop due to Ip tracking.

Figure 5 .

 5 Figure 5. Total plasma current evolution and corresponding loop voltage.

  Tracking of the safety factor profile. Dashed line: reference q value; full line: obtained q value.

  Controller action and ICRH disturbance.

Figure 6 .

 6 Figure 6. Safety factor profile tracking and radio-frequency antenna parameter evolution.

  Tracking of the total plasma current.

  Resulting V loop due to Ip tracking.

Figure 7 .

 7 Figure 7. Total plasma current evolution and corresponding loop voltage.

  Tracking of the safety factor profile. Dashed line: reference q value; full line: closed-loop q value; dash-dotted line: open-loop q value.

  Controller action and ICRH disturbance.

Figure 8 .

 8 Figure 8. Safety factor profile tracking and radio-frequency antenna parameter evolution.

  4.2.3. Second case: Independent I p control, large variations of N , temperature profile disturbed by ICRH heating. The purpose of this second test case is to show the versatility of the controller with respect to different available control parameters.

Acknowlegments This work was carried out within the framework of the European Fusion Development Agreement and the French Research Federation for Fusion Studies. It is supported by the European Communities under the contract of Association between EURATOM and CEA. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The research leading to these results has also received funding from the European Union Seventh Framework Programme [FP7/2007-2013] under grant agreement n o 257462 HYCON2 Network of excellence.

Appendix A.

Proof of Theorem 3.1: Choosing an adequate function f (r), for example the one proposed in [START_REF] Bribiesca Argomedo | D1-input-to-state stability of a time-varying nonhomogeneous diffusive equation subject to boundary disturbances[END_REF], satisfying condition (i) of Theorem 3.2 and with f (1) + f ′ (1) < 0 we have that, along the solution of ( 7)-( 9):

where α 1 > 0. This equation assumes the total current density, defined as in [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics[END_REF]:

to be zero on the last closed magnetic surface. It can be relaxed to assume only uniform boundedness and Lipschitz-continuity in time with some modifications, as presented in [START_REF] Bribiesca Argomedo | D1-input-to-state stability of a time-varying nonhomogeneous diffusive equation subject to boundary disturbances[END_REF], albeit at the expense of more conservative bounds. Differentiating V 2 along the solution of (11), we get:

Considering a control ṼΩ of the form Kζ, with K a row vector with three elements, the previous equation is equivalent to:

If condition (ii) from Theorem 3.2 is satisfied, and from the choice of f , it implies that ζ remains bounded if Plh is bounded, and in particular that, for some positive constant γ:

Recalling the definition of the boundary condition [START_REF] Moreau | Plasma control in Tore Supra[END_REF], and from the definition of C, the derivative of V along the solution of the coupled system is given by:

and hence:

which in turn implies that:

where α . = min {α 1 , α 2 }. This completes the proof.

Appendix B. Brief Introduction to Lyapunov Functions and Stability of Dynamical Systems

This section is intended to provide an intuitive overview of the use of Lyapunov functions for stability analysis of dynamical systems and is provided for pedagogical reasons only. For more precise results on Lyapunov theory for infinite dimensional systems, the reader can refer to [START_REF] Luo | Stability and Stabilization of Infinite Dimensional Systems with Applications[END_REF] or [START_REF] Curtain | An Introduction to Infinite-Dimensional Linear Systems Theory[END_REF]. For the use of strict Lyapunov functions in infinite dimensional systems, see [START_REF] Mazenc | Strict Lyapunov functions for semilinear parabolic partial differential equations[END_REF]. Consider an autonomous (i.e. its evolution depending only on internal variables) finite-dimensional dynamical system with state vector x(t) ∈ R n for all t ∈ [0, T ), whose evolution is given by:

where F : R n → R n is a (possibly nonlinear) Lipschitz function and ẋ denotes the time-derivative of x. The origin of system (B.1) is defined as a Globally Asymptotically Stable equilibrium if for every initial state x 0 ∈ R n , x(t) < +∞ for all t ∈ [0, T ) and x(t) → 0 as t → +∞, where • is a norm in R n . It is a Globally Exponentially Stable equilibrium if for every initial state x 0 ∈ R n , x(t) ≤ ce -αt x 0 for all t ∈ [0, T ), for some positive constants c, α.

A continuously differentiable candidate Lyapunov function V (x) for the system (B.1) has to verify several technical conditions, in particular [START_REF] Khalil | Nonlinear Systems[END_REF]:

is bounded above and below by β( x ) and δ( x ), respectively, two smooth increasing functions such that β(0) = δ(0) = 0.

These conditions make a candidate Lyapunov function analogous to a potential or energy function in a physical system. If it can be shown that V .

= ∂V ∂x • ẋ = ∂V ∂x •f (x) ≤ 0 for all x ∈ R n with the equality occurring only for x = 0, then V is called a global Lyapunov function. The existence of a global Lyapunov function guarantees the global asymptotic stability of the origin of system (B.1).

If the Lyapunov function further satisfies, for some positive constant c, V ≤ -cV (x) for all x ∈ R n , it is called a strict Lyapunov function [START_REF] Sontag | Input to state stability: Basic concepts and results[END_REF]. The existence of a global strict Lyapunov function equivalent to the chosen norm in R n guarantees the global exponential stability of the origin of system (B.1).


Alexandre Soro

Mathieu Cunche

Jérôme Lacan
email: jerome.lacan@isae.fr

Vincent Roca
email: vincent.roca@inria.fr

Erasure Codes with a Banded Structure for Hybrid Iterative-ML Decoding

This paper presents new FEC codes for the erasure channel, LDPC-Band, that have been designed so as to optimize a hybrid iterative-Maximum Likelihood (ML) decoding. Indeed, these codes feature simultaneously a sparse parity check matrix, which allows an efficient use of iterative LDPC decoding, and a generator matrix with a band structure, which allows fast ML decoding on the erasure channel. The combination of these two decoding algorithms leads to erasure codes achieving a very good trade-off between complexity and erasure correction capability.

I. INTRODUCTION AND RELATED WORKS

For the transmission of data packets on erasure channels, linear binary FEC codes often offer the best compromise between fast encoding/decoding operations and a good level of erasure recovery capability. For example, random binary codes have a correction capability very close to channel capacity [START_REF] Studholme | Random matrices and codes for the erasure channel[END_REF]. Unfortunately, the decoding complexity of random codes is often prohibitive because the decoding algorithm, named Maximum-Likelihood (ML) decoding, basically consists in solving a linear system, which can be done by inverting the associated matrix. This matrix represents the relations between the set of received symbols and the set of missing symbols, in our representation a row contains the coefficient used to build an encoding symbol. If k denotes the dimension of the code the inversion has a complexity of O(k 2 ) row operations. In the remaining of the paper the complexity will be evaluated in terms of row operation. A row operation includes the sum of two rows of the matrix and the sum of the corresponding symbols. It is importance to notice that in pratical applications the size of the symbols can be up to several hundred bytes.

To reduce this complexity, Studholme and Blake [START_REF]Windowed erasure codes[END_REF] showed that similar erasure capability can be obtained when the non zero entries of the generator matrix are located in a band of length 2

√ k and where each column contains 2 log(k) nonzero elements. With this improvement, the complexity of the decoding is reduced to O(k 3/2 ) row operations.

The class of Fountain codes, like LT [START_REF] Luby | LT codes[END_REF] or Raptor [START_REF] Shokrollahi | Raptor codes[END_REF] codes can also obtain very good performance in terms of erasure correction capability. An estimation of the performance achieved by Raptor codes with ML decoding is provided in [START_REF] Luby | Reliable multimedia download delivery in cellular broadcast networks[END_REF]. However, despite the description of decoding algorithm provided in [START_REF] Luby | Raptor Forward Error Correction Scheme for Object Delivery[END_REF], the complexity needed to achieve this level of performance is not very clear.

LDPC codes are another class of binary codes providing good level of decoding performance with extremely fast encoding/decoding algorithms [7][8]. Indeed, the classical iterative decoding algorithm for the erasure channel, based on the work of Zyablov [START_REF] Zyablov | Decoding complexity of low-density codes for transmission in a channel with erasures[END_REF], has a linear decoding complexity. The drawback of this algorithm is that it does not reach the performance of ML decoding. Moreover, the sparsity of the matrices also reduces the ML performance compared to full random matrices.

Recently, two independent works [10][11] proposed a hybrid iterative-ML decoding algorithm, where ML decoding is only used when the iterative decoding does not succeed to decode a received codeword.

The new codes we introduce in this paper also rely on a hybrid decoding. Our goal is to build these codes in such a way that both the parity check matrix and the associated generator matrix have good properties for the iterative and ML decoding schemes. More precisely, the parity check matrix should be sparse in order to yield good performance with the iterative decoding. At the same time, the generator matrix should have a banded structure in order to reduce the computational complexity of the ML decoding scheme, as explained in [START_REF]Windowed erasure codes[END_REF]. Thanks to a polynomial representation of the rows and the columns of the generator and parity check matrices, we introduce a method that enabled us to build such codes.

The idea of window based encoding in LDPC codes has been proposed by Haken, Luby and al. in patent [START_REF] Haken | On demand encoding with a window[END_REF]. However, the goal of this patent is only to minimize memory access during encoding, by localizing the access in a window that slides over the input file. Independantly of whether our proposal falls into the scope of this patent or not, we see that, from a purely scientific point of view, the goal of [START_REF] Haken | On demand encoding with a window[END_REF] completely departs from the approach discussed in the current paper, as well as the theoretic tools we introduce to achieve our goals.

The paper is organized as follows: we detail the polynomial approach for the construction of our LDPC-Band codes in Section II. In Section III, we analyse the erasure correction capability of the proposed code and we provide simulation results. Finally, in Section IV we evaluate the complexity of the hybrid decoding. Then we conclude.

II. CONSTRUCTION OF THE CODE

A. General structure of the code

Let G be the systematic binary generator matrix of the banded code: G = (Id|M )

Let H be its associated parity check matrix:

H = (A|U )
Let k be the dimension code and n the length of the code. It

follows that G is a k × n matrix, A is a (n -k) × k matrix and U is a (n -k) × (n -k) binary matrix.
Let us define U as a lower triangular Toeplitz matrix defined by its first column (1, u 1 , . . . , u n-k-1 ). Since all the diagonal elements are equal to 1, U is a full rank matrix. In addition, let us consider the associated polynomial u

(x) = 1 + n-k-1 k=1
u k x k . The coefficient of U are the following:

{u i,j } =    u i-j if i > j 0 if i < j 1 if i = j U =                1 0 0 . . . 0 0 u 1 1 0 u 2 u 1 1 . . . . . . . . . u 2 u 1 . . . u i . . . u 2 . . . 0 . . . u i . . . 1 0 u n-k-1 . . . u i u 1 1               
The relations between the generator matrix G and its associated parity check matrix H : G × H T = 0 give:

M = (U -1 A) T U M T = A
Let B be the band width of the band matrix M studied here, which will have the following form:

M =    m 0,0 m 0,1 . . . m 0,B-1 0 . . . 0 0 m 1,0 m 1,1 . . . m 1,B-1 0 . . . . . . 0 . . . . . . . . . . . . 0    For i = 0, . . . , k -1, let us define the polynomial m i (x) = B-1
j=0 m i,j x j corresponding to the i th row of M . Since A = U M T , it can be shown that, for i = 0, . . . , k -1, the i th column of A is characterized by the polynomial a i (x) such that a i (x) = u(x)m i (x) (see Figure 1). This result is important as it allows to construct banded generator matrices from particular distributions in the parity check matrix. Let us define the Hamming weight of a polynomial as the number of monomials. Then, the number of non-zero elements of a column of A is the weight of a i (x), noted W (a i ).

H =                a 0,0 u 0 a 0,1 a 1,0 u 1 u 0 a 0,2 a 1,1 a 2,0 u 2 u 1 u 0 a 0,3 a 1,2 a 2,1 u 3 u 2 u 1 . . . a 1,3 a 2,2 . . . . . . u 3 u 2 . . . . . . a 2,3 . . . u 3 . . . . . .               

B. Design of the matrices

In order to build a code supporting hybrid iterative/ML decoding, a first constraint is to optimize the iterative decoding on the parity check matrix. This decoding is very sensitive upon the degree of the symbol and check nodes, i.e. the weight of the columns and the rows of the parity check matrix. Thanks to the polynomial approach, these parameters can be controlled. Indeed the weight of the columns and rows of U is defined by the number of monomials in u(x), and the weight of the column of A is defined by the number of monomials in the different a i (x). We will see later how the row weight can be controlled. The second constraint concerns the ML decoding. As explained before, the generator matrix of this code must be band of width B in order to reduce the ML decoding complexity. We do not impose any other constraint on the generator matrix as we assume that the matrix is sufficiently random to allow good correction capabilities under ML decoding.

In order to support an efficient iterative decoding, the columns of the parity check matrix must have a small weight. In particular we must choose a small weight polynomial for u(x), the polynomial defining U the right side of the parity check matrix.

The main point of the process is to find polynomials m i (x) of degree close or equal to B, such that a i (x) have few monomials, where a i (x) = m i (x)u(x). Let us call these polynomials candidate polynomials. Finding these candidate polynomials can be achieved with an exhaustive search, in advance, once the various code parameters (n, k, B) are known. The details of how this search is done does not impact the decoding efficiency and therefore are not described in the reminder of this paper.

Once a polynomial u(x) and a set of candidate polynomials have been selected, the weight of the columns of parity check matrix is fixed. But the weight of the rows can still be modified to fit a specific value. Indeed by permuting two polynomial in M , or by exchanging a polynomial with another candidate polynomial of same weight, we can modify the weight of several rows without changing the weight of the columns.

C. Optimizations

We have now defined a way to implement banded codes for both decoding techniques. However this implementation is limited to a number of encoding symbols of exactly n = 2k+B. Furthermore, these codes turn out to be quite inefficient for the ML decoding in terms of erasure recovery capability, because the first non-systematic symbols only protect a small part of the source symbols. This is a negative side effect of the approach. We now explain how to solve this problem.

Given a polynomial u(x) and a bandwidth B, the set of candidate polynomials usualy contains several elements. Therefore it is possible to choose different polynomials and to interlace them to build the matrix M , and the matrix A inherits from a corresponding interlacement in its polynomial. This technique has two main benefits. The first one is an improvement of the decoding performances for the ML decoding, with limited impacts on the iterative decoding. Indeed, many polynomials interlaced produce more variety in the generator matrix and avoid diagonals and regularities that damage the performances of the ML decoding. The second benefit is the suppression of the side effect mentionned above. The idea is to suppress the B 2 first and last columns of M . As it breaks the polynomial relations between A, U and M because of the first and last B 2 lines of M that are now truncated, we replace these lines with candidate polynomials which are of degree B 2 at maximum, in order to keep the matrix banded, and keep the constraints upon A and U . As a consequence, the encoding packets at the edge of the non-systematic part of the generator matrix are now connected to more source symbols. This also suppress the influence of the band size on the code rate, as the number of symbols produced is now n = 2k, meaning a code rate strictly equal to 1/2. The last point is to handle different code rates. To that goal, we define the f parameter which corresponds to the column shift between consecutive lines of M . For instance, for a ratio of r = 1/2, each successive line is shifted by 1 position from the previous one, i.e. f = 1. In order to handle other ratios, we need to take f = 1/r -1. Moreover, in this case, the first and last B 2f lines are of degree B 2 . If f is not an integer, we should use a family of integer offsets {f i } such as the average value of this set is equal to f .

To summarize, in this section we have defined a way to create erasure band codes that handle efficiently (both in terms of erasure recovery capability and algorithmic complexity) both iterative (thanks to the sparseness of the parity check matrix) and ML decoding (thanks to its banded structure), for different code rates.

III. CORRECTION CAPABILITY ANALYSIS

A. Theoretical Analysis

Influence of the band on the ML performances: For the ML performances on the generator matrix G, the crucial point is to determine a band size that ensures both a quick decoding and good correction capabilities. With the first constraint, the band size used should be as narrow as possible. However, this size should also be sufficient in order to avoid an explosion of the decoding overhead. It is important to notice that the band size corresponds to the maximum degree of the associated binary polynomials m i (x).

Studholme and Blake [START_REF]Windowed erasure codes[END_REF][1] have studied the influence of the window size on binary random matrices constrainted in a window. The conclusion of their work is that for random binary matrices constrained in a window, its size has to be kept close to 2 √ k. Moreover, in this window, the number of non-zero elements must be at least 2 log k. These values guarantees that the probability of being full rank is close to the one of pure random matrices. As our approach uses specific configurations of banded matrices, this value should be considered as a minimum. The band size in our case has to be kept as low as possible and for long dimensions, as close as possible to this value, because it will have a great impact on decoding speed. The size of the band will then be a tradeoff between a good ML decoding speed and a strong ML decoding performance.

Density evolution analysis: A density evolution approach [START_REF] Richardson | Design of capacity approaching irregular low density parity check codes[END_REF] for the design of this code has been studied, but have faced many obstacles. If we choose to consider a infinite matrix with an infinite band width, the density evolution will give the same distributions as for a standard Repeat-Accumulate code. But building a finite length matrix having a such distribution can be difficult as the rows distribution strongly depends on the choice of the polynomials and their permutations. Furthermore as the band is small, a strong irregularity in the joint degree distribution [START_REF] Kasai | Detailed representation of irregular ldpc code ensembles and density evolution[END_REF] may appear.

On the other hand if we choose an infinite matrix with a finite band width, the hypothesis of infinite length cycle does not hold anymore.

B. Simulation Results

In this section, we compare the correction capability of our codes with standard LDPC-Staircase codes. First of all, it is worth noticing that we will use band codes with multiple polynomials interlaced. Indeed, experiments have illustrated the fact that using only one polynomial degrades greatly the decoding performances and will not be discussed here.

The crucial point is now to select a set of polynomials that will respect the conditions in the band, but also that have great performances in both decoding. In this study, we have tried to keep as close as possible of iterative adapted distributions. However, because of the additional conditions of banded matrices, an evolution density approach is not straightforward. It is also worth noticing that the chosen set of polynomials is always a compromise between performance of iterative decoding and ML decoding. This allows the user to adapt the polynomials to the environment of the channel. Experiments have also shown that the way the polynomials are interlaced and their position has only a little impact on the ML decoding performances, even for a random interlacing. The results presented in Table I and Table II are based on a regular interlacing for B = 100 and B = 200 for a code rate of 1/2. We study here, the average overhead required in order to successfully decode the received symbols.

Decoding scheme Iterative

Maximum-Likelihood LDPC-Staircase N 1 = 5 As expected, the band size is a crucial factor for the decoding performances. A band size of 100 is less performant but will have a greater decoding speed. In order to have efficient systems, we can see that the band size has to be large compared to the Windowed Erasure Codes. These results also show that LDPC-Band codes are really close to standard LDPC-Staircase codes, while being more constrained.

IV. COMPLEXITY ANALYSIS

A. Theoretical Analysis

A key point in data transmission is the encoding and decoding speed of the code. In this study we do not include the creation of the code in the complexity analysis. Indeed, the generation of the generator and parity check matrices is straightforward and can be done out of line. On the decoder side, both decoding algorithms apply. The first one is an iterative decoding on the parity check matrix. This decoding is fast and has a linear complexity on the dimension of the code O(k). When the iterative decoding fails, a ML decoding is used on the generator matrix. The complexity of this algorithm is in the general case O(k 2 ) row operations. However the decoding of LDPC-Band codes benefits from the band structure. Thanks to this structure, an optimized LU decomposition of the matrix [START_REF] Golub | Matrix Computations[END_REF] leads to a complexity of O(kB) row operations. This means that for k = 2000 source symbols and a band size of B = 200 symbols, the theoretical speed gain compared to a classical ML decoding will be about one order of magnitude. The complexity obtained here has to be compared with the complexity of O(k 3/2 ) row operations of the Studholme and Blake approach.

However, in practice, the size of the matrices that need to be inverted is lower than k. There are three reasons: (1) whenever a source symbol is received, the corresponding row in the generator matrix is removed. Then (2), only the columns of the non-systematic part that are received are used. It means that for a code rate r the system that has to be inverted is only in the order of (1r)k × (1r)k. Finally (3) the iterative decoding may have rebuilt some missing symbols, thereby reducing the system size.

B. Simulation Results

We carried out several tests to assess the computation benefits of our proposal. These tests were obtained on a 4xIntel Xeon5120 @ 1.86GHz/4 GB RAM/Linux PC. We compared the LDPC-Band codes with a band width of 200 with two other codes:

1) the LDPC-Staircase codes, regular repeat accumulate codes standardized by the IETF [START_REF] Roca | Low Density Parity Check (LDPC) Staircase and Triangle Forward Error Correction (FEC) Schemes[END_REF]. Following the optimisation for the hybrid iterative/ML decoding proposed in [START_REF] Cunche | Optimizing the error recovery capabilities of ldpc-staircase codes featuring a gaussian elimination decoding scheme[END_REF], the degree of the source symbol nodes is set to 5. 2) the random windowed codes proposed in [START_REF]Windowed erasure codes[END_REF]. These codes can be seen as a non systematic LDGM codes with 2log(k) elements per column. Because it is a non-systematic code, iterative decoding cannot be used (unlike standard LDGM codes). Figure 3 shows the average decoding bitrate as a function of the loss probability of the channel for the various code. The LDPC-Band and the LDPC-Staircase performances can be divided into two parts. In the region where the loss probability is low, the iterative decoder is sufficient to recover the losses and the bitrate is high. In the region of high loss probability the ML decoder is needed and the bitrate is reduced. On the opposite, the decoding speed of Windowed code does not depend on the loss probability, since it can only be decoded with a ML decoder.

We see that when the iterative decoding is used, the LDPC-Staircase and the LDPC-Band present comparable decoding speed (with a slight advantage to LDPC-Staircase codes), and they are both significantly faster than Windowed codes (that need ML decoding).

When ML decoding is needed, LDPC-Band codes outperform the other two codes: they are roughly 2 times faster than the Windowed codes, and roughly 4 times faster than LDPC-Staircase codes. Let us remind that the ML decoding of the LDPC-Band and the LDPC-Staircase is equivalent to solve a linear system of size (1r)k × (1r)k, whereas the ML decoding of the windowed erasure codes requires to solve a system of size k ×k (being non systematic, no source symbols is received and the system still has k unknown variables).

Table III present the average ML decoding time for the three compared codes as a function of the number of source symbols k. The LDPC-Band is always faster than the two other codes. Furthermore the ML decoding speed falls more quickly with the two other codes than with the LDPC-Band. The decoding speed with an iterative decoder are shown in Table IV

V. CONCLUSION

In this paper we proposed a flexible scheme that allows the construction of generator and parity check matrices and their efficient hybrid decoding. The iterative decoding on the parity check matrix provides a fast way to recover the source symbols, as decoding complexity is linear. When the environment is too harsh and the loss rate is close to the recovery capability of the codes, a Maximum-Likelihood decoding is applied on the structured generator matrix. Thanks to the band structure of this matrix, the ML decoding complexity is reduced to O(kB) row operations. This property allows the hybrid decoder to sustain high decoding speeds, even in high loss environments.

Thanks to the polynomial representation of the matrices, we presented a practical way of building the matrix of such codes. Our results shows that our LDPC-band codes match the performances of LDPC-staircase codes in terms of iterative decoding speed, while being much faster than these codes when ML decoding is required. Furthermore the LDPC-band codes exhibit erasure recovery capabilities close to standard repeat-accumulate codes using iterative decoding, and are close to the channel capacity when decoded with a ML scheme. In future works, a density evolution study should enable us to further improve the erasure correction capabilites of the iterative decoding scheme, without degrading the ML decoding capability.

  This work was supported by the French ANR grant No 2006 TCOM 019 (CAPRI-FEC project).

Fig. 1 .

 1 Fig. 1. The parity check matrix of the code with a band structure

Fig. 2 .

 2 Fig. 2. Design of the banded matrix of the non-systematic part with side effect removed

Fig. 3 .

 3 Fig. 3. Bitrate W.R.T loss probability, k=2000 code rate=1/2 and symbol size=1024 bytes

TABLE IV AVERAGE

 IV . The LDPC-Band is slower than the LDPC-Staircase but remains very close to it under iterative decoding. ITERATIVE DECODING TIME AS A FUNCTION OF THE OBJECT SIZE, FOR A CODE RATE=1/2 AND SYMBOL SIZE = 1024 BYTES.

	k	1,000	2,000	4,000
	LDPC-Band	326 Mbps	235 Mbps 150 Mbps
	LDPC-Staircase 125 Mbps	60 Mbps	30 Mbps
	Windowed Erasure	220 Mbps	120 Mbps	68 Mbps
		TABLE III
	AVERAGE ML DECODING TIME AS A FUNCTION OF THE OBJECT SIZE, FOR
	A CODE RATE=1/2 AND SYMBOL SIZE = 1024 BYTES.
	k	1,000	2,000	4,000
	LDPC-Band 1100 Mbps	1050 Mbps	900 Mbps
	LDPC-Staircase 1300 Mbps	1200 Mbps	1000 Mbps


Matthieu Lapeyre
email: matthieu.lapeyre@inria.fr

Pierre Rouanet

Pierre Yves Oudeyer

Poppy: a New Bio-Inspired Humanoid Robot Platform for Biped Locomotion and Physical Human-Robot Interaction

Introduction

Humanoid robots are predicted to play a key role in our everyday lives in the future. Yet, several challenges need to be addressed before this becomes a reality. In particular, biped robots need to be able to locomote robustly in human environments, which includes the ability to keep stability when unpredictable physical contact with humans happens. At the same time, these robots need to be capable of rich and safe social and physical interaction with humans, and to adapt to the behaviour and preferences of each particular user.

We should not only try to solve these challenges through artificial cognitive intelligence but also through body intelligence [START_REF] Pfeifer | How the Body Shapes the Way We Think: A New View of Intelligence[END_REF]. On one hand, a way to permit robots to adapt their behaviors to unknown environments is to provide them with learning algorithms based on social guidance [START_REF] Billard | Robot programming by demonstration[END_REF], or on autonomous self-exploration mechanisms [START_REF] Baranes | Active learning of inverse models with intrinsically motivated goal exploration in robots[END_REF][START_REF] Lapeyre | Maturational constraints for motor learning in high-dimensions: the case of biped walking[END_REF]. On the other hand, a part of the computation needed for such adaptation could also be done through the intrinsic mechanics and electronics of the robot, thus providing effective and hyper-responsive reactions while simplifying the algorithms of the different behaviors. Actually, body intelligence, realizing morphological computation, has been argued to facilitate and guide considerably the learning and development of sensorimotor skills [START_REF] Oudeyer | Intrinsically Motivated Learning of Real-World Sensorimotor Skills with Developmental Constraints[END_REF].

No current platform (Nao [START_REF] Gouaillier | The nao humanoid: a combination of performance and affordability[END_REF], Darwin Op [START_REF] Ha | Development of open humanoid platform darwinop[END_REF], Nimbro Op [START_REF] Schwarz | Nimbro-op humanoid teensize open platform[END_REF], HRP-2, ...) does offer both a bio-inspired morphology optimized for both walking and human physical interaction capability (safe, compliant, playful). To address these challenges, while exploring further approaches to the one elaborated for the Acroban robot [START_REF] Ly | Bio-inspired vertebral column, compliance and semi-passive dynamics in a lightweight humanoid robot[END_REF], we have designed a new bio-inspired humanoid robotic platform, called Poppy, which provides some of the software and hardware features needed to explore physical and social interaction together with biped locomotion for personal robots. It presents the following main features:

• Design inspired from the study of the anatomy of the human body and its bio-mechanic, with an emphasis on leg structure;

• Dynamic and reactive: we try to keep the weight of the robot as low as possible (optimized geometry of the pieces and smaller motors);

• Human interaction: screen for communication and compliant physical interaction;

• Practical platform: low cost, ease of use and easy to reproduce through rapid prototyping techniques;

2 Poppy

Poppy (Figure 1) is a humanoid robot, 84cm tall, and 3 kg. It has a large sensorimotor space including 25 Robotix motors (including MX-28 permitting dynamic compliance control), force sensors under its feet and some extra sensors in the head: 2 HD-wide angle-cameras, a stereo-micro and an inertial central unit (IMU 9DoF) plus a large LCD Screen (4 inch) for visual communication (e.g. emotions, instructions or debug).

The poppy morphology is designed based on important aspects of the actual human body. This inspiration is expressed in the whole structure (e.g. the limb proportions) and in particular in the trunk and legs. The weight of all limbs has been minimized thanks to sophisticated girder like design, which was made possible by the use of 3D printing techniques (all limbs were 3D printed).

Poppy uses the bio-inspired trunk system introduced by Acroban [START_REF] Ly | Bio-inspired vertebral column, compliance and semi-passive dynamics in a lightweight humanoid robot[END_REF]. These five motors allow it to reproduce the main DOFs of the human spine [START_REF] Ceccato | Comparison of trunk activity during gait initiation and walking in humans[END_REF]. This feature allows the integration of more natural and fluid motion while improving the user experience during physical interactions. In addition, the spine plays a fundamental role in bipedal walking and postural balance by actively participating in the balancing of the robot.

The legs were designed to increase the stability and agility of the robot during the biped walking by combining bio-inspired, semi-passive, lightweight and mechanicalcomputation features.The architecture of the hips and thighs of Poppy uses biomechanical principles existing in humans. The human femur is actually slightly bent at an angle of about 6 degrees. In addition, the implantation of the femoral head in the hip is on the side. This results in a reduction of the lateral hip movement needed to move the center of gravity from one foot to another and a decrease in the lateral falling speed. In the case of Poppy, the inclination of its thighs by an angle of 6 degrees causes a gain of performance of more than 30% for the two above mentioned points. Another example is Poppy's feet. Poppy has the particularity of having small feet compared to standard humanoids. It has humanly proportioned feet (ie about 15% of its total size). It is also equipped with toe joint actuated by compliance (see Figure 2.a). We believe that these are two key features to obtain a human-like and efficient walking gait. However, that raises problems regarding balance because the support polygon is reduced. We included pressure sensors under each foot in order to get accurate feedback of the current state of the robot (see Figure 2.b).

Open questions

In our current work, we explore the combination of both a bio-inspired body and bio-inspired learning algorithms. We are currently working on experiments involving Poppy to perform skill learning. First we aim at achieving an effective postural balance using the articulated spine, the feet pressure sensors and the IMU. Then, when this is learnt and used as a building block, we will experiment on the learning of biped walking usingonline and developmental learning algorithms such as [START_REF] Lapeyre | Maturational constraints for motor learning in high-dimensions: the case of biped walking[END_REF] and [START_REF] Degris | Off-policy actor-critic[END_REF]. We are expecting to clearly reduce the learning time needed and increase the quality of the learned skills thanks to the bio-inspired morphology of Poppy.

Figure 1 :Figure 2 :

 12 Figure 1: a. Global view of the Poppy platform. b. Zoom on legs design

We are also studying social physical interactions with non-expert users. We plan to conduct user studies to evaluate how playful physical interactions and emotions could improve learning in robotics. We think that the poppy platform could be very suitable for such studies.