# Towards Explainable Metaheuristic: Mining Surrogate Fitness Models for
Importance of Variables

Manjinder Singh<EMAIL_ADDRESS>0000-0003-4720-3473 University of
StirlingStirlingUK ,  Alexander E.I. Brownlee<EMAIL_ADDRESS>0000-0003-2892-5059 University of StirlingStirlingUK and  David Cairns
<EMAIL_ADDRESS>0000-0002-0246-3821 University of StirlingStirlingUK

(2022)

###### Abstract.

Metaheuristic search algorithms look for solutions that either maximise or
minimise a set of objectives, such as cost or performance. However most real-
world optimisation problems consist of nonlinear problems with complex
constraints and conflicting objectives.

The process by which a GA arrives at a solution remains largely unexplained to
the end-user. A poorly understood solution will dent the confidence a user has
in the arrived at solution. We propose that investigation of the variables
that strongly influence solution quality and their relationship would be a
step toward providing an explanation of the near-optimal solution presented by
a metaheuristic.

Through the use of four benchmark problems we use the population data
generated by a Genetic Algorithm (GA) to train a surrogate model, and
investigate the learning of the search space by the surrogate model. We
compare what the surrogate has learned after being trained on population data
generated after the first generation and contrast this with a surrogate model
trained on the population data from all generations.

We show that the surrogate model picks out key characteristics of the problem
as it is trained on population data from each generation. Through mining the
surrogate model we can build a picture of the learning process of a GA, and
thus an explanation of the solution presented by the GA. The aim being to
build trust and confidence in the end-user about the solution presented by the
GA, and encourage adoption of the model.

genetic algorithms, explainability, interpretable, surrogate model, fitness
function, optimization

††journalyear: 2022††copyright: acmlicensed††conference: Genetic and
Evolutionary Computation Conference Companion; July 9–13, 2022; Boston, MA,
USA††booktitle: Genetic and Evolutionary Computation Conference Companion
(GECCO ’22 Companion), July 9–13, 2022, Boston, MA, USA††price: 15.00††doi:
10.1145/3520304.3533966††isbn: 978-1-4503-9268-6/22/07††conference: GECCO ’22:
The Genetic and Evolutionary Computation Conference; July 09–13, 2022; Boston,
MA††ccs: Theory of computation Models of learning††ccs: Theory of computation
Theory of randomized search heuristics

##  1\. Introduction

The growing ubiquity of machine learning (ML) and Artificial Intelligence (AI)
models (Olhede and Wolfe, 2018) to automate decision making in our day-to-day
life has brought with it ethics concerns(Yapo and Weiss, 2018), and a lack of
trust from those users directly impacted by the automated decision making
(Toreini et al., 2020; Aoki, 2020). The interest in explainability is not a
new one (Swartout and Moore, 1993), however the resurgence in interest is
being driven by users demanding an explanation of why a particular decision
was reached. This demand for an explanation has also been enshrined in
European Union (EU) law when the right to explanation was included in the
General Data Protection Regulation (GDPR) (Goodman and Flaxman, 2017) as a
recognition of the rise in importance of AI systems in automated decision
making.

The majority of AI models in use today can be thought of as ”black-boxes”,
with many complex layers consisting of nonlinear transformations (LeCun et
al., 2015). It can be difficult for end users to understand the decision
making process and hence, already, trust in the final output is eroded. This
has led to decisions being made that were biased (Lambrecht and Tucker, 2019),
and some that have led to real harm being done to people due to a lack of
transparency and explainability of these models (Rudin, 2019).

A large area within AI is metaheuristic search algorithms applied to
optimisation problems. These algorithms look for solutions that either
maximise or minimise a set of objectives, such as cost or performance. However
most real-world optimisation problems consist of nonlinear problems with
complex constraints and conflicting objectives. These algorithms find
applications in many areas, including optimisation of transportation for
positive environmental impact (Brownlee et al., 2018), and healthcare (Ochoa
et al., 2020).

Metaheuristic algorithms are not problem dependent, which means we can use
them for a variety of optimisation problems. Many of these algorithms are
inspired by nature. A Genetic Algorithm (GA) takes its inspiration from
evolution. To find an optimal or near-optimal solution for a problem, the GA
is used to optimise a problem by initially creating many random solutions,
populations, and then iteratively updating these solutions through selection,
mutation and crossover, to reach a good enough solution or the most optimal
solution possible. A GA performs well in combinatorial optimisation where we
have a large search space (Anderson and Ferris, 1994). However, despite recent
advances in theoretical understanding of GAs (Malan, 2021; Lehre and Oliveto,
2017), the process by which a GA arrives at a solution remains largely
unexplained to the end-user. Herein lies the problem, firstly a poorly
understood solution will dent the confidence a user has in the arrived at
solution. Secondly during our problem definition we may have excluded some
critical criteria. This may be of particular relevance in a scheduling
optimisation problem, where we may focus on only a distance or financial cost
constrain and so fail to take into account personal preferences and/or
convenience for the end-user. Thirdly, a GA follows a random process, which
can lead to noise in the presented solution (Di Pietro et al., 2004). In this
case it would be useful to know the characteristics of this solution, and be
able to tease out this noise without impacting the solution quality.
Metaheuristic algorithms in general are adept at finding shortcuts in a
problem definition (Lehman et al., 2020). Therefore it is useful to know
whether the solution presented genuinely solves the problem or if the
algorithm merely found a loophole in the problem definition. As demonstrated,
for the Vehicle Routing Problem (VRP)(Arnold and Sörensen, 2019), being able
to explain the characteristics that constitute a near-optimal or sub-optimal
solution also goes a long way towards designing a more efficient algorithm.

Taking all of the above points into consideration we can see that there are
many points at which we can lose a user’s trust when designing a system. This
loss of trust would ultimately lead to lack of adoption of the system, and
hence, wasted resources and suspicion of any systems recommended in the
future.

In this paper we will outline an approach to identify the cardinal
characteristics of metaheuristic derived solutions. We do this for a series of
binary-encoded benchmark problems, by way of identifying which variables have
the greatest influence on solution quality. Our approach will use surrogate
fitness functions (Brownlee and Wright, 2015; Jin, 2011). Usually these are
applied where the fitness function is costly; in parallel to the optimisation
process we train a computationally cheap model, and the majority of calls to
the costly fitness function are replaced by a call to the surrogate model. An
additional benefit of this process is that the surrogate is an explicit
representation of what the metaheuristic has learned about the search space
and more crucially what comprises a good solution. In our study, we take a
high fitness solution for each of the benchmark problems and probe it using
the surrogate model to identify which variables are of greatest importance and
to test their impact on fitness. The present paper extends the approach
initially proposed in (Wallace et al., 2021), by investigating the insights
that might be gained at different stages of the algorithm’s run, rather than
just the final population.

In Section 2 we will review related work around trust and explanation of
metaheuristic optimisation, before going on to focus on the role of surrogate
models in optimisation problems, and our use of them to explain metaheuristic
algorithms. In Section 3 we describe our approach to mining the surrogate
model to explain solution quality. In Section 4 we outline our approach to
mining the surrogate model to extract information about variable importance of
an optimal solution. In Section 5 we demonstrate our methodology using four
benchmark problems, and in Section 6 we present our results and outline our
conclusions and future work in Section 8.

##  2\. Related Work

Explainability in AI is not a new research topic, however, research conducted
in the early 90’s (Swartout and Moore, 1993) tapered off for many decades. The
recent growth and interest in this area is gaining momentum, mainly driven by
a lack of trust and uptake by users of AI systems. The resurgence in this
topic follows the thinking that making AI systems more transparent and
explainable will engender end-user trust in the AI systems and hopefully lead
to greater uptake and use of these systems (Miller, 2017).

There is, however, a gap in research centered around explainability of
metaheuristic algorithms. The most relevant research within this area being
innovization outlined in a paper by Deb et al. (Deb et al., 2014) (Deb and
Srinivasan, 2006), which proposed “innovization” to generate problem based
knowledge alongside the normally generated near-optimal solutions, through
identification of common principles among Pareto-optimal solutions for multi-
objective optimisation problems. The impetus behind this being that most
optimisation techniques adopted are used to surrender a single or small
selection of optimal solutions, their research built on this to allow these
optimisation techniques to find additional problem based knowledge in parallel
to the generated optimal solution(s) via the innovization operation. By taking
a generated set of high-performing trade-off solutions and identifying common
principles hidden within them. The idea is that the main common thread amongst
this set of high-performing solutions will represent properties that ensure
Pareto-Optimality and by extension are valuable properties related to the
problem globally.

More recently proposed by Urquhart et al. (Urquhart et al., 2019) is an
application of Map-Elites to increase trust in metaheuristic algorithms. This
paper aims to address the criticism that end-users have no role in the
construction of the end solution. The authors applied Map-Elites, which allows
for construction of a high-performance solution which is mapped onto a set of
solutions defined by the user, such as cost or time. The solutions being
generated by mutation and recombination, with each solution being assigned a
bin within the solution space, with solutions that are assigned to an already
occupied bin being either rejected or accepted based on a comparison of their
relative fitness, with higher ranked fitness solutions being accepted, so only
the highest fitness being assigned and retained within any particular bin.
Ultimately the end-user is responsible for choosing the preferred solution,
this ensures that input from the user is taken into account and they have a
sense of ownership and a greater level of trust in the proposed solution. Map-
Elites therefore provide a way to filter the solution space, and provides a
set of solutions for the user, from which they can select the one most
applicable to them and their needs. This has the benefit of increasing trust
in the solution selected, because the user is provided an opening to the
process and able to have some measure of influence as to what constitutes a
good solution.

Gaier et al. (Gaier et al., 2017) proposed a hybrid approach by using Map-
Elites alongside a surrogate model to add efficiency to the Map-Elites
process. They proposed reducing the need for the large number of checks
normally required for Map-Elites. The proposed solution, Surrogate-assisted
illumination (SAIL), aims to achieve this by way of integrating an
approximation model (surrogate) alongside intelligent sampling of the fitness
function. As with Map-Elites the search space is partitioned into bins each of
which holds a map with a different layout of feature values. Firstly a
surrogate is constructed based on an initial population of possible solutions
including their fitness scores. Map-Elites is then used to produce solutions
to maximise the fitness function, and generating an acquisition map.
Thereafter, new solutions are sampled from this map and additional
observations are used to iteratively improve the model. With the aim of
looping through this process to generate increasingly better solutions with
higher fitness functions. The performance predictions then being used by Map-
Elites in place of the original fitness function to generate a prediction map
of near optimal representations.

In contrast to the above approaches we do not aim to summarise or present a
set of solutions, rather, we aim to explain a single solution chosen by the
metaheuristic algorithm.

##  3\. Surrogate Models

Evolutionary Algorithms (EA) aim to minimise or maximise a given function
f(X). This is done by initially generating random solutions and evaluating
them against a fitness function, and then generating new solutions, usually
biased towards solutions with higher fitness values. The fitness function is
often the most costly operation in the optimisation process, and to reduce
this cost we can replace calls to the fitness function with calls to a much
cheaper model. This cheaper model, surrogate (Brownlee et al., 2015; Jin,
2011, 2005; Wallace et al., 2021), also represents an explicit model of the
population. We propose mining this model to capture the sensitivity of the
fitness function to the problem variables. The premise being that the
surrogate model is biased by the population of the EA, and therefore contains
another view of the algorithms understanding of the problem.

As in (Wallace et al., 2021), in the work presented here, the EA alternates
between using the surrogate to evaluate solutions and the true fitness
function. The surrogate model is re-trained each time new evaluations are
carried out by the true fitness function, it is re-trained on the current and
all previous populations at alternating intervals. So, we alternate between
using the true fitness function and the surrogate at regular intervals. The
surrogate is implemented by a Support Vector Regression (SVR) model, provided
by scikit-learn (version 1.0.2). We use the default parameter settings for the
SVR model. The model features are simply the problem variables X and the
target for prediction is fitness f(X). We hypothesise that the surrogate model
will retain some of the key properties from the original fitness function,
such as regions of high fitness. The surrogate model is able to evaluate the
solution at a much cheaper cost compared to the real fitness function. The two
benefits of this approach being that firstly the surrogate model can evaluate
the potential solutions in a more timely manner and, secondly can give us an
insight into the algorithms understanding of the search space.

