Proceedings of ITC/ISHW2007

Automatic Analysis of Multichannel Time Series Data Applied to
MHD Fluctuations
D. G. Pretty1) , B. D. Blackwell1) , F. Detering1) , M. Hegland3), J. Howard1) , D. Oliver1), M. J. Hole2)
and J. H. Harris4)
1)

Plasma Research Laboratory, andDepartment of Theoretical Physics, Research School of Physical Sciences and Engineering, and3)
Mathematical Sciences Institute, all of the Australian National University, ACT, 0200, AUSTRALIA
4)Oak Ridge National Laboratory, Tn, USA
2)

We present a data mining technique for the analysis of multichannel oscillatory timeseries data and show
an application using poloidal arrays of magnetic sensors installed in the H-1 heliac. The procedure is highly
automated, and scales well to large datasets. In a preprocessing step, the timeseries data is split into short time
segments to provide time resolution, and each segment is represented by a singular value decomposition (SVD).
By comparing power spectra of the temporal singular vectors, singular values are grouped into subsets which
define fluctuation structures.
Thresholds for the normalised energy of the fluctuation structure and the normalised entropy of the SVD areused to filter the dataset. We assume that distinct classes of fluctuations are localised in the space of phase diﬀerences (n, n+1) between each pair of nearest neighbour channels. An expectation maximisation (EM) clustering
algorithm is used to locate the distinct classes of fluctuations, and a cluster tree mapping is used to visualise
the results. Diﬀerent classes of fluctuations in H-1 distinguished by this procedure are shown to be associatedwith MHD activity around separate resonant surfaces, with corresponding toroidal and poloidal mode numbers.
Equally interesting are some clusters that don’t exhibit this behaviour.
Keywords: Data mining, plasma physics, Mirnov oscillations, magnetic fluctuations

Motivation

Rotational transform parameter scans have been undertaken in the H-1NF heliac [1, 2, 3]. Using 28 Mirnov coils
over 92 discharges, with each shot incrementally changing
the magnetic geometry, a large set of timeseries data wasproduced.
The motivation for the present work has been to findan algorithm which, with minimal human interaction, cangroup together all similar fluctuations over a large number of shots. Here we present our method, showing that it
can be used to discover interesting spectral features and tomap classes of fluctuations to any known parameter of thedataset, e.g.: magnetic geometry.

0.4
0.3

φ = 44.3°
κh = 0Coil 20

Coil 1
Coil 2
Coil 3

0.2
z [m]

1

0.1

Helical
Coil

0

PFC

-0.1
-0.2
-0.3
0.6

0.7

0.8

0.9
1
R [m]

1.1

1.2

1.3

Fig. 1 One of two poloidal Mirnov coil arrays installed in the
H-1 heliac

2

Data

Two toroidally-separated poloidal Mirnov coil arrays wereused in the experimental campaign described here, one array is shown in figure 1. An additional linear array of 5coils is also installed. In all, 28 Mirnov coils are used inthis dataset.
The experimental parameter scanned is κh , the ratio ofcurrent in the helical winding coil to that in the toroidal andpoloidal field coils. To a good approximation, κh scales linearly with rotational transform on axis (ι0 ). Due to the preauthor’s e-mail: david.pretty@anu.edu.au

cisely controllable coil power supplies, very reproducible
plasmas can be formed with an accuracy in vacuum fieldrotational transform of 1 part in 1000. The range of configurations used here is shown in figure 2, from a monotonic
profile with ι0 = 1.12, ιa = 1.28 (κh = 0) to a reverse-shear
profile with ι0 = 1.45, ιa = 1.46, ιmin = 1.41 (κh = 1.07).
The time evolution of Mirnov spectra and lineaveraged electron density for shot typical of this dataset are
shown in figure 3. The RF power is essentially constantat 60 kW, producing peak density of n̄e = 1018 m−3 . The

Proceedings of ITC/ISHW2007
κh=1.07κh=1.00κh=0.80

PNChenergy of a singular value (pk = a2k /E, E = k=1
ak for s.v.
ak and NCh < N s where NCh is the number of channels andN s is the number of samples in ∆t). Noisy short time segments can be filtered with a threshold value of normalisedP
entropy (H = − k pk log pk / log NCh ).
A mode can be be described by several singular values, for example a rotating mode will have two orthogonal bases (i.e. sine and cosine topos (spatial basis vector) with chronos (temporal basis vector) also with π/2phase diﬀerence). We assume that singular vectors whose
chronos have similar power spectra belong to the samemode, grouping together sets of singular values with normalised cross-power above a threshold value γ:

κh=0.60

z [m]

0.2
0
-0.2
1.0

1.2
R [m]

1.4

1.0

1.2
R [m]

1.4

1.0

1.2
R [m]

1.4

z [m]

1.45
10/7

0
-0.2

κh=1.077/5

1.40

κh=1.00κh=0.20

κh=0.80

1.35

0.2
z [m]

4/3
κh=0.60

