Gray-box modeling of resistive wall modes with vacuum-plasma separation and optimal control design for EXTRAP T2R

Fusion Engineering and Design 121 (2017) 245–255

Contents lists available at ScienceDirect

Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes

Gray-box modeling of resistive wall modes with vacuum-plasma separation and optimal control design for EXTRAP T2R A.C. Setiadi a,∗ , P.R. Brunsell a , F. Villone b , S. Mastrostefano b , L. Frassinetti a a b

Department of Fusion Plasma Physics, School of Electrical Engineering, KTH Royal Institute of Technology, SE-10044 Stockholm, Sweden DIEI, Universita di Cassino e del Lazio Meridionale, Italy

h i g h l i g h t s • Gray-box model of RWM that combines finite element model and system identification. • Formulation of optimal model based control for RWM that unify different control logic. • Experimental validation of the gray-box model and controller in EXTRAP T2R RFP.

a r t i c l e

i n f o

Article history: Received 17 January 2017 Received in revised form 6 June 2017 Accepted 11 July 2017 Keywords: Resistive wall mode Reversed-field pinch CARIDDI System identification Optimal feedback control

a b s t r a c t This paper presents a graybox methodology to model the resistive wall mode instability by combining first principle approach and system identification technique. In particular we propose a separate vacuum and plasma modeling with cascade interconnection. The shell is modeled using CARIDDI code which solves the 3D integral formulation of eddy current problem, whereas the plasma response is obtained empirically by system identification. Furthermore the resulting model is used to design an optimal feedback control. The model and feedback control is validated experimentally in EXTRAP T2R reversed-field pinch, where RWMs stabilization and non-axisymmetric mode sustainment is considered. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Successful control of magnetohydrodynamics (MHD) instabilities is important for the operation of a fusion device. The resistive wall mode (RWM) is of concern for the advanced tokamak operation scenario, in which the MHD mode stability of the configuration relies on the presence of a conducting boundary (wall) nearby the plasma. The RWM is a global, kink-like MHD mode, which grows slowly on the time scale of the magnetic field diffusion time of the wall. The origin of the RWM is a wall-stabilized ideal MHD mode, but there are also non-MHD effects that are important for the stability of the RWM, such as plasma rotation in combination with dissipation processes and kinetic effects, mainly resonances with particle drift motions [1]. RWMs are also expected for the Reversed Field Pinch (RFP) devices operated with a thin (resistive) shell. Current driven, nonresonant ideal MHD modes with poloidal mode m = 1 and a wide

∗ Corresponding author. E-mail address: [email protected] (A.C. Setiadi). http://dx.doi.org/10.1016/j.fusengdes.2017.07.011 0920-3796/© 2017 Elsevier B.V. All rights reserved.

range of toroidal mode number n are typically observed experimentally. Active magnetic feedback control of multiple RWMs has been successful on thin shell RFP devices such as RFX-mod and EXTRAP T2R. In addition to the non-resonant RWMs, the EXTRAP T2R plasma exhibits a range of internal tearing modes. These modes are typically observed at saturated amplitudes, rotating toroidally with the plasma fluid. Control of internal modes, in the sense of preventing the appearance of locked modes by sustaining the mode rotation by active feedback control is in principle possible, and has been demonstrated in several tokamak experiments [2,3]. However, in EXTRAP T2R, the typical mode rotation frequencies are much higher than the inverse shell time, and since the active coils are outside the shell, an in-phase feedback control of these modes by high frequency rotating fields would be prevented by the magnetic field screening effect of the shell. The active (low frequency) RWM feedback control in EXTRAP T2R is observed to help sustain the internal mode and plasma rotation, an effect which is attributed to the dynamic suppression of error fields resulting from the feedback control. Since the RWM grows on a relatively slow time scale, active magnetic feedback control is a possibility for achieving reliable

246

A.C. Setiadi et al. / Fusion Engineering and Design 121 (2017) 245–255

advanced tokamak operation. It has been shown in [4–6] that stabilization of RWMs allows the tokamak to reach a normalized beta close to the ideal wall limit. Furthermore, the RWMs, with growth rates determined by the eddy current decay rates, will have to be stabilized if the fusion device needs to operate in steadystate. The reversed-field pinch (RFP) devices such as RFX-mod and EXTRAP T2R are more prone to the current-driven RWMs, since the typical safety factor in RFP is lower compared to the tokamak. There are varieties of controllers that have been proposed to stabilize the RWM both in tokamak and RFP, from the conventional Proportional-Integral-Derivative [7–10] controller to the advanced state-space based controller [5,6,11–14]. The use of an advanced model based controller have been reported to be beneficial for suppressing the RWM [15], owing to the ability of the controller to filter noise, discriminating transient event and optimizing the feedback actuation level [15]. The use of model based controller requires the availability of a model of the plant. The model usually can be obtained by first principle (“white-box” approach) or empirically using system identification technique (“black-box” approach). Several different models are available in interpreting the plasma response that is affected by the RWM. The MARS-F code [16], which includes the toroidal plasma flow, is very well suited in investigating the RWM effects to the plasma response. However, MARS-F is rather limited in capturing the interaction of the plasma with the surrounding 3D geometry. Several finite element codes are available such as CarMa [17,18], and VALEN [19] that are able to model the 3D walls surrounding the plasma and have been proposed for several fusion devices [5,6,9]. The use of these finite element codes usually leads to a complex model which might be an issue when designing a controller. An approximation technique to reduce the complexity of the model is often necessary [5,9,12,20] before using the model for real time controller implementation. A model can also be obtained empirically using system identification technique. The system identification have been extensively used in EXTRAP T2R to obtain RWM model [21–23]. This paper presents a novel RWM modeling approach combining the first principle model and system identification. By proposing a vacuum-plasma cascaded structure the two approaches can be combined to produce a vacuum and plasma model that captures the RWMs dynamics. Dividing the coil-shell-plasma system in two parts is motivated by the fact that it enables more flexibility in the modelling: The coil-shell transfer function can readily be obtained by general FEM codes that includes the 3-D structures of the conducting wall and coil structures (for example the CARIDDI code). In order to obtain the plasma-shell transfer function, taking account of 3-D features of the wall structure, one would have to rely on advanced, specialized numerical codes that combine plasma stability (e.g. MARS, DCON [24] codes) calculation and FEM modeling (e.g. CARIDDI, VALEN codes). However, other possibilities also exist, for example, both models could be obtained by so called empirical modelling, utilizing experimental system identification methods. A combination of numerical and experimental modeling methods is another option. In this work, FEM modelling for the shell is used in combination with empirical modeling by system identification for the plasma. An advantage with this combination of numerical and empirical modeling is that the empirical model to some extent can compensate for inaccuracies in the numerical model, when both are combined to model the plasma case. The modeling approach used in this paper combined the CARIDDI code with the system identification method. Both of these tools are generic, thus the proposed modeling approach is not limited to RFP configuration. Furthermore, an optimal state space controller is designed using the proposed model. Both model and controller is experimentally validated in EXTRAP T2R. Furthermore the relation of the controller with the more well-known RWM con-

