Greedy permanent magnet optimization

Now that the permanent magnet optimization problem has been formulated as sparse regression, we show that the fully binary problem is equivalent to a quadratic program with quadratic equality constraints (QPQC), the binary, grid-aligned problem is equivalent to the quadratic knapsack problem with multiple knapsack constraints (MdQKP), and the binary, single-orientation problem is a binary quadratic program (BQP). These formulations provide a route for the future application of more sophisticated greedy algorithms and performance bounds to permanent magnet optimization.

We then provide a simple greedy algorithm to solve the optimization. Our implementation is computationally fast, intuitive, easy-to-use, and open-source. It a-priori produces sparse, binary solutions that are independent of the initial guess, and additionally illustrates that greedy algorithms can perform comparably to or even better than the state-of-the-art. We propose some likely reasons for the unusual effectiveness of greedy algorithms in the context of permanent magnet optimization in section 3.6. In the present work, there are essentially no hyperparameters, no local approximations, and no multi-stage operations to perform. All of the present work's methodology and results can be found in the SIMSOPT code [3]. As in our recent work, the entirety of the permanent magnet pipeline, i.e. the geometry, optimization tools, post-processing, etc is contained in this single open-source tool.

In order to facilitate an analysis of the truly discrete version of (1), we normalize the m _i by their maximum strengths $\boldsymbol{m}_i \to \boldsymbol{m}_i / m_i^{\text{max}}$, $\boldsymbol{A} \to\boldsymbol{A} \boldsymbol{m}_i^{\text{max}}$, and replace the l₀ norm in (1) with a restriction to $m_i\in\{-1, 0, 1\}$ on the optimization variables,

$$\begin{align} \mathop{\mathrm{arg\;min} }_{\boldsymbol{m}\in\{-1, 0, 1\}}&\frac{1}{2}\|\boldsymbol{A} \boldsymbol{m} - \boldsymbol{b}\|^2_2 + \lambda\|\boldsymbol{B}\boldsymbol{m}\|_2^2,\nonumber\\ \|\boldsymbol{m}_i\|_2^2 &\leqslant 1. \end{align} \tag{ 5 }$$

Here $\boldsymbol{B} = \text{diag}(\boldsymbol{m}_i^{\text{max}}) \in \mathbb{R}^{3D\times 3D}$ and $\boldsymbol{m}_i^{\text{max}} \in \mathbb{R}^{3D}$ is shorthand for the vector formed by triplets of the maximum strengths of each dipole, $[m_1^{\text{max}}, m_1^{\text{max}}, m_1^{\text{max}}, m_2^{\text{max}}, \ldots, m_D^{\text{max}}]$. Note that the grid-aligned property is enforced by the combination of the quadratic constraints and the discrete variables. For the remainder of the present work, all of the optimization variables m , A , etc will refer exclusively to the normalized quantities used in (5).

Despite the reduced form of (5), we have yet to connect it to the larger scientific literature detailing algorithms and performance bounds for constrained, binary minimization of the MSE. To do so, we transform the optimization in (5) to a variant of the well-studied quadratic knapsack problem (QKP) [28] with multiple constraints, MdQKP [29]. First, redefine the variables for optimization to $\boldsymbol{m} = [m^x_1, m^{-x}_1, m^y_1, m_1^{-y}, \ldots, m^{-z}_D]$, so that the problem can be restricted to binary variables. We can encode the grid-aligned property by trading in the quadratic constraints for affine constraints. After reshaping the A and B matrices, (5) becomes,

$$\begin{align} \mathop{\mathrm{arg\;min} }_{\boldsymbol{m}\in\{0, 1\}}&\left(\boldsymbol{m}^T\left[\frac{1}{2}\boldsymbol{A}^T\boldsymbol{A} + \lambda \boldsymbol{B}^T\boldsymbol{B}\right]\boldsymbol{m} - \boldsymbol{b}^T\boldsymbol{A}\boldsymbol{m}\right),\nonumber\\ \sum_{j=1}^{6D}C_{ij}m_{j} &\leqslant 1, \quad \forall i \in \{1, \ldots, D\}. \end{align} \tag{ 6 }$$

Here a choice for C_ij encodes a single constraint per dipole,

