# Control Methods for the Optimization of Plasma Scenarios in a Tokamak

Jacques Blum, Cédric Boulbe, Blaise Faugeras, Holger Heumann

► To cite this version:

Jacques Blum, Cédric Boulbe, Blaise Faugeras, Holger Heumann. Control Methods for the Optimization of Plasma Scenarios in a Tokamak. 27th IFIP Conference on System Modeling and Optimization (CSMO), Jun 2015, Sophia Antipolis, France. pp.1-20, 10.1007/978-3-319-55795-3_1. hal-01626908

HAL Id: hal-01626908

https://inria.hal.science/hal-01626908

Submitted on 31 Oct 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Control methods for the optimization of plasma
cenarios in a tokamak

Jacques Blum, Cédric Boulbe, Blaise Faugeras, and Holger Heumann

Laboratoire J.-A. Dieudonné,
Université de Nice Sophia-Antipolis
Parc Valrose 06108 Nice Cedex 02, France
and
INRIA Castor project
jacques.blum@unice.fr

**Abstract.** 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 cost-function 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.

**Keywords:** plasma equilibrium, optimization, tokamak, magnetohydrodynamics, optimal control

# 1 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:
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic representation of a tokamak</figcaption></figure>

- 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 cost-function 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,...).

## 2 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_\alpha(\mathbf{x}, \mathbf{v}, t)$ (with $\alpha = e$ for electrons and $\alpha = i$ for ions) where $\mathbf{x}$ is the point position and $\mathbf{v}$ the particles velocity. For a collisional plasma the kinetic equations are based on the Fokker-Planck equation

$$ \frac{\partial f_{\alpha}}{\partial t} + (\mathbf{v} \cdot \nabla_{\mathbf{x}}) f_{\alpha} + \frac{\mathbf{F}_{\alpha}}{m_{\alpha}} \cdot \nabla_{\mathbf{v}} f_{\alpha} = C_{\alpha}, \qquad (1) $$

where $m_\alpha$ is the mass of the particles, $\mathbf{F}_\alpha$ the force applied to these particles and $C_\alpha$ 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

$$ n_{\alpha}(\mathbf{x}, t) = \int f_{\alpha}(\mathbf{x}, \mathbf{w}, t) d\mathbf{w}, $$
the fluid velocity by

$$ \mathbf{u}_{\alpha}(\mathbf{x}, t) = \frac{1}{n_{\alpha}} \int f_{\alpha}(\mathbf{x}, \mathbf{w}, t) \mathbf{w} d\mathbf{w}, $$

and the pressure tensor

$$ P_{\alpha}(\mathbf{x}, t) = m_{\alpha} \int f_{\alpha}(\mathbf{x}, \mathbf{w}, t)(\mathbf{w} - \mathbf{u}_{\alpha})(\mathbf{w} - \mathbf{u}_{\alpha})d\mathbf{w}, $$

which under the isotropic assumption becomes

$$ p_{\alpha}(\mathbf{x}, t) = \frac{m_{\alpha}}{3} \int f_{\alpha}(\mathbf{x}, \mathbf{w}, t)(\mathbf{w} - \mathbf{u}_{\alpha})^{2} d\mathbf{w}. $$

Multiplying equation (1) by a test function $\phi(\mathbf{w})$ and integrating over the space of velocities leads to the fluid equations. The first moment (corresponding to $\phi = 1$) gives the equation for the density of particles:

$$ \frac{\partial n_{\alpha}}{\partial t} + \nabla \cdot \int f_{\alpha} \mathbf{w} d\mathbf{w} - \frac{1}{m_{\alpha}} \int \frac{\partial \mathbf{F}_{\alpha}}{\partial \mathbf{w}} f_{\alpha} d\mathbf{w} = 0. $$

Since for electromagnetic forces $\frac{\partial \mathbf{F}_{\alpha}}{\partial \mathbf{w}} = 0$, and since collisions do not change the number of particles one obtains:

$$ \frac{\partial n_{\alpha}}{\partial t} + \nabla \cdot (n_{\alpha} \mathbf{u}_{\alpha}) = 0. $$

The second moment is obtained by taking $\phi = m_{\alpha} \mathbf{w}$ which leads to the momentum equation

$$ m_{\alpha} \frac{\partial}{\partial t} (n_{\alpha} \mathbf{u}_{\alpha}) + m_{\alpha} \nabla_{x} \cdot \int f_{\alpha} \mathbf{w} \mathbf{w} d\mathbf{w} - \int \nabla_{\mathbf{w}} \cdot (\mathbf{F}_{\alpha} \cdot \mathbf{w}) f_{\alpha} d\mathbf{w} = \int m_{\alpha} \mathbf{w} C_{\alpha} d\mathbf{w}, $$

where we have set $\mathbf{w} = (\mathbf{w} - \mathbf{u}_{\alpha}) + \mathbf{u}_{\alpha}$. Using the equation for conservation of the density one gets

$$ m_{\alpha} n_{\alpha} \left( \frac{\partial \mathbf{u}_{\alpha}}{\partial t} + \mathbf{u}_{\alpha} \cdot \nabla \mathbf{u}_{\alpha} \right) = - \nabla \cdot P_{\alpha} + n_{\alpha} \overline{\mathbf{F}}_{\alpha} + \mathbf{R}_{\alpha}, $$

with $\overline{\mathbf{F}}_{\alpha} = Ze(\mathbf{E} + \mathbf{u}_{\alpha} \times \mathbf{B})$ where $Ze$ is the charge of particles and $\mathbf{R}_{\alpha}$ 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.

The single fluid magnetohydrodynamic equations are derived by defining the mass density

$$ \begin{aligned} mn &= m_e n_e + m_i n_i \\ &= m_e Z n_i + m_i n_i \approx m_i n_i \end{aligned} $$
the velocity of the fluid

$$ \boldsymbol{u} = \frac{m_e n_e \boldsymbol{u}_e + m_i n_i \boldsymbol{u}_i}{\rho} \approx \boldsymbol{u}_i, $$

the current density

