H Mameche
email: hamza.mameche@gipsa-lab.fr.

E Witrant

C Prieur

Nonlinear PDE-Based

come

I. INTRODUCTION

The overall objective of controlling the plasma in tokamaks is to steer it towards a desired operation point defining a plasma scenario. A baseline for high-performance scenarios for tokamaks is H-mode (High confinement mode), where the plasma is strongly heated, exceeding a threshold above which the transport of plasma energy in the edge area is reduced and results in a transport barrier. In this mode, the energy confinement time is significantly enhanced, typically by a factor of 2 or more [START_REF] Keilhacker | H-mode confinement in tokamaks[END_REF]. To reach and maintain these high-performance scenarios, plasma profiles control plays a fundamental role. These profiles are: magnetic radial profiles (such as the poloidal magnetic flux Ψ(x), the safety factor q(x) or its inverse ι(x), where x is the location along the small plasma radius), and kinetic profiles, such as electrons and ions temperature and density. These are spatially distributed profiles with two-time scales coupled nonlinear dynamics, hence the difficulty of the control problem. Different approaches have been made to tackle these difficulties. In [START_REF] Moreau | A two-time-scale dynamicmodel approach for magnetic and kinetic profile control in advanced tokamak scenarios on JET[END_REF] and following works, a two-time scale linearized data-driven model was built based on singular perturbations theory, and was used to control the poloidal magnetic flux Ψ(x) and the safety factor q(x) in [START_REF] Moreau | Integrated magnetic and kinetic control of advanced tokamak plasmas on DIII-D based on data-driven models[END_REF]. Other works use first-principles-driven models that capture the dominant and relevant dynamics to synthesize the controller as in [START_REF] Witrant | A control-oriented model of the current profile in tokamak plasma[END_REF] [START_REF] Felici | Real-time control of tokamak plasmas: from control of physics to physics-based control[END_REF]. Taking into account the spatially distributed nature of the dynamics, [START_REF] Boyer | First-principles-driven model-based current profile control for the DIII-D tokamak via LQI optimal control[END_REF] used a spatially discretized model for current profile control, while other works used infinite dimensional theory to synthesize the control algorithms as in [START_REF] Bribiesca Argomedo | Lyapunov-based distributed control of the safety-factor profile in a tokamak plasma[END_REF].

While some previous works emphasized on the coupling between the dynamics and used a linearized model as in [START_REF] Mavkov | Distributed control of coupled inhomogeneous diffusion in tokamak plasmas[END_REF] and the experimental validation in [START_REF] Mavkov | Experimental validation of a Lyapunov-based controller for the plasma safety factor and plasma pressure in the TCV tokamak[END_REF], we are interested in the nonlinearity of the dynamics. In this tokamak control perspective, control of nonlinear PDEs is a recent and challenging topic. An interesting way to tackle the problem is to separate a "slow" finite set of eigenvalues of the system operator that capture the dominant dynamics, and use it as a basis to synthesize the finite-dimensional controller [START_REF] Christofides | Nonlinear and robust control of PDE systems: Methods and applications to transport-reaction processes[END_REF]. But in the case of nonlinear parabolic PDEs, a precise approximation of the dynamics may lead to a large number of modes that should be included [START_REF] Christofides | Control of nonlinear distributed process systems: Recent developments and challenges[END_REF], and this is what motivates us to choose the infinite-dimensional framework to tackle the problem. As in [START_REF] Bribiesca Argomedo | Lyapunov-based distributed control of the safety-factor profile in a tokamak plasma[END_REF], we use a Lyapunov control function approach to synthesize the controller that increases the convergence rate of the closed-loop system. This paper is organized as follows: in Section II, we present the dynamical model for the electron temperature profile including the transport model modification to represent the H-mode pedestal. In Section III, we use the direct Lyapunov method for the stability analysis of the PDE, then we use the Sum of Squares framework to verify the positivity of the differential matrix inequalities emerging from the Lyapunov analysis to check the nominal stability of the system. We then propose a control algorithm to increase the convergence rate of the closed-loop system. Finally, in Section IV we use the RAPTOR simulator [START_REF] Felici | Real-time control of tokamak plasmas: from control of physics to physics-based control[END_REF] to evaluate the control strategy including the input constraints.

