Array magnetics modal analysis for the DIII-D tokamak based on localized time-series modelling

The feature-extraction process can be subdivided into three main steps: (i) time-series modelling, (ii) combined frequency decomposition and mode-selection and (iii) shape estimation. Step (iii) is only necessary if step (i) is computed on a reduced-dimension version of the full set of magnetic probes. If the full set of magnetic probes is used in step (i), then step (ii) would also achieve (iii) as a direct by-product. The reason for using a dimension reduction of the signal is computational complexity as will become clear in section 2.1.3 below. The process is illustrated in the flowchart in figure 2. In the next three sections steps (i) through (iii) will be detailed.

Zoom In Zoom Out Reset image size

2.1.1. Time-series modelling.

Let $\mathcal{Y}_N=\{{\bit y}(k)\}_{k=1}^{N}$ be a block of magnetic probe data where ${\bit y}(k)\in\mathbb{R}^{m\times 1}$ is an m-dimensional vector sample indexed by discrete time k. The components of y(k) are distinct magnetic probe samples at various toroidal and poloidal locations, at time k. For this block of time-series data a discrete-time state-space model of the following form

$$\begin{equation} {\bit x}(k+1) = A_r {\bit x}(k) + K_r {\bit e}(k) \end{equation} \tag{ 1a }$$

$$\begin{equation} {\bit y}(k) = C_r{\bit x}(k) + {\bit e}(k) \end{equation} \tag{ 1b }$$

is to be estimated. Signal model (1a) and (1b) has a state ${\bit x}(k)\in\mathbb{R}^{r\times 1}$ and the estimation procedure assumes that e(k) is a zero-mean sequence of white-noise random variables with some unknown covariance. Note that the time-series state-space dimension is r. In what follows, the method parameters denoted by f, p are respectively called the future and past horizons. The estimation starts by forming a data matrix D_fp with a certain shift structure from the block $\mathcal{Y}_N$:

$$\begin{eqnarray} \fl D_{fp}=\left[\begin{array}{@{}cccc@{}} {\bit y}(1) & {\bit y}(2) & \ldots & {\bit y}(N-p-f+1) \\ {\bit y}(2) & {\bit y}(3) & \ldots & {\bit y}(N-p-f+2) \\ \vdots & \vdots & \ddots & \vdots \\ {\bit y}(p+f) & {\bit y}(p+f+1) & \ldots & {\bit y}(N) \end{array}\right].\nonumber\\ \end{eqnarray} \tag{ 2 }$$

Now the first mp rows of D_fp are denoted by D_p and the last mf rows are denoted by D_f. D_fp is thus partitioned as

$$\begin{equation} D_{fp}=\left[\begin{array}{@{}c@{}} D_p \\ D_f \end{array}\right] \in \mathbb{R}^{(mp+mf) \times (N-p-f+1)} \end{equation} \tag{ 3 }$$

and the idea is to find a linear predictor of D_f given D_p. A basic SSI method goes as follows [24].

$$\begin{equation} R_{fp} = D_f D_p^T \in \mathbb{R}^{mf \times mp} \end{equation} \tag{ 4a }$$

$$\begin{equation} R_{fp} = U\Sigma V^T \end{equation} \tag{ 4b }$$

$$\begin{equation} X_r = \Sigma_r^{1/2} V_r^T D_p. \end{equation} \tag{ 4c }$$

Here Σ is a diagonal matrix, Σ_r = Σ(1 : r, 1 : r), and V_r = V (:, 1 : r). At this point the columns of X_r could be used as the state sequence x(k) in equations (1a) and (1b). This means linear regression can be used to obtain the matrices (A_r, K_r, C_r). But if only the modes and array shapes of equations (1a) and (1b) are of interest, only (A_r, C_r) are required. An alternative method is therefore to use the same SVD as in equation (4b) but instead calculating

$$\begin{equation} O_r = U_r\Sigma_r^{1/2} \end{equation} \tag{ 5 }$$

where U_r = U (:, 1 : r). Here O_r can be interpreted as the 'extended observability matrix' [24, 25] which has a particular shift structure that allows direct extraction of (A_r, C_r):

$$\begin{equation} D_1 A_r = D_2 \end{equation} \tag{ 6a }$$

$$\begin{equation} \widehat{C_r} = O_r(1:m,:) \end{equation} \tag{ 6b }$$

$$\begin{equation} \widehat{A_r} = \left(D_1^{\rm T}D_1\right)^{-1}D_1^{\rm T}D_2 \end{equation} \tag{ 6c }$$

where

$$\begin{equation} D_1 = O_r(1:(mf-m),:) \end{equation} \tag{ 7a }$$

$$\begin{equation} D_2 = O_r((m+1):mf,:) . \end{equation} \tag{ 7b }$$

In eigspec, method (5) through (7b) is implemented as the default calculation.

2.1.2. Modal analysis and selection.

Time-series model (1a) and (1b) can be represented in modal coordinates q = W⁻¹x as

$$\begin{equation} {\bit q}(k+1) = \tilde{A}_r {\bit q}(k) + \tilde{K}_r {\bit e}(k) \end{equation} \tag{ 8a }$$

$$\begin{equation} {\bit y}(k) = \tilde{C}_r{\bit q}(k) + {\bit e}(k) \end{equation} \tag{ 8b }$$

