domain: math phys

text:
Jacques Blum

Holger Heumann

AND Eric Nardon

Xiao Song

Automating the Design of Tokamak Experiment Scenarios

Keywords: optimal control, PDE-constrained optimization, plasma equilibrium, tokamak Automatisation contrôle optimal, optimisation sous contrainte EDP, équilibre du plasma, tokamak

The real-time control of plasma position, shape and current in a tokamak has to be ensured by a number of electrical circuits consisting of voltage suppliers and axisymmetric coils. Finding good target voltages/currents for the control systems is a very laborious, non-trivial task due to non-linear effects of plasma evolution. We introduce here an optimal control formulation to tackle this task and present in detail the main ingredients for finding numerical solutions: the finite element discretization, accurate linearizations and Sequential Quadratic Programming. Case studies for the tokamaks WEST and HL-2M highlight the flexibility and broad scope of the proposed optimal control formulation.

Introduction

Nuclear fusion is a highly exothermic reaction in which two light atomic nuclei fuse to form a heavier nucleus. The peaceful use of such reactions for energy production on earth is a multinational research e↵ort with high impact on the long-term perspective of energy production and consumption. The most promising technology to achieve this goal is currently the tokamak, a torus shaped reactor that uses strong magnetic fields to confine plasma and to achieve the extreme conditions to start the fusion reaction.

The mathematical modeling of nuclear fusion physics in such complicated devices is very diverse (see [17] for an detailed introduction, or [START_REF] Fasoli | Computational challenges in magnetic-confinement fusion physics[END_REF] and [START_REF] Ongena | Magnetic-confinement fusion[END_REF] for a recent overview). We have an extreme range of spatial and temporal scales. Spatial scales range from 10 5 meters for the electron orbit radius to 10 2 meters which is the length of the magnetic field lines in between the reactor walls. The temporal scales range from 10 11 seconds for one period of the electron orbit to 10 2 seconds for the resistive di↵usion time, that describes plasma evolution at the timescale of resistive di↵usion. We are dealing here with multi-scale problems and strong anisotropy. So, the modeling in nuclear fusion ranges from high dimensional kinetic models to completely di↵erent kind of fluid models with various levels of model reductions steps that reduce the high dimensionality but increase the complexity of the models. The reduced-order models are very non-standard and include to a very high degree empirical intuition.

Nevertheless, the operation of thermonuclear fusion experiments relies heavily on such reduced-order models as first principles approaches are far too expensive to address the many interesting and important questions in tokamak engineering: identification of plasma characteristics in real-time, fast pre-shot simulations within virtual tokamak frameworks, fast post-shot simulations within the same framework augmented by measurement data, development of experimental scenarios, design of new devices, improvement of heat load mitigation or assessing novel strategies to improve stability, to name just a few.

In this work we will give a contribution to the development of experimental scenarios. Besides many others an experimental scenario in a tokamak requires prescribed time series of voltages applied in the so-called poloidal field coils that have an axisymmetric geometry and control form and shape of the plasma inside the tokamak. The basis for designing such voltage time series is the numerical solution of the free-boundary equilibrium evolution problem (FBEE). The FBEE problem, as it will be detailed in Section 2, combines the classical Grad-Shafranov equation, describing the magnetohydrodynamic equilibrium in axisymmetric systems, with circuit and induction equations for the poloidal field coils and other passive structures inside the tokamak. For a given time series of voltages Ṽ (t) = (V 1 (t) . . . V n (t)) of n suppliers it allows to predict the evolution the plasma and to check whether a specific choice of voltages can ensure a certain scenario. For example, it can be checked if the voltages enable the transition from limiter to divertor configuration and how much current is accumulated in passive structures and coils, and whether this violates certain technological constraints or not. Hence, the design of voltage trajectories can be regarded as the inverse FBEE problem: We have a certain desired evolution of the plasma in a tokamak in mind and ask to find the voltage trajectories that ensure this. Solving inverse problems is challenging in general.

Here, the choice of correct voltages in the FBEE problem is by no means trivial, and can be a very laborious task. The huge variety of realizations of voltage times series and the non-linear nature of the problem make it very di cult to guess the correct voltages. Moreover, a plain trial and error approach is not possible as the plasma equilibria of interest are often unstable, and small variations of the voltages Ṽ (t) can lead to entirely di↵erent results.

We therefore introduce an optimal control approach to scenario design. As the primal unknown of the FBEE problem is the poloidal magnetic flux we encode the design goal in an objective functional C[ (t), Ṽ (t)] that is large when the design goal is violated and small otherwise. We define a convex regularization functional R[ Ṽ (t)] and formulate a constrained optimization problem [START_REF] Albanese | On the solution of the magnetic flux equation in an infinite domain[END_REF]

min (t), Ṽ (t) C[ (t), Ṽ (t)] + R[ Ṽ (t)]
subject to (t), Ṽ (t) verify the FBEE problem (see Section 2).

Candidates for objective and regularization functionals are

(2) C[ (t), Ṽ (t)] = 1 2 Z T 0 w(t) N desi X i=1 (r i (t), z i (t), t) (r 0 (t), z 0 (t), t) 2 dt and (3) R[ Ṽ (t)] = 1 2 Z T 0 n X j=1 n X i=1 D ij V i (t)V j (t)dt ,
with the weights w(t) non-negative and D positive definite. This choice of the objective functional forces (r, z, t) to be constant on a prescribed set of N desi + 1 points (r i (t), z i (t)) at each instant t and hence can be used to encode a certain desired evolution of the plasma. Other choices of objective functionals could include penalization of induced currents, voltages in suppliers or loop voltages or any other design goal that can be quantified in terms of the poloidal flux and the voltages V i (t). The regularization term R[V (t)] enables to work with a well-posed problem, i.e. a problem that is stable to perturbations on the data. The optimal control formulation (1) alone is only the basis for a powerful tool for scenario design. To put this into practice it is of paramount importance to specify appropriate discretization methods, leading to non-linear algebraic finitedimensional problems, that are then solved by appropriate iteration schemes. Simple gradient descent (see Appendix A), for example, will su↵er from the same di culties as a trial-and-error approach, as it requires to solve repetitively the nonlinear FBEE problem for a varying choice of voltages V i (t). For that reason, a sequential quadratic programming (SQP) approach (see Appendix A) is here the most appropriate choice. Moreover, it is important to formulate SQP for the discretized version of (1) in order to avoid convergence problems due to inaccurate derivatives. The discretization of (1) and the SQP formulation are outlined in Sections 3. 1-3.3. This work does not claim that (2) and ( 3) are the exclusive and optimal choices for scenario design but should highlight how di↵erent choices of objective functions in (1) can be useful to tackle di↵erent tasks in scenario design. We are convinced that robust and easily accessible implementations of (1) will become invaluable tools for designing and improving control for tokamaks. The optimal solutions ⇤ (t) and Ṽ⇤ (t) of the optimal control problem (1) are not only solutions of the FBEE, but for a reasonable choice of the objective functional, the solutions are also close to the design goals. So, if these design goals strongly penalize unstable behavior (e.g. vertical displacement event, strong induced currents), then we have found a voltage times series Ṽ⇤ (t) that provides a plasma equilibrium evolution, even though the equilibrium might pass through physical unstable states. So with this optimal control approach (1) we can separate the problem of finding good feedforward control from the problem of designing e↵ective feedback control as we do not need to include feedback control into our simulation tool to avoid physical instabilities during the simulation. We will illustrate this line of thought with extensive numerical examples in detail in Section 5.

