A Riccati solution for the ideal MHD plasma response with applications to real-time stability control

Active feedback control of plasma stability is essential to achieving high performance in advanced tokamaks. Tokamak experiments with active control systems for tearing modes,² edge-localized modes (ELMs),³ and unstable divertor topologies⁴ have demonstrated success in mitigating such instabilities. Passive, closed-loop control systems⁵ and linear gain controllers⁶ have also demonstrated improved equilibrium and unstable mode control in reversed-field pinch (RFP) and tokamak experiments.

Due to their computational complexity, however, state-of-the-art measures of plasma stability are largely absent from such active feedback methods. Tokamaks are therefore limited to controlling many instabilities after they are observed—e.g., with “catch-and-subdue” strategies, as described in Ref. 2. These methods constrain the real-time capability of stability control systems, and do little to address vulnerability to unplanned-for ideal MHD instabilities (which may grow on timescales as fast as the Alfvén time $τA∼μs$⁠).

While the growth of unstable MHD modes may be quite rapid, the evolution of stable tokamak equilibria is considerably slower. Stable tokamak magnetic field geometries evolve on an inductive timescale $τL/R≳1s≫τA$⁠. This separation of timescales creates a window of opportunity to prevent MHD plasma instabilities before they start. That is, although the instabilities may be too fast to control, a control system operating at a timescale $τcontrol≪τL/R$ might be fast enough to maintain a tokamak's operation in stable equilibria—steering it clear of evolutions through unstable equilibria which spawn uncontrollable MHD instabilities. For example, a tokamak control system that can measure δW in near real-time could be tasked with maintaining a desired “distance-to-unstable-equilibrium”—i.e., $δW>δWmin>0$—as a buffer to the formation of MHD instabilities that accompany unstable equilibria. (δW will be defined at the beginning of Sec. II.)

Any such pre-emptive active feedback MHD control system would require (i) diagnosis and fitting of the plasma equilibrium; (ii) analysis of the stability characteristics of that equilibrium; and (iii) an active controller to steer the plasma away from its stability boundaries, all within a timescale appropriate to the magnetic equilibrium evolution. This work treats only the second of these three control system components—namely, developing the capability for real-time stability analysis.

In particular, in this paper we focus on the real-time analysis of ideal MHD tokamak δW stability. (Additional work on real-time resistive MHD stability analysis is also underway, but lies beyond the scope of this paper.) We demonstrate the computational viability of a near real-time ideal MHD stability analysis by reformulating the δW variational method in a setting more familiar to control theory. Using a standard Hamilton-Jacobi analysis, the plasma response matrix is shown to be the solution of an ideal MHD matrix Riccati differential equation (MRDE). Identifying the ideal MHD δW problem as a Riccati problem brings to bear a range of flexible solution methods which have long been developed in the control theoretic literature.^7–12 Such methods facilitate the near real-time solution of ideal MHD stability problems, and may also enable stability studies for equilibria requiring stiff, high mode-number descriptions (e.g., in stellarator and tokamak edge plasmas).

The δW Riccati formulation not only suggests improved numerical methods for rapid and robust stability analysis, but also highlights an intuitive physical interpretation of the plasma response matrix, seemingly unemphasized in the existing plasma stability literature. In particular, we show that the plasma response matrix is a bilinear form mapping plasma perturbations to their effect on the plasma energy—δW.

We demonstrate the efficacy of Riccati solution methods with a practical example, by adapting Glasser's “Direct Criterion of Newcomb” (DCON) code¹ for use in near real-time. We replace DCON's serial integration of its Euler-Lagrange equations (ELEs) with a domain-decomposed integration of state transition matrices, effecting the code's parallelization and achieving a run-time much less than ITER's MHD evolution time. Although we demonstrate the applicability of Riccati methods to DCON, a non-inertial 2D ideal MHD code, we note that other MHD stability codes, which may incorporate resistive, inertial, 3D, and edge effects, could also benefit from such techniques. The ELITE,^13,14 PEST,^15,16 MARG2D,¹⁷ and resistive DCON¹⁸ stability codes, for example, also solve stiff, spectrally decomposed linear ELEs, and may benefit from the methods described here.