where Ã_r is diagonalized (A_r = WÃ_rW⁻¹), $\tilde{K}_r=W^{-1}K_r$ and $\tilde{C}_r=C_rW$. The modal subsystems are collected from (i) real-valued eigenvalues and (ii) complex-conjugate pairs of eigenvalues. The eigenvalues λ_i are found on the diagonal of Ã_r.

In eigspec, the ith mode shape vector is, nominally, the complex vector found as a column of $\tilde{C}_r$, associated with a complex-conjugate mode i. It is only the actual shape vector if no random projection technique was used. See next section 2.1.3 for the details on this. Regardless of whether random projection was used or not, the next important problem to solve is to determine which linear system modes to keep and which to discard. To do so, a coherence-like quantity is introduced that evaluates the similarity of two complex-valued vectors v, w. This quantity is denoted as MAC(v, w) (acronym for modal assurance criteria) in order to connect with the modal analysis literature where its original usage is to compare finite-element computational modal shapes with experimentally obtained shapes [28, 29]. In eigspec, the MAC is used both for obtaining non-spurious features from time-series data, and for clustering and comparison of different non-spurious shape vectors as will be seen in section 2.2.

The MAC is defined as

$$\begin{equation} {\rm MAC}\left({\bit v},{\bit w}\right)=\frac{({\bit v}^{\dagger}{\bit w})({\bit w}^{\dagger}{\bit v})}{({\bit v}^{\dagger}{\bit v})({\bit w}^{\dagger}{\bit w})}=\frac{|{\bit v}^{\dagger}{\bit w}|^2}{\|{\bit v}\|_2^2\|{\bit w}\|_2^2} \end{equation} \tag{ 9 }$$

where (·)^† denotes complex-conjugate transpose and $\|{\bit v}\|_2^2={\bit v}^{\dagger}{\bit v}$. MAC(v, w) is real-valued in [0, 1]. It holds MAC(r_v exp(ιθ_v)v, r_w exp(ιθ_w)w) = MAC(v, w) (the MAC 'metric' is independent of both phase and amplitude). Modal shape similarity should imply that the MAC is close to unity, and MAC(v, v) = 1 for v ≠ 0. eigspec uses a specific procedure denoted Order-MAC to remove spurious modes in the frequency extraction process. It works as follows.

Apply SSI described in section 2.1.1 to obtain the pair time-series model matrices $(A_{r_1},C_{r_1})$ and $(A_{r_2},C_{r_2})$ with orders r₁ < r₂ (i.e. two state-space models are estimated for the same block data). Extract the modal shape vectors from both models (r₁ and r₂) as outlined above in this section 2.1.2. For each mode in the smaller model with dimension r₁, find its best match (using the MAC value [9)] shape vector in the larger model with dimension r₂. If the MAC value of the best match is higher than a threshold value close to 1, say 0.998, then add the mode of the smaller model to a shortlist of modes. If the MAC value is smaller than the threshold, ignore the mode.

The above description mentions a single threshold value to 'filter out' spurious modes. In addition, eigspec can use a second threshold value to check that the eigenvalue λ_i of the mode is nearly invariant with respect to r₁ and r₂. Specifically, the value $\max(0,1-|\lambda_i^{(r_1)}-\lambda_i^{(r_2)}|/|\lambda_i^{(r_1)}|)$ has to exceed the second threshold for a mode to be shortlisted. Order-MAC is based on observations in OMA [27, 30] that spurious modes can be strongly dependent on the estimation parameters (here the subspace dimension) but that physical modes are often insensitive to the same parameter.

2.1.3. Random projection and shape-vector estimation.

eigspec uses the SSI and MAC algorithms outlined above in a particular way to significantly speed up the calculations. It is observed that multiplying the channel data 'spatially' by a constant matrix does not alter the frequency content of a block of data. This can be exploited to reduce the channel dimension of the analysis, which is important since the computational complexity of SSI grows rapidly with the channel dimension m and the f, p horizons. In fact, the time required to numerically calculate the SVD of the mf × mp matrix R_fp in the basic SSI procedure (4b) above is proportional to m³p²(f + p) using a typical algorithm [12]. The cubic dependence on the array dimension m implies that the SVD calculation would be 8 times quicker if the array sample was 'compressed' to half its original size, say l = m/2. It turns out to be quite effective to use a random projection matrix which is constructed independently of the actual data in the block Y. The idea of using a random projection comes from recent developments in compressive sampling [31]. The Mirnov array in DIII-D has about m ≈ 40 channels (depending on how many of them are working properly and whether all of them are used in the analysis) but it has been found that using about l ≈ 10 channels (from random projection) essentially gives the same results (but is nominally about 64 times faster).

Specifically, if ${\bit y}(k)\in\mathbb{R}^{m\times 1}$ is a magnetic probe array sample at time k, and

$$\begin{equation} Y = \left[\begin{array}{@{}ccc@{}}{\bit y}(1)&\ldots&{\bit y}(N)\end{array}\right] \end{equation} \tag{ 10 }$$

is the block data then introduce a random projection matrix $\Psi\in\mathbb{R}^{l\times m}$ with l ⩽ m and form the smaller block data