1.30

0

9/7
5/4

1.25

-0.2

κh=0.40κh=0.00κh=0.20

6/5

1.20

0.2
κh=0.00z [m]Rotational transformκh=0.40
0.2

1.50

7/6
8/7
9/8

1.15

0
-0.2

1.10
0.00

0.05

0.10

0.15

0.20

1.0

r

1.2
R [m]

1.4

200
150
100
50
0
0

10

20

30

40

50
100

18

ne[×10 m−3]

1.4
1.2

80

1.0

RF

0.8
0.6
0.4

60
40

ne

20

0.2
0.0
0

RF [kW]

Frequency [kHz]Fig. 2 Rotational transform profiles used in this campaign.
Poincaré plots shown for the various magnetic configurations.

10

20

30

40

50

0
60

Time [ms]Fig. 3 Typical shot for this campaign (at κh = 1.0), showing
Mirnov spectra in the top panel and line-averaged electron density in the lower panel. Heating power (60 kW
RF) is essentially constant throughout the discharge.

Mirnov spectra shows multiple co-existing modes, with thehigher frequency features ( f > 50 kHz) showing Alfvénic
density scaling ( f ∝ n−1/2).
e3Data Pre-processing

Data mining procedures generally have 3 stages: firstly thedata needs to be pre-processed into a format suitable forthe main algorithm; secondly the main algorithm (neuralnetwork, clustering, association rules, etc) is applied to thedata, finally the results are visualised and interpreted. The
pre-processing used here involves separating out diﬀerentmode components from the timeseries data and mappingthem to a phase space in which the main clustering algorithm is used to distinguish diﬀerent classes of fluctuation.
For each shot, we split the timeseries data into shorttime segments, in this case we use ∆t = 1 ms. For each
short time segment, we take the singular value decomposition (SVD) [4] of all Mirnov channels. Noise is filteredout of the system by placing a threshold on the normalised

γa,b =

G(a, b)2, G(a, b) = h|F (a)F ∗ (b)|i, (1)G(a, a)G(b, b)

where F is the Fourier transform, and h. . .i represents the
spectral average. An examination of a randomly selected
subset of SVDs showed that a threshold of γ = 0.7 is suitable for the present dataset. Each group of singular values
defines a fluctuation structure.
For each fluctuation structure, we take the inverse
SVD using only the allocated subset of singular values (others are set to zero) to return timeseries data foreach coil representing the given fluctuation. From thesetimeseries we take the phase diﬀerences between nearest
neighbour coils to produce coordinates in NCh -dimensional
phase space (“∆φ-space”) in which the clustering algorithm operates.

4

Clustering

We aim to discover any underlying lower-dimensionalmodel of the dataset, i.e.: to group data points into classes
or modes which can then be mapped back to other properties, such as magnetic geometry. We assume that distinct
modes of fluctuation will be localised in the ∆φ-space defined by the nearest neighbour phase-diﬀerences. This ∆φspace localisation can be easily understood in the simplified case of poloidally equispaced Mirnov arrays in cylindrical geometry – here a mode with poloidal mode number
m will be localised in the region of ∆φ(i, i + 1) = 2πm/N M ,where N M is the number of Mirnov coils in the array.
Clearly the heliac geometry is not so simple, however theclustering algorithm can find arbitrary phase structures andthere is no need to interpret the phase structure before clustering.
We use the expectation maximisation (EM) clustering algorithm which finds the most likely values of latentvariables in a probabilistic model [6]. Here we model thedata by NCl NCh -dimensional Gaussian distributions in ∆φspace, with mean and standard deviation for each cluster
as the set of latent variables. To ensure the results are notbiased by the assumption of Gaussian cluster shape, we
compare the EM clustering results with results from an ag-

Proceedings of ITC/ISHW2007
AH:A

AH:B

(520)

(169)

AH:C

(164)

AH:D

(104)

EM:J

10 (95)

EM:K

10 (352)

EM:L

10 (422)

EM:M

10 (156)

EM:N

10 (57)

EM:O

10 (128)

AH:(Remainder)
(884)

AH:E

EM:B

EM:A

AH:F

(50)

AH:G

(40)

(36)

AH:H

(33)

10 (162)

EM:F

10 (168)

2 (339)1 (2000)

EM:G

EM:C

0.0

10 (308)

EM:I

10 (152)

8 (96)

0.2

0.4

0.6

κh

7/5

4/3

5/4

1 (2000)

6/5

7/6<r> [m]

120
100
80
60
40
20
0

EM:A

EM:H

3 (1309)

5 (433)

EM:D

Freq. [kHz]

EM:E

0.8

1.0

Fig. 4 A Cluster tree representation of the κh scan data. The figure in the bottom left corner contains the whole dataset; its upper panel
shows the fluctuation structures mapped to f and κh , the numbers 1 (2000) at the top right are the tree level, NCl , and cluster
population respectively. The lower panel shows the contours of low-order rational surface within the plasma minor radius (y-axis)
for the configurations. For clarity, only a subset of clusters within the tree have their contents displayed and EM:G has been
displaced to prevent overlap. Vertical parent-child distance is proportional to the distance between cluster means, while line
thickness is inversely proportional to the Gaussian width of the cluster. Several AH clusters are also shown for comparison.