The main variables definitions are given in Table I.

II. SYSTEM DESCRIPTION AND CONTROL PROBLEM

A. Electron temperature dynamics

In this study we focus on the electron temperature T e dynamics, modeled by a diffusion equation under the infinite cylinder geometry hypothesis (the transport is symmetric in the toroidal and the poloidal directions). The only space variable is the cylinder radius ρ and x = ρ/a is the normalized radial variable (a is the small plasma radius). The transport dynamics is modeled as:

3 2 ∂ (n e T e ) ∂t = 1 a 2 1 x ∂ ∂ x (xn e χ e ∂ T e ∂ x ) -P sinks + P sources , (1)
with boundary and initial conditions: ∂ T e ∂ x (0,t) = 0, T e (1,t) = T e,edge (t), ∀t ≥ t 0 and T e (x,t 0 ) = T 0 (x), ∀x ∈ [0, 1].

In this model, n e (x,t) is the electron density, χ e (x,t) is the electron heat diffusivity, T e,edge (t) is the edge temperature and P sources (x,t) = P OH (x,t) + P aux (x,t) is the supplied heating power density. P OH (x,t) is the power density due to ohmic effect and P aux (x,t) is the auxiliary heating sources. P sinks (x,t) represents the lost power density such as electronion equipartition losses P ei (x,t) and radiation losses, and is neglected as in [START_REF] Felici | Real-time control of tokamak plasmas: from control of physics to physics-based control[END_REF].

B. Electron heat diffusivity model

Because of the complexity of the heat transport dynamics, there is no fully analytic model for the heat diffusivity model. Instead, semi-empirical models have been proposed and tested on experimental data. We choose to use the modified bohm/gyrobohm model in [START_REF] Pianroj | Simulations of H-mode plasmas in tokamak using a complete core-edge modeling in the BALDUR code[END_REF] for χ e :

               χ e = χ e c × f s , χ e c = (2χ Be + χ gBe ) f s , χ Be = 4 × 10 -5 R ∇(n e T e )
n e B φ 0 q 2 T e,0.8 -T e,1.0 T e,1.0 ,

χ gBe = 5 × 10 -6 √ T e ∇T e B 2 φ 0 , f s = 1 1+k ω E×B γ IT G 2 × 1 max(1,(s-s thres ) 2 ) (2)
where, χ e c is the classical bohm/gyrobohm model, χ Be is the Bohm diffusivity, χ gBe is the gyro-Bohm diffusivity. B φ 0 is the toroidal magnetic field at the center, R is the major radius, q is the safety factor and T e,0.8 represents the electron temperature at x = 0.8 and T e,1.0 at x = 1. The term ((T e,0.8 -T e,1.0 )/T e,1.0 ) tackles the non-local dependence by representing the phenomena in which the diffusivity increases when the edge temperature is decreased, and vice versa. s is the magnetic shear, ω E×B is the flow shearing rate, γ IT G is the growth rate of ion temperature gradient (ITG), and k a coefficient.

In order to represent the pedestal behaviour corresponding to the H-mode, the classical Bohm/gyro-Bohm model χ e c is then multiplied by f s : the suppression function for the electron thermal diffusivity. In f s , the first term represents the flow shearing rate through ω E×B and the reduction of the turbulence growth rate through the growth rate of ion temperature gradient γ IT G [START_REF] Sugihara | Simulation studies on H-mode pedestal behavioru during type-I ELMs under various plasma conditions[END_REF]. The second term reduces the transport only in the region where the magnetic shear s exceeds a specified threshold s thres . To simplify the incorporation of this modification to the diffusivity model, the suppression function f s was taken as polynomial approximation of the experimental results obtained in [START_REF] Pianroj | Simulations of H-mode plasmas in tokamak using a complete core-edge modeling in the BALDUR code[END_REF].

C. Control problem

