A constrained control strategy for the shape control in thermonuclear fusion tokamaks

Automatica 49 (2013) 169–177

Contents lists available at SciVerse ScienceDirect

Automatica journal homepage: www.elsevier.com/locate/automatica

Brief paper

A constrained control strategy for the shape control in thermonuclear fusion tokamaks✩ Massimiliano Mattei a,1 , Carmelo Vincenzo Labate b , Domenico Famularo c a

DIAM, Seconda Università degli Studi di Napoli and CREATE consortium, Italy

b

DIMET, Università degli Studi ‘‘Mediterranea’’ di Reggio Calabria and CREATE consortium, Italy

c

DEIS, Università degli Studi della Calabria, Italy

article

info

Article history: Received 3 February 2011 Received in revised form 31 May 2012 Accepted 18 July 2012 Available online 1 October 2012 Keywords: Model based control Reference Governor Predictive control Nuclear fusion Tokamaks

abstract The paper deals with the application of the so-called Reference (or Command) Governor constrained control strategy to the shape control of plasmas in thermonuclear fusion reactors with the main scope of optimizing tokamak operations also in conditions very close to the operating envelope limits. A primal inner loop controlling the plasma–wall distance is first designed; the Reference Governor device is then tuned to modify, whenever necessary, the reference signals to the inner loop, on the basis of constraints due to voltage saturations on the power supply converters, limitations of currents in the active control coils, minimum clearance between the plasma surface and the vacuum chamber wall, maximum induced magnetic fields and forces on coils. As usual in model predictive paradigms, the reference signal modification is accomplished through an on-line optimization procedure which embodies plasma model forecasts computed along a finite time virtual receding horizon. The ITER (International Thermonuclear Experimental Reactor) tokamak is assumed as the case study. Numerical simulations are carried out on a finite elements nonlinear model taking into account induced currents in the passive structures. The proposed application shows how almost a hundred constraints can be managed on-line by the Reference Governor. © 2012 Elsevier Ltd. All rights reserved.

1. Introduction In tokamak devices, Wesson (1997), plasma control has gained relevance due to increasing performance demand. To this extent modern tokamaks like ITER are designed to obtain plasmas which are significantly elongated and vertically unstable, placed as close as possible to the metallic facing components. This ensures a good passive stabilization, due to the eddy currents induced in these metallic structures, and an efficient use of the available volume. Contact between plasma and chamber walls is always a major concern during operations. The desired plasma–wall clearance is obtained by regulating currents in a number of PF (Poloidal Field) coils surrounding the plasma ring. Fig. 1 shows a poloidal section of ITER: this is a 500 MW fusion power reactor under construction in Cadarache (France).

✩ The material in this paper was partially presented at the European Control Conference 2009 (ECC’09), August 23–26, 2009, Budapest, Hungary. This paper was recommended for publication in revised form by Associate Editor Huaguang Zhang under the direction of Editor Toshiharu Sugie. E-mail addresses: [email protected] (M. Mattei), [email protected] (C.V. Labate), [email protected] (D. Famularo). 1 Tel.: +39 081 5010406; fax: +39 0817683171.

0005-1098/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2012.09.004

Voltages applied to the 12 PF active coils are generated by power supplies driven by a feedback control system. The feedback action can be divided in two parts: position control guarantees vertical stability of plasma that would be naturally unstable if elongated, whereas shape control guarantees the control of plasma geometrical parameters such as elongation, triangularity, and plasma–wall distance. Usually both plasma position and shape controllers have quite a simple structure and are mainly based on multi-loop Proportional Integral Derivative (PID) actions. In order to maximize the tokamak performance, the distance between the plasma boundary and the vessel at some specific points (gaps) is usually controlled together with the vertical speed to ensure stability. Multiple-input–multiple-output control approaches have been deeply investigated in the last decade to improve shape control performances, Ariola et al. (1999), Crisanti et al. (2003), and Humphreys et al. (2005). The most successful are linear model based approaches which however do not take into account the possibility that some physical variables of interest can exceed their operating limits. In Ambrosino, Ariola, Pironti, and Walker (2001), Varano et al. (2010), and Vitelli et al. (2010), some attempts to deal with control inputs saturation are described. However, during tokamak operations, not only voltages and currents, but also induced

170

M. Mattei et al. / Automatica 49 (2013) 169–177

