Closed-loop MIMO ARX estimation of concurrent external plasma response eigenmodes in magnetic confinement fusion Erik Olofsson1 , Cristian R. Rojas2 , H˚akan Hjalmarsson2, Per Brunsell1 and James R. Drake1 Abstract— Detailed experimental MIMO models of plasma stability behaviour are becoming increasingly important in magnetic confinement fusion (MCF) energy research as an assortment of magnetohydrodynamic (MHD) instabilities develop when fusion performance is pushed. Some of these problems could perhaps be handled by magnetic feedback. We here show a practical method for experimental closed-loop multiinput multi-output (MIMO) characterisation of the macroscopic stability of toroidal MCF devices. It is demonstrated, by application to the MCF experiment EXTRAP T2R, that MHD eigenmodes can be detected using the workhorse MIMO autoregressive exogeneous (ARX) model structure. This is achieved by a simple “dry-wet” two-stage experiment sequence designed to utilise most ARX parameters towards plasma-specific information. Plausibly, the presented methodology could significantly improve highly-desired magnetic feedback accuracy in MCF.
I. INTRODUCTION Forecasts on the ever-increasing energy demand of future earthwide population pose a serious challenge to engineering sciences. Add to this the stepped-up efforts to battle global warming and it follows that the sustainable, virtually zerocarbon footprint, low-radiation, and bulk power output technologies endeavored in thermonuclear fusion reactor research [1], [2] can be well worth both scientific commitment and financial investment. Currently, the most promising concept, known as the tokamak, a russian invention from the 1960s, belong to the major subfield known as magnetic confinement fusion (MCF). MCF strives to negotiate an extremely hot (∼ 108 K) ionized gas to stay put in a torus-like volume by strong (∼ 1 T) and cleverly designed magnetic fields. Not surprisingly, the recently commissioned multi-billion-dollar project ITER (international thermonuclear experimental reactor) is a tokamak. Nominally it should output a handsome 500 MW in transient bursts by fusing two particular hydrogen isotope nuclei: deuterium and tritium. This reaction, abbreviated DT, essentially puts less demand on the plasma density (n) times confinement time (τE ) product (nτE ) than other considered fusion reactions. Deuterium can be extracted from seawater, and tritium can be manufactured by using spare energetic neutrons (emanating from the plasma itself, recall that particles without charge are not deflected by magnetic fields) in a breeding reaction with lithium. It turns out that these tokamaks are increasingly susceptible to magnetohydrodynamic (MHD) instabilities [3], [4], [5] as e.g. the plasma pressure is increased in them. Higher 1 Fusion Plasma Physics (Association EURATOM-VR) School of Electrical Engineering (EES), Royal Institute of Technology (KTH) 2 ACCESS Linnaeus Center KTH/EES, Stockholm, Sweden Corresponding author email:
[email protected]
pressure means more fusion power. A particular dangerous instability is the resistive-wall mode (RWM) [6] that “grows out” of the containing vessel on the eddy current decay timescale of the surrounding passive structures. If unstopped, the RWM will seriously limit reactor performance. This work mainly considers the multiple RWMs residing in an MCF device known as the reversed-field pinch (RFP). In RFPs it has been shown that RWMs can be stabilised by magnetic feedback control using a large array of saddle coils that generate counteracting magnetic fields [7]. We here present a methodology to automate simultaneous detection of experimental MHD RWM-type eigenmodes (stable and unstable) in RFPs [8] from dither-injected feedback stabilised plasma confinement operation data [9]. We use a pragmatic MIMO formulation, the prediction error method (PEM) [10], and numerically tractable and robust computations. Not only can the results conveyed here be useful for tuning, adapting and retooling RWM feedback. Tokamak researchers are currently investigating various approaches to produce resonant magnetic perturbations (RMPs) [11] close to the boundary of the plasma, in order to mitigate and/or trigger the edge-localized mode (ELM). ELMs are indeed a major concern for ITER. More precise in situ plasma response identification could assist in these studies by providing statistically linearised non-trivial nonaxisymmetric magnetic field dynamics in an environment of passive structures external to the vacuum-plasma system. Previous studies on automated experimental identification of MHD external eigenmodes in RFPs [12], [13] have depended on quite strong assumptions to achieve single-input single-output (SISO) programming. MCF device-specific algebra have also been required to construct parameterisations for a Kalman-based PEM. This study, in contrast, is generalised so that almost no assumptions at all on the geometry of the actual MCF device are neccessary. The methods are easy to implement in any capable linear algebra software [14]; in practice a highly important property since MCF experiments typically have an upmarket number of magnetic field pickup and actuation coils (the present study MIMO-analyses 64outputs and 64-inputs). Time-domain MIMO methods have been employed in e.g. structural mechanics [15] to identify vibrational modes. Higher-order ARXs (in conjunction with model reduction) are often a practical alternative in industrial MIMO applications [16]. The presentation is marshalled as follows. The RFP plant EXTRAP T2R is introduced in section II. Section III recalls some prerequisites from system identification and subsequently describes how the signal processing is applied. Next,
section IV describes how the estimated systems can be interpreted and visualised, and sections V and VI respectively document experimental-data results from T2R and provide a concluding commentary. II. MCF PLANT: EXTRAP T2R RFP
happens, a slightly modified ∼ 1 kW power audio amplifier), and s the Laplace argument. w(t) ∈ Rm is a carefully designed injected vector dither signal [13]. However, any switching pseudorandom sequence, binary or multilevel, would quite plausibly also work [12]. Real-time available signals at T2R are digitised at a rate of 10 kHz. Hence the sample time is Ts = 1 × 10−4 s = 100 µs. From hereon, we will chiefly consider discrete signals and systems; k will denote the integer time-index aligned with some real time t. A unit increment k + 1 implies a real time increment t + Ts . The discrete incarnation of (1) will have a delay composed of computational lags1 and finite bandwidth effects from Ci (s). T2R is proudly MIMO: n = m = 64. III. TWO-STAGE MIMO ARX SIGNAL PROCESSING A. Basics of multivariable autoregressive estimation Consider the MIMO ARX(na , nb ) structure
Fig. 1. Cut-away drawing of EXTRAP T2R RFP; showing vacuum vessel (grey), two slotted layers of thin resistive shells (yellow-orange), active saddle coils (red) and sensor saddle coils (blue).
A. The basics: RFPs and T2R EXTRAP T2R, figure 1, is a RFP of moderate size. Its major and minor radii are respectively R = 1.24 m and a = 0.183 m. The aspect ratio R/a = 6.8 is relatively high for an MCF device. RFPs have a poloidal magnetic field strength of the same order as the toroidal magnetic field. For tokamaks, the toroidal magnetic field is much stronger than the poloidal. The plasma current in T2R is around 100 kA. A large tokamak sports a few MAs. RFP plasmas self-organise in a sense of minimised energy subject to a particular global topological constraint [17]. A problem with the RFP is the short energy confinement time τE due to an internally stochastic magnetic field. T2R is well outfitted for investigations of feedback control of plasma MHD; as alluded by figure 1. RFPs and tokamaks have many MHD phenomena in common; such as the RWM. B. Closed-loop stabilisation with injected vector dithering Plant data are collected from closed-loop operation with injected dithering [12], [13]. Let u(t) ∈ Rm be a vector of currents in the actuator coils and y(t) ∈ Rn a vector of timeintegrated induced voltages from the sensor coils. See figure 1. Thus, as follows from Faraday’s law, the ith component of y(t) is a measure of the magnetic flux linking the ith coil at time t. The components of u(t) and y(t) are sorted such that equal indices correspond to coinciding coils. Let (·)i denote the ith component of a vector, then we can model, in continuous time, the decentralised feedback employed for our experiments by (u)i = Ci (s) {−Fi (s) (y)i + (w)i }
(1)
where Fi (s) is a proportional-integral-derivative (PID) system, Ci (s) the dynamics of the ith active coil driver (as it
˜ ˜ A(q)y(k) = B(q)u(k) + e(k) (2) Pna −1 −j−1 ˜ ˜ , and B(q) = where A(q) = I − j=0 Aj q Pnb −1 −j−1 n×n n×m . Aj ∈ R , Bj ∈ R are realj=0 Bj q valued matrices, and q −1 the unit delay operator such that q −1 s(k) = s(k − 1) for arbitrary signal vector s. System output y(k) ∈ Rn , input u(k) ∈ Rm , and Gaussian noise e(k) ∈ Rn . For estimation of (2) we define the linear predictor ˆ (k + 1 |k ) = y
nX a −1
Aj y(k − j) +
nX b −1
Bj u(k − j)
(3)
j=0
j=0
the prediction error ˆ(k) = y(k) − y ˆ (k |k − 1 ) e
(4)
and the prediction error covariance matrix approximation N −1
Reˆeˆ
1 X ˆ(k)ˆ e eT (k) N k=0 T ˆ(k)ˆ ≈ E e e (k)
=
where E [·] is the expectation operator. Introduce the stacked data vector y(k) y(k − 1) .. . y(k − na + 1) z(k) = u(k) u(k − 1) .. . u(k − nb + 1)
(5)
(6)
1 Analog/digital- and digital/analog conversions plus central processor time.
and the auxilliary covariance approximations in the spirit of (5) Ryy Ryz Rzz
= = =
N −1 1 X y(k)yT (k) N
(7a)
1 N
y(k + 1)zT (k)
(7b)
z(k)zT (k)
(7c)
1 N
k=0 N −1 X k=0 N −1 X
C. State-space stacking
k=0
It is easily shown that (5) can be written T + GRzz GT Reˆˆe = Ryy − Ryz GT − GRyz
(8)
(9)
belongs to Rn×(nna +mnb ) . Furthermore, it is straightforward to verify that the covariance matrix norm V (G) = tr (Reˆeˆ )
(10)
ˆ that with tr (·) denoting matrix trace, has a minimizer G solves the normal equations ∂gr,c V = 0 for all r and c, with gr,c being matrix elements of G and ∂(·) the partial ˆ = derivative operator. In summary, our ARX estimate G arg minG V (G) solves the equation GRzz = Ryz
(11)
which represents a sparse linear system with n(nna + mnb ) equations and the same number of unknowns2. B. A dry-wet system identification recipe 1) Signal preconditioning: All signal data are first scaled such that the amplitudes of both y and u are of unity order for better numerical condition of the covariance approximations. In this application we have a particular low-frequency drift that is not of interest for the phenomena we wish to detect. Neither are fluctuations of too high frequency. For this reason, a digital butterworth type bandpass prefilter [19] is applied to all components of y and u independently. 2) Dry system: passive structures and vacuum: Denote data recorded from the vacuum system by D0 = {y0 (k), u0 (k)}k=0...N −1 . We use this vector input-output ˆ 0 from (11). data to directly estimate the ARX predictor G ˆ 0 z0 (k), where z0 (k) is ˆ 0 (k + 1 |k ) = G Thus, y0 (k + 1) ≈ y formed by dataset D0 according to (6). 3) Wet system: passive structures and plasma: Denote data recorded from the plasma system by D1 = {y1 (k), u1 (k)}k=0...N −1 . In order to not refit precious ARX degrees of freedom to known passive vacuum behaviour, we ˆ 0 . Define a residual output y1−0 (k) filter out the effect of G by subtracting the dry prediction. ˆ 0 z1 (k − 1) y1−0 (k) = y1 (k) − G
The MIMO ARX structure (2), with e(k) = 0, can naively be cast in the state-space form x(k + 1) =
where the bundled matrix coefficient G = (A0 . . . Ana −1 B0 . . . Bnb −1 )
The dataset D1−0 = {y1−0 (k), u1 (k)}k=0...N −1 is then supposedly consisting mainly of new information specific to the wet system. Dataset D1−0 then define the residual system ˆ 1−0 . Note that the summed output of systems G ˆ0, estimate G ˆ G1−0 , in parallel is analogous to what we should get for ˆ 1 based on D1 . However by separating we have a system G ˆ (·) are equally sized) system of double dimension (if all G where one half solely focuses on the plasma contribution to the output; y1 (k) = y0 (k) + y1−0 (k).
(12)
2 Numerical solution of (11) based on QR-decomposition [14] is a snug matlab [18] one-liner G=(Rzz’\Ryz’)’; or, even more deceptively: G=Ryz/Rzz;
y(k)
=
ˆ ˆ A(G)x(k) + B(G)u(k) ˆ C(G)x(k)
(13a) (13b)
where the stapled state
y(k) y(k − 1) .. .
y(k − na + 1) x(k) = u(k − 1) u(k − 2) .. . u(k − nb + 1)
(14)
is of dimension Rnna +(m−1)nb and the calligraphic metaˆ of (9) matrices in (13) are block-built from estimates G A0 A1 . . . Ana −1 B1 . . . Bnb −1 I 0 0 0 A= 0 0 0 0 0 0 I 0 (15) B0 0 B= (16) 0 I 0 C = In×n 0n×(na (n−1)+(nb −1)m) (17) In matrices (15) and (16) the matrix partitions are supposed to be inflated by the corresponding sub-blocks in a selfexplanatory way. Matrix (17) is fully explicit, and turns out ˆ in this representation. to be independent of G However (15), (16) and (17) in general do not constitute a minimal realisation [20] of the matrix polynomial system (2). We would like to have the smallest-sized system matrices that can represent (2). It is well-known that this is achieved by removal of unobservable and uncontrollable states [20]. From now on, assume this has been done for (13) and the resulting minimal state dimension is p < nna + (m − 1)nb ; A ∈ Rp×p , B ∈ Rp×m and C ∈ Rn×p . Once reduced, we have structural knowledge of neither the matrices nor state in (13).
In section IV, we will primarily be interested in the eigenvectors and eigenvalues of A of (13), and how they are projected by C.
The continuous-time growth-rate of the jth eigenmode is γj
1 ln λj Ts ι Re λj 1 ln |λj | + arccos Ts Ts |λj |
=
D. A word on order selection and model reduction
=
Although the size of the estimated systems are quite hefty in the present application we have not considered pursuit of more parsimonious representations. The focus of the study is to detect, if at all, MHD physics eigenmodes through MIMO ARX analysis. Reduction of overmodeling is a topic for future development. It can furthermore be imagined that state obfuscation through e.g. balanced truncation and/or some other successfully applied identification method might have been utilised to produce the final estimated MIMO system (13). IV. GEOMETRIC INTERPRETATION AND STATE DEOBFUSCATION A. Basic problem Reconsider the state space model (13). Assume it models the plant time-series data accurately. Assume we have no clue what x might represent. Yet, in our case, we can relate a natural mode of (13) to MHD physics by the spatial pattern it exhibits on the output vector y. Two ways to “decrypt” the states of (13) will be considered. The first approach is most fundamental: sorting, visualisation and labeling of eigenvectors to A in terms of input-output amplification, stability, and output geometry. The second method is rightout procrustean: we seek a state transformation that renders the state vector approximately intelligible, in terms of output geometry.
with Ts the sample time. Stability of eigenmode j is then equivalent to Re γj < 0. C. Output-projection of eigenvectors Recall system (13). Let vj ∈ Cp be an eigenvector to A ∈ Rp×p with eigenvalue λj ∈ Cp . It is possible to label this eigenmode j by studying its projection on the output array y. From (18) we get yj (k) = ∼
Cxj (k) C {aj cos(kωj ) − bj sin(kωj )}
¯j = y
xj (k) = Cj rjk {aj cos(kωj + Dj ) − bj sin(kωj + Dj )} (18) q where rj = α2j + βj2 , λj = αj + ιβj , ι2 = −1, vj = aj + ιbj , ωj = arccos(αj /rj ), and Cj , Dj some constants. A real-valued λj represent the solution xj (k) = Cj λkj vj
(19)
A mode, whether (18) or (19) is unstable if |λj | > 1, and oscillating whenever Im λj = βj 6= 0 (18). It is expected, in this application, that the unstable and nearly marginally stable eigenmodes are mainly of type (19); in practice this could imply βj /αj 1 for type (18).
(21)
where we have set Dj = 0 and discarded the (possibly unstable) amplification Cj rjk . A useful strategy to summarise the meaning of the large identified black-box model (13) in terms of physical eigenmodes is to sample periodograms of (21) spatially. This yields Fourier coefficients of the periodic geometric pattern the eigenmode represents. Fourier modes are commonly used in MHD research. For T2R we can produce a coloured map compactly showing the eigenmode structure by splashing out sampled periodograms horizontally on a vertical position given by the respective eigenvalue magnitudes. We also normalise the eigenvector in a way that is independent of the particular state-space realisation. Define the scaled output
B. Real and complex-conjugate natural modes An autonomous discrete-time system [21] x(k + 1) = Ax(k) where A ∈ Rp×p has the explicit solution x(k) = P k j cj λj vj where (λj , vj ) is an eigenvalue-eigenvector pair to A, and cj some constant. λj can either be real-valued or otherwise can be paired with another complex-conjugate eigenvalue λi = λ?j . It will also hold: vi = vj? , i.e. the corresponding eigenvectors are componentwise complexconjugate. A single-multiplicity conjugate pair (λj , λ?j ) represents a general real-valued solution
(20)
1 yj ||yj ||2
(22)
where yj = Cvj , and ||·||2 the standard Euclidean norm. Welch-like periodograms based on (22) are then manufactured by 2 1 X ˆ l Pˆj (n) = (23a) Yj (n) Nl l
Yˆjl (n) =
N −1 X i=0
e−2πιin/N
n
¯ jl y
¯ jl −ι y i
¯ jl y
=
T −1/2 Cvjl vjl C T Cvjl
vjl
=
aj cos(ϕl ) − bj sin(ϕl )
o (23b) i+N
(23c) (23d)
where N = 32, (·)i denotes the ith component of a vector, ϕl is a random real scalar uniformly distributed on [0, 2π], and vj = aj + ιbj the jth eigenvector to A. Equations (23a) and (23b) are the power spectrum periodogram and the discrete Fourier transform (DFT) respectively, along the spatial domain (as defined by an enumeration of the sensor coils). Equation (23c) originate from (22). The average over l in (23a) is thus in the sense of randomly sampling phases of (21). One-dimensional DFT-based (23) is appropiate for T2R. A tokamak would need a generalised form.
D. Procrustes-type state transformation
CT = Cq
(24)
1.02
1
0.98
0.96
|λ|
Yet again reconsider (13). Suppose we would like to synthesise a state transformation that results in a state with, at least approximately, physical meaning. The geometric idea, of subsection IV-C, of using the sensor array projection to label MHD eigenmodes can be recycled. For MCF experiments it is possible to construct an appropiately dimensioned measurement matrix Cq that maps a stacked vector q of ordered coefficients in some spectral basis set (typically Fourier) to the output y. Thus y = Cx = Cq q, and the meaning of each component of q is known. We now seek a nonsingular matrix T ∈ Rp×p that maps T x = q. We would like to have the property
0.94
at least approximately. Indeed, (24) has a particular wellknown easy-to-compute approximate well-conditioned solution from the orthonormal procrustes [22] problem
0.92
T1 = arg min ||CT − Cq ||F s.t. T T T = I
0.9
T
(25)
with ||·||F betokening the Frobenius norm [23]. The unique solution is T1 = U V T where U ΣV T is the singular value decomposition of C T Cq . On the other hand (24) is solved exactly [24] by T 2 = C † Cq + I − C † C Z (26) †
where Z is an arbitrary matrix and (·) the Moore-Penrose pseudoinverse [23]. A simple method for “completing” (26) is to put Z = T1 . The applicability of an orthonormal transformation presupposes that the matrices are scaled properly. For T2R Cq is basically a nonsquare inverse DFT-type matrix [12], [13]. It is reasonable to assume that Cq has full rank regardless of the particular MCF device it is assembled for since otherwise it would imply that the magnetic field sensors have been distributed in a grossly untoward fashion. V. EXPERIMENTAL RESULTS Data from dithered T2R experiments, processed by ARX na = nb = 3 after signal pretreatment using a 25 − 3000 Hz bandpass 4th order digital Butterworth filter, are shown in figures 2 and 3. These periodograms seem stable in the sense that higher order ARXs look similar. Selection of bandpass filter frequencies with some care also coax lowerorder ARXs to produce comparable images. Figure 2 shows ˆ 1−0 . Figure a plasma stability-periodogram rendered from G ˆ0. 3 depicts a vacuum stability-periodogram rendered from G The sort-of insular unstable peak in figure 2 correspond to the expected spectral peak at n = −11 for T2R [7], in accordance with MHD/RWM stability theory [8]. A sample of this eigenmode in real space is provided by figure 4, showing a predominantly n = −11 undulation. Its continuous-time growth rate is ∼ 89 s−1 ; of the order of the inverse shell time constant τw ≈ 10 ms of T2R. The overall impression from figure 2 is a slightly less unstable T2R than anticipated. Vacuum-mode magnetic field diffusion periodogram, figure 3, shows decaying time-constants with increasing |n|, for the
0.88 −15
−10
−5
0
n
5
10
15
Fig. 2. T2R RFP external plasma response eigenmodes from generic MIMO ˆ 1−0 . n is the toroidal Fourier MHD mode number ARX linear estimation; G [17]. The poloidal mode number is, implicitly, m = 1.
“lower” branch. The broad-band structure above is probably a consequence of spatially localised field penetration due to slots in the resistive shell. Automatic characterisation of such effects is a strong sales argument for the present method; since not including them in plasma reponse estimation could possibly yield a geometric bias.
VI. CONCLUSIONS It was shown, that a full-tilt MIMO ARX direct closedloop system identification approach, can detect eigenmodes, that are related to MHD stability theory of the plasma response for the RFP MCF machine T2R. This is a promising and practical result. MIMO ARXs are easy to use and numerically reliable. It remains, of course, to properly develop quality control of the estimates obtained from the methodology followed here. Experimental nonaxisymmetric plasma response dynamics and stability in complex geometries (slotted shells, ports, holes and openings for diagnostics and microwave heating antennae) is an important, and perhaps unavoidable, topic for any aspiring high-performance thermonuclear fusion reactor. It is imaginable that these MIMO ideas could be developed into various magnetic control and monitoring applications in future reactors, such as ITER.
R EFERENCES
1.02
1
0.98
|λ|
0.96
0.94
0.92
0.9
0.88 −15
−10
−5
0
n
5
10
15
Fig. 3. T2R RFP passive vacuum shell magnetic field diffusion eigenmodes; ˆ0. G
3
2
θ
1
0
−1
−2
−3 −3
−2
−1
0
1
2
3
φ
Fig. 4. A real-space sample of the most unstable perturbation pattern of ˆ 1−0 smoothed over the poloidal θ and toroidal φ angles of the plasma G surface.
[1] J. Li, J. Zhang, and X. Duan, “Magnetic fusion development for global warming suppression,” Nuclear Fusion, vol. 50, no. 1, p. 014005 (5pp), 2010. [Online]. Available: http://stacks.iop.org/0029-5515/50/014005 [2] M. Shimada et al., “Progress in the ITER physics basis, chapter 1: Overview and summary,” Nuclear Fusion, vol. 47, no. 6, pp. S1–S17, 2007. [Online]. Available: http://stacks.iop.org/0029-5515/47/S1 [3] G. Bateman, MHD instabilitites. MIT Press, 1978. [4] W. M. Stacey, Fusion plasma physics. Wiley-WCH, 2005. [5] M. Walker et al., “Emerging applications in tokamak plasma control: Control solutions for next generation tokamaks,” IEEE Control System Magazine, vol. 26, pp. 35–63, 2006. [6] Y. Liu, A. Bondeson, M. Chu, J.-Y. Favez, Y. Gribov, M. Gryaznevich, T. Hender, D. Howell, R. L. Haye, J. Lister, P. de Vries, and E. J. Contributors, “Feedback and rotational stabilization of resistive wall modes in ITER,” Nuclear Fusion, vol. 45, no. 9, pp. 1131–1139, 2005. [Online]. Available: http://stacks.iop.org/0029-5515/45/1131 [7] P. R. Brunsell et al., “Feedback stabilization of multiple resistive wall modes,” Physical Review Letters, vol. 93, no. 22, p. 225001, 2004. [Online]. Available: http://link.aps.org/abstract/PRL/v93/e225001 [8] R. Paccagnella, D. Gregoratto, and A. Bondeson, “Feedback control of resistive wall modes in reversed field pinches,” Nuclear Fusion, vol. 42, no. 9, pp. 1102–1109, 2002. [Online]. Available: http://stacks.iop.org/0029-5515/42/1102 [9] U. Forssell and L. Ljung, “Closed-loop identification revisited,” Automatica, vol. 35, pp. 1215–1241, july 1999. [Online]. Available: http://dx.doi.org/10.1016/S0005-1098(99)00022-9 [10] L. Ljung, System Identification - theory for the user, 2nd ed. PrenticeHall, 1999. [11] T. Evans et al., “Edge stability and transport control with resonant magnetic perturbations in collisionless tokamak plasmas,” Nature Physics, vol. 2, pp. 419–423, 2006. [Online]. Available: http://dx.doi.org/10.1038/nphys312 [12] E. Olofsson, P. Brunsell, and J. Drake, “Closed-loop direct parametric identification of magnetohydrodynamic normal modes spectra in EXTRAP-T2R reversed-field pinch,” Proceedings of the IEEE conference on Control Applications (CCA) & Intelligent Control (ISIC) 2009, pp. 1449–1454, July 2009. [13] E. Olofsson, H. Hjalmarsson, C. R. Rojas, P. Brunsell, and J. R. Drake, “Vector dither experiment design and direct parametric identification of reversed-field pinch normal modes,” Proceedings of the 48th IEEE Conference on Decision and Control CDC/CCC 2009, pp. 1348 –1353, dec. 2009. [14] L. N. Trefethen and D. Bau, Numerical Linear Algebra. Society for Industrial and Applied Mathematics (SIAM), 1997. [15] A. Florakis, S. Fassois, and F. Hemez, “Critical assessment of MIMO structural identification methods: ARX and ARMAX schemes,” Proceedings of the 16th international modal analysis conference, 1998. [16] M. Leskens, L. V. Kessel, and P. V. den Hof, “MIMO closed-loop identification of an MSW incinerator,” Control engineering practice, no. 10, pp. 315–326, 2002. [17] J. Freidberg, Ideal Magnetohydrodynamics. Plenum Press, 1987. [18] “Mathworks Inc.” http://www.mathworks.com/matlab, March 2010. [19] J. F. MacGregor and D. T. Fogal, “Closed-loop identification: the role of the noise model and prefilters,” Journal of Process Control, vol. 5, no. 3, pp. 163 – 171, 1995. [20] T. Kailath, Linear Systems. Prentice-Hall, 1980. [21] O. Galor, Discrete dynamical systems. Springer Verlag, 2007. [22] J. Gower and G. Dijksterhuis, Procrustes problems. Oxford university press, 2004. [23] R. A. Horn and C. R. Johnson, Matrix Analysis. Cambridge University Press, 1990. [24] R. Skelton, T. Iwasaki, and K. Grigoriadis, A unified algebraic approach to linear control design. Taylor & Francis, 1998.