$$ \begin{aligned} \boldsymbol{j} &= -en_e \boldsymbol{u}_e + Z e n_i \boldsymbol{u}_i, \\ &= en_e (\boldsymbol{u}_i - \boldsymbol{u}_e), \end{aligned} $$

and the scalar pressure

$$ p = n_e k T_e + n_i k T_i, $$

where $k$ is the Boltzmann constant. The Maxwell equations need to be added since we are in the presence of a magnetic field **B** and of an electric field **E**. Finally the resistive MHD equations for a single fluid [2] read:

$$ \left\{ \begin{aligned} & \frac{\partial n}{\partial t} + \nabla \cdot (n \boldsymbol{u}) = s && \text{(Conservation of particles)} \\ & mn \left( \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \cdot \nabla \boldsymbol{u} \right) + \nabla p = \boldsymbol{j} \times \boldsymbol{B} && \text{(Conservation of momentum)} \\ & \frac{3}{2} \left( \frac{\partial p}{\partial t} + \boldsymbol{u} \cdot \nabla p \right) + \frac{5}{2} p \nabla \cdot \boldsymbol{u} + \nabla Q = s' && \text{(Conservation of particle energy)} \\ & \nabla \times \boldsymbol{E} = -\frac{\partial \boldsymbol{B}}{\partial t} && \text{(Faraday's law)} \\ & \nabla \cdot \boldsymbol{B} = 0 && \text{(Conservation of $\boldsymbol{B}$)} \\ & \boldsymbol{E} + \boldsymbol{u} \times \boldsymbol{B} = \eta \boldsymbol{j} && \text{(Ohm's law)} \\ & \nabla \times \boldsymbol{H} = \boldsymbol{j} && \text{(Ampere's law)} \\ & \boldsymbol{B} = \mu \boldsymbol{H} && \text{(Magnetic permeability)} \\ & p = nkT && \text{(Law of perfect gases)} \end{aligned} \right. \quad (2) $$

where $n$ denotes the density of the particles, $m$ their mass, $\boldsymbol{u}$ their mean velocity, $p$ their pressure, $T$ their temperature, $Q$ the heat flux, $\eta$ the resistivity tensor, $s$ and $s'$ the source terms and $k$ the Boltzmann constant.

### 3 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 $\tau_A$ is

$$ \tau_A = \frac{a(\mu_0 mn)^{1/2}}{B_0}, $$

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

$$ \tau_n = \frac{a^2}{D}, $$
where $D$ is the particle diffusion coefficient. Likewise, the time constants for diffusion of heat of the electrons and of the ions are

$$ \tau_e = \frac{n_e a^2}{K_e}, $$

$$ \tau_i = \frac{n_i a^2}{K_i}, $$

where $n_e, n_i$ are the density of electrons and ions, respectively, and $K_e, K_i$ are their thermal conductivities. These constants $\tau_n, \tau_e, \tau_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

$$ \tau_r = \frac{\mu_0 a^2}{\eta}, $$

and is of the order of a second.

If a global time constant for plasma diffusion is defined by

$$ \tau_p = \inf(\tau_n, \tau_e, \tau_i, \tau_r), $$

we note that

$$ \tau_A \ll \tau_p. $$

On the diffusion time-scale $\tau_p$ the term $(\frac{\partial u}{\partial t} + u \nabla u)$ is small compared with $\nabla p$ (see [3, 4]) and the equilibrium equation

$$ \nabla p = j \times B \qquad (3) $$

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):

$$ \left\{ \begin{array}{l} \nabla \cdot B = 0, \\ \nabla \times \left( \frac{B}{\mu} \right) = j, \end{array} \right. \qquad (4) $$

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 $\nabla p$ due the kinetic pressure $p$ is equal to the Lorentz force of the magnetic pressure $j \times B$. We deduce immediately from (3) that

$$ B \cdot \nabla p = 0, \qquad (5) $$

and

$$ j \cdot \nabla p = 0. \qquad (6) $$
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, \phi, 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 $\phi$. The magnetic field can be decomposed as $B = B_p+B_\phi$, where $B_p = (B_r, B_z)$ is the poloidal component and $B_\phi$ is the toroidal component. From equation (4) one can define the poloidal flux $\psi(r, z)$ such that

$$ \left\{ \begin{aligned} B_r &= -\frac{1}{r} \frac{\partial \psi}{\partial z}, \\ B_z &= \frac{1}{r} \frac{\partial \psi}{\partial r}. \end{aligned} \right. \qquad (7) $$

Concerning the toroidal component $B_\phi$ we define $f$ by

$$ B_{\phi} = \frac{f}{r} e_{\phi}, \qquad (8) $$

where $e_\phi$ is the unit vector in the toroidal direction, and $f$ is the diamagnetic function. The magnetic field can be written as:

$$ \left\{ \begin{aligned} B &= B_p + B_\phi, \\ B_p &= \frac{1}{r} [\nabla \psi \times e_\phi], \\ B_\phi &= \frac{f}{r} e_\phi. \end{aligned} \right. \qquad (9) $$

According to (9), in an axisymmetric configuration the magnetic surfaces are generated by the rotation of the flux lines $\psi = const.$ around the axis $r = 0$ of the torus.

From (9) and the second relation of (4) we obtain the following expression for $j$:

$$ \left\{ \begin{aligned} j &= j_p + j_\phi, \\ j_p &= \frac{1}{r} [\nabla(\frac{f}{\mu}) \times e_\phi], \\ j_\phi &= (-\Delta^* \psi) e_\phi, \end{aligned} \right. \qquad (10) $$

where $j_p$ and $j_\phi$ are the poloidal and toroidal components respectively of $j$, and the operator $\Delta^*$ is defined by

$$ \Delta^* = \partial_r \left( \frac{1}{\mu r} \partial_r \cdot \right) + \partial_z \left( \frac{1}{\mu r} \partial_z \cdot \right) = \nabla \left( \frac{1}{\mu r} \nabla \cdot \right). \qquad (11) $$
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) \in \Omega_{\infty} = (0, \infty) \times (-\infty, \infty)$ for the poloidal flux $\psi$:

$$-\Delta^*\psi = j_{\phi}. \qquad (12)$$

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Schematic representation of the poloidal plane of a tokamak. $\Omega_p$ is the plasma domain, $\Omega_L$ is the limiter domain accessible to the plasma, $\Omega_{c_i}$ represent poloidal field coils, $\Omega_{ps}$ the passive structures and $\Omega_{Fe}$ the ferromagnetic structures.</figcaption></figure>

The toroidal component of the current density $j_{\phi}$ 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_{\phi}$ are described below:

* $\Omega_L$ is the domain accessible to the plasma. Its boundary is the limiter $\partial\Omega_L$.

* $\Omega_p$ is the plasma domain where relation (5) implies that $\nabla p$ and $\nabla \psi$ are co-linear, and therefore $p$ is constant on each magnetic surface. This can be denoted by

$$p = p(\psi). \qquad (13)$$

Relation (6) combined with the expression (10) implies that $\nabla f$ and $\nabla p$ are co-linear, and therefore $f$ is likewise constant on each magnetic surface

$$f = f(\psi). \qquad (14)$$

The equilibrium relation (3) combined with the expression (9) and (10) for **B** and **j** implies that:

$$\nabla p = -\frac{\Delta^*\psi}{r}\nabla\psi - \frac{f}{\mu_0 r^2}\nabla f, \qquad (15)$$
which leads to the so-called Grad-Shafranov equilibrium equation:

$$-\Delta^* \psi = r p'(\psi) + \frac{1}{\mu_0 r} (f f')(\psi). \qquad (16)$$

Here $\mu$ is equal to the magnetic permeability $\mu_0$ of the vacuum and $\Delta^*$ 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(\psi)$ and $f(\psi)$ which are not directly measured inside the plasma. The plasma domain is unknown, $\Omega_p = \Omega_p(\psi)$, and this is a free boundary problem. This domain is defined by its boundary which is the largest closed $\psi$ iso-contour contained within the limiter $\Omega_L$. The plasma can either be limited if this iso-contour is tangent to the limiter $\partial \Omega_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]).

