2.1. The $\epsilon$-constraint method

Multi-objective optimization problems can be transformed into a series of scalar optimization problems via scalarization methods. One such scalarization method, the $\boldsymbol {\epsilon }$-constraint method (Ehrgott Reference Ehrgott2005), finds (weakly) Pareto optimal points by minimizing a single objective subject to upper bound inequality constraints on all other objectives, with an upper bound parameter $\boldsymbol {\epsilon }\in \mathbb {R}^{{n_{\boldsymbol {F}}}-1}$. To find points on the Pareto front, the $\boldsymbol {\epsilon }$-constraint method solves problems of the form

(2.2)\begin{equation} \left.\begin{array}{c} \min_{\boldsymbol{x}\in\boldsymbol{\varOmega}} \quad F_j(\boldsymbol{x})\\ \text{such that}\quad F_1(\boldsymbol{x}), \ldots,F_{j-1}(\boldsymbol{x}), F_{j+1}(\boldsymbol{x}), \ldots F_m(\boldsymbol{x}) \le \boldsymbol{\epsilon} \end{array}\right\},\end{equation}

where the user selects the index $j$ and vector $\boldsymbol {\epsilon }$. The $\boldsymbol {\epsilon }$-constraint method theoretically guarantees that it can be used to find any point on the Pareto front: if $\boldsymbol {x}$ is a (local) solution to (2.2), then it is weakly (locally) efficient. Furthermore, $\boldsymbol {x}$ is efficient if and only if it is a solution to (2.2) for all choices of $j\in \{1,\ldots,{n_{\boldsymbol {F}}}\}$ (Ehrgott Reference Ehrgott2005).

The $\boldsymbol {\epsilon }$-constraint method is particularly relevant to stellarator optimization since we often have bounds on the range of the objectives, which allows us to easily select the upper bound parameter, $\boldsymbol {\epsilon }$. For example, devices often have an aspect ratio between $3$ and $10$, a coil must be at least as long as the minor circumference of the plasma, and it should not be longer than a few times that. By selecting $\boldsymbol {\epsilon }$, targeted areas of the Pareto front can be searched. Thisis more interpretable than using other scalarization techniques, such as a weighted sum of objectives, which requires selecting abstract weights that do not have any clear relation to the Pareto front.

2.2. Continuation methods