Fig. 1. Block diagram of interconnection between amplifier, vacuum response and plasma response.

trol scheme such as the mode control and explicit mode control [25] is discussed. The paper is organized as follows. Section 2 focuses on describing the modeling approach. Section 3 shows the validation result for the obtained model. Section 4 focuses on the control design methodology. Section 5 discusses the experimental result of the feedback control. Several aspects of the feedback are discussed such as RWM stabilization, non-axisymmetric mode sustainment and robustness of the feedback control under different plasma parameters. Finally Section 6 concludes the paper. 2. Identification of EXTRAP T2R The EXTRAP T2R device has a major radius R = 1.24 m, and a minor radius of a = 183 mm [26]. The conducting boundary of the EXTRAP T2R device consists of three separate layers. The first (innermost) layer is a continuous, thin-walled bellows, stainless steel vacuum vessel extending over the minor radius r = 187–196 mm. The second and third layer of the conducting boundary is made of two thin (0.5 mm thick) layers of copper centered at r = 197.3 and r = 198.3 mm respectively. Each of the two layers has one poloidal and one toroidal gap. Passive shielding of the magnetic field error due to the narrow gaps is accomplished by having the gaps of the two layers displaced relative each other, both poloidally (inboard vs outboard) and toroidally (about half a turn). An array of 128 sensor flux loops at 32 toroidal and 4 poloidal positions measures the radial magnetic field component at the conducting boundary. The flux loops are installed in the interspace between the vessel and the shell. The flux loops are pairwise connected (outboard-inboard, and top-bottom respectively) to form 64 sensor coils. Analog time-integration of the flux loop signals is used for obtaining a set of sensor signals proportional to the radial magnetic fields at the sensor coil positions. There is an array of 128 active saddle coils outside the shell structure at r = 238 mm, with the angular positions (coil centres) coinciding with the flux loops. Similarly as for the sensor flux loops, the active saddle coils are pairwise connected to form 64 active coils. Each active coil is driven by (audio type) power amplifier. The availability of an LTI model of EXTRAP T2R will ease the feedback control design process and real – time implementation of such feedback control. From previous experience, the use of LTI model seems sufficient to capture the dynamics of RWMs in EXTRAP T2R for some small variation of the equilibrium condition [21,22]. The equilibrium of RFP is experimentally characterised by the pinch parameter  = B,a /B  and reversal parameter F = B,a /B . Here B is the magnetic field where the subscript ,  denote the toroidal and poloidal components, the subscript a denotes the evaluation at minor radius and .  is the volume average operator. To model the RWMs in EXTRAP T2R, there are three separate subsystems that need to be considered, namely the power supply for the active coil GA , the shell GV and the plasma RWM dynamics GP . The interconnection of these subsystems is depicted in Fig. 1 and the variables descriptions are given in Table 1. GA is the amplifiercoil transfer function, describing the response of the amplifier with the control coil load. The input variable of GA u indicates the vector

A.C. Setiadi et al. / Fusion Engineering and Design 121 (2017) 245–255 Table 1 Summary of variables for the cascade interconnected model. Variables

Description

u uV y yV yP GA GV GP

Voltage input to the amplifier Current input to the shell Measured radial magnetic field Radial magnetic field generated by the shell Radial magnetic field generated by the plasma Power amplifier transfer function Coil-shell transfer function Plasma transfer function

of 64 control voltages applied to set of the power amplifiers driving the control coils. The output of GA is the vector of 64 current in the coils (uV ). The function describes both the electrical characteristics of the amplifier (gain, bandwidth) and the coil (resistance, inductance). Note that in Fig. 1 we proposed a cascaded structure for the shell and the plasma. This structure is proposed based on the idea that the total magnetic field is a combination from sources that are internal in the plasma and sources that are external to the plasma [21,27]. Furthermore we combined two approaches, “white-box” and “black-box” in obtaining the model for EXTRAP T2R which we will elaborate later. The GV is the coil-shell transfer function, describing the response of the coil and the shell. The input to GV is the vector of 64 control coil current (uV ). The output of GV is the vector of 64 radial magnetic field at the shell radius in a vacuum case without plasma (yV ). The function is determined by the geometrical properties of the coil (width, length, position) and the shell characteristics (thickness, conductivity). GP is the plasma transfer function, describing the change of the vacuum field in a plasma case due to the plasma RWM response. The input to GP is the vector of vacuum radial magnetic field at the shell radius (yV ). The output is the contribution to the radial magnetic field at the same shell radius when plasma is present (yP ). The response is determined by both the plasma equilibrium and the shell characteristics. In EXTRAP T2R it is found that the radial field diffusion time is different in horizontal and vertical direction [21]. In general the radial field diffuses faster for horizontal direction than for vertical direction. This can be attributed to the horizontal gaps of the EXTRAP T2R shell. Furthermore, for RFP, it is known that the growth rate of RWM depends on the plasma equilibrium. For instance, experiments and ideal MHD calculation shows that by increasing the pinch parameter the growth rates of internal RWM decrease while at the same time increase the external RWM growth rates [28,29]. The total radial magnetic field at the position of the shell in the presence of plasma is summed fields, y = yV + yP . Note that, the vacuum magnetic field (yV ) is not the same as the external field produced by the coil as input to a combined “shell-plasma” system. The field yV includes also contribution from eddy current in the shell. Secondly, the plasma correction (yP ) is also in part due to shell eddy current. Thus, the shell characteristics are represented in both functions GV and GP . The “white-box” approach or first principle modelling is used to model the shell. The electrodynamics feature of the shell is modeled by the 3D eddy current code CARIDDI [30]. Note that the CARIDDI code can be coupled with plasma equilibrium solver code to incorporate the RWM dynamics [17,18,31], however in this work we propose to use system identification in combination with the CARIDDI. The motivation behind this approach is related to the modelling error. The discrepancies between the CARIDDI model and measured data might be due to some un-modeled dynamics and some noise. Furthermore, the reduced order model will introduce an approximation error. Since the system identification is a data-driven method which aims to minimize the error between model and measured data, combining both approaches might be