* $\Omega_{Fe}$ represents the ferromagnetic structures. They do not carry any current, $j_{\phi} = 0$ but the magnetic permeability $\mu$ is not constant and depends on the magnetic field:

$$\mu = \mu_{Fe} \left( \frac{|\nabla \psi|^2}{r^2} \right). \qquad (17)$$

* Domains $\Omega_{c_i}$ 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

$$j_{\phi} = \frac{n_i V_i}{R_i |\Omega_{c_i}|} - \frac{2\pi n_i^2}{R_i |\Omega_{c_i}|^2} \int_{\Omega_{c_i}} \dot{\psi} ds, \qquad (18)$$

where $n_i$ is the number of windings in the coil, $|\Omega_{c_i}|$ its section area, $R_i$ its resistance and $\dot{\psi}$ is the time derivative of $\psi$,

* $\Omega_{ps}$ represents passive structures where the current density can be written as

$$j_{\phi} = -\frac{\sigma}{r} \dot{\psi}, \qquad (19)$$

where $\sigma$ is the conductivity.

In summary we are seeking for the poloidal flux $\psi(t)$ that is a solution of (12) with $j_{\phi}$ given by (16), (18) and (19) and verifies boundary conditions

$$\psi(0, z) = 0 \quad \text{and} \quad \lim_{\|(r,z)\| \to +\infty} \psi(r, z) = 0.$$

## 4 Weak Formulation and Discretization

We chose a semi-circle $\Gamma$ of radius $\rho_{\Gamma}$ surrounding the iron domain $\Omega_{Fe}$, the coil domains $\Omega_{c_i}$ and the passive structures domain $\Omega_{ps}$. The truncated domain, we
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Example of a plasma whose boundary is defined by the contact with limiter (left) or by the presence of an X-point (right).</figcaption></figure>

use for our computations, is the domain $\Omega$ having the boundary $\partial\Omega = \Gamma \cup \Gamma_0$, where $\Gamma_0 := \{(0, z), z_{min} \le z \le z_{max}\}$. The weak formulation for $\psi(t)$ uses the following Sobolev space:

$$H := \left\{ \psi : \Omega \to \mathbb{R}, \|\psi\| < \infty, \left\| \frac{|\nabla \psi|}{r} \right\| < \infty, \psi|_{\Gamma_0} = 0 \right\} \cap C^0(\bar{\Omega}),$$

with

$$\|\psi\|^2 = \int_{\Omega} \psi^2 r drdz.$$

It reads as: Given $V(t) = \{V_i(t)\}_{i=1}^N$ find $\psi(t) \in H$ such that for all $\xi \in H$

$$A(\psi(t), \xi) - J_p(\psi(t), \xi) + j^{\text{ps}}(\dot{\psi}(t), \xi) + j^{\text{c}}(\dot{\psi}(t), \xi) + c(\psi(t), \xi) = \ell(V(t), \xi), \quad (20)$$

where

$$A(\psi, \xi) := \int_{\Omega} \frac{1}{\mu(\psi)r} \nabla\psi \cdot \nabla\xi \, drdz,$$

$$J_p(\psi, \xi) := \int_{\Omega_p(\psi)} \left( r S_{p'}(\psi_N) + \frac{1}{\mu_0 r} S_{ff'}(\psi_N) \right) \xi \, drdz,$$

$$\ell(V(t), \xi) := \sum_{i=1}^{N} \frac{n_i}{R_i |\Omega_{c_i}|} V_i(t) \int_{\Omega_{c_i}} \xi \, drdz, \qquad (21)$$

$$j^{\text{ps}}(\psi, \xi) := \int_{\Omega_{\text{ps}}} \frac{\sigma}{r} \psi \xi \, drdz,$$

$$j^{\text{c}}(\psi, \xi) := \sum_{i=1}^{N_i} \frac{2\pi n_i^2}{R_i |\Omega_{c_i}|^2} \int_{\Omega_{c_i}} \psi \, drdz \int_{\Omega_{c_i}} \xi \, drdz,$$
and