$$\begin{align} 0 &\leqslant (m_{i}^{x} + m_{i}^{-x}) + (m_{i}^y + m_{i}^{-y}) + (m_{i}^z + m_{i}^{-z}) \leqslant 1. \end{align} \tag{ 7 }$$

Along with $m_i \in \{0, 1\}$ for all i, the constraints are sufficient for grid-aligned dipoles, and (6) is now in the form of a MdQKP.

The largest MdQKP problem addressed by CPLEX's mixed-integer quadratic program solver in 2012 from Wang et al [30] was of size D = 800 with 15 constraints and it took a few hours to compute on a personal computer. In the permanent magnet optimization scenario with $D \sim 10^5$ and ∼10⁶ constraints, we need a much simpler algorithm.

Before we get to a simpler algorithm, consider the scenario in which the dipoles are a-priori desired to be in a particular coordinate direction, as in the minor radial direction in the MUSE example in section 3. We can then use the optimization variables $\boldsymbol{m} = [m^r_1, m_1^{-r}, m^{r}_2, \ldots, m^{-r}_D]$ and (6) becomes

$$\begin{align} \mathop{\mathrm{arg\;min} }_{\boldsymbol{m}\in\{0, 1\}}&\left(\boldsymbol{m}^T\left[\frac{1}{2}\boldsymbol{A}^T\boldsymbol{A} + \lambda \boldsymbol{B}^T\boldsymbol{B}\right]\boldsymbol{m} - \boldsymbol{b}^T\boldsymbol{A}\boldsymbol{m}\right). \end{align} \tag{ 8 }$$

The reduced constraint $0 \leqslant m_{i}^{r} + m_{i}^{-r} \leqslant 1$ is not required because f_m prevents $m_{i}^{r} = m_{i}^{-r} = 1$. Assuming B has been defined appropriately, i.e. so that $f_m \propto (m^r_i)^2 + (m^{-r}_i)^2$, then the objective can be reduced by changing both $m_i^r$ and $m_i^{-r}$ from 1 to 0. This new solution leaves the f_B term unaffected, while the f_m term is reduced.

Now (8) is the unconstrained quadratic knapsack problem, which is to say it is the unconstrained BQP. Despite the significant reductions, this problem is NP-hard and exact methods of solution have been limited to problems of a few hundred variables. Even the heuristic and greedy methods reported in Kochenberger et al [31] were limited to $D \lesssim 10^4$ or so. Part of the reason for this limit is that often quite sophisticated heuristic methods are used, such as methods incorporating tabu search [32], local search schemes [33], and so on [34]. Similarly, variants with tabu search [35, 36] and genetic algorithms [37] are frequent choices for heuristically solving variants of QKP, but these algorithms are also typically applied to problems with a few hundred optimization variables.

In contrast, we start with a very simple and intuitive greedy algorithm, add some basic backtracking, and already can compete with or even outperform the state-of-the-art algorithms for binary permanent magnet optimization with $D \sim 10^5$. Similar algorithms have been used with some success in the BQP literature for decades [38, 39], but have not been applied for permanent magnet optimization. Reasons for the unusual effectiveness of these simple algorithms in binary permanent magnet optimization are explored in section 3.6. Lastly, the significant literature on heuristic and greedy algorithms for BQPs and QKPs provides a road map for improvements in future algorithms.

Similar in spirit to the binary matching pursuit algorithm in Wen and Li [40], we devise a greedy, binary algorithm for reducing the mean-squared error (MSE) in permanent magnet optimization. This algorithm may be used with objectives other than the MSE, but minimizing other objectives does not typically provide guarantees on the smallness of the MSE. This is problematic, because the MSE in permanent magnet optimization must be greatly reduced for a solution to be useful. A stage-two stellarator optimization must match the normal magnetic field of a stage-1-optimized plasma boundary very closely in order to avoid producing significant deviations from the stage-1 flux surfaces and other optimized physical quantities such as the fast-particle confinement.