247

useful to compensate the error due to un-modeled dynamics and approximation error. Furthermore, typical formulation of system identification directly treats the influence of the noise in the system which is missing in the CARIDDI formulation. 2.1. Power supply identification The first subsystem to be identified is the power supply for the active coils GA . The power supply used in EXTRAP T2R is an array of audio amplifiers. The frequency response of each channel can be approximated by the transfer function kA /( A s + 1). The time constant  A includes both amplifier and coil, and is mainly due to the coil L/R time. For an array of 64 uniform amplifiers, a state space model can be formed:

GA =

⎧ 1 ⎪ x˙ (t) = − I64 xA (t) + I64 u(t) ⎪ ⎪ A A  ⎪ ⎪   ⎨ BA AA

(1)

kA ⎪ uV (t) = xA (t) ⎪ ⎪ A ⎪ ⎪  ⎩ CA

with I denotes the identity matrix. The gain kA and time constant  A can be estimated from open-loop experiment [32]. The coil inductance is affected by the coil proximity to the shell and the resulting modification of the coil time constant is implicitly included in GA through the empirical determination of  A using vacuum measurements. The separation of the amplifier, coil and shell system into two blocks for amplifier-coil (GA ) and coil-shell (GV ) respectively is motivated by the possibility to measure the currents in the active coils. 2.2. Shell model 2.2.1. CARIDDI model of EXTRAP T2R CARIDDI is a finite element code that solves the integral formulation of eddy current problem. Detailed formulation of the CARIDDI code can be seen in [17,30], only the main points are summarized in this paper. The integral formulation of eddy current problem is given by:



J.wdV + Vc

∂ ∂t





∇ V.wdV = 0

A.wdV + Vc

(2)

Vc

Here Vc is the three-dimensional conducting domain, A is the magnetic vector potential which can also be expressed in current density J by the Biot-Savart integral; and V is the electric scalar potential. Additionally,  is the resistivity and w is a suitable weighting function. Electric vector potential T is introduced to guarantee the solenoidality of the current density (J = ∇ × T). The conducting domain Vc is then discretized and the current density is expanded in terms of the edge elements Nk (J = k Ik ∇ × Nk ). By using Galerkin method to solve (2), the following relation is obtained: L

dI V (t) duV (t) + RIV (t) + Mc =0 dt dt

(3)

yV (t) = Qx IV (t) + Qc uV (t) with R, L, Mc , Qx , Qc are the appropriate resistance and inductance matrices [17], whereas IV (t) represents current in the 3D structure and uV represents the current in the coils. Since the matrices are static, the resulting CARIDDI equation (3) belong to LTI system that can be represented in state space form. One possible state space realization can be obtained by introducing a state vari-

248

A.C. Setiadi et al. / Fusion Engineering and Design 121 (2017) 245–255

able xV (t) = LIV (t) + Mc uV (t). The state space representation of the shell is given by:

⎧ x˙ V (t) = (−RL−1 )xV (t) + (RL−1 Mc )uV (t) ⎪ ⎪     ⎨

GV

AV

BV

yV (t) = (Qx L−1 )xV (t) + (Qc − Qx L−1 Mc )uV (t) ⎪ ⎪     ⎩ CV

(4)

DV

The internal state variable in the vacuum state space model xV has the physical interpretation as the magnetic flux linked by the elementary currents loops corresponding to the edge-element basis functions used for the spatial expansion of the current density on the finite elements mesh of the conductor. The vessel is modelled in CARIDDI as a thin toroidal layer with different resistivity in the toroidal and poloidal direction, simulating the anisotropic effect produced by the bellows structure. In order to avoid overlapping of the thin closely spaced layers of the conducting boundary, many elements are required in the angular directions. The 3D finite element model of the vessel is constructed of one element in the radial, 25 elements in the poloidal direction and 60 elements in the toroidal direction, resulting in 1500 conducting elements in total. The 3D mesh of EXTRAP T2R is shown in Fig. 2. Each copper layer is modelled similarly as the vessel but the model consists of somewhat less elements since part of the underlying toroidal grid at the gaps is not used. The CARIDDI code generates a vector space with a dimension related to the number of edges in the finite element mesh [30]. This is also referred to as the degrees of freedom of the model since the actual dimension corresponds to the number of independent loops that can be formed by suitably combining the edges, taking also into account the topological connection degree of the conducting domain. The dimension of present three-layer model of the EXTRAP T2R conducting boundary is 4285. For simplicity, this complex conducting boundary is referred to as the “shell” in the following. 2.2.2. Model reduction The complexity of the model can be an issue if we aim to embed the model into our feedback control algorithm. In such case we need to be able to solve our model in real-time (sampling time in EXTRAP T2R is 0.1 ms). For EXTRAP T2R, it is not feasible to use the full order CARIDDI model in real-time due to the computation limitation. For this reason we employ several model reduction techniques to obtain the optimal low order approximation of the CARIDDI model. Several model reduction techniques have been considered for EXTRAP T2R and the result is presented in Ref. [20]. It is reported in Ref. [20] that the optimal Hankel approximation technique performed better compared to other methods when approximating the shell. Recall that the Hankel norm of a LTI system G is given by

∞

||G||H = sup  00 u ∈ L2

y(t)2 dt

u(t)2 dt −∞

(5)

The optimal Hankel approximation aims to find a low-order approximation of the shell GVr which minimizes the Hankel norm of the error system ||GV (t) − GVr (t)||H . The analytical solution of optimal Hankel approximation for state space have been solved in [33] and the error bound of the approximation can be expressed as a function of the discarded singular values of the original system. A reduced order model of order 118 have been obtained which sufficiently describes the dynamics of the complex CARIDDI model [20]. 2.3. Identification of RWM dynamics The second approach that is used is the “black-box” approach commonly called system identification. System identification is an