$$\begin{equation} \Psi Y = Z = \left[\begin{array}{@{}ccc@{}}{\bit z}(1)&\ldots&{\bit z}(N)\end{array}\right] . \end{equation} \tag{ 11 }$$

The random projection matrix Ψ is constructed such that all its rows are orthogonal. Specifically, ΨΨ^T = I_l where I_l is the l-dimensional identity matrix. In the analysis code eigspec, the lower dimensional signal z(k) is analysed instead of the actual full-array sample y(k). This however obfuscates the shape vectors. To solve this problem, SSI with Order-MAC is used only to obtain a shortlist of frequencies $\Omega=\{\omega_j\}_{j=1}^{d \lesssim r_1/2}$ for each block of projected data Z. The frequency shortlist Ω is then used to approximately factorize the full data Y as a low-rank matrix as follows:

$$\begin{equation} Y = \sum_{i=1}^{d} {\bit d}_i{\bit s}_i^{\rm T} + E = DS^{\rm T} + E \end{equation} \tag{ 12 }$$

where d_i is the ith column of D, s_i is the ith column of S, and E is a residual matrix. The columns of S are prescribed from Ω: each ω_j contributes two columns to S: s_i(k) = cos(ω_jk), s_i+1(k) = sin(ω_jk). Finally, the columns of D are found by least-squares estimation

$$\begin{equation} \widehat{D}=(YS)(S^{\rm T}S)^{-1} \end{equation} \tag{ 13 }$$

and the columns $\widehat{{\bit d}}_i$, $\widehat{{\bit d}}_{i+1}$ are used as the shape vector for the corresponding frequency ω_j which was associated with the creation of the columns s_i and s_i+1. The use of expression (13) minimizes the Frobenius norm of the residual matrix $E=Y-\widehat{D}S^T$ in equation (12).

The final stage of a typical modal extraction is to calculate a root-mean-square (RMS) amplitude for the mode with frequency ω_j. This is done as

$$\begin{equation} {\rm RMS}_j = \sqrt{(\widehat{{\bit d}}_i^{\rm T}\widehat{{\bit d}}_i + \widehat{{\bit d}}_{i+1}^{\rm T}\widehat{{\bit d}}_{i+1})/m}. \end{equation} \tag{ 14 }$$

When the raw Mirnov data is being analysed (time derivative of the poloidal magnetic field component) the amplitude (14) need to be divided by ω_j in order to obtain the overall RMS mode amplitude for mode j.

The output from the feature extraction phase is a small set of modes for each block Y of time-series data. A full-array analysis would repeat this block analysis over the entire tokamak shot dataset by repeatedly shifting the position of the analysis window by a user-determined increment and aggregating the results. The basic visualization of the results would be to scatter the extracted frequencies versus the block-centre time and similarly for RMS amplitudes in other plots. With each extracted frequency ω_j there is also an associated shape vector defined as

$$\begin{equation} {\bit v}_j = \widehat{{\bit d}}_i + \iota \widehat{{\bit d}}_{i+1} \end{equation} \tag{ 15 }$$

using the notation introduced for equation (13). This opens up several other options for visualization, exploration and classification. Section 2.2.1 outlines a clustering technique that can group the entire set of modes found in the shot into classes based on the shape vectors (15) only. The grouping algorithm does not need to know the actual coordinates of the magnetics array probes, but its performance depends implicitly on them (since they determine how separable the modes are). Section 2.2.2 explains a non-parametric regression technique that exploits the Mirnov array probe coordinates in a particular way so that a periodic modal shape can be extrapolated without any reference to poloidal or toroidal period numbers.

2.2.1. Spectral clustering.

2.2.2. Gaussian process regression.

Shape vector (15) is regarded as a spatially discretized sample of a continuous periodic plasma eigenmode shape at a particular time. Gaussian process regression (GPR) [34, 35] turns out to be useful for estimating this continuous mode shape. GPR is a technique where the regularizing properties are removed from a direct parameterization and placed in a kernel instead. A kernel can be thought of as a description of the covariance structure of a random field. The kernel

$$\begin{equation} k({\bit x},{\bit x}^{\prime})=k(\theta,\theta^{\prime},\phi,\phi^{\prime})= k_p(\theta,\theta^{\prime})k_t(\phi,\phi^{\prime}) \end{equation} \tag{ 16 }$$

where x = (θ, φ) and x' = (θ', φ') and

$$\begin{equation} k_p(\theta,\theta^{\prime}) = \exp\left( -2\sin^2((\theta-\theta^{\prime})/2)/\sigma^2_{\theta} \right) \end{equation} \tag{ 17a }$$

$$\begin{equation} k_t(\phi,\phi^{\prime}) = \exp\left( -2\sin^2((\phi-\phi^{\prime})/2)/\sigma^2_{\phi} \right) \end{equation} \tag{ 17b }$$

can be used to model 2π × 2π-periodic patterns on (θ, φ) with two length scales σ_θ, σ_φ [34].

Let x_i = (θ_i, φ_i) be the ith 'training point' with the known target value $f_i\in\mathbb{R}$. Let x be an arbitrary coordinate where an unknown value of a continuous field is to be evaluated, given the discrete training points (i.e. the real and imaginary components of the shape vector (15)). Given the 'hyper-parameters' σ_θ, σ_φ and the noise level σ₁ the kernel (16) is then used as follows.