Even though we can get in principle solutions for any kind of objective functional C[ (t), Ṽ (t)], we would like to stress that the whole power of the approach is hidden in an appropriate choice of the objective functional: an objective functional encoding two competing goals will lead to disappointing solutions and an objective functional not penalizing unstable behavior can lead to unstable optimal solutions. But clearly, deciding whether a cost functional encodes competing goals or whether it penalizes su ciently strong unstable behavior is not easy and requires a good knowledge in both the physical background and the mathematics of optimal control. Optimal control is not an auto pilot, but is capable to leverage considerably the tokamak physicists intuition.

To illustrate the power of the proposed optimal control approach to scenario design we developed an extension of the MATLAB/Octave library FEEQS.M. This library is based on the methods for axisymmetric free boundary plasma equilibria that where described in [START_REF] Blum | Numerical simulation and optimal control in plasma physics[END_REF] and [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF] and utilizes in large parts vectorization. Therefore, the running time is comparable to C/C++ implementations. FEEQS.M is publicly available 1 and a forthcoming release will contain the methods introduced 1 http://www-sop.inria.fr/members/Holger.Heumann/Software.html in this work. The Section 4 contains validation and verification results for these new developments and the Section 5 presents applications for the tokamaks WEST [START_REF] Bourdelle | West physics basis[END_REF] in Cadarache, France and HL-2M [START_REF] Li | The component development status of HL-2M tokamak[END_REF] in Chengdu, China. We can also refer the related [START_REF] Nardon | Magnetic configuration and plasma start-up in the west tokamak[END_REF], where the same approach was used to identify for WEST configurations with regions of low magnetic fields as it is required for a successful start-up.

Free-Boundary Plasma Equilibrium Evolution

The essential equations for describing plasma equilibrium in a tokamak are force balance, the solenoidal condition and Ampère's law that read respectively [START_REF] Bonnans | Numerical Optimization: theoretical and practical aspects : with 26 figures[END_REF] grad

p = J ⇥ B , div B = 0 , curl 1 µ B = J ,
where p is the plasma kinetic pressure, B is the magnetic field, J is the current density and µ the magnetic permeability. In the quasi-static approximation these equations are augmented by Faraday's law in all other conducting structures, and by Ohm's laws in coils and passive structures. Hence for the resistive timescale the plasma is in equilibrium and (4) holds at each instant of time [START_REF] Grad | Classical di↵usion in a tokamak[END_REF].