Continuation methods (Hillermeier et al. Reference Hillermeier2001; Schtze, Dell'Aere & Dellnitz Reference Schütze, Dell'Aere and Dellnitz2005; Gkaragkounis et al. Reference Gkaragkounis, Papoutsis-Kiachagias, Asouti and Giannakoglou2018; Peitz, Ober-Blöbaum & Dellnitz Reference Peitz, Ober-Blöbaum and Dellnitz2019; Vasilopoulos et al. Reference Vasilopoulos, Asouti, Giannakoglou and Meyer2021) in MOO explore the Pareto front locally around a given efficient point. By using optimality conditions and the implicit function theorem (Krantz & Parks Reference Krantz and Parks2002), continuation methods build local models of the Pareto front, which they use to estimate nearby efficient points. Local exploration of the Pareto front with continuation methods comes at a relatively low computational cost when compared to restarting a scalarization solver, such as the $\boldsymbol {\epsilon }$-constraint method, from scratch. However, when combined, continuation methods and scalarization solvers make a powerful pair: the scalarization methods find efficient points over distinct parts of the Pareto front and the low-cost continuation method is used to fill in the space between points. In this section, we discuss a predictor–corrector continuation method based off of the $\boldsymbol {\epsilon }$-constraint method, which first predicts an estimate of an efficient point using a Taylor expansion and subsequently uses the $\boldsymbol {\epsilon }$-constraint method to correct the predictions, see Algorithm . While other predictor–corrector continuation schemes (Schütze et al. Reference Schütze, Dell'Aere and Dellnitz2005; Peitz et al. Reference Peitz, Ober-Blöbaum and Dellnitz2019; Vasilopoulos et al. Reference Vasilopoulos, Asouti, Giannakoglou and Meyer2021) write the expansion in terms of weights, or target points, we expand the Pareto front in terms of $\boldsymbol {\epsilon }$ because it is easily interpreted and allows for the $\boldsymbol {\epsilon }$-constraint method to be used in the corrector step.

Algorithm 1: ϵ-Constraint predictor–corrector method

The key idea of the predictor step is to notice that the $\boldsymbol {\epsilon }$-constraint method allows us to parametrically describe weakly efficient points in terms of $\boldsymbol {\epsilon }$, i.e. we can write weakly efficient points as $\boldsymbol {x}(\boldsymbol {\epsilon })$. In fact, as we will show, under reasonable conditions, $\boldsymbol {x}(\boldsymbol {\epsilon })$ is a continuously differentiable map. The predictor step leverages this parametric representation to approximate the efficient set by a Taylor expansion. Given an efficient point $\boldsymbol {x}$, which is a solution to (2.2) with parameter $\boldsymbol {\epsilon }$, the predictor step locally approximates the efficient set by

(2.3)\begin{equation} \boldsymbol{x}(\boldsymbol{\epsilon}+ {\rm \Delta} \boldsymbol{\epsilon}) \approx \boldsymbol{x}(\boldsymbol{\epsilon}) + \nabla_{\boldsymbol{\epsilon}}\boldsymbol{x}(\boldsymbol{\epsilon})^T{\rm \Delta} \boldsymbol{\epsilon}, \end{equation}

where ${\rm \Delta} \boldsymbol {\epsilon }$ is a small change in $\boldsymbol {\epsilon }$, and $\nabla _{\boldsymbol {\epsilon }} \boldsymbol {x}$ is the Jacobian of $\boldsymbol {x}(\boldsymbol {\epsilon })$ with respect to $\boldsymbol {\epsilon }$. Though less accurate, a zeroth-order Taylor expansion could be used in place of the first-order expansion. A zeroth-order expansion simplifies the predictor–corrector method, Algorithm , by using the last known efficient point $\boldsymbol {x}(\boldsymbol {\epsilon })$ as a prediction of a nearby efficient point $\boldsymbol {x}(\boldsymbol {\epsilon } + {\rm \Delta} \boldsymbol {\epsilon })$, and is equivalent to warm starting the $\boldsymbol {\epsilon }$-constraint method with $\boldsymbol {x}(\boldsymbol {\epsilon })$. The zeroth-order predictor is useful when the Jacobian $\nabla _{\boldsymbol {\epsilon }}\boldsymbol {x}$ is intractable to approximate or the cost of running the $\boldsymbol {\epsilon }$-constraint method is low. However, if $\nabla _{\boldsymbol {\epsilon }}\boldsymbol {x}$ can be approximated, then the first-order expansion will give more accurate predictions of $\boldsymbol {x}(\boldsymbol {\epsilon }+{\rm \Delta} \boldsymbol {\epsilon })$ than the zeroth-order expansion, and should be used. A visualization of the first-order predictor corrector method is shown in figure 3. In the remainder of this section, we will discuss conditions which determine when the expansion, (2.3), exists and a method of computing the Jacobian $\nabla _{\boldsymbol {\epsilon }}\boldsymbol {x}$.

Figure 3. Visualization of the predictor–corrector continuation method.

The existence of the expansion, (2.3), and the computation of the Jacobian $\nabla _{\boldsymbol {\epsilon }}\boldsymbol {x}$ rely on applying the implicit function theorem (Krantz & Parks Reference Krantz and Parks2002) to the first order necessary conditions for optimality of (2.2), the Karush–Kuhn–Tucker (KKT) conditions (Nocedal & Wright Reference Nocedal and Wright1999). The KKT conditions are necessary conditions for a point $\boldsymbol {x}$ to be an optimal solution to an inequality constrained optimization problem, such as the $\boldsymbol {\epsilon }$-constraint problem (2.2). In the case of the $\boldsymbol {\epsilon }$-constraint method, points satisfying the KKT conditions for (2.2)are also weakly locally efficient.

Suppose that $\boldsymbol {x}$ is an optimal solution to (2.2), and $\mathcal {A}\subseteq \{1,\ldots,j-1,j+1,\ldots {n_{\boldsymbol {F}}}\}$ is the set of $m$ active constraints, i.e. the set of constraints $i$ for which $F_i(\boldsymbol {x})=\epsilon _i$. The KKT conditions with the linear independence constraint qualification (LICQ) (Nocedal & Wright Reference Nocedal and Wright1999) state that there exist $\boldsymbol {\lambda }\in \mathbb {R}^{m}$ such that $\boldsymbol {\lambda } \ge 0$ and

(2.4)\begin{gather} \boldsymbol{\nabla} F_j(\boldsymbol{x}) + \sum_{i\in\mathcal{A}} \lambda_i\boldsymbol{\nabla} F_i(\boldsymbol{x}) = 0, \end{gather}

(2.5)\begin{gather}\lambda_i(F_i(\boldsymbol{x}) - \boldsymbol{\epsilon}_i) = 0 \quad \forall \,i\in\mathcal{A}, \end{gather}

(2.6)\begin{gather}F_i(\boldsymbol{x}) \le \boldsymbol{\epsilon} \quad \forall \,i\neq j. \end{gather}

Notice that (2.4) and (2.5) show that if a constraint is active and $\lambda _i >0$, then the solution of the $\boldsymbol {\epsilon }$-constraint method is dependent on $\epsilon _i$ in the sense that a small change in $\epsilon _i$ could change the solution to the system of equations. In this case, it is reasonable to write $\boldsymbol {x}$ as a function of $\boldsymbol {\epsilon }_{\mathcal {A}} = (\epsilon _i)_{i\in \mathcal {A}}$ to stress the dependence of solutions on the parameters, $\boldsymbol {x}(\boldsymbol {\epsilon }_{\mathcal {A}})$. The implicit function theorem gives us the mathematical leeway to do so.

Informally, the implicit function theorem states that if there exists a point $\boldsymbol {x}^*$, with Lagrange multipliers $\boldsymbol {\lambda }^*$ and constraint parameters $\boldsymbol {\epsilon }_{\mathcal {A}}^*$, which satisfy the KKT conditions, and if the Jacobian, $\boldsymbol {D}(\boldsymbol {x}^*,\boldsymbol {\lambda }^*)$, of the left-hand side of (2.4) and (2.5) with respect to $(\boldsymbol {x},\boldsymbol {\lambda })$ is invertible, then $\boldsymbol {x}(\boldsymbol {\epsilon }_{\mathcal {A}}), \boldsymbol {\lambda }(\boldsymbol {\epsilon }_{\mathcal {A}})$ are continuously differentiable functions of $\boldsymbol {\epsilon }_{\mathcal {A}}$ around $\boldsymbol {\epsilon }_{\mathcal {A}}^*$. Taking the form of a continuously differentiable function, the space of solutions to the $\boldsymbol {\epsilon }$-constraint problem (2.2) can be approximated with the first-order Taylor expansion, (2.3). The Jacobian $\nabla _{\boldsymbol {\epsilon }}\boldsymbol {x}$ can be computed by solving the system

(2.7)\begin{equation} \boldsymbol{D}(\boldsymbol{x}^*,\boldsymbol{\lambda}^*) \begin{bmatrix} \nabla_{\boldsymbol{\epsilon}}\boldsymbol{x} \\ \nabla_{\boldsymbol{\epsilon}}\boldsymbol{\lambda} \end{bmatrix} = \begin{bmatrix} \boldsymbol{0} \\ \boldsymbol{\lambda}^* \end{bmatrix}.\end{equation}

The system (2.7) is derived by treating $\boldsymbol {x}$ and $\boldsymbol {\lambda }$ as a function of $\boldsymbol {\epsilon }$ and differentiating (2.4) and (2.5) with respect to $\boldsymbol {\epsilon }$. Importantly, $\boldsymbol {D}$ depends on second derivatives of the objective functions. If $\boldsymbol {D}$ is invertible, then the implicit function theorem guarantees existence of the expansion (2.3), and a unique solution exists to (2.7) (Krantz & Parks Reference Krantz and Parks2002).

Invertibility of the Jacobian matrix $\boldsymbol {D}$ occurs naturally when $\boldsymbol {x},\boldsymbol {\lambda }$ satisfy the KKT conditions, and $\boldsymbol {x}$ is a strict local minimum: $\nabla ^2 F_j(\boldsymbol {x}) + \sum _{i\neq j}\lambda _i \nabla ^2 F_i(\boldsymbol {x})$ is positive definite over the orthogonal subspace to $\{\boldsymbol {\nabla } F_i\}$ for all $i\in \mathcal {A}$ such that $\lambda _i>0$. However, the expansion does not exist if there are no active constraints: if all $F_i(\boldsymbol {x}) < \epsilon _i$, then a small variation in any $\boldsymbol {\epsilon }_i$ should not affect the solution. For example, in figure 1, the predictor–corrector method cannot be used beyond the tips of either branch of the efficient set since the KKT conditions cease to hold. Furthermore, if some constraint $i$ is active but the associated Lagrange multiplier $\lambda _i =0$, then a small increase in $\epsilon _i$ will not change the solution to (2.2), but a small decrease would. In this case, $\lambda _i$ is only sub-differentiable (Clarke Reference Clarke1990) with respect to $\epsilon _i$. While an expansion can be derived in this case, we ignore it for simplicity, as this case only appears at the edges of the efficient set or at saddle points of $F_j$. We are primarily concerned with the case when all active constraints have strictly positive Lagrange multipliers, since in this case, there is a trade-off to be made between objectives: reduction in $F_j$ comes at the cost of increase in $\epsilon _i$.

Now that we have established when a local expansion of the efficient set exists, and how to compute it, we can discuss the predictor–corrector method. The predictor step of the predictor–corrector method takes as input an efficient point $\boldsymbol {x}$ and a step size ${\rm \Delta} \boldsymbol {\epsilon }$. Subsequently, we compute the Lagrange multiplier $\boldsymbol {\lambda }$. If some components of $\boldsymbol {\lambda }$ are strictly positive, then we compute $\boldsymbol {D}(\boldsymbol {x},\boldsymbol {\lambda })$, solve (2.7) for $\nabla _{\boldsymbol {\epsilon }}\boldsymbol {x}$, then predict a weakly efficient point via the expansion, (2.3).

When a simulation is used to compute $\boldsymbol {F}$, it may not be possible to evaluate $\boldsymbol {F}$ at the predicted point since the simulation can ‘fail’. For example, the Variational Moments Equilibrium Code (VMEC) may fail to converge if the predicted point represents an ideal MHD equilibrium with self-intersecting flux surfaces. To proceed with the corrector step, the predicted point must evaluate under $\boldsymbol {F}$. To ensure that such a point can be found, the predicted point can be modified using a backtracking line search (Nocedal & Wright Reference Nocedal and Wright1999). The backtracking procedure iteratively reduces ${\rm \Delta} \boldsymbol {\epsilon }$ to $\gamma {\rm \Delta} \boldsymbol {\epsilon }$ for $\gamma \in (0,1)$ and predicts the value of $\boldsymbol {x}(\boldsymbol {\epsilon }+{\rm \Delta} \boldsymbol {\epsilon })$, until $\boldsymbol {F}$ can be evaluated at the predicted point $\boldsymbol {x}(\boldsymbol {\epsilon }+{\rm \Delta} \boldsymbol {\epsilon })$ or ${\rm \Delta} \boldsymbol {\epsilon }$ is smaller than some tolerance $\tau _{\varDelta }$.

For small enough ${\rm \Delta} \boldsymbol {\epsilon }$, the prediction will typically be a very good estimate of an efficient point. Nonetheless, it is helpful to finely resolve the point, particularly as this allows us to compute the subsequent predictor step accurately. The corrector step corrects the prediction by simply solving the $\boldsymbol {\epsilon }$-constraint problem, (2.2), with new constraint parameters $\boldsymbol {\epsilon } + {\rm \Delta} \boldsymbol {\epsilon }$. This method is detailed in Algorithm .

2.2.1. Computation of the Lagrange multiplier

In practice, we rarely have a point $\boldsymbol {x}$ that exactly satisfies any constraint with equality. Typically, an inequality constraint may approximately be satisfied with equality, i.e. $\overline {F}^j(\boldsymbol {x}) = \boldsymbol {\epsilon } - \nu$ for a small value of $\nu > 0$. In this case, theoretically the constraint is not active and so the Lagrange multipliers are equal to zero. However, the numerics suggest that the constraint is active and the Lagrange multiplier should not be zero. In practice, we determine the active set $\mathcal {A}$ as the set of all constraints with $\nu _i< \tau$ for some small tolerance $\tau$. The Lagrange multiplier can then be computed by solving the KKT conditions for a perturbed version of (2.2) where the constraint right hand side, $\epsilon _i$, is shifted to $\epsilon _i - \nu$ for all constraints in $\mathcal {A}$. In this way, the constraints hold with exact equality, since $\overline {F}^j(\boldsymbol {x}) = \boldsymbol {\epsilon } - \nu$, and we are justified in computing the potentially non-zero Lagrange multipliers. The Lagrange multipliers are then estimated as the solution to

(2.8)\begin{equation} \min_{\boldsymbol{\lambda}\geq 0} \left\|\boldsymbol{\nabla} F_j(\boldsymbol{x}) + \sum_{i\in\mathcal{A}}\lambda_i\boldsymbol{\nabla} F_i(\boldsymbol{x})\right\|^2.\end{equation}

Importantly, this procedure affects the prediction step in that the expansion is taken around $\boldsymbol {\epsilon } + \nu$. However, this is accounted for in step 1 of Algorithm .

3.1. Problem 1: the aspect ratio and quasi-symmetry trade-off

We seek a plasma boundary shape, parametrized by $\boldsymbol {x}$, of an ideal MHD equilibrium for a quasi-helical (QH) stellarator configuration that has minimal aspect ratio $A(\boldsymbol {x})$ and deviation from quasi-symmetry $Q_{M,N}(\boldsymbol {x})$,

(3.1)\begin{equation} \min_{A(\boldsymbol{x})\in[A_l,A_u]} (A(\boldsymbol{x}),Q_{M,N}(\boldsymbol{x})).\end{equation}

Including bound constraints on the aspect ratio $A_l=3, A_u=10$, restricts our decision space to a realistic range of configurations. For convenience, we collect the ${n_{\boldsymbol {F}}}=2$ objectives into the vector $\boldsymbol {F}(\boldsymbol {x}) :=(A(\boldsymbol {x}),Q_{M,N}(\boldsymbol {x}))$.

The violation of quasi-symmetry is defined as the quasi-symmetry-ratio-residual objective employed by Landreman & Paul (Reference Landreman and Paul2022). The objective measures the departure from quasi-symmetry throughout the plasma volume as the sum of flux surface averages across surfaces $(s_1,\ldots,s_{n_s})$,

(3.2)\begin{equation} Q_{M,N}(\boldsymbol{x}) = \sum_{s_j} \left\langle \left( \frac{1}{B^3}[(N-\iota M)\boldsymbol{B}\times\boldsymbol{\nabla} B\boldsymbol{\cdot}\boldsymbol{\nabla}\psi - (MG-NI)\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B]\right)^2\right\rangle.\end{equation}

The helicity parameters, $M,N$, determine the type of quasi-symmetry that $Q_{M,N}$ measures: $M=1,N=0$ being quasi-axisymmetry, and $M=1, N=\pm k{n_{\textrm {fp}}}$ with non-zero integers $k$ being quasi-helical symmetry (Landreman & Paul Reference Landreman and Paul2022). For our experiments, the helicity parameters were set to $M=1,N=-{n_{\textrm {fp}}}$ for the ${n_{\textrm {fp}}}=4$ field period magnetic field. Discretization of the flux surface average over the poloidal and toroidal angles results in an objective with sum-of-squares structure, $Q_{M,N}(\boldsymbol {x}) = \sum _{i=1}^{{n_{Q}}} Q_{M,N,i}^2(\boldsymbol {x})$, where $Q_{M,N,i}$ measures the violation of quasi-symmetry at a point throughout the volume. For a given plasma boundary, $\boldsymbol {x}$, we compute the aspect ratio using the definition in the VMEC code,

(3.3a–c)\begin{equation} A(\boldsymbol{x}) = R_{\text{major}}(\boldsymbol{x})/R_{\text{minor}}(\boldsymbol{x}) \quad R_{\text{minor}} = \sqrt{\bar{C}_A(\boldsymbol{x})/{\rm \pi}} \quad R_{\text{major}(\boldsymbol{x})} = \frac{V(\boldsymbol{x})}{2{\rm \pi}^2R_{\text{minor}}^2(\boldsymbol{x})}, \end{equation}

where $\bar {C}_A$ is the average cross-sectional area of the surface and $V$ is the volume enclosed by the surface.

The decision variables, $\boldsymbol {x}$, are Fourier amplitudes that describe the shape of the plasma boundary. The plasma boundary is represented in terms of the standard cylindrical coordinates $(R,\phi,Z)$ where $R,Z$ are parametrized as a Fourier series in the poloidal and toroidal angles $\theta$ and $\phi$,

(3.4)\begin{equation} \left.\begin{array}{c} \displaystyle R(\theta,\phi) = \sum_{n=0}^{{n_{\rm sm}}} R_{0,n}\cos(- n_{\text{fp}}n\phi) + \sum_{m =1}^{{n_{\rm sm}}} \sum_{n={-}{n_{\rm sm}}}^{{n_{\rm sm}}} R_{m,n}\cos(m\theta - n_{\text{fp}}n\phi),\\ \displaystyle Z(\theta,\phi) = \sum_{n=1}^{{n_{\rm sm}}} Z_{0,n}\sin(- n_{\text{fp}}n\phi) + \sum_{m =1}^{{n_{\rm sm}}} \sum_{n={-}{n_{\rm sm}}}^{{n_{\rm sm}}} Z_{m,n}\sin(m\theta - n_{\text{fp}}n\phi). \end{array}\right\}\end{equation}

