Constrained stellarator coil curvature optimization with FOCUS

Finding less complicated coils that have adequately low field errors is a crucial step in stellarator development. One coil metric that is of high importance is the maximum curvature of the coil centreline, or coil single filament. Conductors cannot be bent below some threshold minimum radius of curvature. High coil curvatures can cause strains to exceed acceptable levels, especially in superconducting coils. We investigate three ways to optimize coil curvature and find that applying penalty functions to the coil curvature solves for coils that have a constrained maximum curvature and low field error. Penalty functions are implemented in FOCUS and coil solutions optimized for an HSX-like ‘plasma boundary’ are presented.


Introduction
Stellarators possess inherent advantages that are important to consider for future fusion power plant designs, including steady-state operation, resiliency to disruptions and theoretically high beta stability limits (Gates et al. 2018). Although stellarators have desirable physics properties, they also require complicated engineering. Generally, stellarators do not require current drive to achieve reasonable rotational transform profiles. To produce non-axisymmetric magnetic fields, non-planar coils are typically used (Merkel 1987;Drevlak 1999;Landreman 2017). Non-planar coils complicate coil design, engineering, assembly and supports. Difficulties arising from non-planar coils led to significant delays during the construction of the W7-X experiment (Riße 2009) and to the cancellation of the NCSX experiment at the Princeton Plasma Physics Laboratory (Neilson et al. 2010).
There is a trade-off between coils that accurately reproduce a desired magnetic field and coils that are simple to manufacture. Stellarator coils need to satisfy engineering constraints such as maximum coil curvature, minimum coil-to-coil separation and minimum coil-to-plasma separation. After such engineering constraints are satisfied one can also optimize for additional considerations such as coil resiliency to displacements, coil spacing for auxiliary systems, coil cost and coil complexity. Coil complexity is left ambiguous here and we do not intend to define a set of quantitative metrics for it. Coil complexity is necessarily dependent on the manufacturing process, whereas the coil maximum curvature is important independent of the manufacturing process. This letter investigates the maximum curvature engineering constraint, just one of many coil considerations.
Previously, both the ONSET and COILOPT codes used curvature objective functions to target coil curvature. The ONSET code uses a penalty function that quadratically penalizes curvatures that exceed some input maximum curvature (see (5) and (6) of Drevlak (1998)). This curvature objective function is not differentiable, which limits the applicability of the objective function to different optimization schemes. Additionally the curvature objective function is dependent on coil lengths. The COILOPT code used an exponentially weighted function to heavily penalize high coil curvatures (see (4) of Strickler, Berry & Hirshman (2002)). One difference between the COILOPT exponential function and our penalty function is that the exponential function penalizes all values of curvature, not just values that are above the maximum. Furthermore, in COILOPT, coils are constrained to lie on a predefined winding surface which greatly limits flexibility.
In this letter we compare three curvature objective functions. In doing so we attempt to elucidate the salient properties of penalty functions that constrain curvatures below some maximum value and objective functions that try to reduce the integrated curvature. While both the ONSET and COILOPT formulations target curvature and are able to reduce coil curvatures, there has not yet been a systematic effort to understand what specific choices of the curvature function yield. This letter provides a first attempt at resolving these questions.
The first two curvature objective functions minimize a coil's average curvature and average curvature squared. The third curvature objective function uses penalty functions to explicitly constrain a coil's maximum curvature. The maximum curvature constraint is of crucial importance since it directly determines whether a coil set can be constructed. For instance, if the minimum radius of curvature is greater than half the small dimension of a coil width, then the coil will self-intersect. Also, smaller radii of curvature directly increase coil strain and can increase the magnetic field strength in the coil, both of which can violate engineering constraints for superconducting coils.
To optimize for a coil's maximum curvature, we will use the single-filament coil optimization code FOCUS (Zhu et al. 2017). Single-filament coils are initial approximations where the coil cross-section is taken to be a point. This coil approximation is appropriate for any type of coil such as low-temperature superconducting, high-temperature superconducting and copper coils. Depending on the size of the coil cross-sections relative to the coil-plasma distances, a multi-filament coil optimization can be done as a subsequent step (Singh et al. 2020). For this letter it is assumed that coil means single-filament coil.
Our motivation and method for constraining maximum coil curvature are described in the rest of this letter. Newly implemented FOCUS objective functions are described in § 2. Coils optimized with these objective functions are shown in § 3. The objective functions and resulting coils are discussed in § 4. Conclusions and future work are discussed in § 5.