The benefit of this approach is that the surrogate model, being an explicit
model of the population, can be mined to gain insights into the learning of
the EA of the search space. This approach has previously been demonstrated
when looking at Estimation of Distribution Algorithms (EDAs), which sample
probabilistic models after first constructing these models to reflect the
distribution of high fitness solutions. It has been demonstrated that
information gained from EDAs can in many cases offer additional insight into
the problem as the solutions presented by an EA.(Brownlee et al., 2010, 2013;
Brownlee, 2016). The usefulness of the information gleaned from a surrogate
model is of course closely coupled to the problem being analysed. The
surrogate model could be inherently interpretable (such as a linear regression
model similar to (Brownlee, 2016) or a decision tree), or we could exploit
existing XAI approaches such as mining the model through probing. We make use
of the latter approach in the present paper.

##  4\. Methodology

To determine the rank and importance of individual variables for each of the
benchmark problems, we will investigate the local sensitivity of the closest-
to-optimal solution (highest fitness) found for each of the problems. We will
take each of these best solutions and evaluate against our surrogate model,
allowing us to determine the change in surrogate fitness due to mutating each
variable in the solution. We will also set the sign of the “importance”
measure to be negative if the corresponding bit in the starting seed solution
was 1: providing an indication of the direction of the relationship between
that variable and the fitness. This gives us a measure of the variables
importance and its correlation with fitness, in the neighbourhood of the high-
fitness solution.

We aim to highlight the variables that impact either negatively or positively
on solution quality. It should be noted that this approach is applied to
benchmark problems for demonstration purposes only, where the simple problem
definition means that the explanations can be compared against known
expectations, and gives a starting point for further research in this area for
real-world problems. The surrogate is trained on data derived from the real
fitness function, hence, it learns the local search space from, and is biased
towards, the best solutions from the metaheuristic. In this way the surrogate
reflects the understanding of the problem as learned by the algorithm.

We will compare the surrogate model trained on only the first generation
against a surrogate trained on all populations over the run of 100
generations. Because the surrogate is biased toward the best solution, we aim
to visualise the surrogates understanding of the problem between the first
generation and over all generations.

A formal definition of our approach is shown in Algorithm 1.

In: $x=(x_{0}\ldots x_{n}),x_{i}=\\{0,1\\}$, near-optimal seed solution found
by GA

In: $S(X)\rightarrow f$, surrogate fitness function to estimate fitness $f$ of
a solution $X$

Out: $C=(c_{0}\ldots c_{n}),c_{i}\in\mathbb{R}$, absolute change to surrogate
fitness for each variable in $x$

$C\leftarrow\emptyset$;

$f_{org}\leftarrow S(x)$ $\triangleright$ surrogate fitness of solution

for for each variable $x_{i}$ do)$i=0$ to $n-1$

$x_{i}\leftarrow(x_{i}+1)\bmod{}2$ $\triangleright$ flip variable $x_{i}$

$\hat{f_{i}}\leftarrow S(x_{i})$ $\triangleright$ surrogate fitness of mutated
solution

if $x_{i}=1$ then $c_{i}=-1*\hat{f_{i}}$ $\triangleright$ Change in surrogate
fitness, optimum should have a 1

else$c_{i}=\hat{f_{i}}$ $\triangleright$ Change in surrogate fitness, optimum
should have a 0

end if

$C\leftarrow\lvert c_{i}\rvert$ $\triangleright$ add to list

end for

Algorithm 1 Probing variables in a solution with respect to the surrogate
fitness function

##  5\. Experiments

We will focus on four well-known bit-string encoded benchmark functions (1D
Checkerboard, 2D Checkerboard, Trap5 and MAXSAT).

The experiment setup is shown in Table 1. Where P is the population size, n
the length of the bit-string, problem-size, maxGen is the maximum number of
generations, mutRate the mutation rate. The Selection operator used is
Tournament with a tournament size of 5. The Crossover operator is Uniform. We
have made little effort to tune the algorithm beyond simple empirical
exploration of these parameters, because the focus of this study is on
explaining the results rather than maximal algorithm performance.

P | n | maxGen | mutRate | Selection | Crossover
---|---|---|---|---|---
100 | 100 | 100 | 0.01 | Tournament(5) | Uniform
Table 1. Experiment Setup

###  5.1. Benchmark Problems

Our experiment focused on four benchmark functions using a bit-string
representation. The four functions were chosen to contrast uni-variate and
multi-variate problems. This will allow us to investigate variations in the
importance attached to each variable and model’s handling of the presence, or
absence, of interactions.

####  5.1.1. 1D Checkerboard

The objective of the 1D-Checkerboard problem (Baluja and Davies, 1997) is to
realise a checkerboard pattern with alternating 1s and 0s. The function scores
the chromosome based on the sum of adjacent variables that do not share the
same value. A formal description of the function is shown in Equation 1.

(1) $f(x)=\sum_{i=0}^{l-2}\begin{Bmatrix}1,&x_{i}\neq x_{i+1}\\\
0,&x_{i}=x_{i+1}\end{Bmatrix}$

####  5.1.2. 2D Checkerboard

The 2D-Checkerboard problem (Baluja and Davies, 1997; Larrañaga, 2002)
introduces bivariate interactions into the problem. While these interactions
are weighted equally in the fitness function, implicitly those in the centre
of the grid have greater impact on fitness. A solution represents the rows of
a s x s grid concatenated into one string, the objective being to realsie a
grid with a checkerboard pattern of alternating 1s and 0s. A formal
description of the function is shown in Equation 2

(2)
$f(x)=4(s-2)^{2}-\sum_{i=2}^{s-1}\sum_{j=2}^{s-1}\begin{Bmatrix}\delta\left(x_{ij},x_{i-1j}\right)+\delta\left(x_{ij},x_{i+1j}\right)\\\
+\delta\left(x_{ij},x_{ij-1}\right)+\delta\left(x_{ij},x_{ij+1}\right)\end{Bmatrix}$

####  5.1.3. Trap5

The Trap-5 problem is designed to be intentionally deceptive (Deb and
Goldberg, 1994), by rewarding steps towards local optima in small groups of
variables in order to force, or deceive, an EA away from the global optimum.
This is of particular concern in algorithms that do not consider interactions
between variables(Goldberg et al., 1992; Brownlee, 2009). The problem
construction consists of bit-strings partitioned into blocks with their
fitness scored separately. A formal description of the function is shown in
Equations 3a and 3b

(3a) $f(x)=\sum_{i=1}^{n/k}trap_{k}({x_{b_{i+1}}}+...+{x_{b_{i+k}}})$

(3b) $trap_{k}(u)=\left\\{f_{\text{high}}\text{ if
}u=k,f_{\text{low}}-u\frac{f_{\text{low}}}{k-1}\text{ otherwise}\right\\}$

####  5.1.4. MAXSAT

The Maximum Satisfiability or MAXSAT problem (Johnson, 1974; Rana and Whitley,
1998; Brownlee et al., 2007) attempts to find a set of values which maximises
the number of satisfied clauses of a fixed predicate logic formula expressed
in conjunctive normal form (CNF). It is known to be NP-complete in its general
form. The MAXSAT is useful for modeling high order interactions, because each
instance of the problem uses a known predefined structure. The problem
construction involves a bit-string, where each bit encodes a predicate
variable in the CNF formula. An individuals fitness is therefore just equal to
the number of satisfied clauses.

###  5.2. Experimental Procedure

For each of the benchmark problems, we started with an initial population of
100 individuals and an Genetic Algorithm (GA) was used to determine a near-
optimal solution. A problem size of 100, population size of 100, tournament
selection with size 5, mutation rate of 0.01, crossover rate of 0.95 and was
run over 100 generations. We used elitism, to guarantee that the single
fittest solution was carried forward for each generation without mutation.

A normal run of a GA would consist of generating an initial population for
which a fitness value is calculated and through selection, crossover and
mutation a second population is generated. This process iterates for a
predetermined number of generations, or until a solution with the possible
fitness, as defined by our fitness function, is found, whichever occurs first.
At the end of this run, the most optimal solution found is presented to the
user. However, the solution presented is devoid of an explanation, as to how
it was found nor why it is considered the fittest solution, other than it has
the best fitness function value, which was determined by the algorithm
designer.

We will attempt to provide an intuitive understanding of the solution
generation process, with the aim of providing an explanation for the end-user
as to how the solution was chosen. We will do this via construction of a
surrogate model, using the population and fitness values for each generation
of the GA run as training data for the surrogate. The aim in this experiment
is not to generate a surrogate to speed up the fitness evaluation process, but
rather to mine the surrogate for importance of variables. We will examine the
surrogate model trained on only the population form the first generation and
one trained on the the population data for all of the runs over 100
generations.

For each problem we used a Support Vector Regression (SVR) model, provided by
scikit-learn (version 1.0.2), as the surrogate. We used the default hyper-
parameters for this model, and made no attempt to tune any of the parameters,
so as to enable reproducible experiments and focus on the learning of the
search space by the surrogate model between the first and all generations of
the GA run.

##  6\. Results and Discussion

###  6.1. Surrogate Trained on First Generation Data

In this section we will consider the results of mining a surrogate fitted to
the first population of the algorithm run. Each plot in Figures 1, 3, 5 and 7
shows the mean contribution to surrogate fitness for each variable. The
positive values indicate that the optimal value was a 0 and negative values
indicate that it was a 1. In Figure 1 we can observe that the model has not
managed to grasp enough of the problem structure, there is a clear partial
chain structure visible, however, it is difficult to see the relative
contribution of each variable due to random noise. We would expect to see an
equal contribution from all the variables. There is a similar picture with
Trap5 (Figure 7) and MAXSAT (Figure 5). Although with the latter there should
be some variation, as variables appear in different clauses and so have
different contributions to overall fitness, there is only a small variation in
the number of clauses a variable appears in, and not as much as suggested by
the bar plot. For 2D checkerboard there is a clearer picture of the problem
structure emerging, suggesting that some variable have twice the importance of
others. This appears to be driven by the suboptimal solutions in the first
generation; as we will see, the model fitted to all generations weights all
variables similarly and shows clearer signs of the checkerboard pattern.

In there results, the model is trained on only the first population, at this
point the GA is just starting an exploration of the search space, and so has
not explored enough of the problem structure to generate meaningful data on
which to train the surrogate model.

Our aim here was to set a starting point to which we can then compare the
results as outlined in Section 6.2

###  6.2. Surrogate Trained on All Generation Data

In this section we will consider the results of mining a surrogate fitted to
all populations of the algorithm run. Each plot in Figures 2, 4, 6 and 8 shows
the mean contribution to surrogate fitness for each variable. The positive
values indicate that the optimal value was a 0 and negative values indicate
that it was a 1.

Comparing these plots to those generated for only the first population data,
we see a marked contrast. We observe that the surrogate model has picked up
key characteristics of the problem; for Trap5 we can now start to see clearer
groupings of variables as we would expect for this problem. However the
surrogate is reflecting that the GA has converged on a local optima, and hence
the clear grouping that has settled on all 0 for variables 85-90.

For the 1D Checkerboard we see a clearer pattern of alternating 1s and 0s and
this pattern is more pronounced, and each variable seems to be contributing
equally to fitness.

For 2D Checkerboard, we see that the variables share a similar relationship to
fitness and fall on or around the middle, and we see some regular groupings of
slightly more dominant variables, which seem mostly clustered around the
central region of our checkerboard.

For MAXSAT, we see similar grouping of clauses and clustering of the variables
that contribute more to overall fitness.

We can note that the explanatory method of using a surrogate model does reveal
that the model has begun to detect the key components of the problem.