Initially, the greedy approach may seem easy to implement but unlikely to generate high-quality solutions. Calculating performance guarantees, e.g. the f_B ratio between the true optimum and greedy solution, for greedy algorithms is not straightforward. However, a greedy solution that produces 100 times larger f_B than the optimum might still be a very useful solution, since permanent magnet grids can sometimes be designed that exhibit solutions satisfying $f_B\sim 10^{-10}$ T² m². A brief review of the relevant mathematical concepts for greedy performance guarantees can be found in appendix. At least in the case of continuous optimization, recent results have shown that, when the MSE is used as the objective function, greedy algorithms retain some reasonable performance guarantees [41, 42]. These guarantees are one justification for using the MSE as an objective for greedy minimization; another justification is simply to show high performance is available using the MSE, as we do in section 3.

Our algorithm is as simple and intuitive as possible. One-by-one, place a maximum-strength permanent magnet (with only one nonzero component in order to be grid-aligned) that maximally reduces f_B . Then add all of the components of this magnet to the list of indices that are off limits to future magnet placement. This process continues until K magnets have been placed. Once a magnet is placed on the grid, it is immutable, so additional algorithm iterations can only add new magnets until the grid has no remaining unfilled locations. We will address this potential issue with our modified algorithm in the following section, which does allow for backtracking during the optimization process to correct for suboptimal magnet placements. The code also allows users to save a solution with K magnets, and use it later as an initial condition for a new optimization that places more than K magnets.

If the permanent magnets are all identical, as is the case for the MUSE grid used in section 3, minimizing $f_B + f_m$ is not required after rescaling because each dipole contributes an identical amount to the regularization f_m . If the grid does not consist of identical elements, the effect of f_m is to preference smaller magnets that contribute most significantly to the solution. Algorithm 1 summarizes the greedy permanent magnet optimization (GPMO) approach.

. 

Algorithm 1. GPMO.
Input: All the physics encoded in fB and $K \equiv$ number of magnets
to place. Initialize $\Gamma^{(1)} = \{1,\ldots,3D\}$, $\boldsymbol{m}^{(1)} = \boldsymbol{0}$.
Output: Solution $\boldsymbol{m}^{(K)}$.
procedure GPMO( A , b )
for $k = 1,\ldots,K$
$i^* = \mathop{\mathrm{arg\;min} }_{i \in \Gamma^{(k)}}f_B(\boldsymbol{m}^{(k)} | m_i = \pm 1) + f_m(\boldsymbol{m}^{(k)} | m_i = \pm 1)$
$m_{i^*}^{(k+1)} = \pm 1$,
$\Gamma^{(k + 1)} = \Gamma^{(k)} - \{i^*, \text{remaining components of the dipole}\}$,
return $\boldsymbol{m}^{(K)}$
end procedure

In words: find the component $i^*$ of m , that when set to ±1, minimizes f_B the most. In order to generate grid-aligned solutions, subtract the chosen $i^*$ index and the remaining components of the same dipole from Γ. Continue until K maximum-strength, grid-aligned dipoles have been placed.

In a greedy approach, we can use the original optimization variables $\boldsymbol{m} = [m^x_1, m^y_1, \ldots, m^{z}_D]$ because after each magnet placement, we can eliminate the related variable components from consideration. The computational bottleneck of iteration k is evaluating $\boldsymbol{A}\boldsymbol{m}$, which can be trivially parallelized with OpenMP. Moreover, the result of the previous matrix-vector product can be stored each iteration, so that in the next iteration only a simple sum is required for checking the contribution of dipole j. For concreteness, define $\Gamma^{(k)}$ as the set of available magnet locations during iteration k and suppose $\Gamma^{(k + 1)} = \Gamma^{(k)} - \{j\}$. Then the contribution of the jth element of m to f_B during iteration $(k+1)$ is,

$$\begin{align} \sum_{l\in\Gamma_c^{(k+1)}} A_{il}m_l^{(k+1)} = A_{ij}m_j + \sum_{l\in\Gamma_c^{(k)}}A_{il}m_l^{(k)}, \end{align} \tag{ 9 }$$

where $\Gamma_c$ denotes the complement of Γ and the right-most term is known from the previous iteration. The worst case computational complexity is ${\sim}{3NDK}$ floating-point operations, parallelized over P OpenMP processes, and the calculation speeds up as K → D. On the Cori KNL machine, placing a very large number ($K = 40\,000$) of magnets for the MUSE stellarator with $D = 75\,460$ potential magnet locations, $N = 64 \times 64$ quadrature points on the plasma boundary, and P = 68 OpenMP processes, takes about 20 min. Given a prescribed permanent magnet grid (and therefore D), the algorithm scales well from the linear dependence on N and K and straightforward parallelization. A computational speed comparison with other algorithms is illustrated in section 3.5.