The number of modes describing the surface ${n_\textrm {sm}} = \{0,1,2,\ldots \}$ can be increased to achieve more intricate boundary representations. Field period symmetry with ${n_{\textrm {fp}}}$ periods and stellarator symmetry have been assumed. The major radius $R_{0,0}$ is held fixed throughout the optimization to fix the scale of the design. The Fourier amplitudes are collected into the decision variable via $\boldsymbol {x} = (R_{0,1},\ldots, Z_{0,0},\ldots )$. The total number of decision variables satisfies ${n_{\boldsymbol {x}}} = 4{n_\textrm {sm}}^2 + 4{n_\textrm {sm}}$.

Numerical experiments were performed using SIMSOPT (Landreman et al. Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021) to handle variables, compute objectives and interface with VMEC which computed the ideal MHD equilibria from the plasma boundary representation. The $\boldsymbol {\epsilon }$-constraint method was used to find points along the Pareto front where quasi-symmetry was set to be the target function for minimization while the aspect ratio was constrained by $\epsilon$. The bound constraints on the aspect ratio made specifying $\epsilon$ straightforward: $\epsilon$ was set to linearly spaced values between $A_l$ and $A_u$. The $\boldsymbol {\epsilon }$-constraint problem (2.2) was solved by reformulating the constrained optimization problem with a quadratic penalty method (Nocedal & Wright Reference Nocedal and Wright1999), see § A. The penalty parameter was increased from an initial value of $1$ by a factor of $10$ at each iteration, and the quadratic penalty sub-problems were solved by applying a Gauss–Newton optimization routine. The optimization was run until the the change in $Q_{M,N}(\boldsymbol {x})$ was of the order of machine precision, $10^{-16}$, or the norm of the gradient reached $10^{-10}$. Upon taking a step with size of the order of $10^{-10}$, the optimization algorithm would terminate due to significant computational noise in the output of VMEC. SIMSOPT was used to compute forward difference gradients via MPI-based concurrent function evaluations. To avoid local minima, the number of Fourier modes, ${n_\textrm {sm}}$, was increased iteratively from $1$ to $5$, reaching ${n_{\boldsymbol {x}}}=120$ variables. The $\epsilon$-constraint problem was solved with the penalty method after each increase of ${n_\textrm {sm}}$.