$$ \mathbf{c}(\psi, \xi) := \frac{1}{\mu_0} \int_{\Gamma} \psi(\mathbf{P}_1) N(\mathbf{P}_1) \xi(\mathbf{P}_1) dS_1 \\ + \frac{1}{2\mu_0} \int_{\Gamma} \int_{\Gamma} (\psi(\mathbf{P}_1) - \psi(\mathbf{P}_2)) M(\mathbf{P}_1, \mathbf{P}_2) (\xi(\mathbf{P}_1) - \xi(\mathbf{P}_2)) dS_1 dS_2, \quad (22) $$

with

$$ M(\mathbf{P}_1, \mathbf{P}_2) = \frac{k_{\mathbf{P}_1, \mathbf{P}_2}}{2\pi(r_1 r_2)^{\frac{3}{2}}} \left( \frac{2 - k_{\mathbf{P}_1, \mathbf{P}_2}^2}{2 - 2k_{\mathbf{P}_1, \mathbf{P}_2}^2} E(k_{\mathbf{P}_1, \mathbf{P}_2}) - K(k_{\mathbf{P}_1, \mathbf{P}_2}) \right), \\ N(\mathbf{P}_1) = \frac{1}{r_1} \left( \frac{1}{\delta_+} + \frac{1}{\delta_-} - \frac{1}{\rho_\Gamma} \right) \text{ and } \delta_\pm = \sqrt{r_1^2 + (\rho_\Gamma \pm z_1)^2}. $$

where $\mathbf{P}_i = (r_i, z_i)$ and $K$ and $E$ the complete elliptic integrals of first and second kind, respectively and

$$ k_{\mathbf{P}_j, \mathbf{P}_k} = \sqrt{\frac{4r_j r_k}{(r_j + r_k)^2 + (z_j - z_k)^2}}. $$