Figure 1. 1D Checkerboard Contribution to Surrogate Fitness per Variable
(First Generation) Figure 2. 1D Checkerboard Contribution to Surrogate Fitness
per Variable (All Generations) Figure 3. 2D Checkerboard Contribution to
Surrogate Fitness per Variable (First Generation) Figure 4. 2D Checkerboard
Contribution to Surrogate Fitness per Variable (All Generations) Figure 5.
MAXSAT Contribution to Surrogate Fitness per Variable (First Generation)
Figure 6. MAXSAT Contribution to Surrogate Fitness per Variable (All
Generations) Figure 7. Trap5 Contribution to Surrogate Fitness per Variable
(First Generation) Figure 8. Trap5 Contribution to Surrogate Fitness per
Variable (All Generations)

##  7\. Acknowledgments

The Data Lab funded project (Grant: 191047)

##  8\. Conclusion

This paper set out to investigate the possibility of mining surrogate models
for explanations about solution quality at different stages of an evolutionary
run. We looked at the surrogate models learning of the search space for four
benchmark problems; 1D Checkerboard, 2D Checkerboard, MAXSAT and Trap5. We did
this through using the population and fitness values for individuals generated
by a GA. We used this population data as input to train an SVR model, which we
mined for a deeper understanding of the solutions generated by the GA for each
generation.

This work has demonstrated a starting point for further exploration and more
importantly the generation of explanations around the near-optimal solution
presented by a GA. We focused on mining the surrogate to understand how the
mean contribution of each variable changed as the GA generated more insights
into the population and problem at hand. In particular, we compared a
surrogate model trained on the initial population against one trained on
population data after 100 generations. These experiments suggest that, for the
purposes of explanation, it is preferable to train the surrogate on all
solutions visited over the course of the GA run. Although, there is also some
evidence to suggest that useful information about the problem can be mined
during the first generation. Further investigation into what is possible at
each stage of the run is thus needed.

Future work will involve analysis of interaction of variables and
investigation of multi-modal problems. As understanding of this approach to
explanations grows we will apply the above method to real-world problems. We
will also further investigate the performance of a surrogate model at
different stages of a GA run; to generate comparisons at different generations
of the GA run and investigate at what point it would be beneficial to switch
from using a costly fitness function to the surrogate model and comparison of
solutions generated between the two scenarios.

We also propose to further extend this work to exploit recent but already well
known techniques for explaining machine learning models, such as permutation
feature importance measurement(Breiman, 2001; Fisher et al., 2018; Wei et al.,
2015) and SHapley Additive exPlanations (SHAP) (Lundberg and Lee, 2017) as a
way to explain variable selection for individual solutions. The work presented
here is a first step towards building a more generalised framework for
explainability within the evolutionary computation (EC) domain.

The ultimate aim of this work is to generate explanations of the solution
generated, by mining the surrogate model, and presenting these explanations to
the end-user to instill more confidence and trust in our model. In turn, this
increased trust in the solutions should lead to greater uptake and increased
benefit from the use of GAs in practice.