To reduce the computational overhead of using the $\epsilon$-constraint method, we used the predictor–corrector method, introduced in § 2.2, to explore the Pareto front between solutions of the $\epsilon$-constraint method. For the high dimensional decision space of interest, the full Hessians $\nabla ^2 Q_{M,N}$ and $\nabla ^2 A$ are too expensive to compute with finite differences. Instead, we approximated the Hessian $\nabla ^2 Q_{M,N}$ with a Gauss–Newton Hessian approximation and use a diagonal second-order central difference approximation to $\nabla ^2 A$.

Figure 4 shows the Pareto front for the problem (3.1). Figure 5 shows three-dimensional renderings of the three solutions highlighted by the square, star and diamond markers in figure 4, as well as contour plots of the field strength in Boozer coordinates. From figure 4, it is clear that low values of the quasi-symmetry violation, $Q_{1,-4}(\boldsymbol {x}) < 10^{-3}$, can be achieved at all aspect ratios in our range. There is a slight trend indicating that quasi-symmetry may be achieved more precisely at higher aspect ratios. Nonetheless, by viewing the contour plots of the magnetic field strength in Boozer coordinates in figure 5, we see that configurations from all parts of the Pareto front have visibly precise quasi-symmetry.

Figure 4. Pareto front of the aspect ratio and quasi-symmetry objectives over the domain $A_l \le A(\boldsymbol {x}) \le A_u$ (black points). The configurations corresponding to the blue star, orange square and green diamond markers are plotted in figure 5. We find that the efficient set undergoes a branch change as the aspect ratio crosses the large gap from $\approx 6.11$ to $\approx 8.5$. The configurations in this range are not Pareto optimal and are shown in grey.