Finally, GPMO can be modified to solve the binary, arbitrary-orientation formulation in (4). The optimal dipole to place during iteration $(k + 1)$ is the dipole that produces the smallest value from solving

$$\begin{align} \min_{\boldsymbol{m}_i}&\left[f_B(\boldsymbol{m}^{(k)} | \boldsymbol{m}_i \neq 0) + f_m(\boldsymbol{m}^{(k)} | \boldsymbol{m}_i \neq 0)\right],\nonumber\\ &\|\boldsymbol{m}_i\|_2^2 = 1. \end{align} \tag{ 10 }$$

Computing this optimal location requires solving $(D - k)$ problems in the form of (10) every iteration. Each problem is a size-three QP with a single quadratic constraint. In the single constraint case, each QP can be transformed into a linear program and efficiently solved [43]. We leave further exploration of this route to future work.

A useful engineering constraint to implement in the GPMO algorithm is the placement of multiple magnets per iteration. Multiple index selections per iteration is now a common variant of the orthogonal matching pursuit and other greedy algorithms in the sparse regression field [44].

This is interesting from an engineering perspective because Algorithm 1 (GPMO) tends to produce solutions that exhibit a number of isolated magnets. This property is disadvantageous from a cost and engineering perspective, and often large portions of the surface of the plasma experiment need to be accessible for diagnostics. Instead of placing magnets one-by-one, Algorithm 1 is modified to place the set of p magnets that most minimizes f_B at each iteration. The set of p magnets can be additionally required to be a continuous set of magnets so no isolated magnets are ever placed. However, there are a combinatorically large number of ways to divide the grid into sets of p magnets, and within each set, an optimization over the p magnets may still be required to best orient the p magnets.

For simplicity, we choose a partitioning in which, at each iteration, a dipole is placed along with its closest available N_adjacent neighbors (some neighbor dipoles may have already been placed in a previous iteration so the algorithm sometimes chooses N_adjacent somewhat more distant neighbors). For convenience, we define N_adjacent by including the dipole in question, so that $N_\text{adjacent} = 1$ refers to using only the selected dipole. All of the neighbor dipoles are used with the same coordinate-aligned orientation as the selected dipole, and the dipole chunk chosen at each iteration is the one that minimizes f_B . This is a sort of coarse-grained approach of setting down large magnets all at once. We denote this algorithm GPMOm and if $N_\text{adjacent} = 1$ GPMOm reduces to the GPMO algorithm. N_adjacent is the only hyperparameter.

There are a number of other ways to dramatically reduce the number of isolated dipoles, such as adding a loss term that penalizes dipole locations with no near neighbors, or requiring the dipoles must be placed next to existing dipoles after some K_min number of dipoles have been placed with the ordinary GPMO algorithm. However, unlike GPMOm, these algorithms cannot guarantee that there will be exactly zero isolated dipoles in the final solution.

Algorithm 1 is a greedy algorithm and therefore will sometimes make poor choices for the magnets. These poor choices are often corrected later in the optimization when many more magnets have been placed. In fact, it can be seen at around $K \sim 10\,000$ in figure 1 that the r-direction-only GPMO solution briefly outperforms the GPMO solution with all $(r, \phi, \theta)$ directions available. The latter algorithm cannot have generated an optimal solution.

Zoom In Zoom Out Reset image size

The ill-conditioned nature of this problem actually facilitates a natural route for error correction; each unhelpful magnet is approximately zeroed out by placing an equal and opposite magnet directly next to it. These close cancellations are somewhat common, especially if many magnets have been placed, $K \sim D$. A natural way to improve these solutions is to perform some backtracking.

After each set of several hundred iterations, the algorithm can prune magnets that are adjacent to another magnet of opposite orientation. The pruning of a close-cancelling pair of dipoles barely changes the f_B metric and is a systematic way to backtrack and correct poor magnet choices made earlier in the algorithm. The backtracking works as follows: N_adjacent neighbors are obtained for each existing dipole, and then a search is made for the nearest oppositely oriented neighbor dipole. If one is found, the neighbor and the existing dipole are both removed from the solution, but the locations are allowed to host magnets in the future.