Construct

$$\begin{eqnarray} \fl K_{11}=\left[\begin{array}{@{}cccc@{}} k({\bit x}_1,{\bit x}_1) & k({\bit x}_1,{\bit x}_2) & \ldots & k({\bit x}_1,{\bit x}_m) \\ k({\bit x}_2,{\bit x}_1) & k({\bit x}_2,{\bit x}_2) & \ldots & k({\bit x}_2,{\bit x}_m) \\ \vdots & \vdots & \ddots & \vdots \\ k({\bit x}_m,{\bit x}_1) & k({\bit x}_m,{\bit x}_2) & \ldots & k({\bit x}_m,{\bit x}_m) \end{array}\right] + \sigma_1^2 I_m\nonumber\\ \end{eqnarray} \tag{ 18 }$$

and

$$\begin{equation} {\bit k}_{12}({\bit x})=\left[\begin{array}{@{}c@{}} k({\bit x}_1,{\bit x}) \\ k({\bit x}_2,{\bit x}) \\ \vdots \\ k({\bit x}_m,{\bit x}) \\ \end{array}\right] \end{equation} \tag{ 19 }$$

for the m training points (components of the shape vector). Set k₂₂ = k(x, x). The predicted noise-free value at the arbitrary coordinate x is then

$$\begin{equation} \widehat{f}({\bit x}) = {\bit k}_{12}({\bit x})^{\rm T} K_{11}^{-1} {\bit f} \end{equation} \tag{ 20 }$$

where f is the vector of target values f_i. The variance of the predicted value at the arbitrary coordinate x is

$$\begin{equation} \widehat{\sigma^2}({\bit x}) = k_{22} -{\bit k}_{12}({\bit x})^T K_{11}^{-1} {\bit k}_{12}({\bit x}). \end{equation} \tag{ 21 }$$

The pattern and shape estimation routines in eigspec use the expressions (20) and (21) to draw smooth estimates of the toroidal and poloidal shapes of an eigenmode. In this case the x_i's are the magnetic probe locations on the tokamak internal surface and the target values f are either the real or the imaginary parts of the shape vector. The smooth toroidal and poloidal shapes are then obtained using (20), with $2\widehat{\sigma}$ error bars from (21), by evaluation at x's along the toroidal and poloidal arrays respectively.

Importantly, there is a consistent method of selecting the set of smoothing hyper-parameters. The logarithm of the likelihood function (log-likelihood) of σ_θ, σ_φ, σ₁, given the data x_i, f_i(i = 1 ... m) is

$$\begin{eqnarray} &&\fl l(\sigma_{\theta},\sigma_{\phi},\sigma_1) = \frac{1}{2}\log\det K_{11} + \frac{1}{2}{\bit f}^{\rm T} K_{11}^{-1} {\bit f} + \frac{m}{2}\log (2\pi)\nonumber\\ \end{eqnarray} \tag{ 22 }$$

and it should be as large as possible. So direct optimization of (22) can be used to regularize the continuous field estimation for the mode shape (20).

After estimating periodic array-aligned poloidal and toroidal shape interpolation on fine grids in the sense above, a 2D pattern is extrapolated as a rank-2 matrix as follows. Let p_re, p_im be the poloidal-array real and imaginary shape estimates. Let t_re, t_im be the toroidal-array real and imaginary shape estimates. The mode pattern is constructed as

$$\begin{eqnarray} \fl P &=& {\rm Re}\left\{({\it p}_{\rm re}^{\rm T}+\iota{\it p}_{\rm im}^{\rm T}) \times ({\it t}_{\rm re}+\iota{\it t}_{\rm im}) \right\} \nonumber\\ \fl &=& {\it p}_{\rm re}^{\rm T}{\it t}_{\rm re} - {\it p}_{\rm im}^{\rm T}{\it t}_{\rm im} \end{eqnarray} \tag{ 23 }$$

and the matrix P is drawn as a contour plot. Pattern estimation (23) is referred to as an extrapolation for the following reason. GPR could be used to evaluate the field over the poloidal–toroidal plane directly. However, since the density of Mirnov db_θ/dt-array probes outside the toroidal and poloidal sub-arrays is basically zero at DIII-D (figure 5), this will not work so well.

The modal analysis performed by eigspec provides a complex-valued eigenvalue λ_i for an array mode (8a) and (8b). In principle there is information on both growth/damping rates and the frequency of a particular mode. It turns out that in typical DIII-D Mirnov array analysis, this extra information (compared to e.g. FFT-based frequency analysis) does not seem to be so easy to exploit. It is nearly always the case that a shortlisted mode has an eigenvalue magnitude |λ_i| very close to unity (this means it is nearly a stationary-amplitude oscillation). Fine details on amplitude dynamics appear to be easier to extract by shifting the block-analysis window by small increments in time instead.

It has also been observed from investigation using basic robust statistics techniques that short outlier-like bursts of Mirnov activity from edge localized modes (ELMs) do not seem to have a significant effect on the performance of the eigspec code. A sub-division of a Mirnov data block into shorter segments and a segment exclusion method based on the minimum covariance determinant (MCD) criterion [36], did not yield a significant change of the modal analysis outputs. A tentative explanation for this is that the ELM signature is comparatively broadband (both in time- and space frequencies). This makes its effect on a line-spectrum code relatively small.