The remainder of this paper is organized as follows: Sec. II demonstrates the equivalence of the ideal MHD Lagrangian variational stability problem to an MRDE, emphasizing a physical intuition for the plasma response matrix. Section III describes the numerical advantage of treating the δW ELE integration as a control theoretic MRDE problem, a problem widely explored in the control theory literature. Section IV presents additional methods used to achieve high performance in the real-time, parallelized adaptation of the DCON code, and characterizes this performance. Section V summarizes and concludes.

The δW stability analysis solves for the least stable perturbations to a plasma equilibrium—in particular, the spatial plasma displacements $Ξ$ that most reduce the plasma's potential energy (releasing, in the process, kinetic energy). This change in potential energy is denoted by the functional $δW[Ξ,Ξ†]$⁠. A plasma equilibrium with δW > 0 is said to be ideal-MHD-stable, while δW < 0 indicates ideal-MHD-instability.

In this section, we present a new, intuitive derivation of the following claim, which introduces the plasma response matrix as a map from perturbations at the plasma edge to their corresponding effect on plasma potential energy:

Claim: The minimum fluid δW, for a stable, axisymmetric plasma equilibrium, is given by

where ψ = 1 indicates the plasma edge, $Ξ(ψ)∈CM$ is a Fourier-decomposed perturbation (spatial displacement) of the plasma on flux surface ψ (truncated at M poloidal mode numbers), and where matrix $P$ is a solution to the MRDE

$WP≡P(1)$ is called the plasma response matrix.

The prime in $P′$ above denotes a derivative with respect to ψ—the flux label (radial) coordinate. The ψ-dependent matrices ${F=F†$⁠, $G=G†$⁠, and $K≠K†}∈CM×M$ describe the ideal MHD couplings between poloidal mode perturbations and the plasma equilibrium, and their resulting effect on the plasma's energy—see a description of their calculation in Ref. 1.

We demonstrate this claim in the following sections by (Sec. II A) reviewing the ideal MHD δW Lagrangian formulation; (Sec. II B) Legendre transforming to a Hamiltonian setting; (Sec. II C) applying Hamilton-Jacobi theory; and (Sec. II D) reformulating the resulting δW analysis as an MRDE Riccati problem.

We begin by restating Eq. (19) of Ref. 1, which neatly captures in matrix form the change in energy δW resulting from the least stable magnetic perturbations of a general axisymmetric toroidal plasma

Here, $0≤ψ≤1$ is a flux surface label (radial) coordinate extending from the magnetic axis to the plasma edge, and plasma perturbation $Ξ∈CM$ is a vector whose entries represent a radial displacement of the plasma in the direction of $ψ̂$⁠, Fourier-decomposed into M modes. (Such modes may be labeled by their toroidal and poloidal mode numbers (n, m), respectively: $Ξ(m,n) exp [i(mθ−nζ)]$⁠.) Although they are Hermitian adjoints, $Ξ$ and $Ξ†$ are taken to be independent dynamical variables, capturing the degrees of freedom in their independent real and imaginary parts. We treat this integral for the perturbed energy as the action integral of the δW Lagrangian, with coordinate ψ acting in lieu of a time parameter. Matrices ${F(ψ),G(ψ), and K(ψ)}$ are as described above.

Since a plasma is unstable to those perturbations which reduce its potential energy, we seek to characterize the perturbations $Ξ$ that minimize δW. A necessary condition at any local extremum of δW is that its variational $δ(δW)$ vanishes for arbitrary variations to the perturbations, ${δΞ,δΞ†}$⁠. It is in this sense that L in Eq. (3) represents the appropriate Lagrangian for stability analysis in our system.