empirical methodology to estimate parameters of a generic model such that the resulting model match the input-output relation as observed from measurement data. To obtain the data for identification a series of experiments were performed. The experiment usually consists of applying a dither signal in the form of pseudorandom signal during feedback controlled plasma operation. A new batch of data was obtained that correspond to a total of 3 s duration of plasma. For each experiment, similar equilibrium parameters  and F evolution are reproduced, this can be seen in Fig. 3. Some data preprocessing is required to decouple the vacuum contribution and plasma contribution from the measured field. It follows from the interconnection that the input to the plasma model is the vacuum magnetic field at the sensor position. The CARIDDI model along with the amplifier model is used to generate a vacuum field contribution given the experimental input data. Additionally the output of the plasma contribution is obtained from the measured experimental data by subtracting with the generated vacuum field contribution. The model parameters estimation procedure follows the optimized Predictor Based Subspace Identification (PBSIDopt)proposed in Ref. [34], where the algorithm aims to estimate a discrete time state space:



GP =

xP (k + 1) = AP xP (k) + BP yv (k) + Ke(k) yP (k) = CP xP (k)

(6)

The matrix K is the Kalman matrix which is inherently estimated by the PBSIDopt algorithm. Furthermore, user can define the order of the model to be estimated. The order of 150 is chosen for EXTRAP T2R based on the observation of the singular value decays of the extended observability matrix [34,35]. 2.4. Vacuum-plasma separation This subsection illustrates the motivation of using the cascaded interconnection shown in Fig. 1. We started our analysis using the CarMa model circuit equation. The coupling of conducting structures and plasma in CarMa relies on defining a coupling surface which is located in between the plasma and the conducting structures. The plasma response to a set of external magnetic field perturbations is computed from MHD equations and represented as static matrix S, in the assumption that the inertia of the plasma is negligible, so that the plasma response is instantaneous. This plasma response is represented as an equivalent surface current Ieq located on the coupling surface that produces the same field as the perturbation outside the coupling surface. Furthermore L denotes the inductance matrix of the shell, R is its resistance matrix and Meq , Mc are the mutual inductance matrices between the coupling surface-shell and actuator coil-shell respectively. The CarMa coupled circuit equation is given by: L

dI eq (t) dI(t) duV (t) + RI(t) + Meq + Mc =0 dt dt dt

(7)

where I is the vector of currents in the 3D surface, Ieq is the vector of current in the coupling surface of the plasma and shell and uV is the vector of current in the actuator coils. The magnetic field (or flux) and the current in the coupling surface for the CarMa model are given by: y(t) = Qx I(t) + Qc uV (t) + Qeq Ieq (t)

(8)

Ieq (t) = S (Qx I(t) + Qc uV (t))

(9)

where Qeq denote the constant matrices that express the effects of the 3D coupling surface current on the field. We assume that the shell current and the output magnetic field are the results of

A.C. Setiadi et al. / Fusion Engineering and Design 121 (2017) 245–255

249

Fig. 2. 3D mesh of EXTRAP T2R in CARIDDI.

Fig. 3. The plasma parameters observed during dithering experiment, shaded area represents the variation of the parameters in the experiments (a) plasma current, (b) pinch parameter  and (c) reversal parameter F.

superposition of the plasma and shell contribution, thus the circuit equation can be partitioned for the vacuum: L

dI V (t) duV (t) + RIV (t) + Mc =0 dt dt

yV (t) = Qx IV (t) + Qc uV (t)

(10) (11)

and consequently for the plasma: L

dI eq (t) dI P (t) + RIP (t) + Meq =0 dt dt

(12)

yP (t) = Qx IP (t) + Qeq Ieq (t)

(13)

Ieq (t) = SyV (t) + SQx xP (t)

(14)

The plasma part can be written as: L∗

dI P (t) dy + RIP (t) + Meq S V = 0 dt dt



yP (t) = Qx + Qeq SQx IP (t) + Qeq SyV (t) ∗

(15) (16)

with L = L + MSQx . From here, several observations can be made. The driving term for the vacuum contribution is the actuator cur-

rent whereas the plasma contribution seems to be driven by the vacuum field. This justifies the proposed cascade interconnection of Fig. 1 and that under the superposition assumptions, the resulting cascade model is consistent with the CarMa model. Thus the isolation of plasma response from measurement can safely be done by subtracting the measured field with the computed vacuum response. In practice, the procedure use the measurements of the magnetic sensors; this practically works provided that these sensors provide a sufficient information about the magnetic field around the plasma up to the desired toroidal and poloidal harmonic content. In this paper, the vacuum contribution is modeled using CARIDDI while the plasma contribution is modeled using system identification. Both approaches will attempt to estimate the system matrices independently. The plasma part of the sensor field yP is here understood as the change of the sensor field y due to the plasma response to a given vacuum field yV (block GP in Fig. 1). The plasma part defined in this way includes also a part which is due to the shell response to the plasma. The plasma contribution of the field yP is obtained experimentally from plasma measurements as

250

A.C. Setiadi et al. / Fusion Engineering and Design 121 (2017) 245–255

yP = y – yV , subtracting out the shell response to the coil field. The part subtracted out is precisely the part that is computed with the CARIDDI code in our case. The present splitting of plasma and vacuum contribution is different from that used in the plasma-shell coupling scheme of the CarMa modelling code [17]. In the CarMa coupling scheme, the plasma response (matrix S in Eq. (9)) is understood as a static response to the external magnetic field, i.e. the magnetic field due to the total currents flowing in the shell. Conversely, in the present approach, the plasma contribution (yP in Eq. (16)) is the dynamical response to the vacuum magnetic field, i.e. the magnetic field due only to vacuum shell currents. The advantage of the present separation of the vacuum and plasma response is that it enables the CARIDDI numerical computation to be combined with empirical (experimental) modelling results in a straight-forward way. The proposed splitting does not fully benefit from passive structure modelling, because we subtract from experimental data only vacuum contribution, but the shell response is still retained in the system to be identified. Although there is no complexity reduction by the proposed splitting, one can argue that the use of system identification in addition to CARIDDI model and model reduction, can be used to compensate the modeling/approximation error in estimating the system by using information obtain from experimental data.