Figure 5. (Left columns) Contour plots of the magnetic field strength in Boozer coordinates $(\theta,\zeta )$ on four flux surfaces $s=0,0.25,0.5,1.0$ for three Pareto optimal configurations denoted by a star (aspect $3.5$), square (aspect $5.6$) and diamond (aspect $8.8$) in figure 4. (Right column) Three-dimensional (3-D) renderings of corresponding Pareto optimal designs. The colour of the 3-D configurations indicates the field strength, where red is stronger.

While it seems that precise quasi-symmetry can be achieved at all aspect ratios considered, it is not clear how well the reduction in quasi-symmetry translates to an improvement in particle confinement. After all, achieving quasi-symmetry is a proxy for the true goal of achieving good confinement. To this end, we computed the fraction of alpha particles lost from $10$ distinct Pareto optimal configurations, with aspect ratios 3, 3.5, 4, 4.5, 4.9, 5.6, 6.0, 8.5, 8.7 and 9. To compute the losses, we scaled the Pareto optimal configurations to the ARIES-CS reactor (Najmabadi et al. Reference Najmabadi, Raffray, Abdel-Khalik, Bromberg, Crosatti, El-Guebaly, Garabedian, Grossman, Henderson and Ibrahim2008) scale ($1.7$ m minor radius and $B_{0,0}(s=0) = 5.7$ T field strength on axis), and traced $5000$ particles born on the $s=1/4$ flux surface as well $5000$ particles born on the $s=1/2$ flux surface until a terminal time of $0.1$ seconds. The reactor scaling and particle count were also used by Bader et al. (Reference Bader, Drevlak, Anderson, Faber, Hegna, Likin, Schmitt and Talmadge2019) and Landreman & Paul (Reference Landreman and Paul2022). Particles were traced according to the vacuum guiding center equations in Boozer coordinates and were deemed lost if they crossed the $s=1$ flux surface. The loss fractions are shown in table 1. Not a single particle born on the $s=1/4$ flux surface was lost from nine of the ten configurations, and the remaining configuration, with aspect ratio $5.6$, lost only $0.22\,\%$ of the alpha particles. A slightly larger fraction of the alpha particles born on the $s=1/2$ flux surface were lost, between $0\,\%$ and $2.6\,\%$for each configuration. The loss fraction of alpha particles born on the $s=1/2$ flux surface increases slightly as the aspect ratio increases up until $A\approx 6$, then drops to approximately zero on the right tail of the Pareto front. It seems that there is no clear trend between the aspect ratio of the Pareto optimal designs and the loss fractions for particles born on either of the two surfaces. Thus, while there is a slight trade-off between quasi-symmetry and aspect ratio, the trade-off between confinement and aspect ratio may be more complex.