FOCUS objective functions
FOCUS is a stellarator coil optimization code that takes advantage of parallel computing. FOCUS solves nonlinear coil optimization problems of the following form: where w i is a user-specified weighting associated with the objective function f i . These objective functions are dependent on and minimized by changing coil geometries and currents. Objective functions in FOCUS are twice differentiable. These objective functions are minimized with gradient-based optimization algorithms, like the conjugate gradient method, or with Hessian-based algorithms, like the Newton method (Zhu et al. 2018). Coils in FOCUS are approximated as space curves and defined using truncated Fourier series as where t ∈ [0, 2π] is our parameterizing variable and N F defines the highest Fourier mode frequency to be used in the optimization. The optimization variables are the Fourier mode amplitudes X c,n , X s,n , Y c,n , Y s,n , Z c,n and Z s,n and coil current for each coil. The primary objective of coil optimization is to satisfy a normal field boundary condition on an input surface. This surface can be arbitrary, but typically is a 'plasma boundary' from a fixed boundary equilibrium. The normal field boundary condition can also be arbitrary, but is typically set by plasma currents and any fixed external current sources. For this letter we use a quadratic, normalized flux objective function: where B v is the magnetic field from coils, n is the input surface unit normal and T Bn is the normal field boundary condition. This objective function will be referred to as a field error objective function. We optimize coils to reconstruct magnetic fields from a zero-beta, fixed boundary equilibrium. This allows us to use the plasma boundary from the equilibrium as our input surface and allows us to set the normal field boundary condition to zero. If the Neumann boundary condition is set to zero, coil currents are not held fixed and the quadratic flux objective function is not normalized by field strength, an additional objective function that optimizes for toroidal flux needs to be included. If it is not included in the optimization, a trivial solution will be solved for where coil currents go to zero, causing the normal field on the boundary to consequently also be zero. For this letter we keep coil currents held fixed, and use the normalized flux objective function (2.3); therefore a toroidal flux objective function is not needed.
'Coil ripple' errors are caused because field strength is large at surface positions close to coils and smaller in between coils. These errors are associated with high toroidal mode numbers corresponding to the number of coils. Coil ripple can be reduced by moving coils farther away from the boundary. Coil ripple errors are visible as varying stripes of minimum and maximum error on the boundary, as in figure 1(a). At larger coil distances it is more difficult to match the general plasma boundary shape even if the field is generated by an external current-carrying surface (Landreman & Boozer 2016). We refer to errors of this type as 'non-ripple errors'.
The inverse of the Biot-Savart law is ill-posed and additional regularization is required in coil optimizations. FOCUS uses coil length as one of the primary regularizations. If there is no penalty on coil length, coils can become arbitrarily long to reduce coil ripple.
plotted on a half-period of the HSX boundary from coils optimized with our penalty curvature objective function (a) and from coils optimized with a quadratic curvature objective function (b).
To obtain the results in this letter we use a quadratic penalty on coil length given as where N c is the total number of coils, L i is the length of the ith coil and L i,0 is the user-specified target length of the ith coil. Including a curvature objective function further regularizes the optimization, but it is still important to include a coil length objective function to ensure that coil solutions are reasonable. Coils optimized with only field error and length objective functions are discussed in § 3. These coils have a maximum curvature that is roughly an order of magnitude larger than that of the Helically Symmetric Experiment (HSX) coils, and therefore motivate the need for a curvature objective function. We now present three candidate curvature objective functions. Including a curvature objective function gives an optimization problem of the following form: min where j = 1, 2, 3 correspond to three possible forms of the curvature objective function. The first curvature objective function (j = 1) is given as with each coil's curvature, κ i , appearing linearly. Alternatively, it may be desirable to penalize areas of high curvature more than areas of low curvature. A second curvature objective function that optimizes for the square of the curvature should accomplish this. This function (j = 2) is given as It is important that we integrate over the coil's arc length |r i | dt and not the parameterizing variable dt. If we were to integrate over the parameterizing variable, the objective function It is also important that we normalize each integral by the corresponding coil length. If we do not do this, coil lengths can decrease to minimize the objective functions despite having regions of higher curvature. The curvature objective functions (2.6) and (2.8) attempt to minimize an integrated function of coil curvature, but they do not constrain it. Therefore we are solving an unconstrained optimization problem. In order to constrain the coil curvature, the third curvature objective function (j = 3) uses a penalty function, making this a constrained optimization problem. We apply a constraint on coil curvature so that it is below some maximum value κ 0 at all points in the set of coils. That is, κ i (t) ≤ κ 0 = 1/R 0 ∀t and ∀i, where R 0 corresponds to the minimum radius of curvature. To apply this constraint we can use the following third curvature objective function which uses a penalty function: (2.10) The third curvature objective function includes a Heaviside function defined in (2.10) which is discontinuous. It is important to notice that at the Heaviside function discontinuity, the term cosh(α(κ i − κ 0 )) − 1 goes to zero which can be seen in figure 2. This preserves the objective function's continuity and moreover its differentiability. This objective function is twice differentiable, f κ,3 ∈ C 2 . FOCUS has multiple optimization algorithms implemented. If the user were to choose one of the quasi-Newton methods, superlinear convergence can be observed since the FOCUS objective functions are twice differentiable (Nocedal & Wright 2006).
The penalty function only penalizes coil locations where coil curvature exceeds the specified maximum curvature, κ 0 . The value of H κ 0 (κ i )(cosh(α(κ i − κ 0 )) − 1) 2 is zero if κ i ≤ κ 0 , positive if κ i > κ 0 and increases exponentially with the constraint violation κ i − κ 0 . To guarantee a constraint on coil curvature, the value of α needs to be iteratively increased to ∞. Here α is called our penalty variable. As α → ∞, the value of f κ,3 → ∞, if any coil has a location where the coil's curvature exceeds the maximum curvature.
This forces all local minima to have coil geometries with curvature less than or equal to the input maximum curvature. To obtain coils using f κ,3 for this letter, the initial value of α is set to 0 and the optimization is run until a local minimum is found. Coils from the initial optimization are then input back into FOCUS with α increased to 0.1 and the optimization is again run until a local minimum is found. This process is repeated with α increasing after each optimization. Coils with a constrained maximum curvature are found after α is increased to 4.