a specific mode number. Furthermore, the statistical consistency of the estimation for the plasma case is evaluated by generating the model with different combination of data sets and evaluating the function (18) similar to procedure performed in [21,22]. The result is displayed in Fig. 5. Fig. 5(a) was obtained from the reduced order CARIDDI model and consistent with past experiences of EXTRAP T2R [21,22] whereas Fig. 5(b) is obtained from the plasma model using system identification. The growth/decay rates obtained resemble the values obtained using ideal MHD calculation in cylindrical geometry [22,28]. 4. Controller design A state space feedback control has been synthesized based on the gray-box model. The control synthesis methodology follows from the standard Linear Quadratic Regulator (LQR) formulation which aims to find an optimal state feedback matrix KLQR such that the feedback actions u(k) =− KLQR x(k) minimized the cost function [36]: V = ∞ x(k) Qx(k) + u(k) Ru(k) + 2x(k) Mu(k) k=0 T

T

T

(19)

3. Model validation

where Q, R and M are positive definite weighting matrices that define the trade-off between state suppression and level of actuation in the cost function. The solution to this optimal control problem is found by solving the matrix equation:

3.1. Input–output agreement

AT SA − S − (AT SB + M)(BT SB + R)

The three subsystems represented by Eqs. (1), (4) and (6) are combined to form the interconnection in Fig. 1. Note that the amplifier model and shell model is first transformed into discrete time state space using first-order method and sampling time of 0.1 ms. The model is then compared with the experimental data to check the agreement with the model. The set of data used for validation is different than the set used for obtaining the model, to avoid any bias in the validation. The agreement is quantified by using Variance-Accounted-For (VAF) function which is defined by:



VAF i =

1−

exp

var yi

− yimodel exp

var(yi

)



× 100%

(17)

where var represents the variance operator. The results are summarized in Fig. 4 where 15 batches of data corresponding to 0.9 s plasma duration with flat-top current are used for the validation. The obtained model shows good agreement with the experimental data, where the worst agreement obtained is around 80%. The worst case agreement seems to be caused by a higher signal-to-noise ratio in that particular channel as can be seen in Fig. 4. 3.2. Modal analysis By performing a modal analysis technique, the growth/decay rate of the relevant toroidal mode number can be extracted from the state space model. The modal analysis technique used in this paper follows from Refs. [22,23] where a bivariate function is evaluated for the shell and plasma model separately:

 †   f ( , n) = log wn C( I − A)−1 B wn 

(18)

where A, B, C are the model matrices and wn is a vector that represent Fourier basis function for a toroidal mode number n and represent a candidate for eigenvalue. Note that the function is evaluated by scanning both the mode number n and the candidate of an eigenvalue . The function in Eq. (18) analogously can be interpreted as a correlation function between a particular Fourier harmonic with each of the candidate eigenvalues of the system . The strongest correlation value indicates the true eigenvalue for

The

state

feedback

−1

matrix

(BT SA + M T )Q = 0 is

constructed

(20) as

−1

KLQR =

(BT SB + R) (BT SA + M T ). The generic LQR formulation can be recast into mode control formulation which is more relevant for RWM control application in EXTRAP T2R since the RWM appears in multiple toroidal mode numbers. A mode control formulation of LQR is done by defining the cost function as: T T¯ V = ∞ y (k) Lyf (k) + u(k) Ru(k); yf (k) = WDFT y(k) k=0 f

(21)

where WDFT is a discrete Fourier Transform matrix that resolves the toroidal mode number of interest. The generic LQR formulation can be cast as the mode control by defining the weighting matrices as: ∗ Q = C T WDFT LW DFT C

R=D

T

M=C

∗ WDFT LDW DFT

T

(22) + R¯

(23)

∗ WDFT LW DFT D

(24)

The control tuning parameter in the mode control formulation of LQR is the diagonal matrix L where each of its diagonal entry represent the weight for a specific toroidal mode number, and R¯ which penalizes the actuator current. To simplify the tuning procedures, in this paper R¯ is set at fixed value of 1. The effects of weight selection will be discussed in Section 5. The following state observer is used which is based on the interconnection of the subsystems:

⎡ xˆ

A (k

+ 1)

⎤

⎣ xˆ V (k + 1) ⎦ xˆ P (k + 1)

⎡A

0

0

⎤ ⎡ xˆ

= ⎣ BV CA

AV

0

⎦ ⎣ xˆ V (k) ⎦

A

⎡ (BP − K)D ⎤ V CA (BP − K)CV BA 0   ⎦ u(k) +⎣ 0 0

AP − KC P

A (k)

xˆ P (k)

⎤

(25)

y(k)

0

K

The state feedback requires the state information which can be estimated from radial magnetic field measurement and voltage given to the amplifier following the relation in Eq. (25). Furthermore based on the state estimation, the output of every subsystem can be reconstructed. This open up possibilities of changing the control

A.C. Setiadi et al. / Fusion Engineering and Design 121 (2017) 245–255

251

Fig. 4. Validation result of the obtained model. (a) Variance-Accounted-For value for each radial magnetic measurement point. (b) Best case comparison of model and experimental data. (c) Worst case comparison of model and experimental data.

Fig. 5. Modal analysis of the shell and plasma. The color represent correlation of a candidate of an eigenvalue ( ) with toroidal mode number (n). The color of the figure represents the correlation value; where the red or yellow color represent strong correlation and blue represent low correlation. Red color in the plot indicates the location of the true eigenvalue for a specific toroidal mode number. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

logic of the feedback controller simply by modifying the weighting matrices. The formulation in Eq. (22) is analogous to the Raw Mode Control (RMC) [25], since the feedback loop aims to suppress the total contribution of the vacuum field and plasma. It is possible to perform Explicit Mode Control (EMC) logic where the vacuum field is subtracted from the measured field [25,8,7] simply by using the output estimation relation:



yˆ V (k) yˆ P (k)



 =

DV CA

CV

0

0

⎡  xˆ A (k) 0 ⎢ ⎣ xˆ V (k) CP

⎤ ⎥ ⎦

(26)

xˆ P (k)(k)

and setting the entries of matrix L that correspond to the estimated vacuum field yˆ V to zero.

Fig. 6. Plasma current dependency with varying state weight L. The case without RWM feedback is shown by the black dash-dotted line.

5. Results from feedback-controlled experiments 5.1. RWM stabilization We implemented the controller in EXTRAP T2R with several L values. In all experiment the weight of the axisymmetric mode is set to 0 to exclude the axisymmetric mode in the RWM feedback loop. The axisymmetric mode in EXTRAP T2R is controlled by the plasma equilibrium controller and it is beneficial to minimise the interference with the RWM controller [37]. Equal weight is set for the rest of the modes. The plasma current is shown in Fig. 6 where a positive correlation of plasma duration can be seen with increasing weight L. The effect of RWM feedback control is clearly observed in Fig. 6. With-