glomerative hierarchical (AH) clustering algorithm. Goodagreement is found between results from both methods.

5

Results

Shown in figure 4 is a cluster tree graph with plots showing
frequency vs. κh of selected clusters, each point plotted ina graph corresponds to a fluctuation structure. The root of
the tree (left side) has the trivial case of NCl = 1, such thatthe whole dataset is contained in a single cluster. The highest level branch shown here (right side) has NCl = 10. The
localisation of the clusters in the f − κh projection shown
here is due to the relation of the phase structure to the magnetic configuration – neither frequency nor κh are includedin the clustering metric.
Vertical displacement at a node in the tree is proportional to the distance in phase space between the parent andchild clusters. Results from AH clustering are also shownfor comparison. The AH clusters show good agreementwith EM clusters at the NCl = 10 end of the tree.
Near the root of the tree we find well defined clusters which have resonant structure in the f − κh projection.
For the EM:B cluster, the frequency minimum occurs at
κh = 0.4, corresponding to the ι = 5/4 rational surface
in the low-shear region of the plasma. The EM:C cluster
has frequency minimum at κh ≃ 0.72, where the ι = 4/3
surface is located in the region of zero-shear. Both modes
show f ∼ |ι − n/m| resonant behavior, as do other clusters

for higher order rational surface configurations. Analysisof the phase structure of these modes has shown that thedominant Fourier components are those expected for the
relevant rational surface, i.e.: (n, m) = (4, 3) for the ι = 4/3
mode [7].
Cluster EM:O has very well defined (n, m) = (0, 0)
structure and also exists in configurations where low orderrational surfaces exist in low-shear regions. Cluster EM:J
appears to be a helical Alfvén eigenmode (HAE) coupled between the (7, 5) and (10, 7) resonance, this δn = 3,
δm = 2 is consistent with the relatively large (n, m) = (3, 2)
Fourier component of the heliac magnetic geometry.

6

Discussion

The clustering occurs only in ∆φ-space, and is unbiasedby the κh and frequency coordinates plotted in cluster tree
figure 4. The clusters can also be mapped to any otherknown plasma properties.
Complications arise in the case of H-1 configurationscans due to the changing shape of magnetic field withκh . The coil coordinates have been mapped to κh -averaged
magnetic angles to account for this.
The process has been kept general enough for it to be
applicable to any set of geometrically ordered timeseries
data. While the simple phase diﬀerence clustering metricworks fine in this case, an alternative metric may be required if timeseries from other diagnostics are included.

Proceedings of ITC/ISHW2007
The scalability depends on the clustering algorithm
used. The EM method scales well with the number of datapoints Nd p (Nd p × NCl ), while the AH method scales poorly(Nd2p ). Variations to the EM clustering algorithm and alternative ways to quantify the results are areas of ongoinginvestigation.
Data mining methods apart from the clustering analysis described here are also being investigated. “Association rule mining” is one such approach. It is designed
to discover temporal patterns in the entire spectrograms
such as shown in figure 3. In this context we search
for rules (frequent patterns) in the form of “if a strong
mode at f ≃ 10 kHz occurs from 10-20ms then a modeat f ≃ 20 kHz from 40-50ms will occur (with e.g.: 50%
probability)”.

7

Acknowledgements

The authors thank the H-1 team for continued supportof experimental operations. This work was performedon the H-1NF National Plasma Fusion Research Facility
established by the Australian Government, and operated
by the Australian National University, with support fromthe Australian Research Council Grant DP0344361 andDP0451960.

References[1] S. M. Hamberger, B. D. Blackwell, L. E. Sharp and D. B.
Shenton, H-1 design and construction. Fusion Technol. 17(1990) 123–130[2] J. H. Harris et al, Fluctuations and stability of plasmas in the
H-1NF heliac. Nucl. Fusion 44 (2004) 279–286[3] B. D. Blackwell, Results from helical axis stellarators, Phys.
Plasmas 8 (2001) 2238–2244[4] T. Dudok de Wit, A.-L. Pecquet, J.-C. Vallet and R. Lima,
The biorthogonal decomposition as a tool for investigating
fluctuations in plasmas. Phys. Plasmas. 1 (1994) 3288–3300[5] A. Dempster, N. Laird and D. Rubin, Maximum likelihoodfrom incomplete data via the EM algorithm. J. Royal Stat.
Soc. 39 (1977) 1–38
[6] W. H. E. Day and H. Edelsbrunner, Eﬃcient Algorithms for
Agglomerative Hierarchical Clustering Methods, J. Classification 1 (1984) 7–24
[7] D. G. Pretty, PhD Thesis, Australian National University(2007)