Under the assumption of perfect axial symmetry, it is convenient to put (4) in a cylindrical coordinate system (r, ', z) and to consider only a meridian section of the tokamak. The primal unknowns are the poloidal magnetic flux = (r, z), the pressure p = p( ) and the poloidal current flux f = f ( ). The poloidal magnetic flux := rA • e ' is the scaled toroidal component of the vector potential A, i.e. B = curl A and e ' the unit vector for '. The poloidal current flux f = rB • e ' is the scaled toroidal component of the magnetic field. We introduce H = [0, 1] ⇥ [ 1, 1], the positive half plane, to denote the meridian plane that contains the cross section of the tokamak device. The geometry of the tokamak determines the various subdomains (see Figure 1):

• the domain F ⇢ H corresponds to those parts that are made of the ferromagnetic circuit,i.e. the iron core and return limbs; • the domain C i ⇢ H corresponds to the M poloidal field coils, where each coil C i has n i wire turns, total resistance R i and cross section |C i |, the coils are part of an electric circuit, that contains also N voltage suppliers; • the domain S ⇢ H corresponds to the passive structures, with conductivity ; • the domain L ⇢ H, bounded by the limiter, corresponds to the domain that is accessible by the plasma.

Then (see [START_REF] Blum | Numerical simulation and optimal control in plasma physics[END_REF] or [START_REF] Jardin | Computational methods in plasma physics[END_REF]), the equilibrium of plasma in a tokamak has to satisfy at each instant the following non-linear initial value problem:

(5)

L (r, z, t) = j(r, (r, z, t), t) in ⌦; (0, z, t) = 0 ; lim k(r,z)k!+1 (r, z, t) = 0 ; (r, z, 0) = 0 (r, z) ,
where L is a non-linear second-order elliptic di↵erential operator (6)

L := @ @r ✓ 1 µ[ ]r @ @r ◆ @ @z ✓ 1 µ[ ]r @ @z ◆ := r • ✓ 1 µ[ ]r r ◆ , r F z C i S P L @L Figure 1.
Left: Geometric description of the poloidal cross section of the tokamak device. Right: Sketch for characteristic plasma shapes during the so-called ramp-up phase. The -isolines are indicated by black lines. In the beginning (first three pictures) the plasma touches the limiter and becomes more and more elongated (limiter configuration) while finally, it moves into the divertor configuration, where the plasma boundary contains an X-point of the poloidal flux. with

(7) µ[ ](r, z) = µ f (|r (r, z)| 2 r 2 ) ( µ 0 in F = µ 0 elsewhere.
Here, r is the 2D gradient in the (r, z)-plane. The current density j is a non-linear function of :

(8) j(r, (t), t) = 8 > > > < > > > : rp 0 ( (t), t) + 1 µ0r ff 0 ( (t), t) in P[ (t)] ; Ii |Ci| in C i ; r @ (t) @t
in S ; 0 e l s e w h e r e , where the plasma section P[ (t)] is the domain bounded by the last closed poloidal magnetic flux line inside L and containing the magnetic axis (r ax , z ax ). The domain L itself is bounded by the limiter and defines the domain that is accessible to the plasma. The magnetic axis is the point (r ax , z ax ) = (r ax [ ], z ax [ ]), where has its global maximum in L. For convenience, we introduce also the coordinates (r bd , z bd ) = (r bd [ ], z bd [ ]) of the point that determines the plasma boundary. (r bd , z bd ) is either a hyperbolic point of (in the case of a divertor configuration) or the contact point with the limiter @L (in the case of a limiter configuration). The di↵erent characteristic shapes of P( ) are illustrated in Figure 1 (right): the boundary of P( ) either touches the boundary of L, the limiter, or the boundary contains one or more hyperbolic points (X-points of ).

The current density j is non-linear in due to the non-linear functions p 0 , ff 0 as well as the definition of the plasma domain P( ).

The total coil current I i in [START_REF] Grad | Classical di↵usion in a tokamak[END_REF] relates to voltages in the suppliers and to self and mutual induction (9)

I i |C i | = n X j=1 R ij V j (t) + m X k=1 S ik Z C k @ (t) @t dr dz , 1  i  m ,
via electric circuit equations. In the case of simple electrical circuits, e.g. circuits where each coil is connected to only one voltage supplier, we have

R ij = ( 0 i 6 = j ni Ri|Ci| i = j and S ij = ( 0 i 6 = j 2⇡n 2 i Ri|Ci| 2 i = j
and we refer to the Appendix B for the general case.

Numerical Methods

There are two di↵erent approaches to arrive at finite dimensional SQP formulations for the constrained optimization problem [START_REF] Albanese | On the solution of the magnetic flux equation in an infinite domain[END_REF]. In the first one, the optimizethen-discretize approach, one computes first the optimality conditions for the continuous optimization problem, and then discretizes them. The second approach, the discretize-then-optimize approach, discretizes directly the optimization problems. This yields a finite dimensional constrained optimization problem for which a SQP formulation follows immediately from the literature [START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Bonnans | Numerical Optimization: theoretical and practical aspects : with 26 figures[END_REF][START_REF] Hinze | Optimization with PDE constraints[END_REF] (see also the Appendix A).

We prefer to work with the discretize-then-optimize approach for the following reasons: the discretize-then-optimize approach yields the exact gradient of the discrete objective function, while the optimize-then-discretize approach yields only an approximation. Both approaches involve approximation, but the optimize-thendiscretize approach does not yield the exact gradient of either the continuous objective in [START_REF] Albanese | On the solution of the magnetic flux equation in an infinite domain[END_REF], or the discretized objective [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF]. Therefore, the optimize-then-discretize approach may result in inconsistent gradients, which, in turn, may cause serious convergence di culties in the optimization process.

Hence, it remains to specify the discretization of non-linear constraints ( 5), ( 8) and ( 9) and the discretization of the cost functionals. Moreover, as we use SQP, we will also have to provide some details on first and second order derivatives of the discretized cost functionals and constraints. The implementation itself is kept flexible, so that cost functionals can be easily changed and modified. Adding new cost functionals encoding new design goals is simple.

For the discretization in time we introduce N T + 1 not necessarily equidistant collocation points T 0  t s  T 1 , 0  s  N T . The discretization in space is based on the finite element method. Further details on the discretization of the constraints and the objective functional of (1) follow in Sections 3.1 and 3.2. The Section 3.3 bundles these two parts and presents the SQP formulation.

3.1. Free-Boundary Plasma Equilibrium Evolution. The current density j in ( 8) is non-linear in due to the non-linear functions p 0 and ff 0 and the definition of the plasma domain P( ). While P( ) is fully determined for a given , the two functions p 0 and ff 0 are not determined by the model ( 5)-( 8). The model ( 5)-( 8) needs to be augmented by the so-called transport and di↵usion equations, that determine p 0 and ff 0 . In this work, we will assume that, up to some scaling coe cient , the functions p 0 and ff 0 are known. But, the domain of p 0 and ff 0 is the interval [ bd , ax ] with the scalar values ax and bd being the flux values at the magnetic axis and at the boundary of the plasma: [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF] ax

[ ] := (r ax [ ], z ax [ ]) , bd [ ] := (r bd [ ], z bd [ ]) .
So, since the domain of p 0 and ff 0 depends on the poloidal magnetic flux itself, it is more practical to supply those profiles as functions of the normalized poloidal flux (r, z):

(11) (r, z) := (r, z) ax [ ] bd [ ] ax [ ] .
These two functions, subsequently termed S p 0 and S ff 0 , have, independently of , a fixed domain [0, 1]. The functions S p 0 and S ff 0 are usually given as polynomial functions or piecewise polynomial functions. The scaling coe cients (t) are determined by fixing the total plasma currents R P( ) j(r, (r, z, t), t) drdz to some given values I p (t):

(t) Z P( ) rS p 0 ( (r, z, t), t) + 1 µ 0 r S ff 0 ( (r, z, t), t) drdz = I p (t) .
Let ⌦ ⇢ H be a su ciently large semi-circle of radius ⇢ , that is centered at the origin and contains the geometry of the tokamak. The boundary @⌦ splits into r=0 := {(r, z) , r = 0} and = @⌦ \ r=0 . We introduce a triangulation ⌦ h of the domain ⌦ that resolves the subdomains L, F, C i , S and use standard H 1 -conforming Lagrangian finite elements with nodal degrees of freedom.

Let b k (r, z) denote the Lagrangian basis functions associated to the vertices of the mesh, e.g. b k vanishes at all mesh vertices except one. Basis functions associated to vertices at r = 0 are excluded from this finite element space X h , as, due to axisymmetry (0, z) = 0. The finite element approximation h of at t = t s is:

(12) s h (r, z) = |X h | X k=1 ks b k (r, z) with ks 2 R, 1  k  |X h |, 0  s  N T .
The domain of the plasma P( h ) of a finite element function h is bounded by a continuous, piecewise straight, closed line. The critical points (r ax ( h ), z ax ( h )) and (r bd ( h ), z bd ( h )) are the coordinates of certain vertices of the mesh. The saddle point of a piecewise linear function h is some vertex (r 0 , z 0 ) with the following property: if (r 1 , z 1 ), (r 2 , z 2 ) . . . (r N , z N ), denote the counterclockwise ordered neighboring vertices the sequence of discrete gradients 0 1 , 0 2 . . . 0 N changes at least four times the sign.

We get the following discretized version (modulo quadrature) of the free-boundary plasma equilibrium evolution: For given evolution Ṽ (t) = (V 1 (t), . . . V m (t)) and I p (t) of the voltages and the total plasma current and for given initial conditions

0 h (r, z) = P |X h | k=1 k0 b k (r, z) we find for 1  s  N T the coe cients ks 2 R, 1  k  |X h | of s h , the values I is 2 R, 1  i  m and s 2 R such that: (13) Z ⌦ 1 rµ[ s h ] r s h • r⇠ drdz s Z P( s h ) rS p 0 ( s h , t s )⇠ + S ff 0 ( s h , t s ) rµ 0 ⇠ drdz m X i=1 I is |C i | Z Ci ⇠drdz + Z S r s h s 1 h t s t s 1 ⇠ drdz + c( s h , ⇠) = 0 8⇠ 2 X h , n X j=1 R ij V j (t s ) + m X k=1 S ik Z C k s h s 1 h t s t s 1 drdz = I is |C i | 1  i  m , s Z P( s h ) rS p 0 ( s h , t s ) + 1 rµ 0 S ff 0 ( s h , t s ) drdz = I s P .
The bilinear form c : V ⇥ V on , accounts for the boundary conditions at infinity [START_REF] Albanese | On the solution of the magnetic flux equation in an infinite domain[END_REF]:

(14) c( , ⇠) := 1 µ 0 Z (P 1 )N (P 1 )⇠(P 1 )dS 1 + 1 2µ 0 Z Z ( (P 1 ) (P 2 ))M (P 1 , P 2 )(⇠(P 1 ) ⇠(P 2 ))dS 1 dS 2 .
with

M (P 1 , P 2 ) = k P1,P2 2⇡(r 1 r 2 ) 3 2 2 k 2 P1,P2 2 2k 2 P1,P2 E(k P1,P2 ) K(k P1,P2 ) ! N (P 1 ) = 1 r 1 ✓ 1 + + 1 1 ⇢ ◆ and ± = q r 2 1 + (z 1 ⇢ ) 2 ,
where P i = (r i , z i ) and K and E are the complete elliptic integrals of first and second kind, respectively and

k Pj ,P k = s 4r j r k (r j + r k ) 2 + (z j z k ) 2 .
We refer to [START_REF] Grandgirard | Modélisation de l'équilibre d'un plasma de tokamak[END_REF]Chapter 2.4] for the details of the derivation. Alternative approaches, that incorporate the boundary conditions at infinity were recently presented in [START_REF] Faugeras | Fem-bem coupling methods for tokamak plasma axisymmetric free-boundary equilibrium computations in unbounded domains[END_REF]. Quadrature rules. The integration over the domain ⌦ is split into a sum of integrals over the triangles T of the mesh and we use the barycenter quadrature rule to approximate these integrals: For a triangle with vertex coordinates a i , a j , a k 2 R 2 the quadrature point is the barycenter b T := 1 3 (a i + a j + a k ) and the quadrature weight ! T is the size |T | of the triangle.

Besides the integrals over ⌦, the weak formulation (13) involves also integrals over the plasma domain P( h ). As the mesh does not resolve the boundary of the plasma domain P( h ), we need to specify also the quadrature rule that is used to approximate integrals over intersections T \ P( h ) of triangles with the plasma domain. We use again barycenter quadrature, but here the quadrature point and weight will depend non-linearly on h , which needs to be taken into account when we use linearizations of ( 13), e.g. in Newton or SQP iterations. The technical details of such linearizations can be found in [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF].

The line integrals over in the definition of c(•, •) are split into line integrals over edges on triangles whose vertices are on . Then, the trapezoidal quadrature rule yields an approximation of su cient accuracy.

3.2.

Cost-and Regularization Functionals. As we already specified the approximation of there is not much left to explain on the discretization of cost functionals that involve the poloidal flux and the voltages.

Here, the same barycentric quadrature rule as in 3.1 is used to approximate area integrals, while time integrals are approximated by the composite trapezoidal rule. This yields expressions that are algebraic in the principal unknowns ks ,

1  k  |X h |, 1  s  N T and I is , 1  i  m, 1  s  N T .
So, computation and implementation of first and second order derivatives of such functionals become a mechanical iterated application of the usual rules of di↵erential calculus.

The discretization of the regularization functionals likewise uses composite trapezoidal rule for time integrals. Moreover, we represent each voltage as a finite series

V i (t) = Nc X j=1 V ij B j (t), 1  i  n ,
where the B j (t) are the basis functions of a polynomial or spline space over [0, T ]. So again this yields expressions that are algebraic in the unknowns

V ij , 1  i  n, 1  j  N c . 3.3. SQP Formulation.
Combining the discretized free-boundary plasma equilibrium evolution [START_REF] Li | The component development status of HL-2M tokamak[END_REF] with discretized objective and regularization functionals we arrive at discretized versions of the optimal control formulation (1) that is of the general form [START_REF] Nocedal | Numerical optimization[END_REF] min [START_REF] Li | The component development status of HL-2M tokamak[END_REF]. The control variable u contains a subset of the expansion coe cients V ij , 1  i  n, 1  j  N c for the voltages. One could for example prescribe the voltage of a few suppliers and treat only the voltages of the remaining suppliers as unknowns. Another possibility would be to treat only certain coe cients as unknown. E.g. if one works with an expansion in a hierarchical basis, it would be beneficial to prescribe the coe cients for low order terms and keep only coe cients of higher order polynomials. So, in general

u 2 R N with 1  N  N c n.
The SQP formulation for the constrained optimization problem (15) involves first and second order derivatives of b(y, u) and J(y, u). But, as we presented in Sections 3.1 and 3.2 explicit expressions for b(y, u) and J(y, u) that are algebraic in u and y, we can also provide explicit expressions for the first and second order derivatives. An inspection of [START_REF] Li | The component development status of HL-2M tokamak[END_REF] shows that the constraint b(y, u) is a ne in the control unknown u, hence has vanishing second order derivatives. Moreover, to avoid the expensive assembling of second order derivatives of b(y, u) we neglect those in the SQP iterations. In the terminology of Newton methods we use rather a quasi SQP method, than an exact SQP method. It is known that such modifications are prone to convergence issues [START_REF] Bonnans | Numerical Optimization: theoretical and practical aspects : with 26 figures[END_REF], but this doesn't seem to be an issue for our specific application. The Algorithm 1 summarizes the quasi SQP formulation, and should be compared to the exact SQP in Algorithm 4 (see Appendix A). We want to stress that the size of Y in Algorithm 4 scales at least linearly (and quadratically if N c = O(N T )) with the number of time steps N T , which causes memory to be the limiting factor for computations with a very large number of time steps. Nevertheless, with the current implementation we can go easily beyond 200 and more time steps, which is su cient for the applications in mind. Moreover, we avoid the introduction of additional parameters due to iterative solver and make use of fast methods for linear systems with multiple righthand sides.

Algorithm 1 SQP (quasi) with direct solver

1: u 1, y 1, y y 0 , u u 0 2: while k uk/kuk > tol, k yk/kyk > tol do 3: ( y, Y) b 1 y (y, u)(b(y, u), b u (y, u))
4:

m 01 J T u (y, u), m 10 Y T J T y (y, u) 5:
M 02 J uu (y, u) u u + u 13: end while

Numerical Validation Tests

We present di↵erent numerical tests, that verify that the implementation of derivatives of the non-linear constraints [START_REF] Li | The component development status of HL-2M tokamak[END_REF] in FEEQS.M is correct. These derivatives are the cornerstone of the proposed optimal control approach and appear through b y (y, u) in SQP (see Algorithm 1).

The following calculations are based on a ITER-like geometry and an equlibrium plasma (see Fig. 2) that corresponds to the currents in the table of Figure 2. The total plasma current is I P = 15 ⇥ 10 6 A and the two profiles S p 0 and S ff 0 in the definition of the plasma current density are ( 16)

S p 0 (y) = r 0 ((1 y) ↵ ) , S ff 0 (y) = (1 )µ 0 r 0 (1 y ↵ ) ,
with the four parameters r 0 = 6.2m , ↵ = 2.0, = 0.5978 and = 1.395.

(17) J( h , b m ) = Z P( h )  rS p 0 ( h ) + 1 rµ 0 S ff 0 ( h , t s ) b m drdz 1  m  |X h | .
This mapping involves the non-linearity due to the free plasma boundary. The linearization is non-standard [START_REF] Heumann | Quasi-static freeboundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF]. In this first test it is verified that for given

h (r, z) = P |X h | k=1 k b k (r, z) and perturbation h (r, z) = P |X h | k=1 k b k (r, z) the imple- mentation of DJ h ( h ) yields: (18) E FD (") := " 1 J h ( h + " h , •) J h ( h , •) DJ h ( h ) • kDJ h ( h ) • k = O(").
The perturbation increment is randomly chosen. In Table 1 we monitor this relative error and observe, as expected, first order convergence.

4.2.

Convergence of the Sensitivities. Another way to validate the implementation of the derivatives appearing in b y (y, u) would consist in solving the nonlinear constraint b(y, u) for y and monitor the convergence. For an accurate implementation of the derivatives the convergence will be quadratic. On the other hand, since in many cases also with inaccurate derivatives one will observe super-linear convergence, this test could fail as indicator for accurate derivatives. The following test is more appropriate: We consider a control u " that is parametrized by " 2 R and suppose we know y 0 that verifies b(y 0 , u 0 ). Then it holds

(19) y 1 " y " = O(" 2 ) ,
where y 1 " is the first Newton iteration for the problem b(y " , u " ) = 0 with initial guess y 0 :

(20) b y (y 0 , u " ) • (y 1 " y 0 ) = b(y 0 , u " )
The result (19) follows from y " = y 0 +y 0 0 "+O(" 2 ) and 0 = b(y 0 , u " )+"b y (y , u " )y 0 0 + O(" 2 ). When the exact derivative b y (y 0 , u " ) in (20) is replaced by some approximation A: A

• (y 1 " y 0 ) = f (u " ) a(y 0 ) .
then y 1 " will be at most first order approximation of y " . The first test of this kind is for the static free-boundary plasma equilibrium (FBE) problem with the following Galerkin formulation: For given currents I 1 , . . . I M : Find h such that (21)

Z ⌦ 1 rµ[ h ] r h • r⇠ drdz Z P( h ) rS p 0 ( h )⇠ + S ff 0 ( h )⇠ rµ 0 drdz = M X i=1 I i |C i | Z Ci ⇠drdz 8⇠ 2 X h Z P( h ) rS p 0 ( h ) + S ff 0 ( h ) rµ 0 drdz = I P .
The discretization and linearization follow as described in Section 3.1 and 3.3.

Then the unperturbed control data u 0 is the data given in the table in Figure 2 and perturbation is a random incremental current for each coil scaled with " = 0.5 0 , . . . 0.5 14 . In Table 1 (4th column) the error E DN (") = ky " y 1 " k is monitored. We observe second order convergence, which shows that we use accurate derivatives. In contrast, the result E CN (") = ky " y 1 " k for a Newton-type iteration method (see 6th column in Table 1) that follows from the discretization of derivatives of the continuous problem [START_REF] Blum | Numerical simulation and optimal control in plasma physics[END_REF][START_REF] Grandgirard | Modélisation de l'équilibre d'un plasma de tokamak[END_REF] yields only first order convergence.

Repeating the same test for the implementation for FBEE problem [START_REF] Li | The component development status of HL-2M tokamak[END_REF], based on the discretization and linarization described in Section 3.1 and 3.3 we observe the expected second order convergence (see Table 1, 8th and 9th column). This reassures that the implementation in FEEQS.M is correct.

Applications for tokamaks WEST and HL-2M

5.1. Fast limiter to divertor transition for the WEST tokamak. In this section we give tests with di↵erent objective functions and regularization terms for the WEST tokamak [START_REF] Bourdelle | West physics basis[END_REF], which is a new tokamak in CEA Cadarache that aims at testing the ITER-like technology for tungsten (W) plasma facing components. Exploring scenarios which can handle fast transitions from limiter to divertor configuration within the engineering limits, such as the maximum and minimum voltages (in Table 2) that can be handled by the coils power supply systems may be a way of reducing plasma pollution by W and will be our goal in the following simulations.

The Figure 3 shows the geometrical setting of WEST. All the results are based on a mesh with 52210 triangles and 27763 vertices. We use N T = 10 time steps with equidistant time step length of 5 ms, the voltage in A coil is fixed to its maximum value of 1400 V in all the following cases to provide energy for flux consumption, and the degree of the polynomial representation of the voltage is 8. The evolution of the total plasma current I p (t) and the parameters ↵(t), (t) and (t) for the current density profiles S p 0 and S ff 0 in ( 16) are based on experimental data of Tore Supra, the predecessor of WEST. 5.1.1. Prescribed level sets at all time steps. As a first approach to the problem, we solve the discretized inverse evolution problem [START_REF] Nocedal | Numerical optimization[END_REF] with J(y, u) := C(y) + R 1 (u), 1. Convergence of the errors E ... and the convergence rate ( log(E...("i+1)) log(E...("i)) log("i+1) log("i)

" i = 1 2 i J in (
): 1.) the finite di↵erence error E FD (18); 2.) the error E DN using the derivatives outlined in Sections 3.1 and 3.3 for the FBE problem (21); 3.) the error E CN using derivatives in [START_REF] Blum | Numerical simulation and optimal control in plasma physics[END_REF][START_REF] Grandgirard | Modélisation de l'équilibre d'un plasma de tokamak[END_REF] for the FBE problem (21); 4.) the error E DN using the derivatives outlined in Sections 3.1 and 3.3 for the FBEE problem [START_REF] Li | The component development status of HL-2M tokamak[END_REF]. Here, C(y) is the discretization of the cost function (2) introduced in Section 1, and has the purpose to prescribe a desired level set for every time step, and the weights are

Coil Maximum Voltage (V) Minimum Voltage (V)

C(y) := 1 2 N T X s=1 N desi X i=1 w s s h (r i , z i ) s h (r 0 , z 0 ) 2 , R 1 (u) := 1 2 w R1 L X i=1 N T X s=1 V i (t s ) 2 .
w i = 1, 1  i  N T .
The weight of the regularization term R 1 (u) is fixed as w R1 = 1 ⇥ 10 11 during this calculation. The results (see Figure 4) show that the prescribed level sets coincide fairly well with level sets of the numerical solution at all times, with a tiny discrepancy when the plasma evolves from limiter to divertor configuration. The evolution of the voltages can be seen in Figure 5, which shows that the voltage limits (the dashed lines) at some time steps are violated.

5.1.2.

Prescribed level sets at all times and constraint penalization term. Secondly, we add two regularization terms to the objective function to penalize the violation of the voltage limits. E.g. we solve now the constrained optimization problem (15) for J(y, u)

:= C(y) + R 1 (u) + R 2 (u) + R 3 (u) with (23) R 2 (u) := 1 2 L X i=1 w R2,i N T X s=1 max(V i (t s ) V i,max , 0) 3 , R 3 (u) := 1 2 L X i=1 w R3,i N T X s=1 max(V i,min V i (t s ), 0) 3 ,
where the weights are w R2,i = w R3,i = 1 ⇥ 10 16 . The level sets of the flux inside and near the plasma are found to be almost the same as in the previous case (shown in Figure 4). But the evolution of the voltages, shown in Figure 6, is slightly di↵erent and almost all voltages remain within their limits.

5.1.3.

Prescribed level sets at start and end. Assuming that level sets are prescribed at every time step is not very relevant in practical applications, where one aims to solve problems that require a few hundred or thousands of time steps. Providing level sets a each time step is a humongous amount of work that moreover confines Voltages (V) the plasma evolution fairly strictly and could prevent from finding the best scenarios. The more relevant task asks to find reasonable trajectories of the flux for prescribed level sets at the beginning and end only. Therefore, we use in this third case the prescribed boundaries only on some of the time steps, meaning we set in the objective function C(y) in ( 22) the weights w 1 , w 2 , w 3 and w 10 to one but the remaining weights w 4 , w 5 , ...w 9 to zero. Introducing

V Dh V Eh V Fh V Fb V Eb V Db
C 0 (y) := 1 2 X s={1,2,3,10} N desi X i=1 s h (r i , z i ) s h (r 0 , z 0 ) 2
for brevity, we solve hence here the optimization problem [START_REF] Nocedal | Numerical optimization[END_REF] with J(y, u)

:= C 0 (y) + R 1 (u) + R 2 (u) + R 3 (u)
, where the weights in the regularization terms R 1 (u), R 2 (u) and R 3 (u) are as in the previous case.

The plasma trajectory of the numerical solution is shown in Figure 7. It is found that the boundaries of the plasma match with the prescribed level sets in the beginning and end. But in between, when there is no prescribed level set, the plasma position and shape is changing drastically. The evolution of the voltages are presented in Figure 8, which shows that the voltages in the time steps without prescribed boundaries are within the bounds and the voltages in the divertor coils are violating the limits only during the last time step.

We obtained a solution that is reasonable for this particular formulation of the optimal control problem for scenario design: As the formulation does not contain any objective for the flux in between the start and the end, we risk in computing solutions with fast changes of plasma position and shape. On the other hand, scenarios with such fast changes are highly undesirable as they are very di cult to control in the experiment. The formulation in the next section is more adequate.

5.1.4.

Prescribed level sets at start and end with penalization term on induced currents. Fast changes of the shape and position of the plasma lead to large induced currents in passive structures S. In order to avoid the fast changes obtained in the previous simulations we add to the objective a new term that penalizes the induced currents Voltages (V) in passive structures S with conductivity . We solve the optimization problem [START_REF] Nocedal | Numerical optimization[END_REF] with J(y, u)

I ps,s (y) := Z S r s+1 h (r, z) s h (r, z) t s+1 t s dr dz 1  s  N T
V Bh V Dh V Eh V Fh V Fb V Eb V Db V Bb
:= C 0 (y) + C 1 (y) + R 1 (u) + R 2 (u) + R 3 (u) and (24) C 1 (y) = 1 2 w C1 N T X s=1 I 2 ps,s (y)
where w C1 = 1 ⇥ 10 6 . The results, shown in Fig. 9 are quite satisfying: the plasma stays near to the equatorial plane even during the time steps when the level set is not prescribed and evolves in a smooth fashion to its final shape. This shows that the penalization term on induced currents allows finding slowly moving solutions even if we prescribed level sets only in the beginning and end. The evolution of the voltages is shown in Figure 10 and we highlight that all the voltages remain within the limits.

These tests show that with di↵erent objective functions and regularization terms, we can get a reliable plasma limiter to divertor transition while respecting the engineering limits (only maximum and minimum voltages where considered in this case). Voltages (V)

V Dh V Eh V Fh V Fb V Eb V Db
Figure 10. Prescribed level sets at start and end with penalization term on induced currents (see 5.1.4): The voltage evolution, dashed lines are the limits.

5.2.

Full scenario computation for HL-2M tokamak. We give in this section a calculation for a full scenario, which includes the plasma current ramp up, flat top and ramp down phases, for the new tokamak of HL-2M [START_REF] Li | The component development status of HL-2M tokamak[END_REF] (see Figure 11 for the cross-section) in China. A major objective of this machine is the exploration of ITER related physical issues. The settings for the computations are the following: The parameters I p , ↵, and are prescribed as in Figure 12. The length of time steps in the ramp up phase (from 0.1s to 2s) and ramp down phase(from 8.1s to 10s) is 0.1s, while in the flat top phase (from 2.3s to 8s) it is 0.3s. In total we have 60 time steps. The objective function is J(y, u) = C(y) + R(u), with C(y) the discretization of (2) and R(u) the discretization of (3). Plasma boundaries are prescribed for every time step and the weights w and D ii are 1 and 1 ⇥ 10 10 and zero otherwise. The computational mesh is divided into 30583 triangles, and the total number of vertices is 15367, and the degree of polynomial to fit the voltage is again fixed as 8. Altogether we end up with a constrained optimization problem that has almost one million unknowns of which approximately a hundred correspond to control unknowns u. We solve this problem in 15 iterations and the computing time is less than 350 s on standard notebook with a 4 ⇥ 2.7 GHz processor and 32 GB memory.

In Figure 13 we show plasma boundaries at some selected time steps. As can be seen, the plasma follows well with reference boundary in the whole time interval. The voltages are shown in Figure 14. This example is a showcase that the presented approach can construct a full operation scenario for tokamaks with reasonable amount of computational power.

Conclusions and Perspectives

The models that describe the evolution of plasma in tokamak devices are highly non-linear. Alone the direct simulation, meaning a simulation that shows only how the plasma will evolve for certain given control inputs is a demanding task. Moreover, the corresponding inverse problem the much more relevant problem when experimenting with tokamaks: what are the control parameters that ensure that the plasma evolves through certain prescribed states? From the engineering point of view it would be beneficial to provide the majority of control as feed-forward control that can be determined a priori, and to minimize the amount of feedback control that needs to be calculated during the operation. Until today, the calculation of the feed-forward control is mainly tackled by ODE models that have the flavor of electrical circuits [START_REF] Ariola | Magnetic Control of Tokamak Plasmas[END_REF]. These ODE models determine the coil voltages that correspond to coil currents of precomputed static equilibria. Particular feedback controllers are then necessary to ensure that the currents in the system are indeed close to the precomputed ones.

In this work we formulated an optimal control problem, that uses the force balance (4) and Maxwell's equations as constraints, to compute directly the optimal feed-forward control. The approach is more consistent, as it uses the full nonlinear model of partial di↵erential equations. Then, in the solution of the optimization problem, it is Newton's method that determines consistently the linearization of this non-linear model.

We presented here the two basic ingredients, the finite element method and sequential quadratic programming, to set up such an optimal control approach for plasma evolution in tokamaks. The finite element method allows for an easy treatment of free plasma boundary and the decaying conditions at infinity. Further, we can, almost automatically, compute linearizations and the adjoint operators that are required for the optimization with sequential quadratic programming. We highlighted that existing implementations of Newton's methods for the direct simulation can be easily extended to an implementation of SQP-type methods for the corresponding inverse or optimal control problems. In the future we will have to augment the non-linear model ( 5)-( 8) with further equations, namely the resistive di↵usion and transport equations, that provide a more complete description.

J uu (u) =Z T (u) ✓ J yy (y(u), u) J yu (y(u), u) J uy (y(u), u) J uu (y(u), u) ◆ Z(u) + Z T (u) ✓ b T yy (y(u), u)p(u) b T yu (y(u), u)p(u) b T uy (y(u), u)p(u) b T uu (y(u), u)p(u) ◆ Z(u) with Z(u) = ✓ b 1 y (y(u), u)b u (y(u), u) 1 ◆ ,
but the requirement of solving repetitively the non-linear constraints b(y(u), u) = 0 remains a big drawback.

Sequential Quadratic Programming (SQP) is an entirely di↵erent approach that avoids this drawback and incorporates at the same time second order information. SQP is one of the most e↵ective methods for non-linear constrained optimization with significant non-linearities in the constraints [START_REF] Nocedal | Numerical optimization[END_REF]Chapter 18]. To motivate SQP we recall that the previous discussion shows that the control u ⇤ is a stationary point of the reduced objective function if and only if there exist states y ⇤ and adjoint states p ⇤ such that (30)