## References

  * (1)
  * Anderson and Ferris (1994) Edward J. Anderson and Michael C. Ferris. 1994.  Genetic Algorithms for Combinatorial Optimization: The Assemble Line Balancing Problem.  _ORSA Journal on Computing_ 6, 2 (May 1994), 161–173.  https://doi.org/10.1287/ijoc.6.2.161
  * Aoki (2020) Naomi Aoki. 2020\.  An experimental study of public trust in AI chatbots in the public sector.  _Government Information Quarterly_ 37, 4 (Oct. 2020), 101490.  https://doi.org/10.1016/j.giq.2020.101490
  * Arnold and Sörensen (2019) Florian Arnold and Kenneth Sörensen. 2019.  What makes a VRP solution good? The generation of problem-specific knowledge for heuristics.  _Computers & Operations Research_ 106 (June 2019), 280–288.  https://doi.org/10.1016/j.cor.2018.02.007
  * Baluja and Davies (1997) Shumeet Baluja and Scott Davies. 1997.  Using Optimal Dependency-Trees for Combinatorial Optimization: Learning the Structure of the Search Space. In _Proceedings of (ICML) International Conference on Machine Learning_. 30 – 38.
  * Breiman (2001) Leo Breiman. 2001\.  Random Forests.  _Machine Learning_ 45, 1 (2001), 5–32.  https://doi.org/10.1023/A:1010933404324
  * Brownlee (2009) Alexander Brownlee. 2009\.  Multivariate Markov networks for fitness modelling in an estimation of distribution algorithm.  (05 2009).
  * Brownlee (2016) Alexander E.I. Brownlee. 2016\.  Mining Markov Network Surrogates for Value-Added Optimisation. In _Proceedings of the 2016 on Genetic and Evolutionary Computation Conference Companion_. ACM, Denver Colorado USA, 1267–1274.  https://doi.org/10.1145/2908961.2931711
  * Brownlee et al. (2013) Alexander E.I. Brownlee, Olivier Regnier-Coudert, John A.W. McCall, Stewart Massie, and Stefan Stulajter. 2013.  An application of a GA with Markov network surrogate to feature selection.  _International Journal of Systems Science_ 44, 11 (Nov. 2013), 2039–2056.  https://doi.org/10.1080/00207721.2012.684449
  * Brownlee et al. (2018) Alexander E.I. Brownlee, Michal Weiszer, Jun Chen, Stefan Ravizza, John R. Woodward, and Edmund K. Burke. 2018.  A fuzzy approach to addressing uncertainty in Airport Ground Movement optimisation.  _Transportation Research Part C: Emerging Technologies_ 92 (July 2018), 150–175.  https://doi.org/10.1016/j.trc.2018.04.020
  * Brownlee et al. (2015) Alexander E.I. Brownlee, John R. Woodward, and Jerry Swan. 2015\.  Metaheuristic Design Pattern: Surrogate Fitness Functions. In _Proceedings of the Companion Publication of the 2015 Annual Conference on Genetic and Evolutionary Computation_. ACM, Madrid Spain, 1261–1264.  https://doi.org/10.1145/2739482.2768499
  * Brownlee and Wright (2015) Alexander E.I. Brownlee and Jonathan A. Wright. 2015.  Constrained, mixed-integer and multi-objective optimisation of building designs by NSGA-II with fitness approximation.  _Applied Soft Computing_ 33 (Aug. 2015), 114–126.  https://doi.org/10.1016/j.asoc.2015.04.010
  * Brownlee et al. (2007) Alexander E. I. Brownlee, John A. W. McCall, and Deryck F. Brown. 2007.  Solving the MAXSAT problem using a multivariate EDA based on Markov networks. In _Proceedings of the 2007 GECCO conference companion on Genetic and evolutionary computation - GECCO ’07_. ACM Press, London, United Kingdom, 2423.  https://doi.org/10.1145/1274000.1274005
  * Brownlee et al. (2010) Alexander E. I. Brownlee, Olivier Regnier-Coudert, John A. W. McCall, and Stewart Massie. 2010\.  Using a Markov network as a surrogate fitness function in a genetic algorithm. In _IEEE Congress on Evolutionary Computation_. IEEE, Barcelona, Spain, 1–8.  https://doi.org/10.1109/CEC.2010.5586548
  * Deb et al. (2014) Kalyanmoy Deb, Sunith Bandaru, David Greiner, António Gaspar-Cunha, and Cem Celal Tutum. 2014\.  An integrated approach to automated innovization for discovering useful design principles: Case studies from engineering.  _Applied Soft Computing_ 15 (Feb. 2014), 42–56.  https://doi.org/10.1016/j.asoc.2013.10.011
  * Deb and Goldberg (1994) Kalyanmoy Deb and David E. Goldberg. 1994.  Sufficient conditions for deceptive and easy binary functions.  _Annals of Mathematics and Artificial Intelligence_ 10, 4 (Dec. 1994), 385–408.  https://doi.org/10.1007/BF01531277
  * Deb and Srinivasan (2006) Kalyanmoy Deb and Aravind Srinivasan. 2006.  Innovization: innovating design principles through optimization. In _Proceedings of the 8th annual conference on Genetic and evolutionary computation - GECCO ’06_. ACM Press, Seattle, Washington, USA, 1629\.  https://doi.org/10.1145/1143997.1144266
  * Di Pietro et al. (2004) A. Di Pietro, L. While, and L. Barone. 2004.  Applying evolutionary algorithms to problems with noisy, time-consuming fitness functions. In _Proceedings of the 2004 Congress on Evolutionary Computation (IEEE Cat. No.04TH8753)_. IEEE, Portland, OR, USA, 1254–1261.  https://doi.org/10.1109/CEC.2004.1331041
  * Fisher et al. (2018) Aaron Fisher, Cynthia Rudin, and Francesca Dominici. 2018\.  All Models are Wrong, but Many are Useful: Learning a Variable’s Importance by Studying an Entire Class of Prediction Models Simultaneously.  (2018).  https://doi.org/10.48550/ARXIV.1801.01489 Publisher: arXiv Version Number: 5.
  * Gaier et al. (2017) Adam Gaier, Alexander Asteroth, and Jean-Baptiste Mouret. 2017\.  Data-Efficient Exploration, Optimization, and Modeling of Diverse Designs through Surrogate-Assisted Illumination.  (2017).  https://doi.org/10.48550/ARXIV.1702.03713 Publisher: arXiv Version Number: 2.
  * Goldberg et al. (1992) David E Goldberg, Kalyanmoy Deb, and Jeffrey Horn. 1992\.  Massive Multimodality, Deception, and Genetic Algorithms.. In _PPSN_ , Vol. 2.
  * Goodman and Flaxman (2017) Bryce Goodman and Seth Flaxman. 2017.  European Union Regulations on Algorithmic Decision-Making and a “Right to Explanation”.  _AI Magazine_ 38, 3 (Oct. 2017), 50–57.  https://doi.org/10.1609/aimag.v38i3.2741
  * Jin (2005) Y. Jin. 2005.  A comprehensive survey of fitness approximation in evolutionary computation.  _Soft Computing_ 9, 1 (Jan. 2005), 3–12.  https://doi.org/10.1007/s00500-003-0328-5
  * Jin (2011) Yaochu Jin. 2011\.  Surrogate-assisted evolutionary computation: Recent advances and future challenges.  _Swarm and Evolutionary Computation_ 1, 2 (June 2011), 61–70.  https://doi.org/10.1016/j.swevo.2011.05.001
  * Johnson (1974) David S. Johnson. 1974\.  Approximation algorithms for combinatorial problems.  _J. Comput. System Sci._ 9, 3 (Dec. 1974), 256–278.  https://doi.org/10.1016/S0022-0000(74)80044-9
  * Lambrecht and Tucker (2019) Anja Lambrecht and Catherine Tucker. 2019.  Algorithmic Bias? An Empirical Study of Apparent Gender-Based Discrimination in the Display of STEM Career Ads.  _Management Science_ 65, 7 (July 2019), 2966–2981.  https://doi.org/10.1287/mnsc.2018.3093
  * Larrañaga (2002) P. Larrañaga. 2002\.  A Review on Estimation of Distribution Algorithms.  In _Estimation of Distribution Algorithms_ , Pedro Larrañaga and Jose A. Lozano (Eds.). Vol. 2. Springer US, Boston, MA, 57–100.  https://doi.org/10.1007/978-1-4615-1539-5_3 Series Title: Genetic Algorithms and Evolutionary Computation.
  * LeCun et al. (2015) Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. 2015\.  Deep learning.  _Nature_ 521, 7553 (May 2015), 436–444.  https://doi.org/10.1038/nature14539
  * Lehman et al. (2020) Joel Lehman, Jeff Clune, Dusan Misevic, Christoph Adami, Lee Altenberg, Julie Beaulieu, Peter J. Bentley, Samuel Bernard, Guillaume Beslon, David M. Bryson, Nick Cheney, Patryk Chrabaszcz, Antoine Cully, Stephane Doncieux, Fred C. Dyer, Kai Olav Ellefsen, Robert Feldt, Stephan Fischer, Stephanie Forrest, Antoine Fŕenoy, Christian Gagńe, Leni Le Goff, Laura M. Grabowski, Babak Hodjat, Frank Hutter, Laurent Keller, Carole Knibbe, Peter Krcah, Richard E. Lenski, Hod Lipson, Robert MacCurdy, Carlos Maestre, Risto Miikkulainen, Sara Mitri, David E. Moriarty, Jean-Baptiste Mouret, Anh Nguyen, Charles Ofria, Marc Parizeau, David Parsons, Robert T. Pennock, William F. Punch, Thomas S. Ray, Marc Schoenauer, Eric Schulte, Karl Sims, Kenneth O. Stanley, François Taddei, Danesh Tarapore, Simon Thibault, Richard Watson, Westley Weimer, and Jason Yosinski. 2020.  The Surprising Creativity of Digital Evolution: A Collection of Anecdotes from the Evolutionary Computation and Artificial Life Research Communities.  _Artificial Life_ 26, 2 (May 2020), 274–306.  https://doi.org/10.1162/artl_a_00319
  * Lehre and Oliveto (2017) Per Kristian Lehre and Pietro S. Oliveto. 2017.  Theoretical Analysis of Stochastic Search Algorithms.  _arXiv:1709.00890 [cs]_ (Sept. 2017).  http://arxiv.org/abs/1709.00890 arXiv: 1709.00890.
  * Lundberg and Lee (2017) Scott Lundberg and Su-In Lee. 2017.  A Unified Approach to Interpreting Model Predictions.  (2017).  https://doi.org/10.48550/ARXIV.1705.07874 Publisher: arXiv Version Number: 2.
  * Malan (2021) Katherine Mary Malan. 2021\.  A Survey of Advances in Landscape Analysis for Optimisation.  _Algorithms_ 14, 2 (Jan. 2021), 40\.  https://doi.org/10.3390/a14020040
  * Miller (2017) Tim Miller. 2017\.  Explanation in Artificial Intelligence: Insights from the Social Sciences.  (2017).  https://doi.org/10.48550/ARXIV.1706.07269 Publisher: arXiv Version Number: 3.
  * Ochoa et al. (2020) Gabriela Ochoa, Lee A. Christie, Alexander E. Brownlee, and Andrew Hoyle. 2020.  Multi-objective evolutionary design of antibiotic treatments.  _Artificial Intelligence in Medicine_ 102 (Jan. 2020), 101759\.  https://doi.org/10.1016/j.artmed.2019.101759
  * Olhede and Wolfe (2018) S. C. Olhede and P. J. Wolfe. 2018.  The growing ubiquity of algorithms in society: implications, impacts and innovations.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 376, 2128 (Sept. 2018), 20170364.  https://doi.org/10.1098/rsta.2017.0364
  * Rana and Whitley (1998) Soraya Rana and Darrell Whitley. 1998.  Genetic algorithm behavior in the MAXSAT domain.  In _Parallel Problem Solving from Nature — PPSN V_ , G. Goos, J. Hartmanis, J. van Leeuwen, Agoston E. Eiben, Thomas Bäck, Marc Schoenauer, and Hans-Paul Schwefel (Eds.). Vol. 1498. Springer Berlin Heidelberg, Berlin, Heidelberg, 785–794.  https://doi.org/10.1007/BFb0056920 Series Title: Lecture Notes in Computer Science.
  * Rudin (2019) Cynthia Rudin. 2019\.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  _Nature Machine Intelligence_ 1, 5 (May 2019), 206–215.  https://doi.org/10.1038/s42256-019-0048-x
  * Swartout and Moore (1993) William R. Swartout and Johanna D. Moore. 1993.  Explanation in Second Generation Expert Systems.  In _Second Generation Expert Systems_ , Jean-Marc David, Jean-Paul Krivine, and Reid Simmons (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 543–585.  https://doi.org/10.1007/978-3-642-77927-5_24
  * Toreini et al. (2020) Ehsan Toreini, Mhairi Aitken, Kovila Coopamootoo, Karen Elliott, Carlos Gonzalez Zelaya, and Aad van Moorsel. 2020.  The relationship between trust in AI and trustworthy machine learning technologies. In _Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency_. ACM, Barcelona Spain, 272–283.  https://doi.org/10.1145/3351095.3372834
  * Urquhart et al. (2019) Neil Urquhart, Michael Guckert, and Simon Powers. 2019\.  Increasing trust in meta-heuristics by using MAP-elites. In _Proceedings of the Genetic and Evolutionary Computation Conference Companion_. ACM, Prague Czech Republic, 1345–1348.  https://doi.org/10.1145/3319619.3326816
  * Wallace et al. (2021) Aidan Wallace, Alexander E. I. Brownlee, and David Cairns. 2021\.  Towards Explaining Metaheuristic Solution Quality by Data Mining Surrogate Fitness Models for Importance of Variables.  In _Artificial Intelligence XXXVIII_ , Max Bramer and Richard Ellis (Eds.). Vol. 13101. Springer International Publishing, Cham, 58–72.  https://doi.org/10.1007/978-3-030-91100-3_5 Series Title: Lecture Notes in Computer Science.
  * Wei et al. (2015) Pengfei Wei, Zhenzhou Lu, and Jingwen Song. 2015.  Variable importance analysis: A comprehensive review.  _Reliability Engineering & System Safety_ 142 (2015), 399–432.  https://doi.org/10.1016/j.ress.2015.05.018
  * Yapo and Weiss (2018) Adrienne Yapo and Joseph Weiss. 2018.  Ethical Implications of Bias in Machine Learning.  https://doi.org/10.24251/HICSS.2018.668


  *

# Stability impacts from the current and pressure profile modifications within
finite sized island

Yuxiang Sun1 & Di Hu1,∗ 1Beihang University, No. 37 Xueyuan Road, Haidian
District, 100191 Beijing, China.  ∗Corresponding author<EMAIL_ADDRESS>
###### Abstract

The stability (or instability) of finite sized magnetic island could play a
significant role in disruption avoidance or disruption mitigation dynamics.
Especially, various current and pressure profile modifications, such as the
current drive and heating caused by electron cyclotron wave, or the radiative
cooling and current expulsion caused by the Shattered Pellet Injection could
be applied within the island to modify its stability, thus change the ensuing
dynamics. In this study, we calculate the mode structure modification caused
by such profile changes within the island using the perturbed equilibrium
approach, thus obtain the change of stability criterion $\Delta^{\prime}$ and
assess the corresponding quasi-linear island stability. The positive helical
current perturbation is found to always stabilize the island, while the
negative one is found to do the opposite, in agreement with previous results.
The pressure bump or hole within the island has a more complicated stability
impact. In the small island regime, its contribution is monotonic, with
pressure bump tends to stabilize the island while pressure hole destabilizes
it. This effect is relatively weak, though, due to the cancellation of the
pressure term’s odd parity contribution in the second derivatives of the mode
structure. In the large island regime, such cancellation is broken due to the
island asymmetry, and the pressure contribution to stability is manifested,
which is non-monotonic. The stability analysis in this paper helps to more
accurately clarify the expected island response in the presence of profile
modifications caused by disruption avoidance or mitigation systems.

Keywords: Tokamak; Magnetic Island; Tearing Mode; MHD Instability
PACS: 52.55.Fa;52.55.-s;52.35.Py;

##  1 Introduction

Due to the nonlinear growth of global magneto-hydrodynamic (MHD)
instabilities, tokamak plasmas would eventually experience disruptions which
not only interrupt the energy output but also cause irreversible damage to the
device[1]. Such damage during the Thermal Quench (TQ) phase of disruptions is
characterized by a large amount of thermal load on the Plasma Facing
Components (PFCs)[1], while the most feared disruptive consequence in the
Current Quench (CQ) phase is the localized deposition of a large amount of
runaway electrons[2]. To avoid such damages and safe guard the continued
operation of future high performance devices, various Disruption Avoidance
Systems (DAS) are designed, such as the Electron Cyclotron Current Drive
(ECCD) or Heating (ECRH) [3, 4, 5, 6], to suppress the instability of the
dominant magnetic islands by modifying the current and pressure profile
within. Facing this, Disruption Mitigation Systems (DMS) such as the Shattered
Pellet Injection (SPI)[7, 11, 12, 9, 8, 10, 13, 16, 14, 15, 17] or other
Massive Material Injection (MMI)[19, 21, 18, 20, 22, 23, 24, 25] techniques
are used to deplete the majority of thermal energy via radiation loss during
the TQ, and raise the density to help with runaway electron suppression during
the CQ phase.

There have been a lot of recent numerical explorations[27, 11, 30, 28, 29, 31,
32, 26, 33, 34, 35] on how the DMS would interplay with global MHD modes. Such
interplay could strongly impact the mitigation efficiency and these studies
provided valuable insight into the disruption mitigation dynamics. However,
current simulations often investigate the scenario of injection into initially
axisymmetric plasmas, and all MHD response is provoked by the injection
itself. In reality, it is most likely that by the time the DMS is triggered,
there is already unstable long wavelength global MHD modes, such as tearing
modes, existing in the target plasma. The interplay between such finite
amplitude mode and the injected materials is of great interest, since it may
have a significant impact on the density transport and the outgoing heat flux
during the mitigation process. Specifically, the islands could be further
destabilized by the helical profile modification within the O-point [28, 34],
prompting premature triggering of the TQ before the fragments penetrate into
the plasma core. Such a scenario is detrimental to the disruption mitigation
efficiency.

While expensive 3D nonlinear MHD codes could be arranged to explore the
aforementioned stability impact from profile modifications caused by either
the DAS and the DMS, quasi-linear stability analysis using helical Grad-
Shafranov equation [36] provides a more agile and cost-effective tool. In our
simulation environment, which was run on a virtual machine on AMD R9 5900HX
cpu, the typical convergence time ranges from $50ms$ to $300ms$, depending on
the strength of profile modification. With better hardware this could be
further improved. Though it may still not be quick enough for real-time
feedback control, it should be possible to obtain a surrogate model for quick
equilibria solution [38, 37], which could potentially be used to satisfy the
demand of the real-time control. This requires dedicated work and is part of
our future work.

The justification for using such a perturbed equilibrium approach for
stability analysis is that the growth time of resistive modes is usually on
the order of millisecond timescale, much longer than the Alfvénic timescale
which is the perpendicular inertia time for magnetized plasma, so that the
plasma can still be considered as a quasi-equilibrium. The advantage of this
approach is its capability of avoiding the singularity in solutions close to
the resonant surface, and accurately describing the impact of profile
modification within the island on the stability criterion $\Delta^{\prime}$
[36].

Considering a dominant finite sized island in a cylindrical analogue of the
tokamak plasma with a toroidal effect multiplier $1-\frac{m^{2}}{n^{2}}$ [39],
the 2D temperature and pressure profile modification with helical symmetry is
described and their impact on the island stability is obtained through the
calculation of $\Delta^{\prime}$ up to the resonant surface. Several different
combinations of the current density and pressure gradient modification are
considered, and it is found that positive helical current perturbation on
island O-point is always strongly stabilizing while the negative one does the
opposite, consistent with previous analytical [40] and numerical results [28,
34]. Furthermore, the pressure gradient modification is found to have a non-
monotonic impact on the stability depending on the island size. It should be
noted that, in the scope of this paper, we only analysed simple perturbation
of current density and pressure gradient without consideration of specific
physical mechanism. In reality, the current and pressure profile modification
could be coupled together. For example, the bootstrap current modification
means certain current profile modification accompanying any given pressure
profile modification. Self-consistent modelling of such profile modifications
is needed for more comprehensive exploration of the island stability, and this
is an immediate future development we will pursue.

The paper is arranged as the following. In Section 2, we introduce our quasi-
linear perturbed equilibrium equations for the mode structure calculation
considering the current and pressure gradient profile modifications. This will
form the governing equation of our numerical solution. Then in Section 3, we
describe our numerical method, set up our simulation and demonstrate the
stability impact of the aforementioned profile modifications. Finally,
conclusions and discussion are given in Section 4.

##  2 Quasi-linear perturbed equilibrium and profile modification for finite
island

The fundamental approach of the perturbed equilibrium method is already
detailed in Ref. [36], but we explicitly write out the essential parts here
for the completeness of the paper. The helical Grad-Shafranov equation for the
cylindrical plasma is [41]:

$\displaystyle\frac{1}{r}\frac{\partial}{\partial
r}\frac{r}{1+k^{2}r^{2}}\frac{\partial\Psi^{*}}{\partial
r}+\frac{1}{r^{2}}\frac{\partial^{2}\Psi^{*}}{\partial\theta^{2}}$
$\displaystyle=$
$\displaystyle-\left(1-\frac{m^{2}}{n^{2}}\right)\bar{p}^{\prime}(\Psi^{*})$
(1)
$\displaystyle+\frac{2k\bar{F}}{(1+k^{2}r^{2})^{2}}-\frac{\bar{F}\bar{F}^{\prime}(\Psi^{*})}{1+k^{2}r^{2}}.$

Where the $k=n/(mR)$, $\bar{F}\equiv B_{z}+krB_{\theta}$. Here $m$ and $n$ are
the wave number of cylindrical plasma, and $B_{z},B_{\theta}$ means the
magnetic induction intensity in $z,\theta$ direction respectively. Here the
toroidal effect multiplier $\left(1-\frac{m^{2}}{n^{2}}\right)$ is added to
mimic the toroidal coupling effect [39]. The helical flux $\Psi^{*}$ is the
magnetic flux through a helical membrane, and we let
$\Psi^{*}=\Psi^{*}_{0}(r)+\psi$ with $\psi$ being the helical flux
perturbation. We take the current density to be ${\bf\bar{j}}\equiv\mu_{0}{\bf
j}$, $\mu_{0}=4\pi\times 10^{-7}H\cdot m^{-1}$. And
$\bar{p}=\mu_{0}p(\bar{\Psi}^{*})$ is the plasma pressure in the helical Grad-
Shafranov equation. And the right hand side can reflect the helical current
term as
$j_{h}=\bar{F}\bar{F}^{\prime}(\Psi^{*})+(1+k^{2}r^{2})\bar{p}^{\prime}(\Psi^{*})$.

With the long wavelength limit $k^{2}r^{2}\ll 1$ and use the 2D version of the
current and pressure profiles, we can rewrite the Eq. (1) as [36]:

$\displaystyle\frac{1}{r}\frac{\partial}{\partial
r}r\frac{\partial\Psi^{*}}{\partial
r}+\frac{1}{r^{2}}\frac{\partial^{2}\Psi^{*}}{\partial\theta^{2}}=-j^{*}(\Psi^{*})+2k\bar{F}-\left(1-\frac{m^{2}}{n^{2}}\right)k^{2}(r^{2}-r^{2}_{s})P(\Psi^{*}),$
(2) $\displaystyle
j^{*}(\Psi^{*})\equiv\left[1+k^{2}r^{2}_{s}\left(1-\frac{m^{2}}{n^{2}}\right)\right]P(\Psi^{*})+T(\Psi^{*}).$
(3)

We take $P(\Psi^{*})\equiv\frac{dp}{d\Psi^{*}}$ as pressure gradient, and
$T(\Psi^{*})\equiv\frac{1}{2}\frac{d}{d\Psi^{*}}(\bar{F})^{2}$. Here the 2D
profile provides a more accurate marginal stability boundary compared with
results obtained from simple 1D perturbed profile, with the 1D profile
predicting instability while the 2D island is actually marginal stable [36].
Combining Eq. (2) and Eq. (3), the helical magnetic flux $\Psi^{*}$ can be
calculated for given boundary conditions once we know the current and pressure
gradient profile as a function of $\Psi^{*}$. We can see from these equations
that the current density perturbation will directly affect the derivative of
the mode structure and thus change the instability criterion
$\Delta^{\prime}$. But the pressure gradient perturbation has to enter the
right hand side together with the $(r^{2}-r_{s}^{2})$ term. In a symmetric
island, the contribution of the pressure gradient should thus be symmetric on
both sides of the resonant surface, so that such a symmetrical island takes no
effect from this term in the calculation of $\Delta^{\prime}$. And only with
the in-out asymmetry of the finite island could this term have an impact on
the stability.

Separating the zeroth and the principle harmonics, we could write down the
perturbed equilibrium equations[36]:

$\displaystyle\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial}{\partial r}\left<\Psi^{*}\right>_{0}\right)$
$\displaystyle=$
$\displaystyle-\left<\left[1+k^{2}r_{s}^{2}\left(1-\frac{m^{2}}{n^{2}}\right)\right]P+T\right>_{0}$
(4)
$\displaystyle-k^{2}\left(r^{2}-r_{s}^{2}\right)\left(1-\frac{m^{2}}{n^{2}}\right)\left<P\right>_{0}+2kB_{z},$
$\displaystyle\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial}{\partial r}\psi\right)$ $\displaystyle=$
$\displaystyle-\left<\left[1+k^{2}r_{s}^{2}\left(1-\frac{m^{2}}{n^{2}}\right)\right]P+T\right>_{1}$
(5)
$\displaystyle-k^{2}\left(r^{2}-r_{s}^{2}\right)\left(1-\frac{m^{2}}{n^{2}}\right)\left<P\right>_{1}+\frac{m^{2}}{r^{2}}\psi.$