2.4.1. Example comparison with cross-spectrum approach.

Approaches based on short-time discrete Fourier transforms (DFTs, FFTs) are limited in frequency resolution by the length of the signal block that is analysed. But the Mirnov array data is rarely a stationary signal, so it becomes important to have the capability to resolve frequencies from quite short segments of data in order to follow transient and time-varying MHD mode activity. An example illustrates the capability of the present method to go below the typical non-dimensional DFT resolution Δω = 2π/N, where N is the number of samples. Two rotating modes are recorded with a signal-to-noise ratio (SNR) of 20. The recording time is N = 200 samples which gives a DFT Δω = 0.0314 rad/sample. The mode frequencies are ω₁ = 0.41 and ω₂ = ω₁ + Δω/2 ≈ 0.4257 (the mode numbers are n₁ = 2 and n₂ = 1 in this case, and the phases are randomized). Two Mirnov probes, selected from a synthetic toroidal array and displaced by Δφ = 33 degrees, are used. Figure 3(a) shows the time traces of these two channels y₁, y₂. As anticipated due to the DFT resolution, the two frequencies are merged into one peak in the cross-power spectrogram of the two probes. This is seen in figure 3(b). Figure 3(c) shows the cross-coherence function for this case. Also displayed in figure 3(b) are two vertical lines. These two lines are the estimated frequencies obtained by the subspace-based method for this segment of data. A zoomed-in version of the peak cross-power spectrum and the two estimated lines are shown in figure 3(d). This example shows that in situations where the FFT-based cross-spectrum approach would have difficulties, the new method can still perform well. Note that toroidal mode-number estimates are also obtained by utilizing knowledge of the geometrical displacement of the probes. The DIII-D tokamak normally records Mirnov channels at a sample rate f_s = 200 kHz. This means that N = 200 samples equals 1 ms block length. And the resolved frequency difference Δω/2 corresponds to 500 Hz. In this particular numerical example the estimated line spectrum frequencies were $\widehat{\omega}_1=0.4094$, $\widehat{\omega}_2=0.4283$ rad/sample (which are close to the true source frequencies ω₁ = 0.41, ω₂ ≈ 0.4257 as defined above).

Zoom In Zoom Out Reset image size

2.4.2. Example comparison with the SVD approach.

Another popular approach to general magnetic probe array analysis is the SVD of the block data matrix (sometimes known as biorthogonal decomposition). The SVD may however not be so well suited for general modal analysis since it imposes constraints on the orthogonality of frequencies with respect to the time frame (as does DFT too) and an orthogonality constraint in a spatial sense. This means that the SVD can fail both in near frequency-degenerate cases as well as in spatially degenerate cases. The reason for this is that the SVD is a low-rank matrix approximation technique in the sense of the Frobenius norm (least squares) [12]. This implies that the chief criteria for SVD is to capture as much signal variation as possible with few terms. In general, this is not the same as extracting sinusoidal components from a data matrix. Figure 4 illustrates a case where the SVD fails to accomplish modal analysis (but seems to achieve data compression) whereas the new method, denoted SSI, succeeds in extracting the modal patterns. In this example the block data length is again N = 200 samples. The normalized frequency of the first mode is ω₁ = 0.41 rad/sample. The second mode frequency is set to ω₂ = ω₁ + Δω/8 (corresponding to a real frequency separation of 125 Hz and a real signal length of 1 ms for the typical sample rate of f_s = 200 kHz). The block data matrix is composed of multi-channel time-series data obtained from a toroidal array with 20 angularly equidistant synthetic probes. The SNR of the data is 20. The toroidal mode numbers are as in the first example in section 2.4.1. Note that the order of the extracted modes is arbitrary. In figure 4, it is clearly seen that the SVD 'modes'(b)–(c) fail to resemble the planar-wave-looking patterns that was used to construct the Mirnov data matrix Y(a). The SSI modes (d) and (e) are very close to the 'true' modes. As in the previous example, the method was not informed that there actually were two modes to be recovered. This is another problem encountered in the use of direct SVD: how many singular vectors to choose and which to pair?

Zoom In Zoom Out Reset image size

Figures 6 and 7 overview the results obtained by running the eigspec analysis on the DIII-D NTM control shot #150792. Figure 6 shows an FFT spectrogram that has been averaged over all available magnetic probes for this shot. Figure 7 shows the eigspec feature extraction results and can be understood to expand the spectrogram picture of figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The purpose of this experiment is to generate a m/n = 2/1 NTM that triggers a 'catch & subdue' control algorithm [38]. The control algorithm monitors the n = 1 RMS amplitude component (which is derived from the toroidal array only). When this amplitude exceeds a preset value, the gyrotrons are powered-up for electron cyclotron current drive (ECCD) on its best estimate for where the q = 2 surface is in the plasma. This type of control requires a trigger on a small NTM signature for it to efficiently remove the NTM. In theory, if the NTM is detected before its amplitude is above a certain critical marginal level, it can be shrunk relatively fast [39, 40]. If the catch is at an amplitude higher than this marginal level, the catch is significantly more power- and time-consuming. It has been observed that prompt distinction of the 2/1 NTM is challenging due to the presence of sawtooth (ST)- and fishbone-related internal kink signatures with the same toroidal period number n = 1. This additional n = 1 signature has an amplitude which is of the order of the desired trigger level for the NTM signature and easily produces false positives. A further complication is that the frequency is also comparable with the expected range of the anticipated NTM frequency.