J T y (y ⇤ , u ⇤ ) + b T y (y ⇤ , u ⇤ )p ⇤ = 0 , J T u (y ⇤ , u ⇤ ) + b T u (y ⇤ , u ⇤ )p ⇤ = 0 , b(y ⇤ , u ⇤ ) = 0 ,
A Newton-type method for solving (30) are iterations of the type

(31) 0 @ H k yy H k yu b T y (y k , u k ) H k uy H k uu b T u (y k , u k ) b y (y k , u k ) b u (y k , u k ) 0 1 A 0 @ y k+1 y k u k+1 u k p k+1 1 A = 0 @ J T y (y k , u k ) J T u (y k , u k ) b(y k , u k ) 1 A with ✓ H k yy H k yu H k uy H k uu ◆ := ✓ H yy (y k , u k , p k ) H yu (y k , u k , p k ) H uy (y k , u k , p k ) H uu (y k , u k , p k ) ◆ := ✓ J yy (y k , u k ) J yu (y k , u k ) J uy (y k , u k ) J uu (y k , u k ) ◆ + ✓ b T yy (y k , u k )p k b T yu (y k , u k )p k b T uy (y k , u k )p k b T uu (y k , u k )p k ◆ ,
which corresponds to the following quadratic optimization problem with linear constraints:

min u k+1 ,y k+1 1 2 ✓ y k+1 y k u k+1 u k ◆ T ✓ H k yy H k yu H k uy H k uu ◆ + ✓ J T y (y k , u k ) J T u (y k , u k ) ◆ ! ✓ y k+1 y k u k+1 u k ◆ s.t b(y k , u k ) + b y (y k , u k )(y k+1 y k ) + b u (y k , u k )(u k+1 u k ) = 0 .
This sequence of quadratic optimization problems is at the origin of the name sequential quadratic programming.

If the linear systems in (31) become too large, it is the common practice to pursue the null space approach to arrive at the SQP formulation with the reduced Hessian. In introducing

(32) Z k = ✓ b 1 y (y k , u k )b u (y k , u k ) 1 ◆ and Y k = ✓ b 1 y (y k , u k ) 0 ◆ we obtain the identity (33) ✓ y k+1 y k u k+1 u k ◆ = Z k (u k+1 u k ) + Y k b(y k , u k )
and find the following linear system for the increment

u k+1 u k (34) M(y k , u k )(u k+1 u k ) = m(y k , u k ) with M(y k , u k ) :=Z T k ✓ H k yy H k yu H k uy H k uu ◆ Z k and m(y k , u k ) :=Z T k ✓✓ J T y (y k , u k ) J T u (y k , u k ) ◆ + ✓ H k yy H k yu H k uy H k uu ◆ Y k b(y k , u k ) ◆ .
It is insightful to compare the expressions involved in the reduced formulation (34) of SQP to the gradient (28) and the Hessian (29) of the reduced objective function: the gradient and Hessian of the reduced objective are equal to m(y, u) and M(y, u) only when the state y and control u verify the non-linear direct problem.