The bilinear form $\mathbf{c}: H \times H \to \mathbb{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 $\mathbf{c}(\cdot, \cdot)$ 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_{ff'}$ in the definition of $J_p(\psi, \xi)$. 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_{ff'}$ as functions of the normalized poloidal flux $\psi_N(r, z)$:

$$ \psi_N(r, z) = \frac{\psi(r, z) - \psi_{\text{ax}}(\psi)}{\psi_{\text{bd}}(\psi) - \psi_{\text{ax}}(\psi)}, \quad (23) $$

where

$$ \begin{aligned} \psi_{\text{ax}}(\psi) &:= \psi(r_{\text{ax}}(\psi), z_{\text{ax}}(\psi)), \\ \psi_{\text{bnd}}(\psi) &:= \psi(r_{\text{bd}}(\psi), z_{\text{bd}}(\psi)) \end{aligned} \quad (24) $$

with $(r_{\text{ax}}(\psi), z_{\text{ax}}(\psi))$ the magnetic axis, where $\psi$ has its global maximum in $\Omega_L$ and $(r_{\text{bnd}}(\psi), z_{\text{bnd}}(\psi))$ the coordinates of the point that determines the plasma boundary. The point $(r_{\text{bnd}}, z_{\text{bnd}})$ is either an X-point of $\psi$ or the contact point with the limiter $\partial\Omega_L$. $S_{p'}$ and $S_{ff'}$, have, independently of $\psi$, a fixed domain $[0, 1]$ and are usually given as (piecewise) polynomial functions. Another frequent a priori model is

$$ S_{p'}(\psi_N) = \lambda \frac{\beta}{r_0} (1 - \psi_N^\alpha)^\gamma, \quad S_{ff'}(\psi_N) = \lambda (1 - \beta) \mu_0 r_0 (1 - \psi_N^\alpha)^\gamma \quad (25) $$

with $r_0$ the major radius of the vacuum chamber and $\alpha, \beta, \gamma \in \mathbb{R}$ given parameters. We refer to [9] for a physical interpretation of these parameters. The parameter $\beta$ is related to the poloidal beta, whereas $\alpha$ and $\gamma$ describe the leakage of the current profile and $\lambda$ 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.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** The different subdomains of the geometry of the tokamak WEST (left) and ITER (right) and triangulations that resolve the geometric details.</figcaption></figure>

In order to prohibit numerical instabilities it is advisable to use implicit time-stepping methods such as implicit Euler, which leads to non-linear finite-dimensional problems. The Newton-type methods for solving such non-linear problems can be based on the Gâteaux derivative

$$D_{\psi}A(\psi, \xi)(\tilde{\psi}) = \int_{\Omega} \frac{1}{\mu(\psi)r} \nabla \tilde{\psi} \cdot \nabla \xi \, drdz \\ - 2 \int_{\Omega_{\text{Fe}}} \frac{\mu'_{\text{Fe}}(\frac{|\nabla \psi|^2}{r^2})}{\mu_{\text{Fe}}^2(\frac{|\nabla \psi|^2}{r^2})} r^3 (\nabla \tilde{\psi} \cdot \nabla \psi)(\nabla \psi \cdot \nabla \xi) \, drdz$$
of A(ψ, ξ) and the Gâteaux derivative

$$
\begin{equation}
\begin{split}
D_{\psi} J_{\mathrm{p}}(\psi, \xi)(\tilde{\psi}) = & \int_{\Omega_{\mathrm{p}}(\psi)} \frac{\partial j_{\mathrm{p}}(r, \psi_{\mathrm{N}}(\psi))}{\partial \psi_{\mathrm{N}}} \frac{\partial \psi_{\mathrm{N}}(\psi)}{\partial \psi} \tilde{\psi} \, \xi \, dr \, dz \\
& - \int_{\Gamma_{\mathrm{p}}(\psi)} j_{\mathrm{p}}(r, 1) |\nabla \psi|^{-1} (\tilde{\psi} - \tilde{\psi}(r_{\mathrm{bd}}(\psi), z_{\mathrm{bd}}(\psi))) \, \xi \, d\Gamma \\
& + \int_{\Omega_{\mathrm{p}}(\psi)} \frac{\partial j_{\mathrm{p}}(r, \psi_{\mathrm{N}}(\psi))}{\partial \psi_{\mathrm{N}}} \frac{\partial \psi_{\mathrm{N}}(\psi)}{\partial \psi_{\mathrm{ax}}} \tilde{\psi}(r_{\mathrm{ax}}(\psi), z_{\mathrm{ax}}(\psi)) \, \xi \, dr \, dz \\
& + \int_{\Omega_{\mathrm{p}}(\psi)} \frac{\partial j_{\mathrm{p}}(r, \psi_{\mathrm{N}}(\psi))}{\partial \psi_{\mathrm{N}}} \frac{\partial \psi_{\mathrm{N}}(\psi)}{\partial \psi_{\mathrm{bd}}} \tilde{\psi}(r_{\mathrm{bd}}(\psi), z_{\mathrm{bd}}(\psi)) \, \xi \, dr \, dz
\end{split}
\tag{26}
\end{equation}
$$

of J_p(ψ, ξ), where Γ_p is the plasma boundary ∂Ω_p and

$$
j_{\mathrm{p}}(r, \psi_{\mathrm{N}}(\psi)) = r S_{p'}(\psi_{\mathrm{N}}(\psi)) + \frac{1}{\mu_0 r} S_{ff'}(\psi_{\mathrm{N}}(\psi)). \quad (27)
$$

The derivation of the linearization D<sub>ψ</sub>J<sub>p</sub>(ψ, ξ)(<span style="text-decoration: overline">ψ</span>) requires to assume that ∇ψ ≠ 0 on ∂Ω<sub>p</sub> and involves shape calculus [12, 13] and the non-trivial derivatives:

$$
D_{\psi} \psi_{\text{ax}}(\psi)(\tilde{\psi}) = \tilde{\psi}(r_{\text{ax}}(\psi), z_{\text{ax}}(\psi)) \text{ and } D_{\psi} \psi_{\text{bd}}(\psi)(\tilde{\psi}) = \tilde{\psi}(r_{\text{bd}}(\psi), z_{\text{bd}}(\psi)).
$$

Clearly, $\nabla \psi \neq 0$ on $\partial \Omega_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(\psi, \xi)$ needs special attention due to the $\psi$-dependent domain of integration. We refer to [14, section 3.2] and [14, section 3.3] for the technical details.

5 The optimal control problem

We intend to determine the voltages $V_i(t)$ applied to the poloidal field circuits so that the plasma boundary $\Gamma_p$ fit to a desired boundary $\Gamma_{\text{desi}}$ during the whole discharge while minimizing a certain energetic cost.

Let $\Gamma_{\text{desi}}(t) \subset \Omega_L$ denote the evolution of a closed line, contained in the domain $\Omega_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 $\Omega_L$. Further let $(r_{\text{desi}}(t), z_{\text{desi}}(t)) \in \Gamma_{\text{desi}}(t)$ and $(r_1(t), z_1(t)), \dots, (r_{N_{\text{desi}}}(t), z_{N_{\text{desi}}}(t)) \in \Gamma_{\text{desi}}(t)$ be $N_{\text{desi}} + 1$ points on that line. We define a quadratic functional $K(\psi)$ that evaluates to zero if $\Gamma_{\text{desi}}(t)$ is an $\psi(t)$-isoline, i.e. if $\psi(t)$ is constant on $\Gamma_{\text{desi}}(t)$:

$$
K(\psi, t) := \frac{1}{2} \left( \sum_{i=1}^{N_{\text{desi}}} \left( \psi(r_i(t), z_i(t)) - \psi(r_{\text{desi}}(t), z_{\text{desi}}(t)) \right)^2 \right). \quad (28)
$$
Another functional, that will serve as regularization, is

$$R(\mathbf{V}(t)) := \sum_{i=1}^{N} \frac{w_i}{2} \mathbf{V}_i^2 \qquad (29)$$

with regularization weights $w_i \ge 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:

$$\min_{\psi(t), \mathbf{V}(t)} \int_{0}^{T} K(\psi(t), t) + R(\mathbf{V}(t)) dt \qquad (30)$$

subject to

$$A(\psi(t), \xi) - J_p(\psi(t), \xi) + j^{\text{ps}}(\dot{\psi}(t), \xi) + j^{\text{c}}(\dot{\psi}(t), \xi) + c(\psi(t), \xi) = \ell(\mathbf{V}(t), \xi) \quad \forall \xi \in H.$$

This minimization problem for transient axisymmetric equilibria extends the minimization problems for static axisymmetric 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 $\phi$ is:

$$\mathcal{L}(\psi(t), \mathbf{V}(t), \phi(t)) = \int_{0}^{T} K(\psi(t), t) + R(\mathbf{V}(t)) dt - \int_{0}^{T} A(\psi(t), \phi(t)) - J_{p}(\psi(t), \phi(t)) + c(\psi(t), \phi(t)) dt - \int_{0}^{T} j^{\text{ps}}(\dot{\psi}(t), \phi(t)) + j^{\text{c}}(\dot{\psi}(t), \phi(t)) - \ell(\mathbf{V}(t), \phi(t)) dt.$$

We can state the first order necessary conditions for optimality under the following three assumptions in the limiter case:

1. $\sup_{\Omega_L} \psi$ is attained at one and only one point $\mathbf{M}_0 = (r_{bd}, z_{bd})$.

2. $\sup_{\Omega_p} \psi$ is attained at one and only one point $\mathbf{M}_1$, which is an interior point of $\Omega_p$ and $\mathbf{M}_1 = (r_{ax}, z_{ax})$. $\psi$ is of class $C^2$ in a neighbourhood of $\mathbf{M}_1$ and the point $\mathbf{M}_1$ is a non-degenerated elliptic point.

3. $\nabla \psi$ vanishes nowhere on $\partial \Omega_p$.

Equivalent necessary conditions can be obtained in the X-point case.

Then necessary conditions for $(\psi(t), \mathbf{V}(t), \phi(t))$ to be a saddle point of $\mathcal{L}$ are obtained, after integrating by parts in time the Lagrangian:

- $\psi(t)$ and $\mathbf{V}(t)$ are solution of the direct problem (20)
– $\psi(t)$ and $\phi(t)$ are solution of the adjoint problem

$$D_{\psi}A(\psi(t), \phi(t))(\xi) - D_{\psi}J_{p}(\psi(t), \phi(t))(\xi) + c(\xi, \phi(t)) - j^{\text{ps}}(\xi, \dot{\phi}(t)) - j^{\text{c}}(\xi, \dot{\phi}(t)) = D_{\psi}K(\psi(t), t)(\xi) \quad \forall \xi \in H \quad (31)$$

with $\phi(T) = 0$ and

$$D_{\psi}K(\psi, t)(\xi) = \sum_{i=1}^{N_{\text{desi}}} (\psi(r_i(t), z_i(t)) - \psi(r_{\text{desi}}(t), z_{\text{desi}}(t))). \\ (\xi(r_i(t), z_i(t)) - \xi(r_{\text{desi}}(t), z_{\text{desi}}(t))).$$

– $V(t)$ and $\phi(t)$ are solution to

$$w_i V_i(t) + \frac{n_i}{R_i |\Omega_{c_i}|} \int_{\Omega_{c_i}} \phi(t) drdz = 0, \quad 1 \le i \le N. \quad (32)$$

The adjoint problem has the following strong formulation:

$$-\Delta^* \phi(t) + 1_{\Omega_{\text{Fe}}} \nabla \cdot \left( 2 \frac{\mu'_{\text{Fe}} (\frac{|\nabla \psi|^2}{r^2})}{\mu^2_{\text{Fe}} (\frac{|\nabla \psi|^2}{r^2})} (\nabla \phi(t) \cdot \nabla \psi) \nabla \psi \right) \\ - 1_{\Omega_{\text{p}}(\psi)} \frac{\partial j_{\text{p}}(r, \psi_{\text{N}}(\psi))}{\partial \psi_{\text{N}}} \frac{\partial \psi_{\text{N}}(\psi)}{\partial \psi} \phi(t) \\ - \delta_{\text{bd}} \int_{\Gamma_{\text{p}}(\psi)} \frac{j_{\text{p}}(r, 1)}{|\nabla \psi|} \phi(t) d\Gamma + (\delta_{\Gamma_{\text{p}}}, \frac{j_{\text{p}}(r, 1)}{|\nabla \psi|} \phi(t)) \\ - \delta_{\text{ax}} \int_{\Omega_{\text{p}}(\psi)} \frac{\partial j_{\text{p}}(r, \psi_{\text{N}}(\psi))}{\partial \psi_{\text{N}}} \frac{\partial \psi_{\text{N}}(\psi)}{\partial \psi_{\text{ax}}} \phi(t) drdz \\ - \delta_{\text{bd}} \int_{\Omega_{\text{p}}(\psi)} \frac{\partial j_{\text{p}}(r, \psi_{\text{N}}(\psi))}{\partial \psi_{\text{N}}} \frac{\partial \psi_{\text{N}}(\psi)}{\partial \psi_{\text{bd}}} \phi(t) drdz \\ - 1_{\Omega_{\text{ps}}} \frac{\sigma}{r} \dot{\phi}(t) - \sum_{i=1}^{N_i} 1_{\Omega_{c_i}} \frac{2\pi n_i^2}{R_i |\Omega_{c_i}|^2} \int_{\Omega_{c_i}} \dot{\phi} drdz \\ = \left( \sum_{i=1}^{N_{\text{desi}}} (\psi(r_i(t), z_i(t), t) - \psi(r_{\text{desi}}(t), z_{\text{desi}}(t), t)) \right) (\delta_{(r_i, z_i)} - \delta_{(r_{\text{desi}}, z_{\text{desi}})})$$

with $\phi(T) = 0$, where $\delta_{\text{ax}}$ and $\delta_{\text{bd}}$ are the Dirac masses at the points $(r_{\text{ax}}, r_{\text{ax}})$ and $(r_{\text{bd}}, r_{\text{bd}})$, respectively. $\delta_{\Gamma_{\text{p}}}$ is the Dirac mass of $\Gamma_{\text{p}}$ with

$$(\delta_{\Gamma_{\text{p}}}, \frac{j_{\text{p}}(r, 1)}{|\nabla \psi|} \phi(t)\xi) = \int_{\Gamma_{\text{p}}} \frac{j_{\text{p}}(r, 1)}{|\nabla \psi|} \phi(t)\xi d\Gamma.$$

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_{\mathbf{u},\mathbf{y}} J(\mathbf{y}, \mathbf{u}) \quad \text{s.t.} \quad \mathbf{B}(\mathbf{y}) = \mathbf{F}(\mathbf{u}), \qquad (33) $$

where $\mathbf{y}$ and $\mathbf{u}$ are the so-called state and control variables. In our setting $\mathbf{y}$ will be the variable that describes the plasma and $\mathbf{u}$ will be the externally applied voltages. We think of $\mathbf{y}$ as the vector of degrees of freedoms describing the space and time evolution of the poloidal flux $\psi$, and $\mathbf{B}(\mathbf{y})$ and $\mathbf{F}(\mathbf{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 non-linear 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

$$ L(\mathbf{y}, \mathbf{u}, \mathbf{p}) = J(\mathbf{y}, \mathbf{u}) + \langle \mathbf{p}, \mathbf{B}(\mathbf{y}) - \mathbf{F}(\mathbf{u}) \rangle, \qquad (34) $$

and the solution of (33) is a stationary point of this Lagrangian:

$$ \begin{aligned} & D_{\mathbf{y}}J(\mathbf{y}, \mathbf{u}) + D_{\mathbf{y}}\mathbf{B}^T(\mathbf{y})\mathbf{p} = 0, \\ & D_{\mathbf{u}}J(\mathbf{y}, \mathbf{u}) - D_{\mathbf{u}}\mathbf{F}^T(\mathbf{u})\mathbf{p} = 0, \\ & \mathbf{B}(\mathbf{y}) - \mathbf{F}(\mathbf{u}) = 0. \end{aligned} \qquad (35) $$

A Newton-type method for solving (35) are iterations of the type

$$ \left( \begin{array}{ccc} \mathbf{H}_{\mathbf{y},\mathbf{y}}^k & \mathbf{H}_{\mathbf{y},\mathbf{u}}^k & D_{\mathbf{y}}\mathbf{B}^T(\mathbf{y}^k) \\ \mathbf{H}_{\mathbf{u},\mathbf{y}}^k & \mathbf{H}_{\mathbf{u},\mathbf{u}}^k & -D_{\mathbf{u}}\mathbf{F}^T(\mathbf{u}^k) \\ D_{\mathbf{y}}\mathbf{B}(\mathbf{y}^k) & -D_{\mathbf{u}}\mathbf{F}(\mathbf{u}^k) & 0 \end{array} \right) \left( \begin{array}{c} \mathbf{y}^{k+1} - \mathbf{y}^k \\ \mathbf{u}^{k+1} - \mathbf{u}^k \\ \mathbf{p}^{k+1} - \mathbf{p}^k \end{array} \right) = - \left( \begin{array}{c} D_{\mathbf{y}}J(\mathbf{y}^k, \mathbf{u}^k) + D_{\mathbf{y}}\mathbf{B}^T(\mathbf{y}^k)\mathbf{p}^k \\ D_{\mathbf{u}}J(\mathbf{y}^k, \mathbf{u}^k) - D_{\mathbf{u}}\mathbf{F}^T(\mathbf{u}^k)\mathbf{p}^k \\ \mathbf{B}(\mathbf{y}^k) - \mathbf{F}(\mathbf{u}^k) \end{array} \right) \qquad (36) $$

with

$$ \left( \begin{array}{cc} \mathbf{H}_{\mathbf{y},\mathbf{y}}^k & \mathbf{H}_{\mathbf{y},\mathbf{u}}^k \\ \mathbf{H}_{\mathbf{u},\mathbf{y}}^k & \mathbf{H}_{\mathbf{u},\mathbf{u}}^k \end{array} \right) = \left( \begin{array}{cc} D_{\mathbf{y},\mathbf{y}}L(\mathbf{y}^k, \mathbf{u}^k, \mathbf{p}^k) & D_{\mathbf{y},\mathbf{u}}L(\mathbf{y}^k, \mathbf{u}^k, \mathbf{p}^k) \\ D_{\mathbf{u},\mathbf{y}}L(\mathbf{y}^k, \mathbf{u}^k, \mathbf{p}^k) & D_{\mathbf{u},\mathbf{u}}L(\mathbf{y}^k, \mathbf{u}^k, \mathbf{p}^k) \end{array} \right). $$

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 $\Delta\mathbf{u}^k := \mathbf{u}^{k+1} - \mathbf{u}^k$:

$$ \mathbf{M}(\mathbf{y}^k, \mathbf{u}^k)\Delta\mathbf{u}^k = -\mathbf{h}(\mathbf{y}^k, \mathbf{u}^k), \qquad (37) $$

where

$$ \mathbf{M}(\mathbf{y}^k, \mathbf{u}^k) := (D_{\mathbf{u}}\mathbf{F}^T(\mathbf{u}^k)D_{\mathbf{y}}\mathbf{B}^{-T}(\mathbf{y}^k)\mathbf{I}_d) \left( \begin{array}{cc} \mathbf{H}_{\mathbf{y},\mathbf{y}}^k & \mathbf{H}_{\mathbf{y},\mathbf{u}}^k \\ \mathbf{H}_{\mathbf{u},\mathbf{y}}^k & \mathbf{H}_{\mathbf{u},\mathbf{u}}^k \end{array} \right) \left( D_{\mathbf{y}}\mathbf{B}^{-1}(\mathbf{y}^k)D_{\mathbf{u}}\mathbf{F}(\mathbf{u}^k) \mathbf{I}_d \right) $$
and

$$ \mathbf{h}(\mathbf{y}^k, \mathbf{u}^k) := D_{\mathbf{u}} J(\mathbf{y}^k, \mathbf{u}^k) + D_{\mathbf{u}} \mathbf{F}^T(\mathbf{u}^k) \lambda^k \\ - (D_{\mathbf{u}} \mathbf{F}^T(\mathbf{u}^k) D_{\mathbf{y}} \mathbf{B}^{-T}(\mathbf{y}^k) \mathbf{H}_{\mathbf{y},\mathbf{y}}^k + \mathbf{H}_{\mathbf{u},\mathbf{y}}^k) D_{\mathbf{y}} \mathbf{B}^{-1}(\mathbf{y}^k) \mathbf{r}(\mathbf{y}^k, \mathbf{u}^k)) $$

with

$$ \lambda^k := D_{\mathbf{y}} \mathbf{B}^{-T}(\mathbf{y}^k) D_{\mathbf{y}} J(\mathbf{y}^k, \mathbf{u}^k), \quad \mathbf{r}(\mathbf{y}^k, \mathbf{u}^k) := \mathbf{B}(\mathbf{y}^k) - \mathbf{F}(\mathbf{u}^k). $$

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_{\mathbf{y}}\mathbf{B}(\mathbf{y}^k)$ and $D_{\mathbf{y}}\mathbf{B}^T(\mathbf{y}^k)$. Alternatively, if we have sufficient memory to store $\mathbf{M}(\cdot, \cdot)$, we can compute $\mathbf{M}(\cdot, \cdot)$ explicitly. Clearly, we never compute neither $D_{\mathbf{y}}\mathbf{B}^{-1}(\mathbf{y}^k)$ nor $D_{\mathbf{y}}\mathbf{B}^{-T}(\mathbf{y}^k)$ explicitly.

Once we know $\Delta\mathbf{u}^k$ we can compute $\mathbf{y}^{k+1}$ and $\mathbf{p}^{k+1}$ by:

$$ \mathbf{y}^{k+1} - \mathbf{y}^k = D_{\mathbf{y}} \mathbf{B}^{-1}(\mathbf{y}^k) D_{\mathbf{u}} \mathbf{F}(\mathbf{u}^k) \Delta \mathbf{u}^k - \mathbf{r}(\mathbf{u}^k, \mathbf{y}^k), \\ \mathbf{p}^{k+1} + \lambda^k = -D_{\mathbf{y}} \mathbf{B}^{-T}(\mathbf{y}^k) (\mathbf{H}_{\mathbf{y},\mathbf{y}}^k (\mathbf{y}^{k+1} - \mathbf{y}^k) + \mathbf{H}_{\mathbf{y},\mathbf{u}}^k (\mathbf{u}^{k+1} - \mathbf{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 $\hat{\mathcal{J}}(\mathbf{u}) := J(\mathbf{y}(\mathbf{u}), \mathbf{u})$, with $\mathbf{B}(\mathbf{y}(\mathbf{u})) = \mathbf{F}(\mathbf{u})$ be the reduced cost function, then we have the following expressions for gradient

$$ D_{\mathbf{u}} \hat{\mathcal{J}}(\mathbf{u}) = D_{\mathbf{u}} J(\mathbf{y}, \mathbf{u}) + D_{\mathbf{u}} \mathbf{F}^T(\mathbf{u}), \lambda $$

and Hessian

$$ D_{\mathbf{u},\mathbf{u}} \hat{\mathcal{J}}(\mathbf{u}) = \mathbf{Z}^T \begin{pmatrix} D_{\mathbf{y},\mathbf{y}} J(\mathbf{y}, \mathbf{u}) & D_{\mathbf{y},\mathbf{u}} J(\mathbf{y}, \mathbf{u}) \\ D_{\mathbf{u},\mathbf{y}} J(\mathbf{y}, \mathbf{u}) & D_{\mathbf{u},\mathbf{u}} J(\mathbf{y}, \mathbf{u}) \end{pmatrix} \mathbf{Z} \\ + \mathbf{Z}^T \begin{pmatrix} -D_{\mathbf{y}}(D_{\mathbf{y}}\mathbf{B}^T(\mathbf{y})\lambda) & 0 \\ 0 & D_{\mathbf{u}}(D_{\mathbf{u}}\mathbf{F}^T(\mathbf{u})\lambda) \end{pmatrix} \mathbf{Z} $$

with

$$ \lambda = D_{\mathbf{y}} \mathbf{B}^{-T}(\mathbf{y}) D_{\mathbf{y}} J(\mathbf{y}, \mathbf{u}) \quad \text{and} \quad \mathbf{Z} = \begin{pmatrix} D_{\mathbf{y}} \mathbf{B}^{-1}(\mathbf{y}) D_{\mathbf{u}} \mathbf{F}(\mathbf{u}) \\ \mathrm{Id} \end{pmatrix} $$

Hence, the reduced gradient $\mathbf{h}(\mathbf{y}^k, \mathbf{u}^k)$ is not the gradient of the reduced cost function, unless the state and control variable $\mathbf{y}^k$ and $\mathbf{u}^k$ verify the equation of state $\mathbf{B}(\mathbf{y}^k) = \mathbf{F}(\mathbf{u}^k)$.
**Preliminary Example** Finally, we would like to show first results for a so-
called 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 im-
portant 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.

<figure><img src="image_6.png" /><figcaption>Fig. 5. The optimal voltages.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 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).</figcaption></figure>
References

1. J.L. Lions. *Optimal control of systems governed by partial differential equations*. Springer, Berlin Heildelberg New York, 1971.

2. S.I. Brangiskii. Consultant Bureau, New York. *Reviews of Plasma Physics*, 1:205–311, 1965.

3. E.K. Maschke and J. Pantuso Sudano. Etude analytique de l'évolution d'un plasma toroidal de type tokamak à section non circulaire (analytic study of the evolution of a toroidal plasma of tokamak type with non-circular section). Technical report, EUR-CEA-FC-668, 1972.

4. D.B. Nelson and H. Grad. Heating and transport in tokamaks of arbitrary shape and beta. Technical report, Oak Ridge Report ORNL/TM-6094, 1978.

5. J. Wesson. *Tokamaks*. Oxford science publications, 2004.

6. R. Albanese, J. Blum, and O. Barbieri. On the solution of the magnetic flux equation in an infinite domain. In *EPS. 8th Europhysics Conference on Computing in Plasma Physics (1986)*, 1986.

7. Virginie Grandgirard. *Modélisation de l'équilibre d'un plasma de tokamak*. PhD thesis, l'Université de Franche-Comté, 1999.

8. G.N. Gatica and G.C. Hsiao. The uncoupling of boundary integral and finite element methods for nonlinear boundary value problems. *J. Math. Anal. Appl.*, 189(2):442–461, 1995.

9. J.L. Luxon and B.B. Brown. Magnetic analysis of non-circular cross-section tokamaks. *Nuclear Fusion*, 22(6):813, 1982.

10. J. Blum, J. Le Foll, and B. Thooris. The self-consistent equilibrium and diffusion code SCED. *Computer Physics Communications*, 24:235–254, 1981.

11. R. Albanese, J. Blum, and O. De Barbieri. Numerical studies of the Next European Torus via the PROTEUS code. In *12th Conf. on Numerical Simulation of Plasmas*, *San Francisco*, 1987.

12. F. Murat and J. Simon. Sur le contrôle par un domaine géométrique. Technical Report 76015, Laboratoire d'Analyse Numérique, Université de Paris 6, 1976.

13. M.C. Delfour and J.-P. Zolésio. *Shapes and geometries*, volume 22 of *Advances in Design and Control*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition, 2011.

14. H. Heumann, J. Blum, C. Boulbe, B. Faugeras, G. Selig, J.-M. Ané, S. Brémond, V. Grandgirard, P. Hertout, and E. Nardon. Quasi-static free-boundary equilibrium of toroidal plasma with cedres++: Computational methods and applications. *Journal of Plasma Physics*, FirstView:1–35, 1 2015.

15. J. Blum. *Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks*. Series in Modern Applied Mathematics. Wiley Gauthier-Villars, Paris, 1989.

16. J. Nocedal and S. J. Wright. *Numerical optimization*. Springer Series in Operations Research and Financial Engineering. Springer, New York, second edition, 2006.