Control methods for the optimization of plasma scenarios in a tokamak

This paper presents the modelling of the evolution of plasma equilibrium in the presence of external poloidal field circuits and passive structures. The optimization of plasma scenarios is formulated as an optimal control problem where the equations for the evolution of the plasma equilibrium are the constraints. The procedure determines the voltages applied to the external circuits that minimize a certain costfunction representing the distance to a desired plasma augmented by an energetic cost of the electrical system. A sequential quadratic programming method is used to solve the minimization of the cost-function and an application to the optimization of a discharge for ITER is shown.


Introduction
A tokamak is an experimental device whose purpose is to confine a plasma (ionized gas) in a magnetic field so as to control the nuclear fusion of atoms of low mass (deuterium, tritium,..) and to produce energy. The magnetic field has two components (see Fig. 1) : a toroidal field created by toroidal field coils, that is necessary for the stability of the plasma, a poloidal field in the section of the torus created by poloidal field coils and by the plasma itself.
The plasma current is obtained by induction from currents in these poloidal field coils. The tokamak thus appears as a transformer whose plasma is the secondary. The currents in the external coils play another role, that of creating and controlling the equilibrium of the plasma. The goal of this paper is to provide a model for the evolution in time of the equilibrium of the plasma and to derive control methods in order to optimize a typical scenario of a discharge of the plasma in a tokamak. There are two approaches for simulating a plasma made of electrons and ions: the microscopic approach based on kinetic equations (Vlasov, Boltzmann, Fokker-Planck) that are 6D (3D in space and 3D in terms of the velocity) and 1D in time.
the macroscopic approach based on magnetohydrodynamics (MHD) equations that are obtained by taking moments of the kinetic equations, and which are 3D in space and 1D in time. The validity of the MHD equations is clearly more restrictive than the one of the kinetic equations. We will present in section 2 the way in which the MHD equations are obtained from the kinetic ones.
At the slow resistive diffusion time-scale, the plasma is in equilibrium at each instant (the kinetic pressure force balances at each point the Lorentz force due to the magnetic field) and hence the plasma follows the so-called quasi-static evolution of the equilibrium. The resistive diffusion in the external passive structures surrounding the plasma and the equations of the circuits of the poloidal field system enable to follow in time this quasi-static evolution. An axisymmetric hypothesis enables to reduce the problem to a 2D p.d.e. formulation, with the Grad-Shafranov equation for the equilibrium of the plasma. The plasma boundary is a free boundary, which is a particular poloidal flux line. It is either the outermost closed flux line inside the limiter, which prevents the plasma from touching the vacuum vessel, or a separatrix (with a hyperbolic X-point), as they are in presence of a poloidal divertor. This equilibrium model will be presented in section 3 of the paper.
In order to solve numerically the set of equations for the poloidal flux, it is necessary to derive the weak formulation of this system and then a finite element method, coupled to Newton iterations for the treatment of the non-linearities, enables to solve the evolution of the equilibrium configuration in a tokamak. This is presented in section 4 of this paper.
A typical discharge in a tokamak is made of several phases : ramp-up of total plasma current, plateau phase (stationary phase), ramp-down. The plasma shape can also move from a small circular plasma (at the beginning of the discharge) to a large elongated one with an X-point. The goal of this work is to determine, thanks to optimal control theory of systems governed by partial differential equations [1], the voltages applied to the poloidal field circuits that achieve at best the desired scenario, by minimizing a certain cost-function which represents the sum of the distance to the desired plasma and of the energetic cost of the electrical system. The introduction of an appropriate lagrangian taking as constraints the equilibrium system of the previous sections and the determination of the corresponding adjoint state enable the computation of the gradient of the costfunction in terms of the adjoint state. The minimization of this cost-function is performed thanks to a SQP (Sequential Quadratic Programming) method. An interesting test-case, solved by using these techniques, will be presented for the ITER (International Thermonuclear Experimental Reactor) tokamak. This is presented in section 5 of this paper. This method has the purpose to replace the empirical methods used commonly to compute the pre-programmed voltages that enable to go from one snapshot to another one. This method can of course be extended to other type of optimization of the scenarios just by modifying the cost-function and the control variables (consumption of flux, desired profile of plasma current density,..).