2. Mathematical modeling of the plasma response Three main subsystems need to be considered to obtain a control oriented model of plasma shape and position evolution in a tokamak: the plasma, the control circuits, and the passive conductors. The physical phenomenon to be controlled is governed by Maxwell’s equations in their quasi-stationary form where the electric field does not vary too rapidly and the current density is divergence free with the constitutive relationships. In axialsymmetric geometry, with cylindrical coordinates (r , ϕ, z), the magnetic field and current density vectors B and J can be expressed in terms of two scalar functions, namely, the poloidal magnetic flux per radians ψ and the poloidal current function f = rBϕ , Bϕ being the toroidal magnetic field. At the time scale of interest for current, position, and shape control, because of the low plasma mass density, inertial effects can be neglected. Hence the plasma momentum balance becomes J × B = ∇ p at equilibrium, that can be rewritten as the well known Grad–Shafranov equation, Wesson (1997): Fig. 1. ITER poloidal section.

∂ ∆ ψ =r ∂r ∗

electromagnetic fields and forces, and shape parameters must stay within prescribed ranges. An important motivation for the present work is that active tokamaks have been designed with large operational margins. This does not happen in ITER, and will not happen in future large machines where the impact of margins on costs grows very rapidly. The need to satisfy input and/or state dependent constraints is in general a relevant problem in control theory and practice. Anti-Windup (AW), Bumpless methods, AW/LQR, and AW/H2, are feedback control methodologies dealing with the presence of input constraints in an indirect manner. Recently, techniques based on invariant sets arguments and predictive control ideas, Kothare, Balakrishnan, and Morari (1996) and Mayne, Rawlings, Rao, and Scokaert (2000), have gained in popularity due to their inherent capability to take the presence of constraints directly into account during the design phase. The control action is computed through the solution of a sequence of optimization problems based on the prediction of the plant state evolution. The objective is to jointly maximize the control performance and enforce the satisfaction of prescribed constraints. The interest towards a practical application of such methodologies has been growing in the last decade due to the availability of fast computing units, Diel, Bock, and Schloder (2005). Some of the predictive control methodologies are mainly devoted to constraints fulfillment leaving control performance satisfaction to traditional regulation frameworks. Such a family of control strategies is known in the literature as the Reference Governor (RG) approach. The RG is a nonlinear device which is added to a pre-compensated plant, designed so as to exhibit stability and good tracking performance in the RG’s absence. At each time instant tk , it computes a modified reference command which, if applied from tk onward, does not produce constraint violations. Such a modified reference command is computed to minimize its distance from the actual desired reference signal, according to an on-line constrained procedure on a receding horizon finite time interval. Many mature assessments of the RG state of the art for linear systems can be found in Bemporad (1998), Bemporad and Mosca (1995), Borrelli, Falcone, Pekar, and Stewart (2009), Casavola, Mosca, and Angeli (2000), Gilbert and Kolmanovsky (2002), Gilbert, Kolmanovsky, and Tin Tan (1995), Kolmanovsky and Sun (2006), and Vahidi, Kolmanovsky, and Stefanopolou (2007).

= −f



1 ∂ψ



µr r ∂ r df

dψ

− µ0 r 2

∂ + ∂z

dp dψ



1 ∂ψ



µr ∂ z

.

(1)

µr = µ/µ0 being the relative magnetic permeability, and µ0 the permeability of free space. To complete modeling with the interaction between plasma and surrounding passive structures and active coils the following partial differential equation problem can be written taking into account the axial-symmetry

 dp df  ∆∗ ψ = −f − µ0 r 2 dψ dψ ∗  ∆∗ ψ = −µ0 rjext ∆ ψ =0

in plasma region in conductors elsewhere

(2)

with the initial and boundary conditions



ψ(r , z , 0) = ψ0 (r , z ), ψ(0, z , t ) = 0, limr 2 +z 2 →∞ ψ(r , z , t ) = 0.

(3)

jext being the toroidal current density in the external conductors and coils. Both jext and ψ depend on space (r and z in axialsymmetric conditions), and time. The above equations are used to calculate the poloidal flux function at time t provided that the plasma boundary can be determined, the toroidal current density in the PF is known and the functions p(ψ) and f (ψ) are defined. The plasma boundary is determined by means of an iterative numerical procedure. Functions p(ψ) and f (ψ) can be expressed in terms of global plasma parameters, for example Ipl , βpol and li , namely plasma current, poloidal beta, and internal inductance. As for the toroidal current density jext , it can be expressed as a linear combination of the PF circuit currents. The dependence of p(ψ) and f (ψ) functions on Ipl , βpol and li is assigned implicitly with a suitable parameterization of such functions. In practice this is obtained by adding constraints to the numerical solution of problem (2) and (3). Therefore, the magnetic flux and the plasma configuration can be determined when prescribing the vector of currents, including poloidal field currents and plasma current along with βp and li . The time evolution of currents is governed by circuit equations driven by voltages in the active coils. Induced eddy currents in the metallic

M. Mattei et al. / Automatica 49 (2013) 169–177

171