Here, $\left<f\right>_{0}$ represents the zeroth harmonic component of
function $f$, while $\left<f\right>_{1}$ represents the principle harmonic
$m,n$. We are only considering a single helicity $m/n$ here, since our
equilibrium approach requires a given symmetric direction which is determined
by the helicity. Although the multi-helicity problem is important in the
triggering of the TQ, the pre-TQ phase often see a single helicity dominant
mode which nonlinearly drives the other helicity modes before the TQ
onset[42]. If we wish to retain the nonlinear coupling by these subdominant
multi-helicity modes in the future, we should retain the principle harmonic of
the nonlinear term involving the other helicity modes in the right hand side
of Eq. (5), this is a potential future development of this work. The above
equations could be re-written as

$\displaystyle\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial}{\partial
r}\left<\Psi^{*}\right>_{0}\right)=-\left<j^{*}\right>_{0}-\left<\delta
j^{*}\right>_{0}+2kB_{z},$ (6)
$\displaystyle\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial}{\partial
r}\psi\right)=-\left<j^{*}\right>_{1}-\left<\delta
j^{*}\right>_{1}+\frac{m^{2}}{r^{2}}\psi.$ (7)

Here we defined

$\displaystyle
j^{*}\left(\Psi^{*}\right)\equiv\left[1+k^{2}r_{s}^{2}\left(1-\frac{m^{2}}{n^{2}}\right)\right]P\left(\Psi^{*}\right)+T\left(\Psi^{*}\right),$
(8) $\displaystyle\delta j^{*}\equiv
k^{2}\left(r^{2}-r_{s}^{2}\right)\left(1-\frac{m^{2}}{n^{2}}\right)P\left(\Psi^{*}\right).$
(9)

With a suitable description of the 2D helical modification of the helical
current density $j^{*}$ and the pressure gradient $P$, all terms on the Right
Hand Side (RHS) of Eq. (6) and Eq. (7) are known, we can then obtain the
perturbed helical flux $\psi$ and the corresponding $\Delta^{\prime}$ from
given boundary condition.

The aforementioned profile modifications can be expressed as the following.
Since both the helical current density and the pressure gradient are still
flux functions, we could first obtain the description of the helical flux near
the resonant surface, then extract the corresponding profiles’ zeroth and
principle harmonics, which will be used in the RHS terms of Eq. (6) and Eq.
(7) respectively. The flux function near the resonant point can be represented
by the following equation[36].

$\displaystyle\Psi^{*}(x,\theta)=\Psi^{*}_{0s}+\frac{1}{2}\Psi^{\prime\prime}_{0s}x^{2}+(\psi_{s}+\psi^{\prime}_{s}x)\cos
m\theta,$ (10)

where $x$ is the radial distance to the resonant surface, and the subscript
‘$s$’ means the value on the resonant surface. The island semi-width _w_ and
asymmetry _s_ are:

$\displaystyle w=2\sqrt{-\frac{\psi_{s}}{\Psi^{\prime\prime}_{0s}}},\qquad
s=\frac{\psi^{\prime}_{s}}{\Psi^{\prime\prime}_{0s}}.$ (11)

We use $x_{Sl}$ and $x_{Sr}$ denote the left and right boundaries of the
separatrix:

$\displaystyle x_{Sl}(\theta)$ $\displaystyle=-x_{S0}-s\cos m\theta,$ (12)
$\displaystyle x_{Sr}(\theta)$ $\displaystyle=x_{S0}-s\cos m\theta,$ (13)
$\displaystyle x_{S0}$ $\displaystyle\equiv|w\cos^{2}{\frac{m\theta}{2}}|.$
(14)

Those equations show that $w$ is the island semi-width, while $s$ induces a
periodic asymmetry for the island along the poloidal direction. It’s obvious
that the asymmetry results in the outward shift of X-point and inward shift of
O-point of the island, causing a corresponding asymmetry in the
$(r^{2}-r_{s}^{2})P$ terms discussed above.

With the island geometry known, we introduce the following flux-coordinate
with the dimension of length squared:

$\displaystyle\chi(\Psi)\equiv\frac{2(\Psi-\psi_{s})}{\Psi^{\prime\prime}_{0}}=x^{2}-x^{2}_{S0}+s^{2}+2xs\cos
m\theta.$ (15)

So that the current density profile can be represented as [36]:

$\displaystyle
j(\chi)=\left\\{\begin{array}[]{ll}j_{0}-j^{\prime}_{l}(\sqrt{\chi+d^{2}_{l}}-d_{l})&+\frac{1}{2}j^{\prime\prime}_{l}\chi,\quad(x<x_{Sl})\\\
j_{0}&+\frac{1}{2}j^{\prime\prime}_{i}\chi,\quad(x_{Sl}<x<x_{Sr})\\\
j_{0}+j^{\prime}_{r}(\sqrt{\chi+d^{2}_{r}}-d_{r})&+\frac{1}{2}j^{\prime\prime}_{r}\chi,\quad(x>x_{Sr})\\\
\end{array}\right.$ (19)

where $w_{l}\equiv w+s=-x_{Sl}(0)$, $w_{r}\equiv w-s=x_{Sr}(0)$, and
$d_{l}\equiv w_{l}+s$, $d_{r}\equiv w_{r}-s$. Like-wise the pressure gradient
could be described as the following:

$\displaystyle
P(\chi)=\left\\{\begin{array}[]{ll}P_{0}-P^{\prime}_{l}(\sqrt{\chi+d^{2}_{l}}-d_{l})&+\frac{1}{2}P^{\prime\prime}_{l}\chi,\quad(x<x_{Sl})\\\
P_{0}&+\frac{1}{2}P^{\prime\prime}_{i}\chi,\quad(x_{Sl}<x<x_{Sr})\\\
P_{0}+P^{\prime}_{r}(\sqrt{\chi+d^{2}_{r}}-d_{r})&+\frac{1}{2}P^{\prime\prime}_{r}\chi,\quad(x>x_{Sr})\\\
\end{array}\right.$ (23)

Here we name the different zones of the MHD by the following way. The area
near the resonant surface $x_{S0}$ is a magnetic island as the magnetic lines
are reconnected to form a shape of island. Then the zones at $-d_{l}<x<-w_{l}$
and $w_{r}<x<d_{r}$ are called buffer zones for they are perturbed by the
magnetic island. And the zones that are outside of the 2D perturbed area are
called unperturbed linear zones.

In Eq. (19) the parameter
$j^{\prime}_{l},j^{\prime\prime}_{l},j^{\prime}_{r},j^{\prime\prime}_{r}$ are
set to match the value and derivative of unperturbed j profile at the edge of
the buffer zone. Besides, thanks to the fact that conductivity in the plasma
is quite high, current perturbations can be limited close to the island
region. So we also choose to keep the total current in quasi-linear zone the
same as the unperturbed results. And the parameter $P^{\prime}_{l}$,
$P^{\prime\prime}_{l}$, $P^{\prime}_{r}$, $P^{\prime\prime}_{r}$ in Eq. (23)
are also iterated to satisfy the matching condition. Here, we consider the
cooling or heating timescale to be fast compared with the global transport
time, hence the pressure profiles outside of the buffer zone would not be
affected significantly. A typical zeroth harmonic result of such profile
modification is shown in Fig. 1, represented by the green solid lines.

In Ref [36], the second-order gradient $j^{\prime\prime}_{i}$ and
$P^{\prime\prime}_{i}$ were set to zero to represent the profile flattening
within the island. As discussed in Section 1, many physical processes could
result in a “bump” or “hole” on the O-point for the corresponding profiles.
For example, the radiative cooling from impurity SPI into the island could
result in both current holes and pressure holes on the O-point [16, 40], while
the electron cyclotron wave injection (like ECRH) could result in bumps within
the island for both profiles [3, 4, 5, 6]. More bizarre cases could exist, for
example during hydrogen isotope injections, where the dilution cooling
resulted in current expulsion from the island interior but didn’t deplete the
pressure due to the low radiation power. Indeed, there may even exist a minor
pressure peak within the island due to the heat transfer into the low
temperature, but high density, island region.

Some examples of profile modifications with selected values of
$j^{\prime\prime}_{i}$ and $P^{\prime\prime}_{i}$ are shown in Fig. 1 for the
zeroth harmonic which is calculated with the setting as semi-width $w=0.08m$
and initial parameters $f_{1}=0.5,f_{2}=0.35$. In Fig. 1(a), the unperturbed
current density profile is shown in purple, a flattened profile is shown in
green and a current hole case is shown in blue. Likewise, the unperturbed,
flattened and pressure hole cases are represented by the same colours in Fig.
1(b). The red vertical dashed lines represent the location of the resonant
surface, the blue dashed lines represent the maximum island width, and the
green dashed lines are the 1D buffer region away from the island which is
impacted by the 2D island profile modification.

|
---|---
(a) | (b)
Figure 1: The current density and pressure of 2/1 mode with setting island
semi-width to $w=0.08m$ and initial parameters to $f_{1}=0.5,f_{2}=0.35$. (a)
Different $j^{\prime\prime}_{i}$ can change the zeroth harmony of current
density greatly. The green lines are the profile of the current density which
is calculated by ideal theory. the red dash lines and blue lines are
calculated by the quasi-linear theory, while the former set
$j^{\prime\prime}_{2}$ to zero and the latter one set $j^{\prime\prime}_{2}$
to 50. (b) As the former, shows the effect of different $P^{\prime\prime}_{i}$
at zeroth harmony, when the general pressure is set as 0.239 MPa. Vertical
lines indicate the different zone of the island.

We wish to remind the readers that these profiles are merely some examples of
possible profile modification caused by material injection or external drives.
The exact shapes of these profile modifications are subjected to detailed
analysis of heating, cooling, current drive and transport processes, which are
beyond the scope of this paper. In the rest of the paper, we will perform
scans on the strength of such profile modifications without delving too much
into the feasibility of such modifications, as our focus is on stability
analysis.

With the above description, we can finally solve the perturbed equilibrium
equation from both the magnetic axis and the plasma boundary, up to both sides
of the resonant surface and obtain the stability criterion

$\displaystyle\Delta^{\prime}=\left.\frac{\psi^{\prime}}{\psi}\right|_{r=r_{s}-\epsilon}^{r=r_{s}+\epsilon}=\frac{\psi^{\prime}(r_{s}+\epsilon)-\psi^{\prime}(r_{s}-\epsilon)}{\psi},\quad\epsilon\rightarrow
0$ (24)

We can solve the mode structure up to the resonant surface in this approach
since we avoid the singularity in the mode structure solution in the perturbed
equilibrium method [36]. We wish to emphasize here that we are calculating
$\Delta^{\prime}$ on both sides of the reconnecting resonant surface, instead
of both sides of the island boundary as conventional wisdom held [43, 40].
Such a $\Delta^{\prime}$ can be used as an indicator of the island stability.
But to obtain the growth rate of the island, one may have to further carry out
the flux-average of Ohm’s law carefully, which is beyond the scope of this
paper.

As a further note, in toroidal geometry, one should also subtract the Pfirsh-
Schluter current contribution [44] to obtain the accurate island stability.
Since we are considering $\Delta^{\prime}$ on both side of the resonant
surface, the flux-average of the Pfirsh-Schluter current on the separatrix
must be used instead of the form given in Ref. [44]. In this study, however,
due to the large major radius we chose, this additional current contribution
is small [45]. Hence we simply neglect this additional contribution here, and
we are satisfied with merely obtaining $\Delta^{\prime}$ as the indication of
island stability. We will take this term into calculation in future researches
to analyse the growth rate for higher inverse aspect-ration tokamak
configuration.

##  3 Stability impact from the current density and pressure profile
modifications

###  3.1 Numerical setup

We consider a simple unperturbed equilibrium profile with the following form:

$\displaystyle\begin{array}[]{rl}j(r)=&j(0)(1-f_{2})\left[(1-f_{1})\left(1-\frac{r^{2}}{a^{2}}\right)+\left(1-\frac{r^{2}}{a^{2}}\right)^{2}\right]+\\\
&j(0)f_{2}\left(1-\frac{r^{2}}{a^{2}}\right)^{3},\\\
p(r)=&p(0)\left(1-\frac{r^{2}}{a^{2}}\right),\end{array}$ (28)

where $a$ is the minor radius. The parameter $f_{1}$ and $f_{2}$ are given
manually to produce an unperturbed profile, upon which the finite island
profile modification would occur. Besides those initial equilibrium profiles,
there are some variables like the island semi-width $w$ and perturbation
parameter $j^{\prime\prime}_{i}$ and $P^{\prime\prime}_{i}$, which we will
vary in our numerical solution to find out their impact on the stability. In
our numerical solution, we also choose some constant parameter values as is
shown in Table 1.

Variable | Value | Physics Meaning
---|---|---
B (T) | 10 | Vacuum external magnetic field
a (m) | 1 | Minor radius
R (m) | 10 | Major radius
Table 1: The initialized parameters for the quasi-linear simulation in this
study.

We are considering cases where there are dominant single helicity modes in the
study. As an example, we investigate the $2/1$ mode and the $3/1$ mode each
with a set of current profiles respectively. The $f_{1}$ and $f_{2}$ will be
given for each case. The safety factor for each current profile can be
obtained by:

$\displaystyle q(r)=\frac{rB_{\varphi}}{R_{0}B_{\theta}}=\frac{2\pi
r^{2}B_{\varphi}}{\mu_{0}I(r)R},$ (29) $\displaystyle
I(r)=2\pi\int_{0}^{r}j(r^{\prime})r^{\prime}dr^{\prime}.$ (30)

The current density, pressure and $q$ profiles for each case are shown in Fig.
2 for the two cases respectively. The initial parameters of 2/1 mode Fig. 2(a)
is $f_{1}=0.5,f_{2}=0.35$, edge safety factor $q(a)=3$ and the initial
parameters of 3/1 mode Fig. 2(b) is $f_{1}=-0.30,f_{2}=-0.20$, edge safety
factor $q(a)=3.95$.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 2: Use Eq.28 to calculate the initial profile used for our simulation,
and use Eq.29 to calculate the safety factor in island. (a) and (b) are
structures of $m=2,n=1$, initial parameter in Eq. 28 $f_{1}=0.5,f_{2}=0.35$,
edge safety factor $q(a)=3$; (c) and (d) are structures of $m=3,n=1$, initial
parameter in Eq. 28 $f_{1}=-0.30,f_{2}=-0.20$, edge safety factor $q(a)=3.95$.

Since Eq. (19) and Eq. (23) are 2D in nature, we have to extract their zeroth
and principle harmonics to be used in the RHS of Eq. (6) and Eq. (7). We begin
with the current density profile. For regions outside of the island but within
the “buffer” region as is marked by the green vertical dashed lines in Fig. 1,
we have:

$\displaystyle j^{(0)}(x)\equiv\frac{1}{2\pi}\int_{-\pi}^{\pi}j(\chi)d\theta,$
(31)

So the current density outside the island becomes:

$\displaystyle
j^{(0)}(x)=\left\\{\begin{array}[]{ll}j_{0}-j^{\prime}_{l}X_{1l}+\frac{1}{2}j^{\prime\prime}_{l}X^{2}_{2}\qquad(x<x_{Sl})\\\
j_{0}+j^{\prime}_{r}X_{1r}+\frac{1}{2}j^{\prime\prime}_{r}X^{2}_{2}\qquad(x>x_{Sr})\\\
\end{array}\right.$ (34)

Here we have defined:

$\displaystyle
X_{1l}=-d_{l}+\frac{2}{\pi}\sqrt{d_{l}^{2}+(x-s)^{2}}E\left(\frac{w^{2}-4sx}{d_{l}^{2}+(x-s)^{2}}\right),$
(35) $\displaystyle
X_{1l}=-d_{r}+\frac{2}{\pi}\sqrt{d_{r}^{2}+(x-s)^{2}}E\left(\frac{w^{2}-4sx}{d_{r}^{2}+(x-s)^{2}}\right),$
(36) $\displaystyle X_{2}^{2}=\frac{1}{2}(2x^{2}+2s^{2}-w^{2}).$ (37)

Here $E(x)$ is the complete elliptic integral of the second kind. As those
equations are functions of $x$, we can calculate the zeroth harmonic of j
profile for any given radius. The perturbed current density in this buffer
region is then calculated by taking the derivative of the zeroth harmonic
current density just like in linear theory.

The profile modifications within the island region as is marked by the blue
dashed lines in Fig. 1 is more complicated due to its 2D nature. We have to
take into account both the “outside island region” and the “inside island
region” expressions in Eq. (19). We define $t\equiv|x-s|/d_{l}$, then Eq. (19)
could be re-written as:

$\displaystyle
j_{-}(t,\theta)=j_{0}-j^{\prime}_{l}d_{l}F_{1}+\frac{1}{2}d^{2}_{l}j^{\prime\prime}_{l}F_{2}+\frac{1}{2}j^{\prime\prime}_{i}F_{3}w^{2},$
(38) $\displaystyle
j_{+}(t,\theta)=j_{0}+j^{\prime}_{r}d_{r}F_{1}+\frac{1}{2}d^{2}_{r}j^{\prime\prime}_{r}F_{2}+\frac{1}{2}j^{\prime\prime}_{i}F_{3}w^{2}.$
(39)

where

$\displaystyle H(t,\theta)\equiv
1,\quad\left(t>\left|\cos\left(\frac{m\theta}{2}\right)\right|\right);$ (40)
$\displaystyle H(t,\theta)\equiv
0,\quad\left(t<\left|\cos\left(\frac{m\theta}{2}\right)\right|\right);\qquad
I\equiv 1-H;$ (41) $\displaystyle
F_{1}(t,\theta)\equiv\left(\sqrt{t^{2}+\sin^{2}{\frac{m\theta}{2}}-1}\right)H(t,\theta);$
(42) $\displaystyle
F_{2}(t,\theta)\equiv\left(t^{2}-\cos^{2}{\frac{m\theta}{2}}\right)H(t,\theta),$
(43) $\displaystyle
F_{3}(t,\theta)\equiv\left(t^{2}-\cos^{2}{\frac{m\theta}{2}}\right)I(t,\theta)$
(44)

It’s obvious that the Eq. (38) and Eq. (39) are functions of both radius $x$
and poloidal angle $\theta$ as mentioned above, to use them in Eq. (6) and Eq.
(7), one have to extract the zeroth and principle Fourier harmonics
respectively.

The pressure gradient profile can be modified like-wise. As we have discussed
above, we consider the pressure increase or decrease within the island to have
a faster timescale compared with that of the global transport, so that the
pressure bump or hole is located within the island and its immediate vicinity.
Once the pressure gradient profile $P(\chi)$ is extracted from Eq. (23), the
perturbed flux could be solved. The pressure profile itself can be obtained by
simply taking the integration from the boundary condition $p(a)=0$:

$\displaystyle dp=Pd\Psi^{*}=Pdx\cdot\frac{d\Psi^{*}}{dx},$ (45)
$\displaystyle\left<p\right>_{0}=\int{\left<P\right>_{0}\left<\frac{d\Psi^{*}}{dx}\right>_{0}dx}+\int{\left<P\right>_{1}^{*}\frac{d\psi}{dx}dx}.$
(46)

Here $\left<P\right>_{1}^{*}$ marks the conjugation of the principle harmonic
of the pressure gradient. The contribution from the
$\left<P\right>_{1}^{*}\frac{d\psi}{dx}$ term is small compared with the
average field term for our parameters considered, however. An example of the
final pressure profile thus obtained has been shown in Fig. 1.

###  3.2 Impact from current density profile modification within the island

(a)  (b) (c)

Figure 3: An example of the 2/1 mode structure for an island semi-wide
$w=0.0801m$, initial parameter $f_{1}=0.5,f_{2}=0.35$ and without pressure
ingredient. (a) The unperturbed current density calculated by a quasi-linear
way with five vertical lines indicating the different zones of the island. (b)
the helical flux $\psi$, and (c) its radial gradient $\partial_{r}\psi$.

In this subsection, we first demonstrate the stability impact from only the
current profile modification, and show that our perturbed equilibrium result
is consistent with previous finite island analyses using constant $\psi$
assumptions on island boundary. Fig.3 shows an example profile we got by our
code which has the island semi-width $w=0.0801m$ and no pressure. Fig.3(a)
presents the current density which is flattened within the magnetic island. As
the figure indicated, the blue dashed lines mark the edge of the magnetic
island, and the green dashed lines divide the area where we need to consider
the 2D perturbations, buffer zones are the areas which between the blue lines
and green lines. While the red dashed line in the middle shows the place of
the resonant surface. Those vertical dashed lines have the same meaning in the
following figures. What’s more, Fig.3(b) and (c) show the calculated magnetic
flux and its first derivative. A notable feature is that the perturbed flux is
not constant-$\psi$ across the island region. The jump in the
$\psi_{r}^{\prime}$ corresponds to the stability criterion $\Delta^{\prime}$.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 4: The current density of 2/1 mode with setting island semi-width to
$w=0.08m$ and initial parameter $f_{1}=0.5,f_{2}=0.35$. The five vertical
lines still indicate the five zones of the island, and the green line in each
figure is the current density for perturbation is zero and the purple line is
for the profile calculated with (a) perturbation $j_{i}^{\prime\prime}$ to
50;(c) set the perturbation $j_{i}^{\prime\prime}$ to -50. (b) and (d) is the
$\partial_{r}\psi$ for the structures (a) and (c). (In figures (a) and (c),
the red lines mark the resonant surface, the blue lines divide the magnetic
island and the green lines are the edge of the buffer zones. In figures (b)
and (d) the black lines mark the resonant surface.)

To demonstrate the effect of the current bump or hole, we show two further
examples in Fig. 4. As we can see, if we give a negative helical current
perturbation to the island, the $\Delta^{\prime}$ increases, which means the
instability grows. On the other hand, if we give a positive perturbation, the
$\Delta^{\prime}$ will decrease to negative and indicates the island becomes
more stable. This is consistent with previous findings in both the ECCD island
stabilizing studies and the cooling island destabilizing studies[16, 40, 3,
4], that a positive helical current perturbation on the island O-point tends
to stabilize the island, while a negative helical current perturbation on the
O-point tends to do the opposite.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 5: Value of $\Delta^{\prime}$ vary by semi-width $w$ and perturbations
$j_{i}^{\prime\prime}$. (a) and (b) Results of 2/1 mode with edge safety
factor $q(a)=3$ and initial parameter $f_{1}=0.5,f_{2}=0.35$; (c) and (d)
Results of 3/1 mode with edge safety factor $q(a)=4.5$ and initial parameter
$f_{1}=0.5,f_{2}=0.35$

A more comprehensive look at the helical current perturbations impact can be
obtained from a parameter scan of island semi-width and current perturbation
amplitude, the result of which is shown in Fig.5, assuming no pressure
contribution. In Fig.5(a) & (c), we show a scan of current perturbation
strength with different fixed island width. First, as the figure indicates,
the stronger the perturbation, the bigger the impact to $\Delta^{\prime}$.
Consistent with the results in Fig. 4, the helical current bump (negative
$j^{\prime\prime}_{i}$) results in a decrease of $\Delta^{\prime}$ and
stabilizes the island, while the helical current hole does the opposite.
Second, the effect of the current bump or hole depends on the island width, so
that for the same $j^{\prime\prime}_{i}$, the wider the island, the stronger
the effect. This is understandable as with a same $j^{\prime\prime}_{i}$ the
more current are perturbed in the larger island, thus the stronger impact.
Such a trend is common for both modes. In Fig.5(b) & (d), these trends are
further confirmed by the contour of $\Delta^{\prime}$ for different island
semi-width $w$ and $j^{\prime\prime}_{i}$. In these figures, the marginal
stable boundary $\Delta^{\prime}=0$ are marked by a solid line, while the
other contours are in dashed lines.

|
---|---
(a) | (b)

(c)

Figure 6: (a) Value of $\Delta^{\prime}$ vary by semi-width and perturbations
$j_{i}^{\prime\prime}$, with given general pressure $p_{0}=0.239$MPa, for 2/1
mode and edge safety factor $q(a)=3$ and initial parameter
$f_{1}=0.5,f_{2}=0.35$. (b) A line of $j_{i}^{\prime\prime}=0$ shows that with
p(0) the $\Delta^{\prime}$ is not monotonous. (c)$\Delta^{\prime}$ changes
with the general pressure and semi-width of the island at no current density
or pressure gradient perturbation.

Adding finite pressure results in the stabilization of small islands due to
the Glasser effect [46] as is shown in Fig. 6, where a simple flattened
pressure profile within the island is considered and we use $m/n=2/1$, island
semi-wide $w=0.0801m$, the initial parameter $f_{1}=0.5,f_{2}=0.35$ to set the
mode. Comparing Fig. 6(a) and Fig.5(b), it is apparent that a strong
stabilizing effect exists for the small island limit, which is caused by the
mode structure modification due to finite pressure gradient near the resonant
surface in toroidal geometry [36]. Such pressure contribution vanishes in the
large island limit, as is shown in Fig. 6(a) & (b), resulting in a linearly
stable but non-linearly unstable scenario. As the pressure increases the
island becomes more and more stable, even with a finite island there comes a
threshold in $p(0)$ beyond which there is no unstable island, as is shown in
Fig. 6(c). Note that here the current perturbation from the neo-classical
effect is not included, otherwise the island would still be unstable from the
decrease of the bootstrap current.

###  3.3 Additional contribution from the pressure profile modification

|
---|---
(a) | (b)
Figure 7: The pressure of the $2/1$ mode with setting island semi-width to 0.0801m and initial parameter $f_{1}=0.5,f_{2}=0.35$. The five vertical lines still indicate the five zones of the island, and the green line in each figure is the pressure gradient for perturbation is zero and the purple line is for the profile calculated with (a) perturbation $P_{i}^{\prime\prime}$ to 50000; (b) set the perturbation $P_{i}^{\prime\prime}$ to -50000.  |
---|---
(a) | (b)
|
(c) | (d)
|
(e) | (f)
Figure 8: Value of $\Delta^{\prime}$ vary by semi-width $w$ and perturbations
$P_{i}^{\prime\prime}$. (a) and (b) are results of 2/1 mode with edge safety
factor $q(a)=3$ and initial parameter $f_{1}=0.5,f_{2}=0.35$; (c) and (d) are
results of 2/1 mode with edge safety factor $q(a)=3$ and initial parameter
$f_{1}=-0.30,f_{2}=-0.20$; (e) and (f) are results of 3/1 mode with edge
safety factor $q(a)=3.95$ and initial parameter $f_{1}=-0.30,f_{2}=-0.20$.

We now consider the stability impact from pressure bumps or holes within the
island. To this purpose, we consider a simple flattened current density
profile within the island but finite $P_{i}^{\prime\prime}$ exists. An example
of such pressure modification is shown in Fig.7 where Fig.7(a) and (b) are the
pressure hole and bump case respectively for the $2/1$ mode case. The
stability criterion as a function for various island semi-width $w$ and
$P^{\prime\prime}_{i}$ is shown in Fig.8 for both the $2/1$ mode case with
$q(a)=3$ and the $3/1$ mode case with $q(a)=4.5$. Here, Fig.8(a) & (b) show
the scan for $2/1$ mode with $f_{1}=0.5$, $f_{2}=0.35$, Fig.8(c) & (d) show
that for $2/1$ mode with $f_{1}=-0.3$, $f_{2}=-0.2$ and Fig.8(e) & (f) show
that for $3/1$ mode with $f_{1}=-0.3$, $f_{2}=-0.2$. In Fig.8(a), (c) & (e),
the pressure modification effect for a series of fixed island semi-width are
shown. In the small island limit, the pressure bump (negative
$P^{\prime\prime}_{i}$) is found to stabilize the island while the pressure
hole (positive $P^{\prime\prime}_{i}$) is found to destabilize the island,
although their impact is weak due to the cancellation between the odd parity
contribution of the $\left(r^{2}-r_{s}^{2}\right)P$ term discussed in Section
2. In the large island limit, the pressure impact becomes non-monotonic. Such
complicated behaviour could be because of the mode structure, thus its
gradient and consequently the island asymmetry, experiencing nonlinear
interplay with the pressure profile modification, as the latter depends on the
former while also affecting the former. Further, the scan of the island size
and the pressure modification for the same cases are shown in Fig.8(b), (d) &
(f). In all cases, a “ridge-like” pattern can be seen going from small island
to large island, so that the finite sized island could be more unstable
compared with linear islands until they tend to saturate. It can also be seen
that the $\Delta^{\prime}$ contour at the small island limit is almost
unaffected by the change of $P^{\prime\prime}_{i}$ while at the large island
limits the effect is more obvious, confirming our previous conjecture that the
island asymmetry at the large island limit tends to produce stronger pressure
gradient effect. Fig.8(b), (d) & (f) also confirm that in the large island
limit, the pressure modification within the island region no longer shows a
monotonic impact on stability. Instead, too large pressure bumps or holes
could both result in a more stable island.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 9: Change of $\Delta^{\prime}$ figure at given current density
perturbation $j_{i}^{\prime\prime}$ in the meantime for 2/1 mode and edge
safety factor $q(a)=3$, initial parameter $f_{1}=0.5,f_{2}=0.35$, for (a) the
$j_{i}^{\prime\prime}=-5$, (b) the $j_{i}^{\prime\prime}=-20$, (c) the
$j_{i}^{\prime\prime}=+10$ and (d) the $j_{i}^{\prime\prime}=+20$.

We now proceed to combine the current and the pressure modifications within
the island. We consider the island stability with several selected
$j_{i}^{\prime\prime}$ and calculate $\Delta^{\prime}$ for various island
semi-width $w$ and pressure modification $P_{i}^{\prime\prime}$. The results
of 2/1 mode were plotted as Fig.9. We wish to remind the readers again that,
although $j^{\prime\prime}_{i}$ is constant for each plot, the current
perturbation would increase in strength as the island grows larger, so that we
will see in Fig.9 the instability could “runaway” for strong negative helical
current perturbation (positive $j^{\prime\prime}_{i}$) as the island grows. In
Fig.9(a) & (b), different values of negative $j^{\prime\prime}_{i}$ are
chosen, corresponding to positive helical current (negative
$j^{\prime\prime}_{i}$) within the island. In those stabilized cases, the
pressure bump (negative $P^{\prime\prime}_{i}$) seems to further stabilize the
island while the pressure hole does the opposite. For the negative helical
current case (positive $j^{\prime\prime}_{i}$) shown in Fig.9(c) & (d), the
pressure bump and hole still show a similar impact as is in Fig.9(a) & (b) for
small island width. For the large island width, however, their influence is
more complicated and non-monotonic, as both strong pressure bumps and holes
could tend to stabilize the island. Such behavior is especially apparent in
the large $P^{\prime\prime}_{i}$, large island regime in Fig.9(d). We wish to
emphasize again that the bootstrap current contribution should be included
within the current modification $j^{\prime\prime}_{i}$ in the above analysis,
so that any pressure contribution discussed here is in addition to its impact
on the current profile.

|
---|---
(a) | (b)
Figure 10: Change of $\Delta^{\prime}$ figure at given current density
perturbation $ji^{\prime\prime}$ in the meantime for 3/1 mode and edge safety
factor $q(a)=4.5$, initial parameter $f_{1}=-0.30,f_{2}=-0.20$, for (a) the
$ji^{\prime\prime}=+10$ and (b) the $ji^{\prime\prime}=+20$.

A similar scan of the $3/1$ mode is plotted in Fig.10. Since the 3/1 mode is
stable enough, we don’t consider applying negative $j_{i}^{\prime\prime}$,
which only decrease $\Delta^{\prime}$. With small island width, the $3/1$ mode
shares the same trend of destabilizing the island with increasing
$P^{\prime\prime}_{i}$ while the same non-monotonic feature appears in the
large island limit, as can be seen by comparing Fig.10(a) & (b) with Fig.9(c)
& (d). Similar “ridge-like” feature in $\Delta^{\prime}$ exists for the small
current hole case ($j^{\prime\prime}_{i}=10$), as can be seen by comparing
Fig.10(a) and Fig.9(c). With a stronger current hole
($j^{\prime\prime}_{i}=20$), the large $P^{\prime\prime}_{i}$ is found to
stabilize the island in the large island limit for the $3/1$ mode, consistent
with the $2/1$ mode trend. We thus anticipate this to be a common feature
that, for a large island, a significant pressure hole would tend to mitigate
the destabilizing effect of a current hole within the island, while it
contributes to instability in the small island cases.

|
---|---
(a) | (b)
Figure 11: Change of $\Delta^{\prime}$ figure at given semi-width in the
meantime for 2/1 mode and edge safety factor $q(a)=3$, initial parameter
$f_{1}=0.5,f_{2}=0.35$, for (a) the semi-width $w=0.0801m$ and (b) the semi-
width $w=0.0253m$.

The above mentioned trend could be more directly shown by considering a scan
of different $j^{\prime\prime}_{i}$ and $P^{\prime\prime}_{i}$ with constant
island width. We choose two different island semi-width $w=0.0235m$ and
$w=0.0801m$, the results are plotted in Fig.11. For the large island case
shown in Fig.11(a), the $P^{\prime\prime}_{i}$ stability impact is once again
seen as non-monotonic, as large positive $P^{\prime\prime}_{i}$ is seen to
mitigate the destabilizing effect of a helical current hole (positive
$j^{\prime\prime}_{i}$), and large negative $P^{\prime\prime}_{i}$ is also
found to stabilize the island. For the smaller island case shown in Fig.11(b),
the $P^{\prime\prime}_{i}$ contribution to island stability is mostly
monotonic, with pressure bump consistently stabilizing the island and pressure
hole consistently destabilizing the island. These different behaviours in
different island width regimes shown in Fig.11 are consistent with what we
discussed above regarding the pressure profile modification stability
influence.

##  4 Conclusion and discussion

The stability analysis of finite size island with current and pressure profile
modifications near the O-point is carried out in this study using a quasi-
linear perturbed equilibrium approach. The helical Grad-Shafranov equation is
considered and the mode structure of the principle harmonic of the perturbed
helical flux is solved up to the resonant surface to obtain the island
stability criterion $\Delta^{\prime}$, where the PS current term from the
curvature and pressure is neglected, because we use a small aspect-ratio. The
profile modification within the island is introduced by considering the 2D
profiles and extracting their zeroth and principle harmonic to feed into the
RHS of the helical Grad-Shafranov equation. Among the RHS terms, the helical
current perturbation contributes consistently to the island’s stability, with
positive helical current perturbation always stabilizing the island while the
negative one does the opposite. These findings are in agreement with the
previous constant $\psi$ analysis. The pressure contribution is more
complicated due to the cancellation of its odd parity contribution in the
second derivative of the perturbed helical flux. In the small island regime,
the pressure modification contribution is monotonic, and the pressure bump is
found to stabilize the island while the pressure hole does the opposite, but
their overall stability impact is weak due to the aforementioned cancellation.
In the large island regime, the pressure impact becomes non-monotonic,
especially when we combine the pressure and current modifications. Both large
pressure bumps and holes are found to stabilize the island in this regime.

Based on our conclusions above, for example, one would expect that the local
heating by electron cyclotron wave within the island would result in an
additional stabilization effect apart from the current drive caused by ECCD.
Another example is the injection of deuterium pellets into the island O-point,
which results in dilution cooling inside the island region. The cooling itself
would result in resistive helical current expulsion from the island O-point,
thus destabilizing. But the low temperature island region would also cause
heat flux to flow into the island O-point, increasing its thermal density and
pressure. Based on our result, one would anticipate the pressure bump caused
by such an injection to mitigate the destabilization caused by the current
expulsion, thus slowing down the island’s growth. These are some examples
derived from our results above.

We would like to emphasize again that the current and pressure perturbation
considered in this study is arbitrary, so long as they can be expressed as a
function of the helical flux. Some immediate future development would be
pursued to obtain a more self-consistent description of the island’s
stability. For example, the bootstrap current modification could be
automatically specified for a given pressure modification, so that the
bootstrap contribution can be taken into account self-consistently[47]. The
same is true for the vortices mode contribution within the island which
impacts the thermal transport, resulting in pressure profile modification thus
bootstrap current modification[48]. Implementing consistent description of the
aforementioned profile modifications will be part of our future works. Another
future development of the current work is to retain the principal harmonic of
the nonlinear terms involving the other helicity modes in the RHS of Eq. (5)
to take into the nonlinear coupling with other islands. The non-circular
shaping effect of the plasma cross-section can be like-wise be added.
Furthermore, the convective term could be added into the force balance
equation to take the finite flow effect into account in the equilibrium
solution.

Overall, the study carried out here provides insight into the island’s
response to various disruption avoidance or mitigation scheme such as current
drive and local heating caused by cyclotron waves, or local cooling and
current expulsion caused by massive material injections. Understandings on the
above dynamics would help to achieve more efficient disruption avoidance and
mitigation in future machines.

Acknowledgments

The authors thank L.E. Zakharov for fruitful discussion and permission to
carry out code development upon his code infrastructure. This work is
supported by the National MCF Energy R&D Program of China under Grant No.
2019YFE03010001 and the National Natural Science Foundation of China under
Grant No. 11905004.

References

## References

  * [1] Lehnen M, Aleynikova K, Aleynikov P B, et al  2015 J. Nucl. Mater. 463 39
  * [2] Fil A, Nardon E, Hoelzl M, et al and JET Contributors  2015 Physics of Plasmas 22 062509
  * [3] Zohm H, Gantenbein G, Giruzzi G, et al  1999 Nuclear Fusion 39 577
  * [4] Giruzzi G and Zabiego M  2001 Fusion Engineering and Design 53 43
  * [5] Tang Weikang, Wang Zheng-Xiong, Wei Lai, et al  2020 Nuclear Fusion 60, 026015
  * [6] Zhang Y, Wang X J, Zhang X D, et al  2021 Nuclear Fusion 61, 096028
  * [7] Shiraki D, Commaux N, Baylor L R, et al  2018 Nuclear Fusion 58 056006
  * [8] Li Y, Chen Z Y, Wei Y N, et al  2018 Rev. Sci. Instrum. 89, 10K116
  * [9] Xu H B, Nie L, Zhu G L, et al  2018 Fusion Sci. Technol. 75 98
  * [10] Baylor L R, Meitner S J, Gebhart T E, et al  2019 Nuclear Fusion 59 066008
  * [11] Raman R, Sweeney R, Moyer R A, et al.  2020 Nuclear Fusion 60 036014
  * [12] Sweeney R, Baylor L, Eidietis N W, et al  2020 Nuclear Fusion 61 066040
  * [13] Park SooHwan, Lee Kunsu, Baylor Larry R, et al  2020 Fusion Engineering and Design 154 111535
  * [14] Sheikh U A, Shiraki D, Sweeney R, et al  2021 Nuclear Fusion 61 126043
  * [15] Park SooHwan, Lee Kunsu, Lee HyunMyung, et al  2021 Fusion Engineering and Design 164 112200
  * [16] Jachmich S, Kruezi U, Lehnen M, Baruzzo M, et al and JET Contributors  2022 Nuclear Fusion 62 026012
  * [17] Jang Juhyeok, Kim Jayhyun, Kim Jaewook, et al  2022 Fusion Engineering and Design 180 113172
  * [18] Pautasso G, Zhang Y, Reiter B, et al  2011 Nuclear Fusion 51 103009
  * [19] Lehnen M, Alonso A, Arnoux G, et al  2011 Nuclear Fusion 51 123010
  * [20] Hollmann E M, Arnoux G, Commaux N, et al  2011 J. Nucl. Mater. 415 S27
  * [21] Reux C, Plyusnin V, Alper B, et al  2015 J. Nucl. Mater. 463 143
  * [22] Duan X R, Liu Yi, Xu M, et al  2017 Nuclear Fusion 57 102013
  * [23] Chen D L, Shen B, Granetz R S et al  2018 Nuclear Fusion 58 036003
  * [24] Ding Yonghua, Chen Zhongyong, Chen Zhipeng et al  2018 Plasma Sci. Technol. 20 125101
  * [25] Hollmann E M, Parks P B, Shiraki D, et al  2019 Phys. Rev. Lett. 122 065001
  * [26] Izzo V A  2020 Nuclear Fusion 60 066023
  * [27] Nardon E, Fil A, Hoelzl M, Huijsmans G, et al and JET contributors  2016 Plasma Phys. Control. Fusion 59 014006
  * [28] Hu D, Nardon E, Lehnen M, et al  2018 Nuclear Fusion 58 126025
  * [29] Kim C C, Liu Y, Parks P B, et al  2019 Physics of Plasmas 26 042510
  * [30] Ferraro N M, Lyons B C, Kim C C, Liu Y Q and Jardin S C  2019 Nuclear Fusion 59 016001
  * [31] Hoelzl M, Hu D, Nardon E and Huijsmans G T A  2020 Physics of Plasmas 27 022510
  * [32] Nardon E, Hu D, Hoelzl M and Bonfiglio D  2020 Nuclear Fusion 60 126040
  * [33] Zafar A, Zhu P, Ali A, Zeng S Y, Li H L  2021 Plasma Sci. Technol. 23 075103
  * [34] Hu D, Nardon E, Hoelzl M et al  2021 Nuclear Fusion 61 026015
  * [35] Zeng S Y, Zhu P, Izzo V A, Li H L and Jiang Z H  2021 Nuclear Fusion 62 026015
  * [36] Hu D and Zakharov L E  2015 Journal of Plasma Physics, 81(6) 515810602
  * [37] Liu Y, Akcay C, Lao L L and Sun X  2022 Nucl. Fusion 62 126067
  * [38] Lao L L, Kruger S, Akcay C, et al  2022 Plasma Phys. Control. Fusion 64 074001
  * [39] Zakharov L E  1978 Nuclear Fusion 18 335
  * [40] White R B, Gates D A and Brennan D P  2015 Physics of Plasmas 22 022514
  * [41] Zakharov L E and Shafranov V D.  1986 Reviews of plasma physics, Volume 11 (New York: Consultants Bureau) p.153
  * [42] Yu Q, Günter S, Lackner K, Gude A and Maraschek M  2000 Nuclear Fusion, 40 2031
  * [43] Rutherford P H  1973 Physics of Fluids 16 1903
  * [44] Kotschenreuther M, Hazeltine R D, and Morrison R J  1985 Physics of Fluids, 28 1985
  * [45] Wesson J, Campbel D J, Connor J W, et al  2011 Tokamaks (4th edition) (New York: Oxford) pp.152
  * [46] Glasser A H, Greene J M and Johnson J L  Resistive instabilities in general toroidal plasma configuration.  1975 Physics of Fluids, 18 875
  * [47] Kikuchi M and Azumi M  1995 Plasma Physics and Controlled Fusion 37 1215
  * [48] Hornsby W A, Peeters A G, Snodin A P, et al  2010 Physics of Plasmas, 17 092301