Noting that L is Hermitian, and assuming our system is stable to all perturbations—(i.e., L is positive definite)—then L is convex, and we are free to Legendre transform this variational problem to its Hamiltonian formalism.¹⁹ (As our method is intended to pre-emptively analyze plasma stability before MHD instabilities arise, this is an entirely natural assumption to make.) Defining $q1≡Ξ$ and $q2†≡Ξ†$⁠, our system's canonical momenta (absorbing 1/2μ₀ for convenience) are

We thus derive the quadratic Hermitian Hamiltonian

When treating the system dynamically, we have been careful to separate $q1$ from $q2†$⁠. Nevertheless, we note that $(q1)†=q2†$ and $(p1†)†=p2$⁠. Throughout the remainder of this paper, we duly omit such subscripts.

The original Lagrangian problem may therefore be reframed; we must find the plasma perturbations q (and their conjugate momenta p) which extremize the action according to the Hamiltonian of Eq. (5),

After an integration by parts, the perturbations which satisfy this variation are found to be those obeying Hamilton's equations of motion (EOM)

which can be expressed compactly in the following 2M degree-of-freedom linear dynamical system, the δW EOM:

The integration by parts of the second term in Eq. (6) yields an additional boundary term, which must vanish at extrema of the action:

In axisymmetric toroidal magnetic systems, regularity conditions require that all modes—except those with $|m|=1$—vanish at the magnetic axis: $Ξ(0)≡q(0)=0$⁠. The (n, m) = (1, 1) mode, on the other hand, finitely displaces the magnetic axis such that $Ξ(1,1)(0)≠0$⁠. (Such a displacement of the axis gives rise to sawteeth oscillations.)

Numerically, however, it is insufficient to choose the corresponding “natural” set of initial conditions: q(0) = 0 for all modes $|m|≠1$⁠, and p(0) = 0 for $|m|=1$⁠. This is because the magnetic axis represents a regular singular point of the ordinary differential equation (ODE), Eq. (8). (As we shall later discuss, a Riccati solution for this ODE with such initial conditions would be ill-conditioned.) As suggested in Ref. 1, a rigorous treatment of the magnetic axis would therefore require an asymptotic expansion of the M regular extremal modes at a small distance from the axis. Such an effort lies beyond the scope of this work, and we shall at present impose q(0) = 0 boundary conditions for all modes. We note that this choice effectively omits extremal $|m|=1$ modes from our analysis.

As for boundary conditions at the plasma edge, we are in general interested in a tokamak equilibrium's stability to both internal plasma perturbations, for which $Ξ(1)=0$⁠, and external perturbations, which have displacements at the plasma edge. Consequently, we do not impose edge boundary conditions a priori, and will be interested in the behavior of all M independent solutions to Eq. (8) consistent with the boundary conditions we have selected on the axis.

We now solve our system using a classical strategy of Hamiltonian theory, canonically transforming to the corresponding Hamilton-Jacobi problem (as described in Ref. 19). In this approach, the problem is transformed by performing a canonical transformation that zeros the Hamiltonian everywhere. Recalling that any canonical transformation $(q,p;H)↦(Q,P;K)$ must preserve the Lagrangian up to a total “time” derivative

we use a type-2 generating function, demanding that

and setting K = 0. This canonical transformation, substituted into Eq. (10), yields three equations (or five, including equivalent adjoints)

Here, $H=H(q,q†,p,p†,ψ)$ is our Hamiltonian, and S is Hamilton's principal function, a solution to

for which we must solve. Substituting Eqs. (12) into Eq. (5), we find that S satisfies the following first-order second-degree partial differential equation (PDE):

Subscripts of S in the equation above denote partial derivatives, e.g., $Sq≡∂S/∂q$⁠.

Before solving this PDE for S (as we will do in Sec. II D), we note that K = 0 implies constant $Q(ψ)≡β$ and $P(ψ)≡α$⁠, because $Q′=∂K/∂P=0=−∂K/∂Q=P′$⁠. Substituting from Eq. (12) for the partial derivatives of S, therefore, the total derivative of S is found to be

Upon comparison with the first line of Eq. (5), we note that we have recovered the well-known result of Hamilton-Jacobi theory