Table 1. Fraction of alpha particles lost for ten Pareto optimal configurations from figure 4. The loss fractions were computed by sampling $5000$ particles on the $s=1/4$ flux surface or $s=1/2$ flux surface and evolving their trajectories by the vacuum guiding centre equations in Boozer coordinates until a terminal time of $0.1$ seconds or the particle breached the $s=1$ surface and was lost.

Spatially, we may expect or hope that the set of desirable configurations forms a compact, connected region. However, we find the contrary – the efficient set is not connected. Qualitatively, the large gap in the Pareto front in figure 4 hints that the efficient set is undergoing a ‘branch change’. Branched curves, or branches, are parametric families of solutions to (2.4)–(2.6) which can be computed by continuation. Two branches can either be connected by a non-smoothness or be entirely disconnected sets, as in the case of figure 1. Branch changes are reflected in the Pareto front via a non-smoothness or gap where the Pareto front transitions from representing points on one branch to representing points on another branch. Quantitatively, we measure the branch change by evaluating the Lagrange multipliers at the left edge of the gap in the Pareto front (roughly $A = 6.11$) in figure 4; a branch change is indicated by $\lambda = 0$, since this implies that the aspect ratio can be increased without causing a reduction in the quasi-symmetry. Indeed, we find that the Lagrange multipliers are approximately zero and hence that the efficient set corresponding to the small piece of the Pareto front, with aspect ratio greater than or equal to $8.5$, is disconnected from the remainder of the efficient set. Within the gap in the Pareto front, configurations still can achieve exceptional levels of quasi-symmetry, as shown by the grey points in figure 4. The increase in quasi-symmetry in the gap may not be a physical phenomena, but may instead be due to poor numerical refinement in VMEC, local optimization or the representation of decision variables. It is interesting to see nonetheless that the set of desirable configurations is disparate.