out feedback control the plasma terminates prematurely at around 15 ms, coincident with an exponentially growing radial magnetic field at the plasma boundary, indicating the growth of a linearly unstable RWM as can be seen in Fig. 7. Feedback control results in a stabilization of the RWM as shown by the suppression of the radial magnetic field growth. The early plasma termination is avoided, and the pulse length is with feedback control in the range 40–80 ms. The plasma pulse termination with feedback control is associated with the appearance of locked tearing modes at low plasma rotation. The plasma rotation is strongly affected by low amplitude radial (error) magnetic fields at the plasma boundary. The feedback control effectively performs both RWM stabilization and dynamical error field correction. The error field correction performed by the

252

A.C. Setiadi et al. / Fusion Engineering and Design 121 (2017) 245–255

Fig. 7. Root mean square of radial field dependency with varying state weight L. The case without RWM feedback is shown by the black dash-dotted line.

Fig. 9. Amplitude of one of the most unstable RWM mode with different control logic. The dashed lines in the figure represent the time when the feedback controller was activated.

Fig. 8. Merit functions dependency with varying state weight L.

feedback control acts to sustain the plasma rotation and thereby lengthen the plasma pulse. The plasma pulse length with feedback control is shown here as an indicator of the performance of the feedback controller in terms of error field suppression for the different cases. Note that, even with a relatively small weight L = 10 the plasma duration is improved when compared to the open loop case and significant improvement achieved when L is increased further. To illustrate the effect of L with the final suppression level of the field we introduce two merit functions: J1 =

15 

rmsn (y¯ f (k))

(27)

var n (y¯ f (k))

(28)

5.2. Mode sustainment

n=−16;n = / 0

J2 =

15 

resent total variation of the non-axisymmetric modes. Judging from J1 and J2 it seems that the optimal value for L lies around 100–150. We also perform the feedback experiment using the EMC logic scheme and the result is shown in Fig. 9. The EMC shot is obtained by removing the vacuum field contribution of −16 ≤ n ≤ −9 and 9 ≤ n ≤ 15. Fig. 9 shows that one of the most unstable RWM n =−10 is still stabilized by the EMC even though the feedback loop neglected the vacuum contribution. In terms of the remaining steady state error, it may appear that the EMC performs slightly worse than the RMC. However, the observed difference is expected since the EMC does not act on the part of the error field which is produced to the coil current. The steady state error is mainly due to intrinsic machine error fields which are compensated, or corrected, by the RWM controller. In order to perform this compensation, the RMC acts on the total error signal, that is the sum of the machine error field and the control field. Indeed, the RMC performs the dual goal of mode stabilization and error field compensation, while the EMC only targets mode stabilization. However the flexibility of the proposed controller to use different control logic can be useful to investigate the interaction of vacuum field to the plasma stability.

n=−16;n = / 0

red where the rmsn is the root-mean-square operator of a toroidal mode n and varn is the variance operator of a toroidal mode n. Furthermore, y¯ f is the radial magnetic field in Fourier domain, that is normalised with B,a . The normalisation is carried out to incorporate the plasma current variation in each shot. Both the merit functions are computed for the duration of plasma current flat-top. The effect of increasing L in J1 and J2 is shown in Fig. 8. The value J1 represent the overall magnetic field suppression by the controller, and can be seen to have positive effect with increasing L although at J1 seems to converge after the gain is higher than 200. However, in the experiment, it is observed that at high L = 500, some high frequency oscillation is appears. This might be an indication of approaching stability limit or the controller overcompensating the high frequency noise. This is captured by function J2 which rep-

The proposed feedback control can be extended to output tracking case, where instead of suppressing all of the modes, the feedback controller aims to sustain arbitrary modes to a specific amplitude and phase. This is done by defining the tracking error as a new state variable xe (k + 1) = xe (k) + (r(k) − y(k)). The controller synthesis process then follows in similar manner as in previous section. The experimental result is shown in Fig. 10, where the dashed line indicates the time when the feedback control is switched on. The black dots in Fig. 10 represent the amplitude of the nonaxisymmetric modes that are not sustained by the feedback. In Fig. 10 two different cases is shown. The first case is where the feedback controller is set to sustain n = 6 at amplitude 0.21 mT and the second case is with mode n = 10 set to the same amplitude. The mode n = 6 is expected to be an unstable mode where the mode n = 10 is expected to be a stable mode. This affects the actuator current waveform required to sustain the modes. Since n = 6 is unstable (or marginally stable) an initial “kick” requires to bring the mode to the desired amplitude but only a small current requires to maintain

A.C. Setiadi et al. / Fusion Engineering and Design 121 (2017) 245–255

253

Fig. 12. The plot of reversal parameter F against pinch parameter .

Fig. 10. Radial field sustainment result (a) toroidal mode n = 6, (b) toroidal mode n = 10, (c) actuator current for mode n = 6 sustainment and (d) actuator current for mode n = 10 sustainment. The dashed lines in the figure represent the time when the feedback controller was activated. The dotted lines in (a) and (b) represent the rest of the toroidal mode number during mode sustainment experiment.

it. Mode n = 10 is a stable mode therefore a higher constant current is required to maintain the field amplitude.

that the growth rates of the RWMs in EXTRAP T2R depends of the equilibrium parameters [28]. To this end, we evaluate the merit functions for the RWMs with internal helicity handedness (same as the RFP equilibrium field internal to the reversal surface)(m = 1, n < 0) and RWMs with external helicity handedness (same the RFP equilibrium field external to the reversal surface) (m = 1, n > 0) separately for multiple equilibrium parameters. J1int =

rmsn (y¯ f (k))

(29)

n=−1

5.3. Robustness against different plasma condition The robustness of the controller was also tested against several different plasma parameters. Fig. 11 illustrate the variation of plasma current, and the experimental equilibrium parameters. It is known that the parameters  and F cannot be changed independently, rather they follow a fixed monotonic relation (F is a monotonic function of )[38]. The F −  curve for the obtained experiments is shown in Fig. 12. Experimentally, it has been shown

−15 

J1ext =

15 

rmsn (y¯ f (k))

(30)

n=1

The result can be seen in Fig. 13. As the  increases from the default parameter, the merit function of the internal modes J1int initially decrease while the external modes J1ext tend to slowly increase. This is consistent with the observation in reference [28]