We may choose the overall constant that can be freely added to any solution of Eq. (14) in order to achieve a physical interpretation of S. In particular, the regularity assumption of our perturbations at the magnetic axis, (as well as the vanishing measure dψ at the axis), encourages us to set S(ψ = 0) = 0. S may then be regarded as the action of our system, which tracks the “cumulative change” in plasma energy due to an extremal perturbation $Ξ(ψ)$ between the magnetic axis and a given flux surface ψ—i.e., on the interval [0, ψ]. The importance of this choice for S is underscored by noting that, at the plasma edge, S(ψ = 1) corresponds to the value of δW

Equation (14) is therefore a dynamical equation for the action and, by extension, for the cumulative energetic cost of a plasma perturbation. When solved at the plasma edge, Eq. (14) implicitly maps the edge perturbation to its corresponding least stable excited modes, and measures the impact of such modes on the plasma potential energy. Indeed, given initial conditions, a solution for the action S to Eq. (14) would yield the corresponding solutions for q and p as well.

We may now solve for S. Let us assume (and we shall prove in a moment) that our dynamical variables obey a linear dependence

for some matrix $P$⁠. (Note matrix $P$—sans serif—is not to be confused with transformed canonical coordinate vector P.) Using the product rule to take the derivative of Eq. (18), and substituting from our dynamical Eq. (8), we derive the ideal MHD δW Riccati formulation

(By the self-adjointness of this ODE, we observe that $P$ is everywhere Hermitian if it is anywhere Hermitian.) Once the linear relationship in the dynamical variables Eq. (18) is assumed, the Riccati Eq. (19) is an equivalent formulation of the dynamical Eq. (8) resulting from the δW analysis. We note the equivalence of the Riccati matrix $P$ to the critical determinant matrix $W$ in Eq. (58) of Ref. 1. Indeed, the Riccati formulation is a natural setting for the study of a Hamiltonian system's stability.

We now observe that, given Eqs. (12) and (18), it is natural to attempt a solution for S of the form

We see that substituting this relation for S into the dynamical action PDE of Eq. (14) immediately reproduces the Riccati ODE in Eq. (19). The validity of the solution Eq. (20) for S is thereby demonstrated.

Despite the complexity of the nonlinear PDE Eq. (14), the Hamiltonian structure of our problem allows us to show that the above solution for S is nearly unique. In particular, the linearity assumption of Eq. (18) and the first of the transformed dynamical Eqs. (12) require $S(q,q†,α,α†,ψ)=q†P(ψ)q+f(α,α†,ψ)$ for some function f. Upon substitution into Eq. (14), then, and after applying Eq. (19) for $P′=∂P/∂ψ$⁠, one discovers that this solution for S must have $∂f/∂ψ=0$ as well. Since the vectors ${α,α†}$ are constant for all ψ, we are left to conclude that, indeed, Eq. (20) is the unique solution of our PDE, up to a constant.

We have only left to show, then, that there exists a matrix $P(ψ)$ satisfying the linear relationship in Eq. (18). We do so by construction, finding $P$ explicitly, as follows. We concisely rewrite Eq. (8) in the form

We then consider the fundamental matrix of solutions $Φ$ for this ODE, a 2M × 2M matrix whose columns form a complete basis for all possible solutions of the system

Note that $Φ0≡Φ(0)$ may in general be any nonsingular matrix, so long as it spans the space of all possible initial conditions for the ODE. (We choose $Φ0=1$ in the definition above merely for its simplicity.). Since the fundamental matrix forms a complete basis, a solution satisfying any boundary condition can be formed from a linear combination of its columns—a feature we now put to use.

Since $Φ0=1$ is nonsingular, $x0=Φ0c$ for some constant coefficient vector c. But by ODE linearity, it must then hold true for all ψ that $x(ψ)=Φ(ψ)c$⁠. Combining these relationships, we find a familiar result from linear ODE theory