All the Pareto optimal configurations found achieve exceptional levels of quasi-symmetry and particle confinement. From this, it is clear that the aspect ratio is not limiting the device performance. However, engineering criteria on coil shape, fabrication tolerances and placement tolerances may exhibit a greater trade-off with particle confinement, since the magnetic field generated by coils may not align well with the target magnetic field computed by ideal MHD. In the next section, we consider the trade-off between coil length and the ability for coils to reproduce a target magnetic field.

3.2. Problem 2: the coil length and quadratic flux trade-off

Given the shape of the last close flux surface, $\mathcal {S}$, the ‘stage-two’ stellarator optimization problem seeks to find magnetic coils with magnetic field $\boldsymbol {B}$, such that the field is orthogonal to the surface normal, $\boldsymbol {B}\boldsymbol {\cdot } \boldsymbol {n}=0$ (Merkel Reference Merkel1987; Zhu et al. Reference Zhu, Hudson, Song and Wan2017). Coils with a tight constraint on their length will not be able to form intricate shapes and reduce the normal component of the magnetic field to zero, whereas coils which are longer can form arbitrary shapes that better minimize $\boldsymbol {B}\boldsymbol {\cdot } \boldsymbol {n}$. Long coils, on the other hand, are undesirable from a financial and engineering standpoint: longer coils are more expensive to build, often have higher curvature and are more difficult to fit into the space around the device. In this experiment, we seek to understand the trade-off between coil length and the ability of coils to reproduce a target magnetic field.

We develop insight into this trade-off by solving the multi-objective problem

(3.5)\begin{equation} \min_{\boldsymbol{x}\in\boldsymbol{\varOmega}} (J_{\boldsymbol{n}}, L),\end{equation}

where $J_{\boldsymbol {n}}$ is the quadratic flux,

(3.6)\begin{equation} J_{\boldsymbol{n}}(\boldsymbol{x}) = \frac{1}{2}\int_{\mathcal{S}}(\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{n})^2 \, {\rm d} s, \end{equation}

$L$ is the total length of the ${n_{C}}$ coils considered,

(3.7)\begin{equation} L = \sum_{i=1}^{n_{C}} L_i, \end{equation}

and the decision variables $\boldsymbol {x}$ represent coil shape parameters and the coil currents. The quadratic flux is a standard metric used to find coils in the stage-two problem (Merkel Reference Merkel1987; Landreman Reference Landreman2017; Zhu et al. Reference Zhu, Hudson, Song and Wan2017; Singh et al. Reference Singh, Kruger, Bader, Zhu, Hudson and Anderson2020; Kruger et al. Reference Kruger, Zhu, Bader, Anderson and Singh2021; Glas et al. Reference Glas, Padidar, Kellison and Bindel2022; Wechsung et al. Reference Wechsung, Giuliani, Landreman, Cerfon and Stadler2022a). To solve the bi-objective coil problem, (3.5), we use the coil optimization framework implemented in SIMSOPT. The modular coils are described by a Fourier representation of the Cartesian coordinates with ${n_\textrm {cm}}$ modes. The Fourier representation of the $x$-coordinate of the $i$th coil is

(3.8)\begin{equation} x^i(t,\boldsymbol{x}) = x_{c,0}^i + \sum_{k=1}^{{n_{\rm cm}}}[x_{c,k}^i\cos(kt) + x_{s,k}^i\sin(kt)], \end{equation}

with analogous forms for the $y$ and $z$ coordinates. Each coil is additionally equipped with a current, though the current of the first coil is held fixed at $I_1 = 10^{5}$ A throughout the optimization. The number of Fourier amplitudes per coil is $3(2{n_\textrm {cm}} + 1)$ making the total number of design variables ${n_{\boldsymbol {x}}} = {n_{C}} 3(2{n_\textrm {cm}} + 1) + {n_{C}} -1$. The coils considered satisfy field-period symmetry, as well as stellarator symmetry. Thus, the quadratic flux is only computed over a half-field-period, the total coil length is only taken over the coils in a half-field-period and the decision variables $\boldsymbol {x}$ only represent the ${n_{C}}$ coils in a half-field-period. The plasma boundary shape $\mathcal {S}$ used in this experiment is that of the ${n_{\textrm {fp}}}=2$ field-period quasi-axisymmetric (QA) configuration from Landreman & Paul (Reference Landreman and Paul2022), which we will refer to as LP-QA. The ${n_\textrm {cm}}=5$ Fourier modes were used to describe the coil shapes.

To avoid ill-posedness of the objectives due to non-uniqueness of the underlying coil parametrization, a penalty on the variation of the arc length was used to encourage the coils to have uniform incremental arc length. The arc length variation objective discretizes each coil into $n_{\textrm {av}}$ segments and computes the variance of the set of arc lengths $\{l_i\}_{i=1}^{n_{\textrm {av}}}$ for each coil, i.e. for a single coil, $V =\textrm {Var}[\{l_i\}_{i=1}^{n_{\textrm {av}}}]$. The total arc length variation over the coil set is

(3.9)\begin{equation} J_{\text{av}}(\boldsymbol{x}) = \sum_{i=1}^{n_{C}} V_i.\end{equation}