Applications
FOCUS is used to solve for coils, targeting an idealized boundary similar to that of the HSX that has been smoothed to remove errors produced by the actual HSX coils. We refer to this smoothed boundary as the HSX boundary. HSX is a four-field-period device with 48 coils. Because of symmetries, there are only six unique coil shapes. The normal field boundary condition is set to zero, T Bn = 0, representing a zero beta equilibrium with no fixed external magnets. We initialize the optimizations with filamentary coils that are roughly equal to the centroids of the HSX finite build coil cross-sections. For the remainder of this letter, the term 'HSX coils' refers to these filamentary coils.
We present four optimizations to compare the different objective functions: one with no curvature optimization and one for each of the three curvature objective functions. Coils in these optimizations have equal currents of 151 kA and currents are held fixed during optimization. We do not change each coil's target length between optimizations, we use the nonlinear conjugate gradient algorithm and strong Wolfe conditions (Nocedal & Wright 2006) and we set Fourier series truncation to be N F = 16 for the coils in each optimization. Target coil lengths are chosen to be slightly larger than the HSX coil lengths. The four optimizations are run until local minima are found.
Choice of objective function weightings will significantly change the solution space and resulting coil solutions. Objective function weightings are given in table 1. The value of coil length objective function weightings is chosen such that the difference between the resulting coil lengths and target coil lengths is of the order of centimetres for each of the four optimization runs. Weighting for the linear curvature objective function was adjusted until the field error objective function value f Bn for the resulting coils was close to the f Bn value for the HSX coils. Weighting for the quadratic curvature objective function is set in the same way. It should be noted that the penalty function variable α is more important than the penalty curvature objective function weighting.
To evaluate each optimization result, we compare against the HSX coils. The HSX coils and the four FOCUS-optimized coil sets can be seen in figure 3. A comparison of the coil curvatures can be seen in figure 4. Important metrics to consider for each coil set are the field error objective function, maximum field error, maximum curvature |B v | from coils with no curvature optimization is plotted on a half-period of the HSX boundary in (a) and |B v | from coils optimized with our penalty curvature objective function is plotted on a half-period of the HSX boundary in (b). High |B v | is plotted in black and low |B v | is plotted in yellow.
and minimum coil-to-coil separation. These values are given in table 2. For the penalty curvature objective function, the value of κ 0 is set equal to the HSX coil maximum curvature and we can see that our penalty curvature objective function has been used to successfully constrain the coil curvature.

No curvature optimization
If N F = 16 and coil curvature is not optimized for, it is evident that FOCUS does not solve for reasonable coils due to high coil curvatures, with radius of curvature in places reduced below 1 cm. One technique to reduce coil curvatures is to decrease the value of N F . Even for low values of N F , coil curvature can still become too high. Using a low value of N F restricts coil shaping and leads to higher field errors.