We denote this algorithm as GPMOb, and there are now in principle two hyperparameters: the maximum distance to consider between two oppositely oriented dipoles to be eliminated (here determined implicitly by choosing how many neighbor dipoles are considered adjacent via the N_adjacent hyperparameter), and how often to perform the backtracking, K_b . However, it was empirically found that the f_B performance was insensitive for values of $K_b \lesssim 500$ since this meant the backtracking was performed frequently enough to prevent large accumulations of poorly placed magnets. In the present work, the value $K_b = 200$ was used for generating all of the GPMOb results.

The greedy approach allows for useful engineering constraints to be straightforwardly incorporated into new permanent magnet configurations. Here, we briefly point out another opportunity with the greedy algorithms. One potentially advantageous property is to generate a configuration that is built from one continuous block of magnetic material. This can be generated greedily by, after placing the first magnet, always placing magnets directly next to existing magnets. With stellarator and field-period symmetries, this means the permanent magnet configuration is simply $2\times N_{fp}$ geometrically identical permanent magnet blocks.

The greedy algorithm performs remarkably well. Figure 1 illustrates that GPMO obtains comparable (but marginally worse) f_B performance with the same number of magnets as the FAMUS and relax-and-split solutions. The GPMO solution with $K \sim 12\,500$ and only r-aligned dipoles looks very similar to the FAMUS solution. But GPMO has no hyperparameters to tune and the solution took ∼5 min to compute. Its simplicity facilitates the exploration of more exotic configurations such as φ-aligned or θ-aligned magnet configurations. A reasonably accurate solution $f_B \sim 6\times 10^{-7}$ T² m² can even be generated with an all-φ-aligned configuration.

Consider the remarkable fact that the GPMO, FAMUS, and relax-and-split solutions qualitatively match. The three algorithms have very different degrees of freedom, define different optimization problems, and solve the problem with disparate iterative approaches. There are two effects that seem like plausible reasons for this convergence between algorithms. First, certain grid locations may be essential for properly minimizing f_B , e.g. locations near highly-shaped magnetic fields on the plasma boundary. Second, the requirement on each of the algorithms that the end solution must consist of maximum strength, binary magnets effectively regularizes much of the permanent magnet optimization space. In other words, there may be a vast number of permanent magnet configurations that minimize f_B to high-performance levels, but far fewer such configurations have binary distributions. The specific reasons for the high performance of GPMO is further remarked upon in section 3.6.

GPMO also shows that after ∼15 000 magnets have been placed on this grid, placing additional magnets has fairly negligible impact on f_B . There appears to be a f_B floor to the MUSE solutions, which we speculate comes from the binary nature of the magnets. The intuition that the solution can always be marginally improved by placing another magnet appears likely to hold only if the magnet directions and magnitudes can vary continuously. The greedy solution with $40\,000$ magnets is illustrated in figure 1 to show that much more complicated permanent magnet configurations are produced by many more φ and θ-aligned dipoles, but these changes have almost no impact on f_B .

Lastly, for all the optimizers, figure 1 exhibits some localized and regular $\boldsymbol{B} \cdot \hat{\boldsymbol{n}}$ residuals on the plasma surface. The two local peaks at φ = 0 and $\phi = \pi$, clearly visible in the FAMUS and GPMO radial-direction-only results, come from the MUSE grid. The grid has diagnostic ports at these locations, and it is challenging to compensate for these holes in the grid with radial dipoles placed nearby.

The plasma surface also exhibit stripes and other regular structures in most of the $\boldsymbol{B} \cdot \hat{\boldsymbol{n}}$ residual plots, but it is hard to conclude about the source of these patterns. We posit that these patterns are once again reflective of the available locations on the MUSE grid, which is a set of nested toroidal shells. A more geometrically complex grid that mimics the plasma surface might provide permanent magnet locations that can better address these errors.