In general, iterative methods, such as the conjugate gradient (CG) methods, are used to solve (34). Within each iteration step two linear systems corresponding to b 1 y (y k , u k ) and b T y (y k , u k ) (see Algorithm 3) need to be inverted. The CG algorithm (Algorithm 3) appears di↵erent than standard formulations, as we update within the CG-iterations not only the control but also the state unknown. This avoids solving one additional non-linear direct problem after each CG call.

Algorithm 3 SQP with CG iterative solver (less memory intensive) The external circuit of poloidal field coils and suppliers is represented by a directed graph with nodes and directed edges between nodes. We assign to each directed edge s = (i, j) between two nodes with index i and j a (directed) voltage V s and current I s . For an edge s representing a coil C we have (35)

I s = Z ss V s + ( FC ( d dt )) s := n s R s V s 2⇡ n 2 s R s 1 C Z C d dt drdz ,
where R s and n s are the total resistance and the wire turns of the coil.

Then we introduce the incidence matrix G that has entry G s,i = 0 when the node with index i is not contained in the edge with index s = (i, j), or entry G s,i = 1( 1) when the node with index i is contained in the edge with index s = (i, j) and induced orientation coincides (coincides not) with the orientation of the edge. Likewise we can introduce oriented polygons, whose boundaries are the edges of the graph and an incidence matrix C for edges and polygons. Then we have CG = 0 and the Kirchho↵ current and voltage laws are:

G T Ĩ = 0 and C Ṽ = 0 ,
where the components of Ṽ and Ĩ are the voltages and currents associated to the edges.

In the following we assume that the edges of the circuit correspond either to a coil or to an external voltage supplier/source. We introduce the subscripts S and C to distinguish between edges corresponding to supplies and edges corresponding to coils, e.g. Well-posedness requires to fix Ũ at one node and to remove the Kircho↵ current equation for the same node. We choose an arbitrary node and fix Ũ to zero at this point. This is equivalent to cancel a column (row) in G (G T ), hence in the subsequent lines G • (G T

• ) always refers to the reduced matrices. After some tedious calculations we find

  , u) s.t. b(y, u) = 0 . The state variable y contains the unknowns ks 1  k  |X h |, 1  s  N T of the poloidal flux and the scaling parameter s , 1  s  N T , hence y 2 R (|X h |+1)N T . Likewise the constraint b(y, u) = 0 corresponds to (|X h | + 1)N T coupled non-linear equations as the unknowns for the currents I is , 1  i  m, 1  s  N T have been eliminated from the formulation

4. 1 .

 1 Derivatives due to the free-boundary. Let DJ h ( h ) 2 R |X h |⇥|X h | denote the derivative of the discretization J h ( h , b m ) (via quadrature from Section 3.1) of the non-linear mapping

Figure 2 .

 2 Figure 2. Left: The plasma (flux lines and flux intensity) that correspond to the currents in the table. Center: The triangulation of the ITER-geometry. Right: Data for coils. Coils U1-U3 and L1-L3 are upper and lower central solinoid coils. Coils P1-P6 are the poloidal field coils.

Figure 3 .

 3 Figure 3. Cross-section of the WEST tokamak, red rectangles are the poloidal coils C i and iron areas F are in purple, the dark grey areas are the vacuum vessels, the divertor (X) coils casings and the passive plates S .