To formulate our control system from the electron temperature dynamics in II-A and II-B, we consider these assumptions: Assumption 1. Because of the slow and small time variations of the q profile compared to the T e one, and in order to decouple the two dynamics in χ Be (2), we choose q as a time-averaged polynomial approximation of a nominal q(x,t) profile simulated in RAPTOR. Similarly, L T e is a time average for ((T e,0.8 -T e,1.0 )/T e,1.0 ) in (2), and n e is taken as a time-fixed, line-averaged density ne which is a classical variable for tokamaks, this assumption is made to focus on the electron temperature dynamics and its nonlinear coupling with the electron heat diffusivity. Assumption 2. To deal with the terms with absolute value in χ Be and χ gBe (2), we assume that ∂ T e ∂ x ≤ 0 almost always and everywhere. This assumption is verified experimentally given that auxiliary power sources target mainly the central region (x ≤ x inv ), and act with a wide angle when deposited on the core region (x 1 ≤ x ≤ x ped ), (x inv ≈ 0.45 and x ped ≈ 0.9 in TCV tokamak H-mode plasmas [START_REF] Kim | Simple predictive electron transport models applied to sawtoothing plasmas[END_REF]). We also assume that T e,edge is very small compared to the temperature in the center of the plasma and we can consider it as zero.

We get the control system described by a nonlinear diffusion PDE:

∂ T e ∂t = A x ∂ ∂ x x B(x) +C(x) √ T e ∂ T e ∂ x 2 +u, ∀x ∈ [0, 1] (3) with boundary conditions, ∂ T e ∂ x (0,t) = 0, T e (1,t) = 0, ∀t ≥ t 0 (4)
where:

A = 2 3a 2 , B(x) = -8×10 -5 RL Te B φ 0 q 2 (x) f s (x), C(x) = -5×10 -6 B 2 φ 0 f s (x).
The control problem is to prove the nominal stability of an arbitrary equilibrium point ( ū(x), Te (x)) of the system (3)-( 4), and to synthesize a control strategy to ensure the tracking of a reference electron temperature profile T e,re f .

III. STABILITY ANALYSIS AND DISTRIBUTED CONTROL

A. Stability of the open-loop system

In this section, we develop a Lyapunov function for the nonlinear diffusion PDE in (3)-( 4) to prove the nominal stability of the equilibrium ( ū(x), Te (x)). To do that, we consider the following Lyapunov function candidate:

V 1 (T e ) = 1 2 1 0 xP T e (x)(T e -Te ) 2 dx (5)
where Te is the steady-state equilibrium resulting from a constant input ū. P T e (x) is a weighting function with P T e (x) > 0 for all x ∈ [0, 1] to ensure the positivity of V 1 (T e ). The Lyapunov function candidate is (xP T e (x))-weighted L 2 ([0, 1]) norm squared. We multiply by x to handle the singularity at x = 0 when we differentiate V 1 (T e ). Using a similar development as in [START_REF] Valmorbida | Stability analysis for a class of partial differential equations via semidefinite programming[END_REF], we present the main result of the paper.

Theorem 1: Suppose that for a given positive number α 1 , there exist a polynomial P T e (x) and a 5 × 5 symmetric polynomial matrix H(x), such that P T e (x) > 0 for all x ∈ [0, 1], H(0) ≥ 0, H 1,1 (1) ≤ 0 and:

F(x) + H(x) ≥ 0, ∀x ∈ [0, 1] (6)
where F(x) + H(x) is defined in [START_REF] Mavkov | Distributed control of coupled inhomogeneous diffusion in tokamak plasmas[END_REF] and its elements are defined in [START_REF] Mavkov | Experimental validation of a Lyapunov-based controller for the plasma safety factor and plasma pressure in the TCV tokamak[END_REF].

Then the time derivative V1 of V 1 defined in ( 5) along the solutions of ( 3)-( 4) verifies:

V1 (T e ) ≤ -α 1 V 1 (T e ) + 1 0 xP T e (x)(T e -Te ) ũ dx (7)
where ũ is defined as ũ = uū.