However, some localized errors are probably an inevitable consequence of using permanent magnets to minimize the global f_B metric. The inverse cubic fall-off of a dipole field fundamentally implies that the magnets can only contribute quite localized fields to the plasma surface. It would be interesting and feasible to penalize additional terms related to, for instance, the smoothness of the $\boldsymbol{B} \cdot \hat{\boldsymbol{n}}$ residual on the plasma surface since this would amount to another linear least-squares term in the optimization. Another way to potentially address this issue is to use the full near-field expressions for permanent magnets. For rectangular cubes, the near-field expressions remain linear in the dipole moments, and therefore the optimization structure remains unchanged [45].

Next, we compare the GPMOm approach against the baseline algorithm. Recall that the the motivation for GPMOm is to produce solutions with very few isolated dipoles. Figure 2 illustrates the results using $N_\text{adjacent} = \{1,2,4,7, 20\}$. The f_B performance decreases fairly rapidly with increasing N_adjacent because we are performing an increasingly coarse-grained optimization. However, if we regard each set of N_adjacent, identically grid-aligned magnets as a single larger permanent magnet, f_B actually improves with increasing N_adjacent, although the f_B performance floor still decreases. Examining the m_r component of the GPMOm permanent magnet solutions shows that increasing N_adjacent does lead to bigger magnetic chunks used in the solution. Interestingly, as N_adjacent increases, so does the number of dipoles that have an equal and oppositely-oriented dipole directly adjacent, further motivating the backtracking approach from section 2.5. This property of the solutions may in fact be responsible for the reduced f_B performance at larger N_adjacent; it is easier to make greedy mistakes with a coarser grid of larger magnets.

Zoom In Zoom Out Reset image size

Here, we show that the greedy algorithm with backtracking actually consistently outperforms the state-of-the-art algorithms on the MUSE grid, and further below we repeat this achievement in section 3.5 on another permanent magnet configuration. Backtracking near cancellations, or even not-so-near cancellations, is a very effective way to generate high-quality solutions from greedy permanent magnet optimization. Systematically increasing N_adjacent will sweep out a valley of solutions, since initially increasing N_adjacent will improve solutions from the backtracking, but the limit $N_\text{adjacent} \to D$ will eventually cause reduced performance, since this limit requires that the algorithm use only three (instead of 6) possible directions. In other words, it will generate configurations with suboptimal directionality, e.g. producing dipoles that only point in the positive r direction, the $-\theta$ direction, and the $-\phi$ direction.

Figure 3 illustrates the results of backtracking with different values of the hyperparameter N_adjacent. Small but important improvements to the MUSE solutions are available with GPMOb, such that, given the same number of magnets, GPMOb outperforms the binary FAMUS and relax-and-split solutions. Even with a very large value of $N_\text{adjacent} = 4000$, GPMOb can find an accurate solution. Furthermore, much larger improvements of several orders of magnitude reduction in f_B appear accessible in other permanent magnet configurations, such as the NCSX example detailed below.

Zoom In Zoom Out Reset image size

Lastly, we provide a comparison of the CPU run times for the FAMUS, relax-and-split, and GPMO algorithms for obtaining a target $f_B \sim 10^{-7}$–10⁻⁸ T² m². The CPU run time is approximated as wall time × the number of CPUs, and allows us to compare run times while controlling for parallelization. We make no claims that any of the algorithms are optimally calculated or could not be sped up significantly. Conclusions about general computational speed are quite difficult to make because of a number of complications detailed below.

Significant variations in algorithm speed can occur simply from hyperparameter changes. For instance, the NCSX example described in section 3.5 is a particularly bad case for the relax-and-split algorithm, which currently relies on a suboptimization using projected gradient descent (PGD); PGD can progress very slowly when the solution lies near the constraint boundary. However, in the nonbinary, convex scenario, adding some Tikhonov regularization to the relax-and-split algorithm moves the optimum into the feasible region and then the algorithm is quite rapid.

For a direct comparison with GPMO, grid-aligned and binary FAMUS and relax-and-split solutions are presented, although generating binary solutions typically requires significant additional iterations for these algorithms. It follows that the run time estimates in table 1 for FAMUS and relax-and-split could be significantly sped up if fully binary solutions are not required. Despite these caveats, the takeaway is that GPMO calculates an accurate, binary, grid-aligned, MUSE solution much faster than the other algorithms. GPMO computes an accurate r-direction-only solution at least 10–15 times faster than FAMUS and an accurate grid-aligned solution at least 100 times faster than relax-and-split.