(The economy of the choice $Φ0=1$ lies in the cancellation of the inverse factor above.) In this way, $Φ(ψ)$ is to be regarded as the state transition matrix—or propagator—of the system, which maps ODE solutions forward in ψ.

Therefore, given the fundamental matrix of solutions $Φ(ψ)$⁠, we may map forward any initial conditions that are imposed at ψ = 0. In particular, we may write Eq. (23) as

for any initial conditions q(0) = q₀, p(0) = p₀.

This expression is sufficient to assert a general linear dependence between q and p in any such linear system of ODEs. After all, given $p0=P0q0$⁠, we see that

To specify this linear dependence for our q₀ = 0 boundary conditions at the axis, however, we must modify this analysis; in particular, we set $q0=0·p0$⁠. Substituting this relation into Eq. (24), after a short manipulation we again find a linear mapping from $q(ψ)$ to $p(ψ)$

This relation is in the desired form of Eq. (18), and satisfies our choice of boundary conditions. The matrix factor in Eq. (26) thus provides a solution to the ideal MHD δW Riccati equation, proving our claim.

We briefly comment here on the numerical difficulty posed by the “natural” boundary conditions for the tokamak, with $q0≠0$ for $|m|=1$ modes, as discussed in Sec. II B. In particular, such “natural” boundary conditions fail to produce a well-conditioned Riccati solution because in such a case there is no linear transformation from the set of all q₀ to the set of all p₀, nor vice versa. In other words, the “natural” boundaries require that the set of M q₀ vectors and the set of M p₀ vectors are both degenerate. An asymptotic expansion of solutions slightly off-axis would repair this degeneracy.

At the plasma edge, the Riccati solution

comprises the plasma response matrix. It may be regarded as the plasma permeability, or the “conjugate momentum” created per unit of plasma displacement—i.e., $WP∼p/q$⁠.

The Riccati formulation clarifies how $WP$ is appropriately viewed as the bilinear form which maps a perturbation q(1) at the plasma surface to its associated energetic cost, δW. This interpretation emphasizes the importance of its eigenvalues in studies of plasma stability. In particular, the system's energy due to a perturbation is expressly given by combining Eqs. (17) and (20)

The foregoing relation appears in Ref. 1 in the context of a proof of Newcomb's criterion. However, recognition of its connection to Hamilton-Jacobi theory and the Riccati equation—and, by extension, the wider control literature—has been absent in the plasma stability literature.

We note that the quadratic Hamiltonian of Eq. (5) immediately renders our problem on a similar footing to linear quadratic regulator (LQR) control problems that feature quadratic cost functions. This is the cause of the appearance of the Riccati equation in each case. Identifying our problem as a Riccati problem is quite profitable: Techniques for solving such problems are well-developed, even when the ODEs involved are stiff, or singular. Such techniques can therefore be leveraged to improve plasma stability codes, as we now show.

Identifying the δW stability problem as an MRDE connects plasma stability calculations to the wider academic literature devoted to solving Riccati problems. The stiffness of ideal MHD ODE solutions is common in Riccati boundary value problems (BVPs), and methods for overcoming it are widely discussed, for example in Refs. 8–11. Solutions for approaching singular Riccati ODEs feature in the literature as well: Ref. 11, for example, describes a solution method for singular Riccati ODEs that integrates through the equation's singularities.

Such literature suggests that in solving for $WP$⁠, it may at times be advantageous not to integrate the Hamiltonian system of Eq. (8), but rather solve for the plasma response matrix via the Riccati Eq. (19), directly. In high mode-number 3D stellarator and tokamak edge calculations, for example, where stiffness can render the Hamiltonian system nearly insoluble, such an alternative could prove decisive.

The parallel implementation of Riccati solutions has also received wide attention in the literature.¹² A common technique for solving matrix Riccati control problems, which immediately renders their solution in parallel, is that of state transition matrices⁷ (introduced in Sec. II D). We briefly explore how state transition matrices afford a parallelizable solution to the ideal MHD Riccati problem, and further mitigate the difficulties caused by stiffness as well as the singularity of solutions at the rational surfaces. We then discuss the parallel adaptation of DCON, using these features, in Sec. IV.