structures can be modeled as currents flowing in circuits driven by voltages generated by electromagnetic induction. Solutions of problem (2) and (3) can be numerically found by means of Finite Element Methods (FEMs), Albanese, Calabró, Mattei, and Villone (2003). A control design oriented procedure to linearize the plasma model response is described in Albanese and Villone (1998). This procedure finally approaches to a linear dynamic plasma-circuit model in the form: Lδ ˙I + Rδ I = δ V + LE δ w ˙ δ y = C δ I + F δw



(4)

where V is the vector of voltages applied to the circuits (zero for passive coils), R is the circuit resistance matrix, L is the matrix of self- and mutual inductances between the plasma, the coils, and T T T the equivalent circuits of the passive structures, I = [Iacti v e Ieddy Ipl ] is the vector of active currents (a linear combination of PF currents), passive (eddy) and plasma currents, and w = [βpol li ]T is a vector of plasma parameters typically assumed as external disturbances; y is a vector of outputs, and δ denotes a variation with respect to a nominal equilibrium condition. When passive metallic structures are discretized with a FEM approach, a high number of corresponding passive circuits and currents has to be accounted for in the model. A model reduction can be achieved by means of a modal decomposition. The unique unstable mode and those modes whose time constants are not much faster than the closed loop characteristic times are usually considered. The controllers discussed in the paper have been designed on reduced order linear models although they have been tested on high dimensional nonlinear models.

Fig. 2. Definition of the controlled gaps, at left, and of the maximum plasma displacement (black markers), at right.