Table 1. Approximate CPU run times (calculated roughly as wall time × number of CPUs to control for different parallelization) for the different algorithms to achieve $f_B \sim 10^{-7}$–10⁻⁸ T² m² using the MUSE permanent magnet manifold.

Algorithm	Time (h)	Binary	Grid-aligned
GPMO	13	$\checkmark$	$\checkmark$
GPMOb	27	$\checkmark$	$\checkmark$
GPMO (r-only)	6	$\checkmark$	$\checkmark$
FAMUS (r-only)	100	$\checkmark$	$\checkmark$
Relax-and-split (convex)	1	×	×
Relax-and-split	1000	$\checkmark$	$\checkmark$

To conclude the results section, we illustrate an example for which no algorithm to date has been able to generate an accurate, binary solution. This list includes the GPMO algorithm, and seems to suggest fairly strong evidence that there are permanent magnet grids that do not support binary solutions that can minimize f_B to tolerable levels. There are a number of potential solutions available by altering the grid: defining grids that are more complex, closer to the plasma, or allow for more permanent magnet volume. Weakly breaking the binary constraint appears to be another way to obtain more accurate solutions. However, we show below that the greedy algorithm with backtracking can significantly improve the GPMO solution to tolerable levels, and we obtain the best binary solution for this grid by over an order of magnitude in f_B .

NCSX was a planned quasi-axisymmetric stellarator that was partially built at the Princeton Plasma Physics Laboratory. It was originally designed with 18 modular coils and 18 planar coils. The equilibrium of interest, C09R00, was scaled to have an on-axis magnetic field strength of 0.5 T, which is the maximum field produced by the existing planar coils. C09R00 also exhibits a three-fold field symmetry, major radius of 1.44 m, minor radius of 0.32 m, and volume-averaged plasma beta $\langle \beta\rangle = 4.09\%$. However, for direct comparison with the FAMUS solution in [17] and the relax-and-split solution in [21], the C09R00 shape is used but with no plasma current, i.e. $\langle \beta\rangle = 0$, and the toroidal field was taken to be perfectly toroidal with no ripple. The coordinate system is cylindrical and the same grid is used for all the algorithms. The grid has a resolution of 14 points radially and 57 344 locations ($\times 2N_{fp}$ for all N_fp field periods), so $\boldsymbol{A} \in \mathbb{R}^{4096 \times 172032}$. FAMUS is used with a level of regularization for a plausible but nonbinary and non-grid-aligned solution $f_B \approx 1.6 \times 10^{-6}$ T² m² and relax-and-split omits Tikhonov regularization with a similar $f_B \approx 1.6 \times 10^{-6}$ T² m² from a weakly nonbinary and cylindrical-coordinate-aligned solution. The relax-and-split algorithm produces both a fully binary and weakly nonbinary solution (both aligned to cylindrical coordinates), but the fully binary solution achieves only $f_B \approx 5 \times 10^{-4}$ T² m². In order to compare the nonbinary solutions with the a-priori binary relax-and-split and greedy solutions, we define the effective volume of the permanent magnet region,

$$\begin{align} V_\text{eff} = \frac{1}{M_0}\sum_{j=1}^D \|\boldsymbol{m}_j\|_2, \end{align} \tag{ 11 }$$

where M₀ is the maximum magnetization of the magnets. FAMUS and relax-and-split both produce $V_\text{eff} \approx 2.3$ m³. To measure how binary a solution is, we define the binary fraction