The magnetohydrodynamic equations
A plasma is a ionized gas composed of ions and electrons. The kinetic equations describe the plasma thanks to a distribution function f α (x, v, t) (with α = e for electrons and α = i for ions) where x is the point position and v the particles velocity. For a collisional plasma the kinetic equations are based on the Fokker-Planck equation where m α is the mass of the particles, F α the force applied to these particles and C α the term due to collisions between particles. This microscopic approach requires the resolution of a partial differential equations in 6 dimensions (space and velocity) plus the time dimension. This is extremely difficult from a computational point of view. Therefore from this equation one derives a macroscopic representation based on the fluid equations in the following way. Let us define the density of particles by the fluid velocity by and the pressure tensor which under the isotropic assumption becomes Multiplying equation (1) by a test function φ(w) and integrating over the space of velocities leads to the fluid equations. The first moment (corresponding to φ = 1) gives the equation for the density of particles: Since for electromagnetic forces ∂F α ∂w = 0, and since collisions do not change the number of particles one obtains: The second moment is obtained by taking φ = m α w which leads to the momentum equation where we have set w = (w − u α ) + u α . Using the equation for conservation of the density one gets where Ze is the charge of particles and R α is the change rate of the momentum due to collisions. The third moment gives the energy equation which needs to be complemented with closing relations on the heat flux. These latter come from a transport model. Finally the resistive MHD equations for a single fluid [2] read: where n denotes the density of the particles, m their mass, u their mean velocity, p their pressure, T their temperature, Q the heat flux, η the resistivity tensor, s and s the source terms and k the Boltzmann constant.