F(x) + H(x) =                       • • • • • • • • • • • • • f 9 + h 1 • • • • • • • • • • • • f 7 + h 10 h 2 • • • • • • • • • • • 0 h 11 h 3 • • • • • • • • • • f 8 + h 21 h 18 h 12 f 7 + h 4 • • • • • • • • • f 5 h 22 0 h 13 h 5 • • • • • • • • 0 0 f 5 0 0 0 | • • • • • • h 28 f 6 h 26 h 23 h 19 f 8 + h 14 h 6 • • • • • • f 1 h 29 0 f 6 h 24 f 5 h 15 h 7 • • • • • 0 f 3 f 1 0 0 f 6 0 f 5 0 • • • • 0 0 0 f 3 0 0 0 0 0 0 • • • h 32 f 4 f 2 3 4 h 31 h 30 f 1 h 27 h 25 h 20 h 16 f 8 + h 8 • • 0 h 33 0 f 4 f 2 h 31 f 3 f 1 2h 27 f 6 h 17 h 9 • 0 0 0 0 0 f 4 f 2 0 f 3 f 1 0 f 6 f 5 0 0 0 0 0 0 0 f 4 0 0 f 3 0 0 0                       (8)
Proof: The time derivative of the Lyapunov function ( 5) along ( 3)-( 4) is:

V1 (T e ) =
V1 (T e ) + α 1 V 1 (T e ) = - 1 0 xAP T e (x) B(x) +C(x) √ T e ∂ T e ∂ x 3 dx - 1 0 xAP T e (x)T e B(x) +C(x) √ T e ∂ T e ∂ x 2 dx + 1 0 xA(P T e (x) Te + P T e (x) Te ) B(x) +C(x) √ T e ∂ T e ∂ x 2 dx + α 1 2 1 0 xP T e (x)(T e -Te ) 2 dx + 1 0 xP T e (x)(T e -Te )( ū + ũ) dx = - 1 0 G(x, √ T e , D 1 T e ) dx + 1 0 xP T e (x)(T e -Te ) ũ dx (10)
At this stage we perform a change of variable τ = √ T e , to ensure the polynomial dependence of G on the dependent variable and its derivatives. We can then write G in terms of τ as:

G(x, D 1 τ) = ξ (D 1 τ) T F(x)ξ (D 1 τ) (11)
where F(x) is defined in ( 8)-( 9) by only taking the f i elements, and:

ξ (D 1 τ) := 1, τ, ∂ τ ∂ x , τ 2 , τ ∂ τ ∂ x , ∂ τ ∂ x 2 , τ 3 , τ 2 ∂ τ ∂ x , τ ∂ τ ∂ x 2 , ∂ τ ∂ x 3 , τ 4 , τ 3 ∂ τ ∂ x , τ 2 ∂ τ ∂ x 2 , τ ∂ τ ∂ x 3 T
Note that:

1 0 ξ (D 1 τ) T F(x)ξ (D 1 τ) dx = 1 0 ξ (D 1 τ) T F(x)ξ (D 1 τ) + d dx µ(τ) T H(x)µ(τ) dx -µ(τ) T H(x)µ(τ) 1 0 (12)
where: µ(τ) := 1, τ, τ 2 , τ 3 , τ 4 T , and H(x) is a 5 × 5 symmetric polynomial matrix. Using ( 12) and since we have that: µ(τ(1)) = [1, 0, 0, 0, 0] T , we rewrite [START_REF] Christofides | Nonlinear and robust control of PDE systems: Methods and applications to transport-reaction processes[END_REF] as: .