$$\begin{align} f_\delta = 1 - \frac{\#\{i | \delta \leqslant \|\boldsymbol{m}_i\|_2 \leqslant 1-\delta\}}{D}. \end{align} \tag{ 12 }$$

FAMUS exhibits $f_{0.01} = 0.57$ and the nonbinary relax-and-split solution exhibits $f_{0.01} = 0.84$.

Although the NCSX grid contains magnets of differing sizes, for simplicity Tikhonov regularization was omitted for the GPMO and GPMOb results presented in this section. All of the algorithm solutions are plotted together in figure 4. Even with many binary magnets, GPMO is not able to get below $f_B \sim 10^{-4}$ T² m², which is not sufficient for an accurate solution. Breaking the binary constraint, as in FAMUS or relax-and-split, allows for significantly improved f_B .

Zoom In Zoom Out Reset image size

However, figure 4 additionally illustrates that the GPMOb results are able to obtain a much improved $f_B \sim 10^{-5}$ using $N_\text{adjacent} = 100$, albeit at somewhat larger values of $V_\text{eff}\approx 3.27$ m⁻³. This is the best binary solution for this configuration by more than an order of magnitude in f_B , and it easily meets the engineering tolerance for NCSX, defined in the PM4STELL collaboration [46, 47] as $b_n \lt 0.5\%$ (their definition of b_n differs slightly with equation (3) but we found empirically that both definitions produce roughly the same value). Below this threshold the plasma boundary and other plasma properties agreed well with the reference equilibrium. The work by Hammond et al [26] performs permanent magnet optimization for a similar NCSX example, using finite β and a set of rectangular permanent magnets, and is able to achieve marginally improved performance for f_B . However, it is hard to conclude about the relative merits of the algorithms since their method achieves improved performance but uses more directional degrees of freedom, is not fully binary, and requires a much larger $V_\text{eff} = 5.74$ m⁻³. Moreover, the volume may be larger primarily because of finite β and the differences in the magnet geometry and orientations.

It is possible that FAMUS and relax-and-split can generate similarly high-performance binary solutions at these larger values of V_eff, but it was the speed and greedy design of GPMOb that facilitated exploring this higher-volume, less sparse part of the parameter space anyways. The illustration of the components of the m solution reveals geometrically interesting magnet structures generated by the backtracking.

Despite the improvements with GPMOb, the single discovered high-quality solution is not particularly sparse. This example provides some evidence that the advantageous engineering properties of all-maximum-magnitude dipoles may not be not worth a likely reduction in f_B accuracy, reduction in sparsity, or increase in magnet volume. Purchasing a small number of weak permanent magnets along with a large number of binary permanent magnets may be a useful way to reduce f_B significantly while still facilitating a magnet configuration that is low-cost and low-complexity.

An explanation is warranted for GPMO's strong performance. Part of the explanation derives from the strong ill-posedness of the problem, which itself derives from the grid satisfying $D \gg 1$. Unlike in many other optimization problems, there are no optimization variables such that, once a variable is determined, the solution is permanently and unrecoverably sub-standard. This capability for recovery comes from several sources: each dipole independently (no coupling between magnets) contributes a very small amount to minimizing f_B , there are no unique and irreplaceable dipole locations (this may not be the case for a permanent magnet grid where some permanent magnets are much larger or much smaller than others, or for a grid with complex or sharp geometrical features), and there are a very large number of degrees of freedom available to either drown out the contribution of a poorly-placed magnet or place an equal and oppositely-oriented magnet directly next to it. So the greedy solutions may use more magnets than necessary, but do not seem to suffer, as in many other scientific fields, from getting inexorably trapped in bad local minima.

Indeed, similar explanations hold for the greater performance of the GPMOb algorithm. Physically, the primary way that solutions can be non-optimal is that a dipole magnetic field is not being fully utilized for f_B minimization because of the presence of a neighbor dipole field that is oppositely oriented. Some amount of oppositely oriented magnets is clearly a benefit for generating complicated magnetic fields on the plasma surface, but systematically pruning close, equal and opposite neighbors seems to work very well for improving the greedy solutions. Moreover, this backtracking is built into the greedy algorithm, i.e. there is no multi-stage process to perform, and adds minor computational cost to the calculation of the solution.

A secondary route for non-optimality appears to be partial cancellations occurring between magnets of different coordinate orientations, which can occur because the magnetic field contribution of a dipole is proportional to a term like $(\boldsymbol{m}_i \cdot \boldsymbol{r}_i)\boldsymbol{r}_i$, where r _i is the spatial location of the ith dipole. The GPMOb algorithm does not address this type of non-optimality, and we speculate here that this may be one reason why the r-direction only algorithm still outperforms the full $(r, \phi, \theta)$ GPMOb algorithm around $K \sim 10\,000$. Future work with more advanced backtracking could address this secondary avenue for non-optimality.