Fig. 11. The plasma parameters observed during robustness test experiment, shaded area represents the variation of the parameters in the experiments (a) plasma current, (b) pinch parameter  and (c) reversal parameter F.

254

A.C. Setiadi et al. / Fusion Engineering and Design 121 (2017) 245–255

the gray-box modelling scheme, together with a suitable model reduction technique can result in a practical control model for implementation in a model-based state space controller. The controller based on the cascade interconnection is flexible, enabling test of different control logic. Experimental results show successful stabilization of RMWs and arbitrary mode sustainment. A test of controller robustness against varying equilibrium parameters was done. Although a detailed comparison of control methods was not the aim of this work, it can be noted that the performance of the controller is comparable to previously tested controllers based on the empirical black-box modelling method. Acknowledgement The authors would like to thank Andrea Polsinelli for his help in producing the CARIDDI model. References Fig. 13. Merit functions dependency of reversal parameter  (a) internal modes and (b) external modes; dashed line represent approximation by 4th order polynomial, dotted line represent the default plasma parameters.

The analysis provided in [29,38], which is based on linear cylindrical MHD model also support this result. In [29,38] it was suggested that the change in the growth rates is due to the change of the vacuum potential energy. Since the helical winding of the internal modes are opposite of the external modes, thus the change in the potential energy affects both modes in an opposite manner. At higher value of  both J1int and J1ext increase with the rate of increase for the external modes noticeably faster. The best magnetic field suppression is observed when operating at slightly higher  than the default parameter at which the model was identified. This can be explained by the RWM growth rates dependence on . In the region of the tested equilibrium parameters, the performance degradation of the feedback controller is considerably smaller than the case with a sub-optimal gain L. In term of performance robustness of the feedback control, in the worst case shot that we encounter, the merit function increased to around 10% from the default region, this value is relatively small compared to the case with sub-optimal L. It is possible that this performance degradation can be compensated by re-tuning the weight L. 6. Conclusion We have presented a methodology to obtain gray-box model of RWMs in EXTRAP T2R by combining the reduced order CARIDDI model with the system identification technique. A general motivation for a gray-box modelling approach is to improve, or refine a white-box model using empirical data. The obtained model shows good input-output agreement and the eigenvalues that are consistent with linear cylindrical model prediction. The improved model is envisioned to perform better than white-box model since it can utilize the empirically obtained data to compensate modeling error which may arise from stochastic noise or limited information on the actual physical system. Furthermore it maintains physical interpretation of the model which might not be straightforward in a black-box model. The present results show that accuracy comparable to the black-box model can be obtained using the gray-box modeling method, with the additional advantage of providing a state space model with physically interpretable internal states. The model can readily be used for model based control design, where a mode control formulation of LQR is used. It has also been demonstrated, in experimental tests at the EXTRAP T2R device that

[1] Z.R. Wang, S.C. Guo, Y.Q. Liu, Drift kinetic effects on the resistive wall mode stability – comparison between reversed field pinches and tokamaks, Phys. Plasmas 19 (7) (2012) 072518, http://dx.doi.org/10.1063/1.4737200. [2] G.A. Navratil, C. Cates, M.E. Mauel, D. Maurer, D. Nadle, E. Taylor, Q. Xiao, W.A. Reass, G.A. Wurden, Active control of 2/1 magnetic islands in a tokamak, Phys. Plasmas 5 (5) (1998) 1855–1863, http://dx.doi.org/10.1063/1.872856. [3] F.A.G. Volpe, M.E. Austin, R.J.L. Haye, J. Lohr, R. Prater, E.J. Strait, A.S. Welander, Advanced techniques for neoclassical tearing mode control in DIII-D, Phys. Plasmas 16 (10) (2009) 102502, http://dx.doi.org/10.1063/1.3232325. [4] A. Garofalo, et al., Stability and control of resistive wall modes in high beta, low rotation DIII-D plasmas, Nucl. Fusion 47 (9) (2007) 1121–1130, http://dx. doi.org/10.1088/0029-5515/47/9/008. [5] O. Katsuro-Hopkins, J. Bialek, D. Maurer, G. Navratil, Enhanced ITER resistive wall mode feedback performance using optimal control techniques, Nucl. Fusion 47 (9) (2007) 1157–1165, http://dx.doi.org/10.1088/0029-5515/47/9/ 012. [6] S.A. Sabbagh, et al., Overview of physics results from the conclusive operation of the National Spherical Torus Experiment, Nucl. Fusion 53 (10) (2013) 104007, http://dx.doi.org/10.1088/0029-5515/53/10/104007. [7] E.J. Strait, et al., Resistive wall stabilization of high-beta plasmas in DIII-D, Nucl. Fusion 43 (6) (2003) 430–440, http://dx.doi.org/10.1088/0029-5515/43/ 6/306. [8] M. Okabayashi, et al., Active feedback stabilization of the resistive wall mode on the DIII-D device, Phys. Plasmas 8 (5) (2001) 2071–2082, http://dx.doi.org/ 10.1063/1.1351823. [9] G. Marchiori, et al., Dynamic simulator of RWM control for fusion devices: modelling and experimental validation on RFX-mod, Nucl. Fusion 52 (2) (2012) 023020, http://dx.doi.org/10.1088/0029-5515/52/2/023020. [10] P. Zanca, et al., Beyond the intelligent shell concept: the clean-mode-control, Nucl. Fusion 47 (11) (2007) 1425–1436, http://dx.doi.org/10.1088/00295515/47/11/004. [11] J.M. Hanson, et al., A digital control system for external magnetohydrodynamic modes in tokamak plasmas, Rev. Sci. Instrum. 80 (4) (2009) 043503, http://dx.doi.org/10.1063/1.3112607. [12] M. Ariola, G. De Tommasi, A. Pironti, F. Villone, Control of resistive wall modes in tokamak plasmas, Control Eng. Pract. 24 (2014) 15–24, http://dx.doi.org/10. 1016/j.conengprac.2013.11.009. [13] A.C. Setiadi, P.R. Brunsell, L. Frassinetti, Implementation of model predictive control for resistive wall mode stabilization on EXTRAP T2R, Plasma Phys. Control. Fusion 57 (10) (2015) 104005, http://dx.doi.org/10.1088/0741-3335/ 57/10/104005. [14] A.C. Setiadi, P.R. Brunsell, L. Frassinetti, Improved model predictive control of resistive wall modes by error field estimator in EXTRAP T2R, Plasma Phys. Control. Fusion 58 (12) (2016) 124002, http://dx.doi.org/10.1088/0741-3335/ 58/12/124002. [15] E.J. Strait, Magnetic control of magnetohydrodynamic instabilities in tokamaks, Phys. Plasmas 22 (2) (2015) 021803, http://dx.doi.org/10.1063/1. 4902126. [16] Y.Q. Liu, et al., Feedback stabilization of nonaxisymmetric resistive wall modes in tokamaks. I. Electromagnetic model, Phys. Plasmas 7 (2000) 3681, http://dx.doi.org/10.1063/1.1287744. [17] A. Portone, et al., Linearly perturbed MHD equilibria and 3D eddy current coupling via the control surface method, Plasma Phys. Control. Fusion 50 (8) (2008) 085004, http://dx.doi.org/10.1088/0741-3335/50/8/085004. [18] F. Villone, et al., Effects of three-dimensional electromagnetic structures on resistive-wall-mode stability of reversed field pinches, Phys. Rev. Lett. 100 (25) (2008) 255005, http://dx.doi.org/10.1103/PhysRevLett.100.255005. [19] J. Bialek, A.H. Boozer, M.E. Mauel, G.A. Navratil, Modeling of active control of external magnetohydrodynamic instabilities, Phys. Plasmas 8 (5) (2001) 2170, http://dx.doi.org/10.1063/1.1362532.