V1 (T e ) + α 1 V 1 (T e ) = - 1 0 ξ (D 1 τ) T (F(x) + H(x))ξ (D 1 τ) + H 1,1 (1)
h 1 = dH 1,1 dx , h 2 = dH 1,2 dx , h 3 = H 1,2 , h 4 = 2 3 dH 1,3 dx , h 5 = H 1,3 , h 6 = 1 2 dH 1,4 dx , h 7 = H 1,4 , h 8 = 2 5 dH 1,5 dx , h 9 = H 1,5 , h 10 = 2 3 dH 1,3 dx + dH 2,2 dx , h 11 = H 1,3 + H 2,2 , h 12 = 1 2 dH 1,4 dx + dH 2,3 dx , h 13 = H 1,4 + 2H 2,3 , h 14 = 2 5 dH 1,5 dx + dH 2,4 dx , h 15 = H 1,5 + 3H 2,4 , h 16 = dH 2,5 dx , h 17 = 4H 2,5 , h 18 = H 1,4 + H 2,3 , h 19 = H 1,5 + H 2,4 , h 20 = H 2,5 , h 21 = 2 5 dH 1,5 dx + dH 3,3 dx , h 22 = H 1,5 + 2H 3,3 , h 23 =
where H(x) is defined in ( 8)-( 9) by only taking the h i elements in [START_REF] Mavkov | Distributed control of coupled inhomogeneous diffusion in tokamak plasmas[END_REF] and set the f i elements to 0.

Then if H(0) ≥ 0, H 1,1 (1) ≤ 0 and F(x) + H ≥ 0, ∀x ∈ [0, 1], we get that:

V1 (T e ) ≤ -α 1 V 1 (T e ) + 1 0
xP T e (x)(T e -Te ) ũ dx Remark 1. Notice that unsimilarly to [START_REF] Valmorbida | Stability analysis for a class of partial differential equations via semidefinite programming[END_REF] where in a system model ∂W ∂t = F(x, D β W ), F is polynomial on its second argument, it is not the case for our model, therefore we had to apply integration by part followed by a change of variable to obtain G(x, D 1 τ) to which we can apply the rest of the procedure.

Remark 2. As a result of the theorem, by setting ũ = 0, we get that the equilibrium ( ū, Te ) is exponentially stable in V 1 .

B. Calculation of the weighting Function

Since F(x) contains P T 2 (x) and its derivative, and H(x) contains continuously differentiable functions and their derivatives, ( 6) is a differential matrix inequality. Using the theorem assumption that P T e (x) and the elements of H(x) are polynomials in x, it is possible to formulate the positivity of F(x) + H(x) in x ∈ [0, 1] as a convex optimization problem in the form of Semidefinite Programming using the following corollary of the Putinar's Positivstellensatz theorem [START_REF] Valmorbida | Semi-definite programming and functional inequalities for distributed parameter systems[END_REF]: Corollary 1.1: If there exists N(x) ∈ Σ 14×14 [x] (Sum Of Squares polynomial matrix of order 14) such that:

F(x) + H(x) -N(x)x(1 -x) ∈ Σ 14×14 [x] (13)
then ( 6) holds.

In this context, the problem of finding P T (x) is formulated as the following feasibility Sum Of Squares problem (SOSP): Find P T e (x), H(x), N(x). Sub ject to :

P T e (x) > 0 in [0, 1], F(x) + H(x) -N(x)x(1 -x) ∈ Σ 14×14 [x], N(x) ∈ Σ 14×14 [x], H(0) ≥ 0, H 1,1 (1) ≤ 0.
This problem is then solved using Yalmip with the Sum Of Squares Module [START_REF] Löfberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF], and the resulting P T e (x) is a decreasing polynomial strictly positive on [0, 1].

C. Distributed control

Based on the previous analysis, and to perform reference tracking control, we define the control strategy while ensuring exponential stability of the closed-loop system as follows:

Corollary 1.2: If the conditions of Theorem 1 are verified, we choose the control input u ctrl = ū + ũ, where ũ is calculated to verify the equality:

1 0 xP T e (x)(T e -Te ) ũ dx = -α 2 V 1 (T e ) (14)
with α 2 > 0, a tuning parameter. With this control we get:

V1 (T e ) ≤ -(α 1 + α 2 )V 1 (T e )
and the system (3)-( 4) with ( 14) is exponentially stable in V 1 with convergence rate α 1 + α 2 . An explicit control law from ( 14) is the proportional controller:

u ctrl = ū -α 2 2
(T e -Te ).