The bi-objective problem was solved by applying the $\boldsymbol {\epsilon }$-constraint method, where quadratic flux was the objective and the total coil length was the constraint. The predictor–corrector method was not used to investigate this trade-off, since the computational expense of solving the coil optimization problem with the $\boldsymbol {\epsilon }$-constraint method was minimal. However, a warm-starting procedure or ‘zeroth-order’ predictor–corrector method, was used in which solutions of the $\boldsymbol {\epsilon }$-constraint method were used as a good initial guess of efficient points with similar values of $\boldsymbol {\epsilon }$. The $\boldsymbol {\epsilon }$-constraint problems were solved with a quadratic penalty method (Nocedal & Wright Reference Nocedal and Wright1999), see § A. The penalty parameter was increased from an initial value of $1$ by a factor of $10$ at each iteration, and the quadratic penalty sub-problems were solved by the L-BFGS-B algorithm (Zhu et al. Reference Zhu, Byrd, Lu and Nocedal1997) until the norm of the gradient reached $10^{-14}$ or $30\ 000$ iterations ellapsed. The arc length variation (3.9) was included as a penalty term with fixed penalty parameter $10^{-4}$.

Figure 6 shows the Pareto front for the bi-objective coil optimization problem (3.5) when ${n_{C}}=4$ coils are used per half-field-period. Figure 6 also shows three-dimensional renderings of the two Pareto optimal coils sets, highlighted as the star and diamond in the plot of the Pareto front. Figure 6 shows that increasing the coil length substantially improves the ability for coils to reduce the quadratic flux up until the average coil length over the minor circumference of the plasma boundary is approximately $6$, after which we numerically find no improvement. In addition, we find that when the average coil length over the minor circumference of the plasma boundary exceeds roughly $6$, the coils become exceedingly complex: coils curvature becomes large and coils begin to pass under one another, competing for space. At these coil lengths, coil curvature and the coil-to-coil separation should also be constrained.

Figure 6. (a) Pareto front for coil length versus quadratic flux. The $x$-axis is non-dimensionalized by dividing by the minor circumference of the surface ($2{\rm \pi}$ times the minor radius). (b) Three-dimensional rendering of the Pareto optimal coil sets, denoted by the blue star (smaller blue coils) and the orange square (larger orange coils).

While figure 6 shows that constraining the total coil length limits the ability for coils to reproduce a magnetic field, it does not show how the limits on the total coil length relate to a loss in quasi-symmetry or loss in confinement. In figure 7, we measure the extent to which the Pareto optimal coil configurations from figure 6 generate quasi-symmetric magnetic fields. To do so, we fit a quadratic flux minimizing (QFM) surface (Dewar, Hudson & Price Reference Dewar, Hudson and Price1994) to the coil-generated magnetic field $\boldsymbol {B}$, evaluate VMEC using the QFM surface as the boundary shape and compute the quasi-symmetric metric $Q_{1,0}$. The QFM surface is constrained to have identical volume to the LP-QA target surface. We find the QFM surface by solving the following problem in SIMSOPT:

(3.10)\begin{gather} \min_{S} \ \frac{\displaystyle \int_S(\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{n})^2 \, {\rm d} s}{\displaystyle \int_S\|\boldsymbol{B}\|^2 \, {\rm d} s}, \end{gather}

(3.11)\begin{gather}\text{such that} \ \text{Volume}(S) = \text{Volume}(\text{LP-QA}). \end{gather}

Figure 7 shows that the quasi-axisymmetry degrades significantly, by four orders of magnitude, as the coil length constraints are tightened. Unlike the longer coil sets, the shorter coil sets generate significant coil ripple, degrading the quasi-symmetry. While it may not be possible to entirely avoid the loss of quasi-symmetry when coil constraints are tightened, the loss may be diminished by using single-stage optimization approaches to design the coils. Part of the loss of quasi-symmetry is inherent to the two-stage approach to coil design: the coils are optimized to reproduce a target magnetic field, but not to generate a quasi-symmetric field. Single-stage optimization approaches (Giuliani et al. Reference Giuliani, Wechsung, Cerfon, Stadler and Landreman2022a,Reference Giuliani, Wechsung, Landreman, Stadler and Cerfonb ; Jorge et al. Reference Jorge, Goodman, Landreman, Rodrigues and Wechsung2023), however, directly optimize coils for quasi-symmetry, ensuring that they achieve optimal quasi-symmetry levels over the space of coils satisfying the coil length constraint. Single-stage approaches would hence improve the unfortunate trade-off between coil length and quasi-symmetry.

Figure 7. (a) Trade-off between quasi-axisymmetry and coil length for the Pareto optimal points shown in figure 6. For each Pareto optimal coil set, the quasi-axisymmetry metric $Q_{1,0}$ was computed by running VMEC with a boundary surfaced obtained by fitting a quadratic flux minimizing surface to the field generated by the coils. (b,c) Contour plots of the coil generated field strength in the Boozer toroidal and poloidal angles $\zeta,\theta$ on the $s=1$ surface of the LP-QA configuration for the coil sets of the blue star and orange square shown in figure 6.