3. The control problem and the primal-inner controller The proposed plasma controller has three main objectives: vertical stabilization, plasma current and shape control. 11 voltages can be independently driven in ITER to make shaping: P1÷P6, CS3U, CS3L, CS2U, CS2L coil currents are controlled separately; CS1L and CS1U coils are connected in series (see Fig. 1). VS (Vertical Stabilization) can be achieved controlling the vertical speed of the plasma centroid through the production of a suitable radial field. In ITER this control action, obtained by driving the (P2-P3)–(P4-P5) coils connected in antiseries with an additional VS power supply, counteracts the vertical instability in the time scale of hundreds of ms. In the proposed application, a SISO (Single Input–Single Output) proportional control is designed to guarantee vertical stabilization assuming the time derivative of the vertical centroid position zpl as a controlled and measured variable. The typical crossover frequency for the VS loop is about 15 rad/s. Plasma current control can be obtained by means of the socalled transformer pattern of active coil currents, i.e. a combination of currents chosen to induce plasma current without significant changes in plasma shape. A purely proportional action is adopted to control Ipl with a response time in the order of 10 s. Shape control is implemented in the form of gap control, i.e. the control of a predefined number of reference plasma–wall distances called gaps. The controlled gaps adopted in the paper, depicted in Fig. 2, are the plasma top gap (#21), two lateral gaps on the equatorial plane (#9 and #30), and the position of the separatrix legs with respect to the divertor plates (LS and RS gaps). Such a choice is quite common in shape control strategies since their control allows changing global plasma shape parameters such as triangularity, volume, and elongation. Gap control is achieved by means of a multivariable PI (Proportional-Integral) action, whose gains are designed with a

Fig. 3. Representation of the pre-compensated system. Primal-inner shape control is designed starting from a plant model including vertical stabilization and plasma current control.

linear quadratic optimal control approach. Such a control action is tuned to operate in the same time scale of the plasma current loop. It is assumed that the reduced order system state is reconstructed via Kalman filtering (Fig. 3). A detailed description of the primal-inner controller structure shown in Fig. 2, which is quite standard in control engineering practice, is outside the scope of the present work. The closed loop system behavior in the presence of a multi-step reference signals on gaps is shown in Fig. 4–5. 4. The reference governor In its most common formulation, Bemporad and Mosca (1995) and Gilbert et al. (1995), and according to the scheme depicted in Fig. 6, the RG approach assumes a discrete time, pre-compensated linear time-invariant model obtained by means of a system linearization around a nominal equilibrium condition: x(tk+1 ) = Φ x(tk ) + Gg (tk ) + Gd d(tk ) y(tk ) = Hy x(tk ) c (tk ) = Hc x(tk ) + Lg (tk ) + Ld d(tk )



(5)

where tk = t0 + kTs , k ∈ Z0+ , t0 and Ts being the initial time instant and the sampling interval respectively; x(tk ) ∈ Rnx is the state vector including plant and primal controller states; g (tk ) ∈ Rnr is the modified reference signal, coinciding with the

172

M. Mattei et al. / Automatica 49 (2013) 169–177

Fig. 4. Inner-primal controller performance: vertical position and speed, controlled gaps compared with reference signals.

Fig. 5. Primal-inner controller closed loop performance: states (δIactiv e and δIeddy ) and command inputs (Vcirc ).

original reference signal r (tk ) ∈ Rnr in the absence of constraints; d(tk ) ∈ Rnd is the disturbance vector (βpol and li in our working belonging example)  to a closed convex and compact set DD = d ∈ Rnd : Ud ≤ h¯ with U ∈ Rnu ×nd , nu ≥ nd a full column rank matrix, and h¯ ∈ Rnu a vector of non-negative constraints (h¯ p ≥ 0, p = 1, . . . , nu ); y(tk ) ∈ Rnr is the controlled output vector (gaps), required to track r (tk ); c (tk ) ∈ C ⊂ Rnc is a vector containing linear combinations of plant states and inputs to be constrained (currents, voltages, fields, forces, gaps), C = {c ∈ Rnc : Tc ≤ f }, with T ∈ Rnt ×nc , nt ≥ nc , a full column rank matrix, and f ∈ Rnt a vector of constraints.

Fig. 6. Reference Governor control scheme. xp and xc represent the plant and the controller states respectively.

M. Mattei et al. / Automatica 49 (2013) 169–177

Under the assumptions that system (5) is asymptotically stable and offset-free (i.e. Hy (Inx − Φ )−1 G = Inr ) the RG design problem can be stated as follows.

173

The obtained constant reference signal is applied to the plant in the [tk , tk+1 [ time interval, and the procedure is repeated at the next time tk+1 on the basis of the new measured state x(tk+1 ). Consequently, if we denote with

RG problem. Find, at each discrete time instant tk , a memoryless command function of the current state and reference signal g (tk ) = g (x(tk ), r (tk )) which is the best approximation of r (tk ), and compatible with the prescribed constraints under all feasible disturbance sequences. In order to achieve a numerical procedure, the role of external disturbances and references w.r.t. the prescribed constraint set must be clarified. The disturbance effect is taken into account via P-difference argument. Starting from the constraints set C the following recursion is obtained:

the set of all constant references ω ∈ W , whose corresponding c-evolutions starting from an initial condition x(tk ), at time tk , satisfy the constraints also during transients, and provided that V (x(tk )) is nonempty, closed and convex for all tk , the RG command is the solution of the following constrained optimization problem:

C0 := C ∼ Ld D ,

where

Ch := Ch−1 ∼ Hc Φ C∞ :=

∞ 

h −1

Gd D ,

h > 0,

(6)

Ci

i =0

where

A ∼ B = a ∈ R n : a + b ∈ A , ∀b ∈ B .





(7)

Sets Ch turn out to be suitable convex restrictions of the original constraints convex set C . Since C is compact and 0nd ∈ D , it is possible to show that, if nonempty, all Ch ’s are compact and satisfy Ch ⊂ Ch−1 . In this case the set C∞ results in a non-empty convex and compact set. When considering the set valued evolution of the state due to the disturbance action xd (th ) it can be proved that



{xd (th )} ∈ ∆h

(8)

{d(th )}hl=0 ∈D

where ∆0 = {0nd }, ∆h = ∆h−1 + Φ h−1 Gd D . Moreover, each ∆h is compact and, since Φ is a stability matrix, the sequence ∆h has an Hausdorff compact limit ∆∞ (see Gilbert and Kolmanovsky (2002)). The set ∆∞ is the smallest closed set which contains the state trajectory xd (th ) generated by all possible disturbance sequences {d(th )} ∈ D . The consequence is that the following implication holds cf (th , x, ω) ∈ Ch , ∀h ⇒ c (th , x, ω, D ) ∈ C

(9)

where cf denotes the disturbance-free constraints evolutions under the action of a constant input ω. The admissible references are characterized by defining the convex and closed set W (assumed nonempty) which characterizes all constant inputs ω ∈ Rnr whose corresponding disturbance free steady-state solutions of Eq. (5), c¯ω = Hc (Inx − Φ )−1 Gω + Lω

(10)

c¯ (tk+h , x(tk ), ω) = Hc

Φ x(tk ) + h





g (tk ) = arg min J (r (tk ), ω)

(12)

(13)

ω∈V (x(tk ))

J (r (tk ), ω) = ∥ω − r (tk )∥2Σ .

(14)

∥x∥2Σ := xT Σ x being a weighted norm with Σ a positive definite symmetric matrix. The minimizer of (13) represents the best approximation of the reference signal r (tk ) which, if constantly applied from tk onwards to (5), would never produce constraints violation. The optimization problem (13) is a QP (Quadratic Programming) problem with an infinite number of constraints where the feasible set is non-void (the equilibrium state/modified command couple (xω , g ) = (0nx , 0nr ) is always a solution of (13)). The implementation of the RG strategy requires then a computable way to overcome the inherent difficulties related to the cardinality of V (x(tk )). Nonetheless, V (x(tk )) can be finitely determined. There exists, in fact, an integer k∗ such that if c¯ (tk+h , x(tk ), ω) ∈ Ch , ∀h ∈ {0, 1, . . . , k∗ } , then c¯ (tk+h , x(tk ), ω) ∈ Ch , ∀h ∈ Z0+ . The horizon length value k∗ can be obtained according to an algorithm which is based on the solution of the following optimization problem: Gk (j) =

max

x∈Rnx ,ω∈W ξ

Tj c¯ (tk , x, ω) − fjk

subject to Tj c¯ (ti , x, ω) ≤ fji ,

i = 0, . . . , k − 1

(15)

where Tj , j = 1, . . . , nt , denotes the j-th row of the matrix (see the definition of vector c) and fji , i = 0, . . . , k − 1, have the following expression:

 0 fj = fj − sup Tj Ld d    d∈D   1 fj = fj0 − sup Tj Hc Gd d d∈D    k−1 k−2  fj = fj − sup Tj Hc Φ k−1−j Gd d.

(16)

d∈D

satisfy the given constraints with a prescribed tolerance level ξ > 0. The RG strategy consists in choosing, at each time step tk , a constant virtual reference ω ∈ W such that the corresponding disturbance-free constrained output evolution starting from the initial measured state x(tk )



V (x(tk )) := ω ∈ W ξ : c¯ (tk+h , x(tk ), ω) ∈ C , ∀h > 0

h −1 

Step 1.; k = 1;

 Φ

h−i−1

Gω

Notice that in Eq. (16) the presence of possible disturbances affecting the constrained output is taken into account via Pdifference arguments. The following algorithm is adopted to derive the constraint horizon length:

+ Lω

(11)

i=0

fulfills the restricted sets Ch , ∀h > k, which take into consideration the disturbance effects, and its distance from the reference r (tk ) is minimal.

Step 2. Find Gk (j), j = 1, . . . , nt , solving Problem (14); Step 3. If Gk (j) ≤ 0, j = 1, . . . , nt ; then, set k∗ = k, and stop; else, k = k + 1, go to Step 2; end. The above algorithm is executed off-line and, once k∗ is computed, the optimization problem (15) is converted into a QP

174

M. Mattei et al. / Automatica 49 (2013) 169–177

problem with a finite number of linear constraints to be solved on-line, that is:

If we also consider the tolerance margin on the constraints fulfillment ξ , we have the following approximation of C∞ :

g (tk ) := arg min J (r (tk ), ω)

aξ C∞ ≈ C∞ (ε) = Ckε ∼ Bε ∼ Bξ



ω∈W ξ

s.t. THc Φ h x(tk ) + T

h−1  i =0



h = 0, . . . , k .

(17)

It can be proved that, Bemporad and Mosca (1995), if the model plant (5) is asymptotically stable, it satisfies the offset-free assumption, and V (x(tk )) is nonempty, then: the minimizer in (17) uniquely exists at each time instant; V (x(tk + h)) is nonempty ∀h > k along the trajectories generated by the RG command (viability property); the constrained output c (tk ) always fulfills constraints ∀ tk and the overall closed loop system is asymptotically stable; whenever r (tk ) ≡ r, g (tk ) monotonically converges in finite time to either r, or its best admissible approximation in the sense of the optimization problem rˆ := arg minω∈W ξ J (r , ω). Then, by the offset-free condition limtk →+∞ y¯ (tk ) = rˆ , y¯ (tk ) being the disturbance free component of the plant output. As for the computational aspects of the proposed technique, it is worth pointing out that the off-line problem (15) is a LP (Linear Programming) problem, which can be solved using standard simplex algorithms, whereas QP problem (17), which has to be solved on-line, can be solved through many solvers embedded in the main commercial optimization and control software packages. As it will be shown in the experiment section QP solvers provide a solution fast enough to be implemented on-line for our application studies. However we cannot theoretically guarantee that a solution is given within a sampling time interval. In the event that, at a given time instant tk , the RG fails in the sense that the modified reference computation time exceeds the sampling interval, it is possible to resort to a non-optimal though feasible, and compatible with constraints, strategy g (tk ) ← g¯ (tk−1 ), g¯ (tk−1 ) denoting the optimal RG modified reference computed at the previous time instant. This strategy is justified from a theoretical point of view by the viability property mentioned in the previous subsection and discussed in Bemporad and Mosca (1995), and Gilbert et al. (1995). The computation of the set W ξ is a key point in the numerical set-up of the RG. The numerical characterization of the set C∞ , which is mandatory to compute W ξ , can be obtained by a a approximating it with a convenient C∞ (ε) such that C∞ (ε) ⊂ a C∞ ⊂ C∞ (ε) + Bε , where Bε represents a ball of radius ε (safety level) centered at the origin. Such a set is computable in a finite number of steps. In fact it can be shown that

C∞ = Ck ∼ Hc Φ i Gd D .

(18)

Matrix Φ is a stability matrix and D is a compact, closed and convex set. This implies the joint of two positive  existence  constants M, λ ∈ (0, 1) such that Φ k 2 ≤ M λk , and of the quantity dmax := maxd∈D ∥d∥2 . As a consequence we have that, for a given ε > 0, there always exists a positive integer such that Hc Φ i Gd D ⊂ Bε ,

∀k > k ε .

(19)

i =k

Once the prescribed tolerance is fixed and M , λ and dmax are determined, thanks to the following inequality dmax σ¯ (Hc )σ¯ (Gd )M

∞ 

λi ≤ ε

(20)

i=kε

the quantity kε can be numerically computed as kε =

ln(ε) + ln(1 − λ) − ln (σ¯ (Hc )σ¯ (Gd )Mdmax ) ln(λ)

(22)

and the approximation (22) can then be used to efficiently compute the set in the following way:

Φ i−h−1 Gω + TLω ≤ f h

∗

∞ 



.

(21)

W ξ = ω ∈ Rnr : THc Inx − Φ

≤ f kε



 − (ξ + ε) TjT Tj ,

−1

Gω + TLω



j = 1, . . . , nt .

(23)

5. Numerical simulations on the ITER tokamak The high performance requested during ITER operations in terms of plasma volume and elongation, and the high costs needed to increase the machine flexibility, led us to consider operating conditions which are very close to physical limitations in terms of currents and voltages, induced forces and maximum fields on coils, and shape parameters. A 15 MA plasma at SOF (Start Of Flat top) has been chosen to evaluate the performance of the RG control strategy. Nominal values for the li and βpol are 0.85 and 0.1 respectively. A total number of 87 active hard constraints is imposed:

• 11 (+11) constraints on maximum (minimum) active currents: (±) 44 kA for CS, 45 kA for P1, P3, P4 and P6, 40 kA for P2 and P5;

• 11 (+11) constraints on maximum (minimum) voltages: 3.0 kV on CS1, 1.5 kV on all the other circuits;

• 11 constraints for maximum field in coils: 13 T on CS, 6 T for P1, 4 T for P2-P4, 5 T for P5, 6 T for P6;

• 38 (+38) constraints on gaps (see Fig. 2); • 2 constraints on the minimum distance between first and second plasma separatrix (20 mm, both on left and right sides);

• 12 constraints on vertical forces induced on CS ± 75 MN. The simulation analysis has been carried out considering ‘‘plasma maneuvers’’, characterized by desired variations of the plasma shape with respect to the nominal equilibrium conditions. The application of the RG strategy has been dictated by the violation of physical limits observed in several constrained outputs, when forcing the pre-compensated system to make some of the simulated plasma maneuvers. The RG procedure has been set in operation on a full nonlinear model of the plasma evolution, implemented in the CREATE-NL free boundary equilibrium code, able to solve PDE problem (2) and (3) via FEM. This is a detailed simulation code experimentally validated and used in the design of controllers for several tokamaks including JET, Albanese et al. (2003). In order to fit model (4) into the standard structure (5) required to apply the RG (as in Fig. 6), the following assumptions are done. (i) Assuming δ I − L−1 LE δw as state variable the dependence of δ˙x on δ w ˙ in (4) disappears. (ii) Assuming as controlled output δ y − F δw the dependence of the controlled output on the disturbance disappears. This requires the measurability of li and βpol which are well estimated on tokamaks. (iii) The active coils and plasma currents are measurable, and a reconstruction of the reduced set of currents in the passive structure is available. A standard Zero Order Hold scheme with a sampling time Ts = 0.01 s has been used for the discretization of the precompensated system. A safety margin of ξ = 5% has been taken into account on every constraint and a virtual horizon k∗ = 92 has been computed solving the LP problem (15) through the off-line procedure presented in Section 4. Then, by solving problem (17) at every time instant tk , reference signals g (tk ), representing the best approximation of

M. Mattei et al. / Automatica 49 (2013) 169–177

175

Fig. 7. Reference signals r (t ) (blue continuous line) compared with the RG outputs g (t ) (red dashed line) (m).

r (tk ) compatible with the prescribed constraints, have been determined on-line. An identity Σ matrix has been assumed to equally weight all the input channels. The following discussion is referred to two simulations (one in the absence and one in the presence of external disturbances) made to reduce top gap 21, and move plasma towards the torus axis (gaps 9 and 30 are increased and reduced respectively), with a fixed position of LS and RS gaps. Fig. 7 shows a comparison between the original reference signals and the reference signals modified by the RG in the disturbance free case. The presence of a stationary deviation between the reference and the output is due to constraints on CS1 and P6 currents, vertical CS coil forces and dsep limitations. Fig. 8 gives an idea of what would happen if forcing the pre-compensated system without the RG action: constraints are violated and systems moves to instability (run had to be stopped around 2.5 s). Fig. 9 shows the evolution of some of the constrained variables in the presence of the RG. The same reference signals were assumed also in the presence of a disturbance on poloidal beta (increase of 0.05 in about 5 s) and internal inductance (increase of 0.15 in the same time). The presence of disturbances is in practice accounted for with a progressive reduction of the constraints set along the virtual time axis. In particular, starting from the initial set C , the prediction uses a set of suitable restrictions Ch , h = 0, . . . , +8, according to recursion (6), where conditions at infinity are approximated by the finite virtual horizon kε computed off-line through (21). Fig. 10 shows the nonlinear evolution of plasma shape through equilibrium snapshots at different time instants. It is worth remarking that, in the absence of RG action, the run had to be stopped around 2.5 s due to the saturation of some of the control variables and a consequent plasma–wall contact. The ‘‘quadprog’’ solver of the Matlab Optimization Toolbox has been used to solve QP problem (17) and the ‘‘cputime’’ routine has been exploited to obtain timing performance on a Dual Core Duo 2.4 GHz processor. Although theoretical results cannot assure that the computation time is compatible with the prescribed sampling period in both the application cases considered the required computational burden was acceptable. An average time response of 7.8 ms, for the

Fig. 8. Closed loop run without RG, nonlinear time evolution of voltages (request to the power supplies), normalized with respect to their saturation values. Red markers highlight the violation of constraints. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

disturbance free case, and 8.1 ms, for the disturbance acting case, was observed, both compatible with the 10 ms sampling time used for the discretization. Nevertheless, as noted in Section 3, in the event that, at a certain step, the current RG reference is not available ‘‘in time’’, as happened in a few snapshots in our simulations, control inputs can be safely recovered by applying the modified reference signal provided by the RG at the previous time step. 6. Conclusions Large and costly nuclear fusion tokamaks as ITER will push the operation towards the limits of their envelope forcing also feedback controllers to face with nonlinear constraints. This paper explored the possibility to apply a Reference Governor constrained control strategy to the shape control of plasmas. Voltage saturations on the power supply converters, currents limitations in the active control coils, minimum clearance between

176

M. Mattei et al. / Automatica 49 (2013) 169–177

Fig. 9. Closed loop simulation with RG, nonlinear time evolution of some constrained outputs normalized to their maximum allowable value (no disturbance case).

Fig. 10. Result of the simulation in the presence of disturbance. Starting initial equilibrium (magenta line), plasma shape evolution in the presence (blue line) and in absence of RG action (green line) are shown. In the latter case, run stops at about 2.5 s because of loss of control. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the plasma surface and the vacuum chamber wall, maximum induced magnetic fields and forces on coils are explicitly taken into account by this model based predictive constrained control technique. Results on ITER 15 MA plasmas appear to be encouraging. The RG allows us to enforce constraints and seems to be tractable from a computational point of view. It is expected that large margins can be obtained with high performance hardware and optimized software.

References Albanese, R., Calabró, G., Mattei, M., & Villone, F. (2003). Plasma response models for current, shape and position control in JET. Fusion Engineering and Design, 66–68, 715–718. Albanese, R., & Villone, F. (1998). The linearized CREATE-L plasma response model for the control of current, position and shape in tokamaks. Nuclear Fusion, 38(5), 723–738. Ambrosino, G., Ariola, M., Pironti, A., & Walker, M. (2001). A control scheme to deal with coil current saturation in a tokamak. IEEE Transactions on Control Systems Technology, 9(6), 831–838.

M. Mattei et al. / Automatica 49 (2013) 169–177 Ariola, M., et al. (1999). A modern plasma controller tested on the TCV Tokamak. Fusion Technology, 36(2), 126–138. Bemporad, A. (1998). Reference governor for constrained nonlinear systems. IEEE Transactions on Automatic Control, 43(3), 415–419. Bemporad, A., & Mosca, E. 1995 Nonlinear predictive reference governor for constrained control systems. In: Proceedings of the 34th IEEE Conference on Decision and Control, New Orleans (LA), pp. 1205–1210. Borrelli, F., Falcone, P., Pekar, J., & Stewart, G. (2009). Reference governor for constrained piecewise affine systems. Journal of Process Control, 19(8), 1229–1237. Casavola, A., Mosca, E., & Angeli, D. (2000). Robust command governors for constrained linear systems. IEEE Transactions on Automatic Control, 45(11), 2071–2077. Crisanti, F., Albanese, R., Ambrosino, G., Ariola, M., Lister, J., Mattei, M., Milani, M., Pironti, A., Sartori, F., & Villone, F. (2003). Upgrade of the present JET shape and vertical stability controller. Fusion Engineering and Design, 66–68, 803–807. Diel, M., Bock, H.G., & Schloder, J.P. 2005 Real-Time iterations for nonlinear optimal feedback control. In: Proceedings of the joint Conference on Decision and Control - European Control Conference, Sevilla (Spain), pp. 5871–5876. Gilbert, E. G., & Kolmanovsky, I. (2002). Nonlinear tracking control in the presence of state and control constraints: A generalized reference governor. Automatica, 38(12), 2063–2073. Gilbert, E. G., Kolmanovsky, I., & Tin Tan, K. (1995). Discrete-Time reference governors and the nonlinear control of systems with state and control constraints. IEEE Transactions on Automatic Control, 5, 478–504. Humphreys, D. A., Deranian, R. D., Ferron, J. R., Jayakumar, R. J., Johnson, R. D., Khayrutdinov, R. R., La Haye, R. J., Leuer, J. A., Makowski, M. A., Penaflor, B. G., Walker, M. L., & Welander, A. S. (2005). High performance integrated plasma control in DIII-D. Fusion Engineering and Design, 74, 666–669. Kolmanovsky, I., & Sun, J. (2006). Parameter governors for discrete-time nonlinear systems with pointwise-in-time state and control constraints. Automatica, 42(5), 841–848. Kothare, M. V., Balakrishnan, V., & Morari, M. (1996). Robust constrained model predictive control using linear matrix inequalities. Automatica, 32(10), 1361–1379. Mayne, D. Q., Rawlings, J. B., Rao, C. V., & Scokaert, P. O. M. (2000). Constrained model predictive control: stability and optimality. Automatica, 36(6), 789–814. Vahidi, A., Kolmanovsky, I., & Stefanopolou, A. (2007). Constrained handling in a fuel cell system: A fast reference governor approach. IEEE Transactions on Control Systems Technology, 15(1), 86–98. Varano, G., Boncagni, L., Galeani, S., Granucci, G., Vitale, V., & Zaccarian, L. 2010 Results on plasma position and elongation regulation at FTU using dynamic input allocation. In: Proceedings on the 47th Conference on Decision and Control, Atlanta (GA), USA, pp. 2753–2758. Vitelli, R., Boncagni, L., Mecocci, F., Podda, S., Vitale, V., & Zaccarian, L. 2010 An antiwindup-based solution for the low current nonlinearity compensation on the FTU horizontal position controller. In: Proceedings on the 47th Conference on Decision and Control, Atlanta (GA), USA, pp. 2735–2740. Wesson, J. (1997). Tokamaks (second ed.). Oxford: Clarendon Press.

177

Massimiliano Mattei received the Laurea degree in aeronautical engineering and the Ph.D. in electronic engineering from the University of Napoli ‘‘Federico II’’ in 1993 and 1997 respectively. Professor Mattei is an expert in the area of modeling and control of plasmas in thermonuclear fusion reactors. He regularly collaborates with JET, ITER, Fusion for Energy, EFDA, ENEA, and many other worldwide institutions in this area. His research interests also include flight control and autonomous vehicles guidance, control and navigation. He is a Professor of Flight Mechanics and Control at the Second University of Napoli where he is also the head of the Department of Industrial and Information Engineering.

Carmelo Vincenzo Labate received the Laurea degree in electronic engineering ant the Ph.D. in automation engineering from the University Mediterranea of Reggio Calabria, Italy, in 2007 and 2011 respectively. His current research interests include predictive constrained control and nuclear fusion, with main applications in modeling and control of tokamak plasmas. Since 2009 he regularly cooperates with CREATE Consortium, ENEA, EFDA and F4E, via several internships as Visiting Researcher at JET, (Culham, UK) and FTU (Frascati, IT) labs. Since 2011 his research activity, carried out as Visiting Researcher at VTT Technical Research Centre of Finland (Tampere, FI), is focused on remote handling of tokamak components.

Domenico Famularo received the Laurea degree in computer engineering from the University of Calabria, Italy, in 1991 and the Ph.D in computational mechanics from the University of Rome, Italy, in 1996. From 1991 to 2000 he was a Research Associate at the University of Calabria. In 1997 he was a visiting Scholar Research at the University of New Mexico (Albuquerque, NM, USA) and in 1999 he covered the same position at the University of Southern California (Los Angeles, CA). He was a Researcher at the Istituto per il Calcolo e le Reti ad Alte Prestazioni (ICAR) - Consiglio Nazionale delle Ricerche (CNR) and since 2005 he is an Associate Professor at the University of Calabria. His current research interests include control under constraints, control reconfiguration for fault tolerant systems and networked control systems.