IV. CONTROL IMPLEMENTATION AND RESULTS

Now we aim to implement the control strategy on the RAPTOR simulator with TCV tokamak settings. RAPTOR is a lightweight code used to simulate simplified nonlinear plasma transport physics, and is also used as a real-time state observer for the TCV tokamak. We first show the shape constraints on the input, then we present the RAPTOR configuration used for the simulation, we formulate an optimization problem to find the engineering parameters to apply as inputs to the system, and finally we present the results of the simulation where we test the performance and robustness of the controller by adding disturbances and time-delays.

A. Input constraints

During the beginning of the discharge, the plasma heating comes from the induced plasma current I p whose time evolution scenario is computed offline and defines the discharge phases (ramp-up, flat-top and rump-down). The plasma current is thus not considered as an input to the system (a feedback loop sets the voltage on the poloidal field coils to generate the desired I p ). On the other hand, the auxiliary heating sources such as Electron Cyclotron Resonance Heating (ECRH), Electron Cyclotron Current Drive (ECCD) and NBI are used for the online distributed control. They are subject to profile shape constraints (see Fig. 1), and their power densities are approximated by weighted Gaussian distributions, for the case of ECRH/ECCD antennas, we have for the ith actuator [START_REF] Felici | Real-time control of tokamak plasmas: from control of physics to physics-based control[END_REF]:

P aux,i (x,t) = P i (t) e z a 0 e z V dx, z = -4(x -x dep,i ) 2 w 2 dep,i (15)
where w dep is the deposit width and x dep is the location of the peak of the deposit. In most cases the x dep and w dep are fixed, and the amplitude of the power density for the i-th actuator P i (t) in ( 15) is used as an input in the model and:

P aux (x,t) = ∑ i P aux,i (x,t) (16)

B. RAPTOR configuration

We use a configuration where the plasma current ramps up from 80 kA to 120 kA and we use two EC antennas as auxiliary sources, with power amplitude P 1 (t) and P 2 (t), The first one is a pure ECRH heating source deposited at x dep = 0 with deposition angle of w dep = 0.35, and the second one is an EC current drive (ECCD) source with x dep = 0.4 and w dep = 0.35. The engineering inputs P 1 (t) and P 2 (t) are limited to 1 MW.

For the heat diffusivity χ e we choose the model introduced in [START_REF] Kim | Simple predictive electron transport models applied to sawtoothing plasmas[END_REF] (equation ( 13)), with the parameters used in the paper to fit experimental measurements of the T e profile in H-mode TCV plasma. In

C. The control algorithm

Taking into consideration the constraints on the actuators, we are not able to implement the control algorithm presented in Corollary 1.2. A practical implementation that we use is to calculate and apply at each time step the actuator inputs u ac = [P 1 , P 2 ] T that solve the squared difference minimization problem of ( 14 where (T e,re f , u re f ) are the reference temperature profile and its corresponding input, that is the imposed equilibrium point of the system in closed loop with the controller [START_REF] Kim | Simple predictive electron transport models applied to sawtoothing plasmas[END_REF].

D. Simulation results

We now present the simulation results with different scenarios in order to test the performance and the robustness of the control algorithm, and we compare the open-loop and the closed-loop responses. To test the behaviour of the controller when changing the operating point we used a 2-stage reference profile.

1) Adding disturbance:

The controller disturbance attenuation is tested by adding a third ECCD source at x dep = 0.2 with w dep = 0.35 and P 3 = 0.1 MW at t = 0.2 s. The results are shown in Fig. 3, where we see the time evolution of the T e profile in various space points, as well as the time evolution of u ac that solves the optimization problem [START_REF] Löfberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF]. We see that the tracking performance in closed loop is significantly improved compared to the open-loop case. The time response in closed loop is further reduced in the region where the actuators are most efficient (x dep = 0 and x dep = 0.4). Before t = 0.2 we see that the optimization routine retrieved the original values P 1 = 0.3 MW and P 2 = 0.3 MW with a short saturation, and after t = 0.2 we can notice that the controller could not totally compensate the effect of the disturbance due to the fact that the disturbance source is not collocated with the input ones. In Fig. 4, we show a comparison of the time evolution of the Lyapunov function in open loop and closed loop. We can see that the practical implementation of the controller succeeded to improve the convergence rate. 2) Adding time-delays: Because of the fast dynamics of the system, the computation and the transportation of the control signal is considered as a time-delay. To take it into account we include a time-delay of 5 ms in the control loop. In Fig. 5 we see the deterioration of the performance in presence of such time-delay. The time-delay induces an overshoot in the closed-loop response, more important in the plasma centre (where there is more actuation) but which stays within reasonable bounds. The proposed control strategy is thus reasonably robust to time-delays.