Viewing our ODE's fundamental matrix of solutions $Φ$ as a control theoretic state transition matrix—or propagator—proves quite useful in solving the δW EOM [Eq. (8)]. In particular, as seen in Eq. (24), $Φ$ is an operator that maps solutions forward in ψ. As a result, any interval of the ODE Eq. (22) can be subdivided, so that $Φ$ is recovered from a multiplication of its subpropagators

where $Φ(ψi,ψi)=1$ is the initial condition for subinterval i. Integration of Eq. (22) across subdomains is thereby reduced to the multiplication of matrices that may each be independently calculated, which affords an integration of the EOM in parallel.

Furthermore, since the ultimate aim of our calculation is $WP=P(1)=Φpp(1)Φqp−1(1)$⁠, we are free to transform the cumulative subpropagator $∀ N≥1$

by any right-multiplied linear operator of the form $RN=(XN0YNAN)$ with $AN$ nonsingular. For example, the product

leaves

invariant. In particular, this enables us to perform Gaussian elimination (via column reduction) on the right-side columns of each cumulative subpropagator, separating orders of magnitude spanned by the stiff solutions of the ODE. This crucial advantage mitigates otherwise catastrophic numerical error in taking the matrix product of subpropagators, which may span many orders of magnitude (⁠$∼O(10M)$⁠, say).

Another advantage of the domain-decomposed integration of the state transition matrix is its suitability for integrating near singular rational surfaces. We note that subpropagators are invertible, and therefore satisfy

This feature admits a convenient reversibility of the direction of integration, which may be leveraged to assist at the singular surfaces, as follows.

As is well known, the ideal MHD δW EOM have regular singularities at rational surfaces, at the magnetic axis, and at the separatrix. However, as was previously noted in Ref. 20, the solutions to the EOM are well-behaved when they are integrated away from the singular surfaces; a solution which asymptotically diverges in the forward direction of integration decays in the reverse direction. Therefore, no numerical instability is created, and the integrated modes retain their linear independence. (Such behavior is illustrated in Fig. 1.)

As a result, integration across a singular surface may be achieved by integrating backward from the left side of the singular surface, and taking the matrix inverse of the resulting subpropagator. This is then multiplied with the forward-integrated subpropagator on the right of the singular surface. Using asymptotic expansions at the singular surfaces as initial conditions renders this integration even more numerically robust by effectively separating the solution subspaces.

It is worth further emphasizing the freedom that subdividing the domain of integration affords. Its adaptability may render it useful for calculations beyond the ideal MHD model, providing a convenient technique for resistive and high toroidal mode number MHD stability problems, for example. The authors' early efforts at adapting resistive MHD $Δ′$ calculations—which are known for their numerical intractability—have demonstrated the method's flexibility and robustness, though a full discussion of this extension lies beyond the scope of the present work.

We chose to apply Riccati methods to the ideal DCON code due to its comparatively rapid solution time and widespread use. Using the techniques described in Sec. III, we adapt the ODE integration of DCON, which comprises ∼90% of the code's total runtime.

In our application of the Riccati solution techniques to a parallel adaptation of DCON, we integrate the ODE of Eq. (22) using the ZVODE²¹ complex adaptive integrator. We parallelize this integration with OpenMP,²² subdividing the interval between the magnetic axis and plasma edge. The increased computational cost of our ODE near ill-behaved surfaces (to wit, the magnetic axis and separatrix, and the rational surfaces between them), requires a careful division of integration intervals.

We find that a modification to a grid-packing algorithm suggested in Ref. 18 works adequately to load-balance the computation. In particular, we numerically fit a form-factor

to estimate the time of integration near a point ψ, located at some distance from each nonanalytical surface $ψsi$⁠. The coefficients ${αi,βi}$ are determined by the surface type of $ψsi$ (axial, rational surface, or edge separatrix). The total time of integration over an interval $[ψ1,ψ2]$ is therefore estimated to be