Linear curvature optimization
Including the linear curvature objective function produces coils with low average curvature. Unfortunately, the maximum curvature of these coils is 74.4 m −1 which is significantly larger than that of the HSX coils that have a maximum curvature of 12.3 m −1 , and both coil sets have roughly equal f Bn values. After close inspection it can be concluded that some regions of the coils have been straightened out while other regions of the coils still have sharp bends. High curvatures are not heavily penalized as long as they exist over short arc lengths.

Quadratic curvature optimization
The quadratic curvature objective function should penalize regions with high curvatures more than does the linear curvature function. Coils generated with the quadratic curvature function have a maximum curvature of 7.62 m −1 which is lower than the HSX coil's maximum curvature of 12.3 m −1 . Again both coil sets have roughly equal values of f Bn . It can be seen from figure 1(b) that the field error distribution does not appear to be dominated by coil ripple. A large portion of the error comes from non-ripple error. For  roughly equal values of f Bn , it appears that the quadratic curvature objective function outperforms the linear curvature objective function with regards to the maximum curvature. Coils optimized with the linear curvature objective function have a maximum curvature that is ≈10 times higher than coils optimized with the quadratic curvature objective function.

Penalty function curvature optimization
The penalty function is designed to only penalize positions on the coils where the curvature exceeds a prescribed maximum, here set to be the HSX coil curvature maximum, κ 0 = 12.3. Coils solved for with the penalty function are shown to have a maximum curvature equal to that of the HSX coils. The value of f Bn is more than an order of magnitude lower than that of the HSX coils. Remarkably, the value of f Bn for these coils is close to the value of f Bn from the FOCUS coils that are not optimized for curvature. Since we do not penalize locations of coils that have curvatures below our specified maximum curvature, the coils have more flexibility to accurately reconstruct the plasma boundary. It can be seen from figure 1(a) that the field error distribution appears to be dominated by coil ripple and the normalized ripple error is roughly 0.25 %. This means that the normalized non-ripple error is even lower than 0.25 % and because of this we have higher confidence that the general shape of the 'plasma boundary' is reconstructed well.

Conclusions
Without a curvature objective function, optimizations in FOCUS can quickly solve for coils that have unreasonably high curvatures. There are techniques that can be used to limit coil curvatures, but typically at the cost of increasing field errors. Because of this, multiple curvature objective functions have been implemented in FOCUS. We see that a quadratic curvature objective function seems to outperform a linear curvature objective function with respect to maximum coil curvature, but does not outperform with respect to maximum field error. Neither result is satisfactory. This motivated a constraint on maximum curvature as opposed to optimizing for curvature at every location of the coil. To apply this constraint we used penalty functions. The resulting coils with the penalty functions outperformed both the linear and quadratic objective functions with regard to matching the target boundary while still remaining below the maximum curvature threshold.
If coil buildability is our primary concern, then a coil's maximum curvature seems to be a more apt metric than the average coil curvature or average squared curvature. We see that constraining a coil's maximum curvature with a penalty function can solve for coils with a very low value of the field error objective function, f Bn . A hyperbolic cosine penalty function was chosen for the optimization mainly due to its exponential penalty and its differentiability.
The penalty objective function achieves the desired behaviour of only penalizing locations of the coils where curvatures exceed a maximum curvature. However, we can see in figure 3(b) that there are small wiggles in the coils produced with the penalty objective function. These deviations may decrease field errors by a small amount, but more importantly, they may also significantly increase 'coil complexity'. Removing these additional wiggles is a topic that can be addressed with additional curvature objective functions that specifically target coil complexity. This is left for a future publication.
In this letter we demonstrate how to use penalty functions to constrain coil curvature. We intend to also use penalty functions to constrain the minimum coil-to-coil separation. Many more useful applications of penalty functions exist outside of coil optimization. For instance, when optimizing for stellarator equilibria, it is often desirable to constrain rotational transform profiles so as to not cross low-order rational values. Penalty functions are easy to implement and do not require additional optimization algorithms for their implementation.
First derivatives of the three curvature objective functions are now given. Here X i is an arbitrary Fourier mode amplitude, X i ∈ {X c,0,i , . . . , X c,N F ,i , X s,1,i , . . . , X s,N F ,i }, that defines the x position of the ith coil. Derivatives for Y i or Z i Fourier mode amplitudes take a similar form. The equation for ∂L i /∂X i is given in (15) in Zhu et al. (2017). The derivatives are