Equilibrium of a plasma in a tokamak
In order to simplify system (2) some characteristic time constants of the plasma need to be defined. The Alfven time constant τ A is where a is the minor radius of the plasma and B 0 is the toroidal magnetic field. It is of the order of a microsecond for present tokamaks. The diffusion time constant of the particle density n is where D is the particle diffusion coefficient. Likewise, the time constants for diffusion of heat of the electrons and of the ions are where n e , n i are the density of electrons and ions, respectively, and K e , K i are their thermal conductivities. These constants τ n , τ e , τ i are of the order of a millisecond on tokamaks currently operating. Finally, the resistive time constant for the diffusion of the current density and magnetic field in the plasma is given by and is of the order of a second. If a global time constant for plasma diffusion is defined by On the diffusion time-scale τ p the term ( ∂u ∂t + u∇u) is small compared with ∇p (see [3,4]) and the equilibrium equation is thus satisfied at every instant in the plasma. Consequently the equations which govern the equilibrium of a plasma in the presence of a magnetic field in a tokamak are on the one hand Maxwell's equations satisfied in the whole of space (including the plasma): and on the other hand the equilibrium equation (3) for the plasma itself. Equation (3) means that the plasma is in equilibrium when the force ∇p due the kinetic pressure p is equal to the Lorentz force of the magnetic pressure j × B. We deduce immediately from (3) that and j · ∇p = 0.
Thus for a plasma in equilibrium the field lines and the current lines lie on isobaric surfaces (p = const.); these surfaces, generated by the field lines, are called magnetic surfaces. In order for them to remain within a bounded volume of space it is necessary that they have a toroidal topology. These surfaces form a family of nested tori. The innermost torus degenerates into a curve which is called the magnetic axis. In a cylindrical coordinate system (r, φ, z) (where r = 0 is the major axis of the torus) the hypothesis of axial symmetry consists in assuming that the magnetic field B is independent of the toroidal angle φ. The magnetic field can be decomposed as Concerning the toroidal component B φ we define f by where e φ is the unit vector in the toroidal direction, and f is the diamagnetic function. The magnetic field can be written as: According to (9), in an axisymmetric configuration the magnetic surfaces are generated by the rotation of the flux lines ψ = const. around the axis r = 0 of the torus. From (9) and the second relation of (4) we obtain the following expression for j: where j p and j φ are the poloidal and toroidal components respectively of j, and the operator ∆ * is defined by Expressions (9) and (10) for B and j are valid in the whole of space since they involve only Maxwell's equations and the hypothesis of axisymmetry. Hence they can be reduced to one equation given in 2 space dimensions in the poloidal plane (r, z) ∈ Ω ∞ = (0, ∞) × (−∞, ∞) for the poloidal flux ψ: The toroidal component of the current density j φ is zero everywhere outside the plasma domain, the poloidal field coils and the passive structures. The different sub-domains of the poloidal plane of a tokamak (see Fig. 2) as well as the corresponding expression for j φ are described below: • Ω L is the domain accessible to the plasma. Its boundary is the limiter ∂Ω L .
• Ω p is the plasma domain where relation (5) implies that ∇p and ∇ψ are co-linear, and therefore p is constant on each magnetic surface. This can be denoted by p = p(ψ).
Relation (6) combined with the expression (10) implies that ∇f and ∇p are co-linear, and therefore f is likewise constant on each magnetic surface The equilibrium relation (3) combined with the expression (9) and (10) for B and j implies that: which leads to the so-called Grad-Shafranov equilibrium equation: Here µ is equal to the magnetic permeability µ 0 of the vacuum and ∆ * is a linear elliptic operator. From (10) it is clear that right-hand side of (16) represents the toroidal component of the plasma current density. It involves functions p(ψ) and f (ψ) which are not directly measured inside the plasma. The plasma domain is unknown, Ω p = Ω p (ψ), and this is a free boundary problem. This domain is defined by its boundary which is the largest closed ψ iso-contour contained within the limiter Ω L . The plasma can either be limited if this iso-contour is tangent to the limiter ∂Ω L (see figure 3, left) or defined by the presence of a saddle-point also called X-point (see figure 3, right). In the later configuration which is obtained in presence of a divertor, the plasma does not touch any physical component and the performances and the confinement of the plasma are improved (see [5]).
• Ω Fe represents the ferromagnetic structures. They do not carry any current, j φ = 0 but the magnetic permeability µ is not constant and depends on the magnetic field: • Domains Ω ci represent the poloidal field coils carrying currents. If we consider that the voltages V i applied to these coils are given, using Faraday and Ohm laws the current density can be written as where n i is the number of windings in the coil, |Ω ci | its section area, R i its resistance andψ is the time derivative of ψ, • Ω ps represents passive structures where the current density can be written as where σ is the conductivity.
In summary we are seeking for the poloidal flux ψ(t) that is a solution of (12) with j φ given by (16)