In this shot, ECCD is triggered by the n = 1 RMS level at t = 2230 ms. As seen in figure 7 this manages to make the 2/1 mode disappear before t = 4000 ms (the 2/1 mode is labelled by cluster index 4, explained below). The ECCD is turned off at t = 5600 ms and the 2/1 mode reappears around t = 6000 ms.

The estimation parameters for the SSI Order-MAC code for this example is as follows. The block length N = 400 corresponds to a real-time block duration of 2 ms since the sample frequency is 200 kHz. The block increment is equal to the block length, so there is no overlap (but all data is used). The random projection dimension is l = 12 and the future and past horizons are respectively f = 10 and p = 20. The Order-MAC orders are r₁ = 12 and r₂ = 20. This means that the maximum number of modes that can be extracted from each block is 6 (at two states per mode). The spurious mode threshold values are both set to 0.998. The post-processing of the extracted modes was done using a 60-NN MAC similarity matrix. The mode extraction took about 30 s and the sparse-method clustering took around 25 s (on a typical laptop computer). For comparison, about 4 s is required to calculate an FFT-based mean spectrogram as shown in figure 6 (on the same hardware, using the same block size and block stride, and with zero-padded FFT of length 2048).

The clustering of the extracted modes from shot #150792 are shown in figures 7(a) and (b). It turns out that eigspec can achieve an automatic separation of the 2/1 NTM (blue symbols, cluster 4) from the ST-related internal kink (red symbols, cluster 1). The dominant classes of modes are represented as patterns in (θ, φ)-space in figure 8. The phases of the patterns are arbitrary. Note that figure 8(b), cluster 2 in figure 7, is a 3/2 pattern that corresponds to the illustration in figure 1 (if it was the only mode in the present shot). The modes labelled as 1 and 4 in figure 7 respectively correspond to the extrapolated patterns of figures 8(a) and (b). It can be seen that these two types of modal patterns show high similarity if only the outboard midplane is considered. This is the case when only the toroidal array is employed. If more poloidal shape information is used then the modal patterns appear to be separable. Specifically, the difference between these modal shapes becomes large near the inboard midplane. The internal kink mode in figure 8(c) even seems to change helicity along the poloidal angle θ. This has previously been dubbed 'phase folding' [3]. It will be further commented in section 3.3 below.

Zoom In Zoom Out Reset image size

The application of the eigspec analysis toolkit on this NTM control experiment suggests that it is possible to develop a classifier system (based on the full-array shape vector (15) extracted by the present algorithm) that can reliably distinguish between the m/n = 2/1 mode and the n = 1 ST-related internal kink mode. It does not seem enough to restrict the analysis to the toroidal array (both n = 1, and both can have comparable frequencies). Detailed poloidal shape information seems to be mandatory to clearly separate these two classes of modal pattern.

The quiescent H-mode (QH-mode) is characterized by H-mode-like edge plasma conditions but without the typical ELM activity [41, 42]. ELMs are believed to be a serious danger to the plasma wall in large machines such as ITER since they result in intense localized heat loads. The QH-mode may offer a regime that avoids the ELMs but otherwise retains the H-mode confinement advantage. In place of the ELMs, edge transport is facilitated by MHD phenomena known as edge harmonic oscillations (EHOs). The EHO seems to exist in edge plasma conditions close to the onset of ELMs. EHOs have been observed to increase edge particle transport without degrading the edge thermal transport barrier. It is thought of as a saturated kink-peeling mode driven by edge current density and edge rotational shear. EHOs have been observed both as narrow-band oscillations and broadband activity. The experimental characterization of the EHO is usually limited to assigning its toroidal period number n and its oscillation frequency f.

The present analysis toolkit can provide an experimental (external magnetics) characterization of the narrow-band EHO phenomenon. This is exemplified here for DIII-D shot #153343. This experiment attempted to demonstrate interaction with the EHO using internal saddle loop coils (excited at several kHz to match the anticipated EHO frequencies). A shape-similarity map of the extracted modes for the timeframe 2500 to 3500 ms is shown in figure 9. The plot is produced by assigning a colour index to each extracted mode based on its shape vector w and a reference shape vector v. The reference shape vector v is located at t ≈ 2932 ms and its frequency is f ≈ 4.28 kHz. The reference shape vector is marked by a red square marker in the figure. The colour index is then simply the MAC value (9) between w and v. The colour index is given the colour according to the colour bar on the right-hand side in figure 9. The estimation parameters in this case were similar to the NTM clustering example above. The block length was N = 500 (2.5 ms at f_s = 200 kHz), the block-step size was 50 samples or 0.25 ms. The past and future horizons were p = 20 and f = 10. The random projection dimension was l = 12. The Order-MAC orders were r₁ = 12, r₂ = 20. The Order-MAC threshold values were set to 0.995.