where γ₀ is fit to the unavoidable initialization time of the integrator on each subinterval. The grid is then chosen to minimize τ for a given number of subintervals.

The importance of grid-packing is emphasized in Fig. 2. It is apparent that without an effective packing scheme, the benefit of parallelization would be quite limited; a single interval might otherwise require runtime comparable to runtimes of the entire parallelized code.

As shown in Fig. 3, the eigenvalues of DCON's plasma response matrix are reproduced with high accuracy using our new parallel methodology. For this ITER-relevant CHEASE-generated²³ EFIT²⁴ equilibrium with two rational surfaces, even without a complete optimization of the parallel DCON code, our initial implementation consistently achieves a fluid ODE integration time of $∼185$ ms. Given a projected ITER confinement time of $τE∼5$ s,²⁵ such computation speed may be sufficient for an ideal MHD control system implementation at ITER. After further scaling, such speed may enable a demonstration of ITER-relevant control schemes on operational machines as well, such as DIII-D, where confinement time is estimated to be $τE∼200$ ms.²⁶ The above results were achieved using 2.4 GHz Intel Broadwell processors, with 28 cores.

As a further test of our new method, we plot in Fig. 4 the ten smallest energy eigenvalues of the plasma response matrix over the course of DIII-D shot #163518.²⁷ This QH-mode shot was chosen because its lack of ELMs allowed for a fair comparison between the two methods—(i.e., one without ELM-induced reconstruction issues). The accurate, rapid analysis of such a shot's time-varying δW—which corresponds to the smallest plotted energy eigenvalue—would be critical to any real-time ideal MHD control system of the type we have described.

One factor ostensibly limiting the scalability of our approach is a tradeoff between the speedup of a parallel integration over a finer subdivision of intervals, on the one hand, and the increased time required to matrix-multiply the resulting subpropagators, on the other. As one might expect, however, this constraint turns out to be quite loose, due to the computational ease of matrix multiplication. As depicted in Fig. 5, the benefit of parallelization is likely to extend to a much larger number of processors.

We have leveraged the δW problem's quadratic Hamiltonian structure to derive an alternate representation of the ideal MHD plasma response as a solution to a Riccati differential equation. In doing so, we have sought to connect ideal MHD stability to the wider control literature and highlighted the plasma response matrix's intuitive interpretation as the bilinear form mapping perturbations to their energetic cost. We have used control theoretic Riccati methods to adapt the generic tokamak δW stability analysis and demonstrated the success of this new methodology by modifying the ideal DCON code to achieve runtimes appropriate for tokamak ideal MHD stability control. In particular, we have shown the numerical advantages of using common MRDE state transition matrix techniques in finding robust solutions for the δW EOM in parallel.

We surmise that the abundant literature on the solutions of Riccati equations may be exploited to simplify ideal MHD δW stability analyses beyond the treatment explored here. As noted earlier, one might, for example, consider integrating the MRDE, Eq. (19) itself, to solve for the plasma response matrix $WP$ in stiff stellarator and tokamak edge calculations. Although the singularities of this equation can be just as virulent as the singularities of the dynamical ODE in Eq. (8)—these nonanalytical features appear unavoidably in the coefficients of the MRDE themselves, after all—asymptotic expansions would likely prove useful in the solution of Eq. (19) at singular surfaces. A further advantage of solving Eq. (19) directly is that where $P$ is ill-behaved, the dual Riccati equation for $S≡P−1$ can be equivalently solved instead:

(This expression is derived from the formula

for the matrix derivative.)

We would like to acknowledge Zhirui Wang, Jong-Kyu Park, Ian Cosden, and Keith Erickson for helpful discussions. This research was supported in part by the United States Department of Energy (DOE) under Contract Nos. DE-AC02–09CH11466, Early Career Research Program: DE-SC0015878, and DE-SC0016201. The digital data for this paper can be found in http://arks.princeton.edu/ark:/88435/dsp01gx41mm508.