V. CONCLUSION

In this work, the stability analysis and control of the electron temperature profile in H-mode tokamak plasmas was addressed, with dynamics described by a first-principle control-oriented model.

The stability of the resulting nonlinear parabolic PDE was studied with the Lyapunov approach. The sum of squares framework was used to compute the weighting functions that ensure the positivity of the resulting integral inequalities.

A control strategy was proposed to ensure a good tracking of the electron temperature profile in closed loop at an increased convergence rate. To evaluate the control strategy, RAPTOR plasma simulator was used and the control input constraints were taken into account to derive the engineering control parameters. The simulation results show a good performance of the controller in tracking the H-mode TCV plasma electron temperature profile. The robustness of the controller was investigated with respect to input disturbances and by adding an time-delay.

1 0 1 0 1 0

 111 xP T e (x) ∂ T e ∂t (T e -Te ) dx = xP T e (x)(T e -Te ) parts we get: V1 (T e ) = AP T e (x)(T e -Te ) x B(x) +C(x) T e (x)(T e -Te )) xP T e (x)(T e -Te )u dx With m ∈ N * , define D m T e := T e , ∂ T e ∂ x , . . . , ∂ m T e ∂ x m T a vector of spacial partial derivatives up to order m. Considering the boundary conditions (4) we get:

1 0

 1 µ(τ(0)) T H(0)µ(τ(0)) + xP T e (x)(T e -Te ) ũ dx f 1 = 4 9 (ABxP T e ), f 2 = 2 3 (ACxP T e ), f 3 = 4 5 (ABxP T e ), f 4 = (ACxP T e ), f 5 = -4 9 (ABx(P T e Te + P T Te )), f 6 = -2 5 (ACx(P T e Te + P T Te )), f 7 = 1 3 (α 1 xP T e Te -xP T e ū), f 8 = -1 10 (α 1 xP T e ), f 9 = (xP T e Te ū -1 2 α 1 xP T e Te 2 )

dH 3, 4 dx

 4 , h 24 = 3H 3,4 , h 25 = dH 3,5 dx , h 26 = 2H 3,4 , h 27 = 2H 3,5 , h 28 = dH 4,4 dx , h 29 = 3H 4,4 , h 30 = dH 4,5 dx , h 31 = 4H 4,5 , h 32 = dH 5,5 dx , h 33 = 4H 5,5 .

2 Fig. 1 .Fig. 2 .

 212 Fig. 1. Auxiliary ECRH (P aux,1 ) and ECCD (P aux,2 ) sources distributions.

): u * ac = arg min u ac 1 0 2 dx ( 17 )

 1217 xP T e (x) (P aux (u ac ) -u re f ) + α 2 2(T e -T e,re f ) subject to: 0 ≤ P 1 ≤ 1MW, 0 ≤ P 2 ≤ 1MW.

Fig. 3 .

 3 Fig.3. T e tracking and time-evolution of P 1,2 when changing the set point and introducing a disturbance at t = 0.2 s.

Fig. 4 .

 4 Fig. 4. Open-loop and closed-loop convergence rates of the Lyapunov function.

Fig. 5 .

 5 Fig.5. T e tracking and time-evolution of P 1,2 when changing the set point and introducing a disturbance at t = 0.2 s in the presence of an input timedelay of 5 ms.

| | | | | | | | | | | | | | | | | | | | | | | | | | | •