Zoom In Zoom Out Reset image size

It can be seen that the similarity measure (9) seems to take a high value for a contiguous trace. This is indeed the evolution of a clear n = 1 component in the Mirnov data. The straight lines with some inclination that can be seen in the figure are actually the 'modes' extracted for the internal coil coupling to the Mirnov array (frequency sweeps starting at 4.5 kHz at t = 2500 ms and coming down to 2.5 kHz at t = 3000 ms, and similarly in the timeframe t = 3000 ms to t = 3500 ms). It turns out here that as this externally applied frequency aligns with the EHO n = 1 trace frequency, a strong modulation results, after t ≈ 3100 ms. In fact both the amplitude (not shown here) and the frequency of the EHO trace are strongly affected. This will be analysed in more detail elsewhere. The point here is that the eigspec toolkit can effectively give an overview of the Mirnov data. Specifically, a modulation effect of the actual EHO shape appears to also exist (i.e. the colour variation of the n = 1 time trace in the modulated part of figure 9).

Figure 10 shows the spatial extrapolation results using the GPR technique for the specific reference shape vector used in figure 9. Figures 10(a) and (b) should be understood in relation to figure 5. The dashed lines in figure 5 are the one-dimensional directions along which figures 10(a) and (b) are evaluated. These evaluations (using the GPR technique) give $2\widehat{\sigma}$ 'error bars' which are indicated in the plots as dashed lines (the area within the dashed curves should really be thought of as a ∼ 95% confidence interval for the Gaussian process, given the hyper-parameters, described in section 2.2.2 and references therein). The markers in the poloidal and toroidal shape plots are the actual shape-vector values. These markers should be thought of as noisy observations. The smooth periodic curves are interpreted as the underlying functions which could have generated the noisy observations. The confidence interval is not for the possible observations but for the underlying smooth periodic function. This implies that the markers can be outside the region between the dashed lines. Figure 10(c) is the combination of the real- and imaginary-part smooth curve pairs in the sense of equation (23).

Zoom In Zoom Out Reset image size

An MHD stability code that provides EHO eigenmodes could be compared to fully empirical modal plots such as those in figure 10. Perhaps the easiest way to achieve this is to transform the theoretical mode to the actual Mirnov array coordinates as a forward problem (including sensitivity studies).

The ST cycle was first observed almost 40 years ago [43]. It starts by a peaking current density due to decreased resistivity with temperature in the plasma core. The peaking is concomitant with a 'precursor' n = 1 magnetics oscillation. The current peaking tends to push the central safety factor q₀ close to 1 and triggers an internal kink or quasi-interchange mode with m = n = 1. An abrupt re-arrangement of the plasma core flattens the temperature and current profiles. After this 'crash' event, the central q₀ tends to be above 1 and the cycle restarts. The details of the nonlinear ST cycle are complicated and is a field of active research. The internal kink and the ST cycle significantly depend on the shaping parameters of the tokamak plasma, such as ellipticity and triangularity [44, 45]. Sometimes there is a 'postcursor' magnetics oscillation that decays after the crash (these are also known as 'successor' oscillations). Sometimes the crash triggers other MHD modes such as NTMs [1].

STs are not believed to pose a direct problem for burning plasma fusion performance. Instead a main concern is that fast ions in hot plasmas appear to stabilize the internal kink mode. This has an effect of producing a slow ST cycle characterized by large ST crashes and a cycle time which is longer than an energy confinement. These 'monster' STs can be more prone to trigger NTMs than smaller STs. One recent approach to destabilize STs (to make the ST more frequent, smaller and less prone to trigger NTMs) is to use ECCD [46, 47].

For the present purposes, in connection with the NTM discrimination problem described in section 3.1, it becomes relevant to study the variability of ST precursor and postcursor modal patterns. Another motivation for studying the evolution of ST precursor Mirnov array signal patterns is to provide detailed empirical observations that can be related to ST cycle modelling. In what follows, an ST cycle with both pre- and postcursor Mirnov oscillations is analysed with eigspec. The plasma is an inboard-limited L-mode oval-shaped low-q discharge, #154887, designed to study external kink stability. At t = 1980 ms the safety factor key values are q₉₅ ≈ 2.1, q_min ≈ 1.0. Figure 11 shows an ST cycle as seen on the magnetics from this shot. A rapid growth in amplitude of the |n| = 1 component, associated with a nearly constant frequency, is observed for ∼10 ms. See figures 11(b) and (c) for the evolution of these quantities. Times t₁ and t₂ belong to this precursor growth phase. The ST crash event happens at t ≈ 1999.5 ms. The postcursor oscillation is observed to decay relatively slowly in amplitude and has a frequency which is higher than the precursor frequency. Times t₃ through t₆ belong to this postcursor decay phase.

Zoom In Zoom Out Reset image size