Figure 4 .Figure 5 .

 45 Figure 4. Prescribed level sets at all time steps (see 5.1.1): the magenta contour is the plasma boundary and the blue lines are the level sets inside the plasma, the red points are the desired boundaries.

Figure 6 .

 6 Figure 6. Prescribed level sets at all time steps and constraint penalization term (see 5.1.2): The voltage evolution, dashed lines are the limits.

Figure 7 .

 7 Figure 7. Prescribed level sets at start and end (see 5.1.3): the magenta contour is the plasma boundary and the blue lines are the level set inside plasma, the red points are the desired boundaries.

Figure 8 .

 8 Figure 8. Prescribed level sets at start and end (see 5.1.3): The voltage evolution, dashed lines are the limits.

Figure 9 .

 9 Figure 9. Prescribed level sets at start and end with penalization term on induced currents (see 5.1.4): the magenta contour is the plasma boundary and the blue lines are the level sets inside the plasma, the red points are the desired boundaries.

Figure 11 .

 11 Figure 11. Cross-section of the HL-2M tokamak, red rectangles are the poloidal coils C i , the grey areas are the two layers vacuum vessel, and the black line is the limiter L.

Figure 12 .

 12 Figure 12. The prescribed I p , and the parameters of the plasma profile for the full scenario simulation of HL-2M.

Figure 13 . 2

 132 Figure 13. Level lines (blue) of the magnetic flux inside the plasma at selected times for the calculation of a full scenario (ramp up, flat top, ramp down) for HL-2M based on 60 timesteps with prescribed boundaries (green dots) at each time step.

2 : 1 y 2 :

 212 while k uk > 0, k yk > 0 do u)(J T y (y, u) + H yy (y, u, p) y)5: r J T u (y, u) + H uy (y, u, p) y + b T u (y, u) p (y, u)b u (y, u)s 9: b b T y (y, u)(H yu (y, u, p)s + H yy (y, u, p)a) 10: e s H uu (y, u, p)s + H uy (y, u, p)a + b T u (y, u)b whileAlternatively, if we have su cient memory to store M(•, •), we can compute M(•, •) explicitly. We never compute neither b 1 y (y k , u k ) nor b T y (y k , u k ) explicitly. This alternative approach is summarized in Algorithm 4.Algorithm 4 SQP with direct solver (memory intensive) while k uk > 0, k yk > 0 do 3:( y, Y) b 1 y (y, u)(b(y, u), b u (y, u)) 4: m J T u (y, u) + Y T J T y (y, u) + H uy (y,u, p) y + Y T H yy (y, u, p) y 5: M H uu (y, u, p) + Y T H yu (y, u, p) + H uy (y, u, p)Y + Y T H yy (y, u, p)Y u)(J T y (y, u) + H yy (y, u, p) y + H yu (y, u, p) u) 11: end while Appendix B. Circuit Equations

  G T Ĩ = G T S ĨS + G T C ĨC = 0. Moreover we introduce a (node) potential Ũ with Ṽ = G Ũ and hence ṼS = G S Ũ and ṼC = G C Ũ . In summary we arrive at: be close with the notation from CEDRES++,

  R = I + ZG C KG T C + ZG C KG T S MG S KG T C Q = G C KG T C + G C KG T S MG S KG T C and S = ZG C KG T S M , T = G C KG T S M , where K 1 = G T C ZG C , and M 1 = G S KG T S .

Table 2 .

 2 Maximum and minimum voltages on the coils systems of WEST.

	A	1400	1400
	Bh	1400	1400
	Dh	2500	2500
	Eh	2500	2500
	F h	2500	2500
	F b	2500	2500
	Eb	2500	2500
	Db	2500	2500
	Bb	1400	1400
	X U	300	300
	X L	300	300
	in which:
	(22)

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Figure [START_REF] Nardon | Magnetic configuration and plasma start-up in the west tokamak[END_REF]. The voltages of the inverse evolution calculation for HL-2M.

Appendix A. Fast Algorithms for Optimal Control Problems

The following section is a short summary on algorithms for general optimal control problems, where we focus on finite dimensional optimal control problems. This simplifies considerably the presentation, and is also more relevant for this work as we always work with discretized versions of (1). While most of the subsequent methods are well-known and can be found in excellent text books such as [START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Bonnans | Numerical Optimization: theoretical and practical aspects : with 26 figures[END_REF][START_REF] Hinze | Optimization with PDE constraints[END_REF], we prefer to include this discussion to keep the presentation self-consistent.

We consider the following generic version of an optimal control problem (25) min

where y 2 R n and u 2 R m are the so-called state and control variables. The constraint b(y, u) = 0 with b(y, u) 2 R n is the discretization (in space and time) of ( 5), ( 8) and [START_REF] Grandgirard | Modélisation de l'équilibre d'un plasma de tokamak[END_REF]. The optimization problem (25) is a constrained optimization problem that can be recast as an unconstrained optimization problem. For this we introduce the

Some control u ⇤ is called the stationary point of the reduced objective function if and only if

An iterative approach to the solution of (26) faces two major challenges. First, in each iteration we have to solve the full non-linear problem b(y(u), u) = 0 for y(u) and second, we also need to compute the sensitivities y u (u). The second challenge


SPADE: A SMALL PARTICLE DETECTION METHOD USING A DICTIONARY OF SHAPES WITHIN THE MARKED POINT PROCESS FRAMEWORK

Nicolas Cedilnik † , Eric Debreuve ¶ , Fabienne De Graeve § , Florence Besse § , Xavier Descombes † Université Côte d'Azur, Team Morpheme † Inria ¶ Centre National de la Recherche Scientifique § Valrose Biology Institute

ABSTRACT

We introduce SPADE (small particle detection), a novel method to detect and characterize objects only a few pixel large on microscopy images, based on the use of a dictionary of shapes (Figure 1). It involves the minimization of an energy function embedding a prior and a data fidelity term to match shapes of the dictionary with patterns in images. The data fidelity term is a comparison between the pixels' intensities within a potential match and the surrounding pixels' intensities. It can be defined in different ways, from simple averages differences to formulas incorporating variance of the intensities. The prior forbids any overlap between the matched shapes. The sensitivity of the detector can be tuned using a minimal threshold on the data fidelity term.

Using synthetic images of increasing blur, i.e., detection difficulty, we showed that, on certain images, SPADE can outperform widely used methods in this classical task of biological image processing (Figure 2). We also show that the data fidelity term must be adapted to the nature of the images in order to get the best results.

We then apply SPADE, in a real-case scenario, on confocal fluorescence microscopy images of D. melanogaster cultured cells where cytoplasmic ribonucleoproteic granules have to be detected and characterized (Figure 3). We demonstrate that our method can be calibrated for high-throughput screening imaging study, generating a F 1 detection score of 0.76 after cross-validation using a trained biologist expertise as ground truth (data not shown).

An open-source python implementation of SPADE has been released on the Python Package Index1 . Different dictonaries of shapes and data fidelity terms could be used for different applications.

Fig. 2 .Fig. 3 .

 23 Fig. 1. SPADE default dictionary of shapes and their contours

https://pypi.python.org/pypi/ small-particle-detection/