A.C. Setiadi et al. / Fusion Engineering and Design 121 (2017) 245–255 [20] A.C. Setiadi, et al., Reduced order modelling of resistive wall modes in EXTRAP T2R reversed-field pinch, in: 43rd EPS Conf. Plasma Phys., Brussels, p. P2.047, 2016. [21] K.E.J. Olofsson, P.R. Brunsell, J.R. Drake, Measurements of the vacuum-plasma response in EXTRAP T2R using generic closed-loop subspace system identification, Fusion Eng. Des. 87 (12) (2012) 1926–1929, http://dx.doi.org/ 10.1016/j.fusengdes.2012.04.001. [22] K.E.J. Olofsson, P.R. Brunsell, J.R. Drake, Experimental modal analysis of resistive wall toroidal pinch plasma dynamics, Nucl. Fusion 53 (7) (2013) 072003, http://dx.doi.org/10.1088/0029-5515/53/7/072003. [23] K.E.J. Olofsson, et al., Subspace identification analysis of RFX and T2R reversed-field pinches, Control Eng. Pract. 21 (7) (2013) 917–929, http://dx. doi.org/10.1016/j.conengprac.2013.03.004. [24] A.H. Glasser, M.S. Chance, Determination of free boundary ideal MHD stability with DCON and VACUUM, APS Meeting Abstracts (1997). [25] P. Zanca, et al., Advanced feedback control of magnetohydrodynamic instabilities: comparison of compensation techniques for radial sensors, Plasma Phys. Control. Fusion 54 (12) (2012) 124018, http://dx.doi.org/10. 1088/0741-3335/54/12/124018. [26] P.R. Brunsell, H. Bergsåker, M. Cecconello, J.R. Drake, R.M. Gravestijn, A. Hedqvist, J.-A. Malmberg, Initial results from the rebuilt EXTRAP T2R RFP device, Plasma Phys. Control. Fusion 43 (2001) 1457–1470, http://dx.doi.org/ 10.1088/0741-3335/43/11/303. [27] V.D. Pustovitov, et al., Decoupling in the problem of tokamak plasma response to asymmetric magnetic perturbations, Plasma Phys. Control. Fusion 50 (10) (2008) 105001, http://dx.doi.org/10.1088/0741-3335/50/10/105001. [28] P.R. Brunsell, J.-A. Malmberg, D. Yadikin, M. Cecconello, Resistive wall modes in the EXTRAP T2R reversed-field pinch, Phys. Plasmas 10 (10) (2003) 3823, http://dx.doi.org/10.1063/1.1604775. [29] Z.R. Wang, et al., Comparison between cylindrical model and experimental observation on the study of resistive wall mode in reversed field pinch plasmas, Phys. Plasmas 17 (5) (2010) 052501, http://dx.doi.org/10.1063/1. 3389229.

255

[30] R. Albanese, G. Rubinacci, Finite element methods for the solution of 3D eddy current problems, Adv. Imaging Electron Phys. 102 (1997) 1–86, http://dx.doi. org/10.1016/S1076-5670(08)70121-6. [31] F. Villone, et al., GPU-accelerated analysis of vertical instabilities in ITER including three-dimensional volumetric conducting structures, Plasma Phys. Control. Fusion 54 (8) (2012) 085003, http://dx.doi.org/10.1088/0741-3335/ 54/8/085003. [32] K.E.J. Olofsson, et al., Synthesis and operation of an FFT-decoupled fixed-order reversed-field pinch plasma control system based on identification data, Plasma Phys. Control. Fusion 52 (10) (2010) 104005, http://dx.doi.org/10. 1088/0741-3335/52/10/104005. [33] K. Glover, All optimal Hankel-norm approximations of linear multivariable systems and their L∞ -error bounds, Int. J. Control 39 (6) (1984) 1115–1193, http://dx.doi.org/10.1080/00207178408933239. [34] A. Chiuso, The role of vector autoregressive modeling in predictor-based subspace identification, Automatica 43 (6) (2007) 1034–1048, http://dx.doi. org/10.1016/j.automatica.2006.12.009. [35] J.W. Van Wingerden, The asymptotic variance of the PBSID opt algorithm, in: IFAC Proc. Vol., Vol. 16, Elsevier, 2012, pp. 1167–1172, http://dx.doi.org/10. 3182/20120711-3-BE-2027.00228. [36] S. Skogestad, I. Postlethwaite, Multivariable Feedback Control: Analysis and Design, John Wiley, 2005. [37] A.C. Setiadi, P.R. Brunsell, L. Frassinetti, Design and operation of fast model predictive controller for stabilization of magnetohydrodynamic modes in a fusion device, Conf. Decis. Control (2015) 7347–7352, http://dx.doi.org/10. 1109/CDC.2015.7403379. [38] Z. Wang, S. Guo, Physical understanding of the instability spectrum and the feedback control of resistive wall modes in reversed field pinch, Nucl. Fusion 51 (5) (2011) 053004, http://dx.doi.org/10.1088/0029-5515/51/5/053004.