An |n| = 2 signal component is seen to grow and decay together with the |n| = 1 component which carries the main array signal amplitude. This signal component can be interpreted as a spatial harmonic of the ST pre- and postcursor since it has twice the frequency of the |n| = 1 signal component (and thus the same toroidal phase velocity). The ST crash also triggers an |n| = 3 signal component. This is interpreted as a 4/3 TM that rapidly decays and disappears. The array analysis code can resolve the evolution of the ST modal patterns present in the overview figure 11. The modal patterns for the |n| = 1 array signal component at times t₁–t₆ are depicted in figure 12. The precursor phase (first row; t₁, t₂) exhibits a strong inboard side pattern that seems to have reversed helicity. This is qualitatively similar to the phase-folded pattern for the ST shape representative in section 3.1. But the present plasma is quite different. The second row (t₃, t₄) exhibits a more kink-like pattern. It appears that the ST crash rapidly changes the inboard-side magnetics array pattern. The third row (t₅, t₆) shows a similar pattern to the second row. But it seems that the outboard side bulge in the modal pattern decays faster than does the inboard side. Plasma #154887 is closer to the inboard side wall than the typical diverted DIII-D plasma. The modal patterns for the |n| = 2 array signal component at times t₁–t₄ are depicted in figure 13. The precursor phase (first row; t₁, t₂) exhibits an inboard side pattern with reversed helicity, analogous to the |n| = 1 pattern. The second row (t₃, t₄) exhibits a more kink-like pattern, also in analogy with the |n| = 1 case. In contrast to the |n| = 1 case, the spatial harmonic does not seem to have a clear decay difference for the inboard and outboard sides. The |n| = 2 array signal is either not present or too weak to be detected by the method with the current estimation parameters at times t₅, t₆. Not shown here is the modal pattern for the triggered |n| = 3 signal component. It has a dominant m/n = 4/3 geometry without any 'phase folding'.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For all the shots analysed so far with the new code eigspec, the precursor phase of an ST cycle is invariably associated with a phase-folded modal pattern. In particular, this is the case for both the H-mode monster STs and ECCD destabilized STs in [47]. A typical phase-folded precursor modal pattern from ECCD-destabilized discharge #145692 is depicted in figure 14. A similar modal pattern is observed for the precursor oscillation in the non-ECCD monster-ST discharge #145861 from the same study [47]. The details of the modal patterns appears to differ slightly depending on the plasma. This is not surprising since axisymmetric shaping parameters and distances to the Mirnov probes from the plasma surface are expected to influence these details. But the overall characteristics of the precursor modal patterns seem to persist. The simplest interpretation of the phase-folded modal pattern is possibly to think of it as a consequence of the internal ideal kink mode that effectively displaces (tilts) the plasma axis in a m/n = 1/1 fashion. The crudest model of this is a single 1/1 current filament. Such a filament could create a phase-folded pattern. The same inclination of the filament is then seen from opposite sides when viewed from the outboard and inboard respectively.

Zoom In Zoom Out Reset image size

The next section 3.4 below shows another (related) class of MHD modes that seems to exhibit phase folding on the Mirnov array.

This final DIII-D analysis example illustrates that eigspec can be used to decompose dynamic and transient events. Figure 15 shows a single fishbone event [48] from DIII-D shot #155279. The shot is an ITER-like plasma with q₉₅ ≈ 5.1, q₀ ≈ 1.0 and β_N ≈ 2.9 (at t = 3205 ms), designed to explore steady-state scenarios (focus on current-profile optimization). The fishbone instability is believed to be excited by energetic particles when q_min is in the vicinity of 1 [49]. It is considered to be a finite-frequency branch of the ideal internal kink mode [50–52]. It can significantly deteriorate the efficiency of beam heating.

Zoom In Zoom Out Reset image size

The estimation parameters for the SSI/Order-MAC algorithms in this case were as follows. The block length was N = 100 (0.5 ms at f_s = 200 kHz), the block-step size was 2 samples (i.e. 10 µs). The past and future horizons were p = 25 and f = 10. The random projection dimension was l = 16. The Order-MAC orders were r₁ = 12, r₂ = 20. The Order-MAC threshold values were set to 0.992. Notice that 0.5 ms is about a tenth of the length of the entire fishbone burst as seen in figure 15(a). This means that the eigspec-analysis is indeed achieving localized time-series modelling for the fishbone event.

Figures 15(b) and (c) indicate that the fishbone in the present plasma is mainly described by an amplitude envelope and a frequency sweep. The frequency starts at about 20 kHz and sweeps down to about 15 kHz. The amplitude rapidly goes up and then down. The amplitude decay is somewhat slower than the growth. The times marked by t₁ and t₂ in figure 15 are used to explore the modal shape pattern. These times respectively exemplify different fishbone phases with amplitude growth and decay. It turns out that the fishbone event analysed here is well described by a single modal pattern. Thus the amplitude envelope can be simply interpreted as a scaling of the modal pattern. The modal pattern for time t₂ is shown in figure 16. It would look the same if it was plotted for t₁ (modulo a toroidal phase shift). It is observed that the fishbone modal pattern in this analysis exhibits phase folding, in analogy with the ST precursor oscillations above in sections 3.3 and 3.1. The phase folding observed for the fishbone in figure 16 suggests that it has a significant m/n = 1/1 internal kink mode character. The frequency of the fishbone oscillation is consistent with the plasma core velocity. Compared to e.g. wavelet spectrograms of the fishbone burst as seen in [51] (experimentally) and [52] (theoretically), the signal 'parametrization' (figure 15) produced by eigspec seems clearer and may be easier to interpret with respect to models.

Zoom In Zoom Out Reset image size