Weak Formulation and Discretization
We chose a semi-circle Γ of radius ρ Γ surrounding the iron domain Ω Fe , the coil domains Ω ci and the passive structures domain Ω ps . The truncated domain, we use for our computations, is the domain Ω having the boundary ∂Ω = Γ ∪ Γ 0 , where Γ 0 := {(0, z), z min ≤ z ≤ z max }. The weak formulation for ψ(t) uses the following Sobolev space: It reads as: Given and c(ψ, ξ) := 1 µ 0 Γ ψ(P 1 )N (P 1 )ξ(P 1 )dS 1 where P i = (r i , z i ) and K and E the complete elliptic integrals of first and second kind, respectively and The bilinear form c : H × H → R is accounting for the boundary conditions at infinity [6]. We refer to [7,Chapter 2.4] for the details of the derivation. The bilinear form c(·, ·) follows basically from the so-called uncoupling procedure in [8] for the usual coupling of boundary integral and finite element methods. As we focus here on the equilibrium problem the two functions p and f f have to be supplied as data, called S p and S f f in the definition of J p (ψ, ξ). While the domain of p and f f depends on the poloidal flux itself, it is more practical to supply those profiles S p and S f f as functions of the normalized poloidal flux ψ N (r, z): where ψ ax (ψ) := ψ(r ax (ψ), z ax (ψ)), ψ bnd (ψ) := ψ(r bd (ψ), z bd (ψ)) (24) with (r ax (ψ), z ax (ψ)) the magnetic axis, where ψ has its global maximum in Ω L and (r bnd (ψ), z bnd (ψ)) the coordinates of the point that determines the plasma boundary. The point (r bnd , z bnd ) is either an X-point of ψ or the contact point with the limiter ∂Ω L . S p and S f f , have, independently of ψ, a fixed domain [0, 1] and are usually given as (piecewise) polynomial functions. Another frequent a priori model is with r 0 the major radius of the vacuum chamber and α, β, γ ∈ R given parameters. We refer to [9] for a physical interpretation of these parameters. The parameter β is related to the poloidal beta, whereas α and γ describe the peakage of the current profile and λ is a normalization factor.
Numerical Methods It is straightforward to combine Galerkin methods in space and time-stepping schemes to get approximation schemes for solving (20) numerically. For the choice of the spatial discretization, the fine details of realistic tokamak sections (see Figure 4) give here favor to finite element spaces based on triangular meshes. Since for many years now the piecewise affine approximations are the standard choice for the stationary free-boundary equilibrium problems [10,11,7], we stay also here with linear Lagrangian finite elements for the discretization in space. Higher order methods are likewise implementable. In order to prohibit numerical instablities it is advisable to use implicit time-stepping methods such as implicit Euler, which leads to non-linear finitedimensional problems. The Newton-type methods for solving such non-linear problems can be based on the Gâteaux derivative of A(ψ, ξ) and the Gâteaux derivative where Γ p is the plasma boundary ∂Ω p and The derivation of the linearization D ψ J p (ψ, ξ)( ψ) requires to assume that ∇ψ = 0 on ∂Ω p and involves shape calculus [12,13] and the non-trivial derivatives: D ψ ψ ax (ψ)(ψ) =ψ(r ax (ψ), z ax (ψ)) and D ψ ψ bd (ψ)(ψ) =ψ(r bd (ψ), z bd (ψ)).
Clearly, ∇ψ = 0 on ∂Ω p will not be true for the nowadays important X-point equilibria. Nevertheless this theoretical difficulty is not very essential for practical computations. In [14] it is pointed out that accurate Newton methods for discretized versions of the weak formulation (20) need to use accurate derivatives for the discretized non-linear operator, which is not necessarily equal to the discretization of the analytical derivatives. Here, the discretization and linearization of J p (ψ, ξ) needs special attention due to the ψ-dependent domain of integration. We refer to [

The optimal control problem
We intend to determine the voltages V i (t) applied to the poloidal field circuits so that the plasma boundary Γ p fit to a desired boundary Γ desi during the whole discharge while minimizing a certain energetic cost. Let Γ desi (t) ⊂ Ω L denote the evolution of a closed line, contained in the domain Ω L that is either smooth and touches the limiter at one point or has at least one corner. The former case prescribes a desired plasma boundary that touches the limiter. The latter case aims at a plasma with X-point that is entirely in the interior of Ω L . Further let (r desi (t), z desi (t)) ∈ Γ desi (t) and (r 1 (t), z 1 (t)), . . . , (r N desi (t), z N desi (t)) ∈ Γ desi (t) be N desi + 1 points on that line. We define a quadratic functional K(ψ) that evaluates to zero if Γ desi (t) is an Another functional, that will serve as regularization, is with regularization weights w i ≥ 0. The regularization functional penalizes the strength of the voltages V i and represents the energetic cost in the coil system. We consider the following minimization problem: subject to This minimization problem for transient axisymmetric equilibria extends the minimization problems for static axisysmmetric equilibria introduced in [15, Chapter II]. Hence, theoretical assertions for (30) such as the first order necessary conditions for optimality follow by similar arguments as those in [15, page 80-84].
The Lagrangian for the optimization problem (30) with Lagrange multiplier φ is: We can state the first order necessary conditions for optimality under the following three assumptions in the limiter case: 1. sup ΩL ψ is attained at one and only one point M 0 = (r bd , z bd ). ψ(t) and φ(t) are solution of the adjoint problem

sup
with φ(T ) = 0 and -V(t) and φ(t) are solution to The adjoint problem has the following strong formulation: with φ(T ) = 0, where δ ax and δ bd are the Dirac masses at the points (r ax , r ax ) and (r bd , r bd ), respectively. δ Γp is the Dirac mass of Γ p with Equation (32) is the Euler equation for the minimization of (30). Equations (20), (31) and (32) constitute the optimality system for problem (30).
Numerical Methods The discretization of our minimization problem (30) builds on the space-time discretization for (20) that we outlined in the previous section. Next, the discrete minimization problem can be recast as the following constrained optimization problem min u,y where y and u are the so-called state and control variables. In our setting y will be the variable that describes the plasma and u will be the externally applied voltages. We think of y as the vector of degrees of freedoms describing the space and time evolution of the poloidal flux ψ, and B(y) and F(u) are the discretizations of the non-linear operators in the variational formulation (20). Sequential Quadratic Programming (SQP) is one of the most effective methods for nonlinear constrained optimization with significant non-linearities in the constraints [16,Chapter 18]. SQP methods find a numerical solution by generating iteration steps that minimize quadratic cost functions subject to linear constraints. The Lagrange function formalism in combination with Newton-type iterations is one approach to derive the SQP-methods: the Lagrangian for (33) is and the solution of (33) is a stationary point of this Lagrangian: If the linear systems in (36) become too large, we are pursuing the null space approach to arrive at the SQP formulation with the reduced Hessian for the increment ∆u k := u k+1 − u k : where We are using iterative methods, e.g. the conjugate gradient methods, to solve (37). Since in our case the number of control variables will be small we can expect convergence within very few iterations. Within each iteration step of the iterative method, we still have to solve the two linear systems corresponding to D y B(y k ) and D y B T (y k ). Alternatively, if we have sufficient memory to store M(·, ·), we can compute M(·, ·) explicitly. Clearly, we never compute neither D y B −1 (y k ) nor D y B −T (y k ) explicitly.
Once we know ∆u k we can compute y k+1 and p k+1 by: y k+1 − y k = D y B −1 (y k )D u F(u k )∆u k − r(u k , y k ), p k+1 + λ k = −D y B −T (y k )(H k y,y (y k+1 − y k ) + H k y,u (u k+1 − u k )). We would like to highlight that the SQP-method relies on proper derivatives of the non-linear operators B and F. In our case F is affine, hence the derivative of B remains the most difficult part. On the other hand these are exactly the same terms that appear in the Newton iterations for the direct problem (20) and we can reuse the methodology presented at the end of Section 4. For practical purposes we do neglect all involved second order derivatives of B.
It is very instrumental to compare the expression involved in the reduced formulation (37) of SQP to the gradient and the Hessian of the reduced cost function, that would appear when using algorithms for unconstrained optimization problems.
Let J(u) := J(y(u), u), with B(y(u)) = F(u) be the reduced cost function, then we have the following expressions for gradient Hence, the reduced gradient h(y k , u k ) is not the gradient of the reduced cost function, unless the state and control variable y k and u k verify the equation of state B(y k ) = F(u k ).
Preliminary Example Finally, we would like to show first results for a socalled ramp-up scenario in an ITER-like tokamak, where the plasma evolves from a small circular to a large elongated plasma. The optimal coil voltages are depicted in Figure 5. Then, if we use those as data to solve the direct problem we verify that the plasma boundary follows indeed the prescribed trajectory (see Figure 6).

Conclusion
The study and the optimization of scenarios is more and more important for the realization of objectives in magnetic confinement controlled fusion and will certainly be crucial for the ITER project. The first results presented in this paper are very encouraging and are the starting point of the development of new tools devoted to the preparation of scenarios of the future devices.   6. Optimal control for a ramp-up scenario: the plasma boundary (green) follows the prescribed boundary (black points), snapshots at t = 0, 2, 6, 10, 20, 30, 40, 45, 50, 54, 58, 60s (from left to right, top to down).