Optimization of modular and helical coils applying genetic algorithm and fully-three-dimensional B-spline curves

A new numerical method for designing the external coils of a stellarator is presented. In this method, the shape of filamentary coils is expressed using fully three-dimensional B-spline curves that are not necessarily constrained on a winding surface. The control points of B-spline curves are optimized together with the coil position and current to minimize an objective function, which is defined using normal field components and engineering constraints. The genetic algorithm is employed to minimize the objective function for arbitrary combinations of modular, helical, and circular poloidal field coils without giving any specific initial guess of coil shapes. A new numerical code genetic optimizer using sequence of points for external coil is developed on the basis of this method, and successfully found optimized modular coils for the stellarators CFQS and Wendelstein 7-X. We also found a specific pattern of helical coil arrangement that can reproduce these optimized stellarators while creating divertor legs outside of the closed magnetic surfaces.


Introduction
Designing the external coils that produce the optimized threedimensional magnetic field with minimized error field while satisfying the engineering feasibility is one of the key steps of the development of stellarator-type fusion reactors. As the * Author to whom any correspondence should be addressed. a See Klinger et al 2019 (https://doi.org/10.1088/1741-4326/ab03a7) for the W7-X Team.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. numerical optimizers [1,2] of the magnetic configuration have been established, numerical schemes and codes to seek for or to optimize the external coils have also been developed.
Most of these codes assume the existence of a toroidal winding surface on which external coils are wound. NESCOIL [3] and REGCOIL [4] optimize a current potential Φ that produces a divergence-free surface current on a winding surface. A set of discrete coils is obtained as a set of contours of Φ expressed by a finite number of Fourier modes. Nonlinear optimization codes that can take into account engineering constraints on coils have also been developed, such as ONSET [5], COILOPT [6], and COILOPT++ [7]. These codes, sometimes utilizing the result of the previous two codes as an initial guess, optimize more directly the coil shape on a two-dimensional toroidal winding surface, together with the shape of the winding surface itself. The practical applications of these codes so far are mostly related to optimization of modular coils, which are eventually represented as current circuits on a single, toroidal winding surface.
On the other hand, using other types of coils as primary magnetic field coils is an option. As represented by some of the early stellarators, heliacs, torsatrons, and heliotrons, the non-axisymmetric toroidal magnetic field is not always realized only by coils on a single winding surface. For example, the heliotron configuration of the large helical device [8] (LHD) uses a pair of helical coils and three pairs of poloidal field (PF) coils, which cannot be expressed by a single current surface. Use of helical coils has several inherent advantages such as the presence of robust divertor legs [9] and easier access to the plasma or inside the vacuum vessel. The codes assuming a single winding surface may be used to design auxiliary coils for a known set of the primary coils in an iterative way. However, they cannot treat such different types of coils simultaneously as primary coils.
There is also another approach to coil designing in which a winding surface is not assumed, which has been implemented in ONSET and the FOCUS code [10,11]. In this approach, the coils are treated as spatial curves that can move and shape in the three-dimensional real space. In the reference [11], the shape of the coils was expressed using Fourier series and numerically optimized using modified Newton method. This approach can optimize a coil configuration composed of different types of coils. An initial guess of coil configuration, which may not be too far from the unknown 'optimized' one, may be required to start a calculation. This has two drawbacks. One is that there is arbitrariness in the arrangement and the shapes of initial coils we give to a completely new configuration whose external coils have never been given. Second, the optimized coils that can be obtained may depend on a given initial coil configuration and may not always be the optimal solution in the global sense.
One of the optimization approaches that may potentially overcome these drawbacks is the genetic algorithm (GA). The GA has been implemented in other coil optimization code and STELLOPT [12,13]. In the GA, a set of trial solutions called 'individuals' constitutes a generation. In the selection process of the GA, individuals that are well adapted to the 'environment' are preferentially selected, and their traits are inherited by the individuals of the next generation. In the context of the coil optimization, the degree of adaptation to the environment may be interpreted as the smallness of the error field components, as well as how well the engineering constraints are met. It should be noted that the GA usually requires more iterations than linear problem and gradient-based nonlinear methods. It is expected that, however, if we introduce a random set of coil configurations as the individuals of the first generation, for a given combination of different coil types, we may obtain optimum shape and arrangement of coils automatically by applying the GA.
In this paper, we apply the GA to perform optimizations of external coils without assuming the existence of a winding surface, nor giving any specific initial guess for coil shape from other codes. To express the coil shape in three-dimensional real space, we use an interpolation of a sequence of spatial points rather than a Fourier expansion. Some of the existing codes, such as ONSET (cubic spline) and COILOPT++ (cubic B-spline) can use this type of coil expression in the two-dimensional space of a winding surface. In this study, we employ a fully three-dimensional cubic B-spline curve. A B-spline curve [14] is drawn as a smooth interpolation of a sequence of spatial points using basis functions. The optimization of a coil configuration is formulated as a combinatorial optimization problem of the real space positions of a sequence of points together with the coil currents. Our goal is to demonstrate that the combination of these two methods is practically useful in designing external coils, not only modular coils but also helical coils.
In section 2, the methodology of the coil optimization using the B-spline curves and GA will be explained, including the basis of B-spline curve and the details of the objective functions used in the optimization. The explained methodology is implemented in the newly developed numerical code, GOSPEL genetic optimizer using sequence of points for external coil (GOSPEL). In section 3, GOSPEL is applied to two different types of optimized stellarators, namely, CFQS [15,16] and Wendelstein 7-X [17,18], which are significantly different in the symmetry of the magnetic field, aspect ratio, and toroidal period, etc. Optimization of both modular and helical coils is performed targeting these magnetic configurations without giving any initial guesses from other codes. The dependency of convergence on the parameter specific to the GA is also investigated. The conclusions of this study are given in section 4.

B-spline expression of coil shape
We approximate the shape of a coil as a spatial curve that closes on itself. A position on the curve is expressed using a parameter 0 t 1 as x(t), where x(0) = x(1).
A (cubic) B-spline spatial curve x(t) is written in terms of a sum of the position vectors of discrete spatial points (control points) weighted with the value of the basis functions. Let us use total N c independent spatial vectors x i (i = 1, 2, . . . , N c ) as the control points. Each of these control points are evenly distributed on the t-axis, as t = t 1 , t 2 , . . . , t N c , where t i = i Δt = i/N c . Using the basis functions b i (t), a closed B-spline curve can be written as follows: where we have added three control points near t = 0 and t = 1 to close the curve, namely, The basis functions are defined as where s i (t) = (t − t i )/Δt + 2. The basis function b i (t) has its peak at t = t i as shown in figure 1, and has been normalized to satisfy for 0 t 1. The derivatives of x are calculated using the derivatives of the basis function. Figure 2 shows a comparison of a closed B-spline curve and an ordinary cubic-spline curve, each defined with 12 control points in a two-dimensional space. Unlike the ordinary spline curve passing through all of the control points, the B-spline curve does not necessarily pass through the control points.
In this study, we assume a stellarator symmetry of the configuration. This can reduce the degree of freedom to 1/(2N t ), or to 1/(4N t ) for a symmetric helical coil, of a fully asymmetric case for a configuration with the number of toroidal period N t . It is useful to separate translational movement and rigid body rotation of whole coil from local modification of coil shape. We introduce the vertical displacement Z n and the toroidal rotation angle φ n of the nth coil as variables to be optimized. These translational movement and rotation are introduced to efficiently explore the possible coil configurations. Actually, these movements do not add new degree of freedom.
For modular coils, the shape of nth coil is expressed using N c control points measured in a local Cartesian coordinates, . . , N c , with its origin set to (R = 0, Z = Z n ) and unit vectorsx =R(φ n ),ŷ =φ(φ n ),ẑ =Ẑ. For helical coils, we define the control points for one toroidal period in cylindrical coordinates, To satisfy the symmetry, two different types of helical coils should be considered: self-symmetric coil and paired-symmetric coil. For the former type of coil, control points should be symmetric and Z n = φ n = 0. The latter type of coil should not always be symmetric itself, but always constitute a symmetric pair. If the nth and (n + 1)th coils constitute a symmetric pair, following relations hold: The shape and position of a PF coil pairs is specified by the radius and the height of the coil, R (n) , Z (n) . Finally, the coil current, I n , completes the set of parameters to be optimized for each of the coils. Figure 3 shows a concept of coil optimization using the GA. In the optimization using the GA, a set of trial coil configurations constitutes a generation. The value of the objective function is evaluated for all of the coil configurations in the  current generation. The coil configurations in the current generation are ranked according to their fitness, or the smallness of the value of the objective function. New coil configurations are generated by crossover of pairs of coil configurations that are randomly selected from the current generation with a probability depending on their rank. Random modification on each of the new coil parameters called mutation is imposed at a prescribed probability. Finally, the coil configurations in the current generation are replaced with the new coil configurations. The specific models used in this study will be explained in the following subsections.

Initial generation.
For modular coils, we generate an initial ensemble of coils by mixing two kinds of initial coils randomly and imposing random modulations to them. The first kind of coils are planar coils that enclose the target surface poloidally. The second kind of coils are poloidally-closed circuits on a toroidal surface obtained by extending the target surface outwards. The latter ones are more 'aligned' to the target surface and closer to typical shape of modular coils.
For helical coils, we also introduce two kinds of initial shape, which are both expressed by a simple formula, where σ R± and σ Z± are uncorrelated random signs thereby randomizing the direction of helical winding. Thus, helical coils whose direction of helical winding is different from that of other helical coils, and/or target surface, can be generated.
The first kind has a fixed major radius R 0 that is equal to that of the target surface. The minor radii, a R and a Z , are randomly modified while their ensemble averages are kept to be slightly larger (∼ 110%) than the averaged minor radii of the target surface. This kind of helical coil usually circulates around the target surface both poloidally and toroidally without significant interference with the target surface. The second kind of coil has random major radius whose ensemble average is the same as the first kind of coil. The minor radii for this kind of coil is also random, but their averages are set to be smaller (∼ 10%) than the minor radii of the target surface. This kind of coil does not always circulate around the magnetic axis, and can behave like a distorted PF coil.
To generate a self-symmetric helical coil, only the first kind of coil is adopted with no vertical shift nor the toroidal rotation (Z n = Z 0 = 0, φ n = δφ = 0). For a symmetric pair, the first and second kinds are randomly adopted at an equal probability, and the vertical shift and toroidal rotation angle are randomly given. Control points with respect to toroidal angles of the initial generation are not randomized, and given by φ (n) For the PF coils, the radius and height are randomly given. For each of the coils, coil current I n are given by imposing random modulation to a reference value of current. Because our objective function explained in the next subsection does not depend on the magnetic strength but the normal vector on the target surface, the choice of this reference value does not affect the calculation. For modular coils, all coil currents are set in the same poloidal direction corresponding to a clockwise toroidal magnetic field. For helical and PF coils, the direction of coil currents are randomly given. To avoid divergence of coil currents and to perform the appropriate mutation of the coil currents, we normalize all the coil currents for each of the coil configurations at every iteration in optimization loop.

Fitness and selection of parents.
We define the fitness of a coil configuration as the smallness of the objective function f = 6 k=1 f k , where f k is the component for different targets. The first component is the significance of field component  normal to the target magnetic surface, where M is the total number of reference points distributed on the target magnetic surface, b j = B j /|B j | is the unit vector in the direction of the magnetic field generated by the coils and n j is a normal vector of the target magnetic surface, respectively. The subscript j denotes these quantities are evaluated at jth reference point. Second component is introduced to restrict the minimum value of local curvature radius, and is defined as where C 2 is a constant parameter, min(r c ) is the minimum value of local curvature radius of the whole coil configuration, and r cmin is the acceptable minimum value of local curvature radius. As long as a coil configuration satisfies min(r c ) > r cmin , coil curvature will not affect the total fitness. The third component corresponds to a target on local torsion, which is defined similarly, where max(τ) is the maximum value of local torsion of a whole coil configuration and τ max is the acceptable maximum value.
Although we include this component in the objective function, putting a limitation on local torsion can be a too strong limitation on the possible coil shape since we assume a single filament. Thus, in this study, we set C 3 = 0. The fourth component is the target on the length of each of the coils, where N coil is the number of independent coils, l (n) c and l (n) cmax are the total length and the acceptable maximum value for nth coil, respectively. The fifth and sixth targets are related to relative distances, and are defined as where min(Δ c-c ) is the minimum value of the coil-to-coil distance, min(Δ c-s ) is the minimum value of the distance between the coil and the target magnetic surface. Δ c−c min and Δ c−s min are the acceptable minimum values. The value of constant factors are set to C 1 = 1 and C 2-6 = 10 in this study.
Each of the coil configurations in the same generation is ranked with respect to the fitness from the best configuration with N rank = 1 to the worst configuration with N rank = N. Here N is the total number of coil configurations included in each of the generations. The coil configuration with the highest fitness or N rank = 1 in the present generation is kept as an 'elite' in the next generation without any changes to the coil parameters as shown in the conceptual diagram of figure 3. This promotes a monotonic decrease of f .
Excluding some of the coil configurations with small fitness from the selection of parents can accelerate the convergence of the minimum value of f . On the other hand, such exclusion generally deteriorates the genetic diversity and may increase the probability of converging to a local minimum at an early phase of optimization. We left the number of configurations to be excluded from the reproduction process, N ex , as one of the free parameters. Dependence of the convergence on N ex will be checked in the next section.
A pair of coil configurations are randomly selected from N − N ex coil configurations as parents with the following (relative) probability: where μ is a parameter controlling dependence of a chance for a coil configuration to create its offspring on N rank . μ = 0 corresponds to an equal probability for all N − N ex configurations. By using larger μ, the probability becomes higher for configurations with higher fitness and convergence is accelerated. In this study, we use the value μ = 0.44. Each of the selected pairs produce a pair of coil configurations according to the crossover law explained in the next subsection. Selection of parents continues until the next generation is complete. We allow the same pairs of parents to appear multiple times in the selection loop, while every coil configuration is prohibited from making a pair with itself because our crossover law only produces a copy of the parent in that case.

Crossover and mutation.
Blend crossover [19] (BLX-α) is employed as the crossover method. This crossover method generates new parameters for child, h ch , from the parameters of the parents, h p1 , h p2 , according to the following probabilistic operation: where σ is a random number ranging from 0 to 1. α 0 is a parameter extending the range of possible value of h ch . The corresponding coordinates of control points, vertical shift, toroidal rotation angle, coil currents, (and the radius and height for the PF coils) of the new coil configuration are generated using this crossover operation from the parental coil configurations. In this study we use a value α = 0.25. Finally, mutation is imposed on the generated new coil configurations. Mutation is a probabilistic process. There is an arbitrariness in what mutation laws are employed, while the choice of mutation law can affect the convergence of f . The following mutation laws are chosen after trial and error. Systematic studies on the effect of the mutation law on the quality of optimization are left to future works.  The mutation laws we employ for modular and helical coil are summarized in the diagram shown in figure 4. In these laws, mutations in the modular coils occur in only one of the modular coils at once, and mutation of a control point and other parameters (Z n , φ n , I n ) do not occur simultaneously. Thus, the modular coils are not drastically modified. When mutation occurs to a control point, random fluctuations ranging from −20% to 20% of the major radius of the target surface are added to each of the x, y, z coordinates of the control point. Similarly, the range of fluctuation to be added is from −10% to 10% of the major radius of the target surface for Z n and from −5% to 5% of the interval between the coils of the evenly-spaced case for φ n . The coil current I n is modified by multiplying a random factor in the range of 0.9 to 1.1 to the original value.
For helical coils, each of the independent coils experience mutation at a probability of 1/2, and the mutation of the control points, Z n , φ n and I n occurs independently. The range of mutations are from −20% to 20% of the minor radius of the target surface for the R and Z coordinates of a control point and from −10% to 10% of the interval between control points of the evenly-spaced case for the φ coordinate of a control point. The range of mutation of Z n is from −40% to 40% of the minor radius of the target surface. For the PF coils, mutation of the radius, height and coil current occur each at an equal probability of 1/4. The range of mutation is from −20% to 20% of the minor radius of the target surface for the radius and height. In the trial and error process, it is found that the radius of a PF coil tends to become smaller in the subsequent generations once it has become smaller than a helical coil. To avoid this, the above mutation probability for the radius (1/4) includes the 1/20 probability of more drastic mutations ranging up to twice the major radius of the target surface. A similar drastic mutation is also applied for the height.
The range of mutation for I n of the PF and helical coils are defined using a reference value, that is the maximum value of I n of the PF and helical coils of each of the coil configurations. The range is from −10% to 10% of this reference value for the PF coils and from −50% to 50% for the helical coils. In addition, when mutations of the coil currents occur to the helical and PF coils, the direction of the coil current is also reversed at a probability of 1/2 (helical) and 2/5 (PF), respectively.

Numerical code: GOSPEL
A new numerical code for stellarator coil optimization, GOSPEL, is newly developed on the basis of the method described above. The code is written in FORTRAN and parallelized using OpenMP. Some of the notable features of GOSPEL are as follows: • Arbitrary number of control points/coil can be used.
• Can be run without giving initial guess for coil shape.
• Restart from previous run with larger/smaller number of control points and different constraint parameters are possible. • Both free curve and curve constrained on a winding surface (not so strict, as explained above) can be assumed.
In the following sections, we will utilize the restart function for efficiently finding a smooth coil configuration with smaller error field by gradually increasing a number of control points from a previous run.

Modular coil optimization
In this section, we apply the developed code, GOSPEL, to two different types of optimized stellarators, namely, CFQS and Wendelstein 7-X. CFQS is a low-aspect ratio (∼4.3) quasiaxisymmetric stellarator with toroidal period 2. Wendelstein 7-X is a quasi-omnigeneous optimized stellarator with relatively large aspect ratio (∼10) and toroidal period 5.

Quasi-axisymmetric configuration of CFQS
To check the practicality of the developed code, we firstly apply GOSPEL to CFQS [15,16]. The poloidal cross sections and the boundary magnetic surface of this configuration are shown in figure 5.
We assume the same number of the coils as the original design [16] and optimize the coil configuration with and without the engineering constraints. First, we perform the calculation without engineering constraints. The number of coil configurations per generation is set to N = 480. The number of coil configurations excluded from the selection process is set to N ex = 80. Figure 6 shows the coil configurations with the highest fitness of several generations for N c = 8. We can see the evolution of coil configuration converging to a typical stellarator-type modular coils as the generation proceeds. The number of control points per coil is gradually increased from N c = 8 to 14. For the best configuration at the 1200th generation with N c = 14, the normal-field component reaches f 1 = 5.85 × 10 −3 . The minimum curvature radius and the coil-tocoil separation distance are 12.4 cm and 7.10 cm, respectively. Figure 7 shows the evolution of the minimum value of f for the number of control points per modular coil N c . We can see f almost converges at N c = 14.
Next, we impose engineering constraints, namely, r cmin = 21.5 cm, min(Δ c-c ) = 18.5 cm, which are from the reference [16]. Also we put constraints on the maximum length of each of the independent coils, l (n) cmax , where n = 1-4. The values of l (n) cmax are also set equal to the length of the original design. The coil configuration obtained without engineering constraints was used as the initial data to restart the calculation. The numbers of control points per coil and coil configurations per generation were set to N c = 16 and N/N ex = 1000/200, respectively. The calculation was truncated at the 4000th generation and we successfully obtained an optimized configuration with f 1 = 9.22 × 10 −3 , min(r c ) = 22.1 cm and min(Δ c-c ) = 18.7 cm. The minimum curvature radius, coil length, and coil current of the four independent coils of the optimized configurations are summarized in table 1. The coil#1 is located at the vertically elongated cross section and the coil#4 at the triangular cross section of the magnetic configuration. We can see that the coil configuration with larger curvature radius and shorter coil length was obtained by imposing the engineering constraints. Figure 8 shows the optimized coils obtained with the engineering constraints. To improve the readability of the threedimensional geometry of coils, these coils are drawn as those with square cross-section. The Poincáre plots of the vacuum magnetic field lines generated by these coil configurations are plotted together with the vacuum equilibrium by VMEC [20] in figure 9. We can see good agreements of these Poincáre plots with that by VMEC. Figure 10 compares the rotational transform profiles as a function of R measured at (Z, φ) = (0, 0).  We can also see relatively good agreements for the rotational transform profiles.
The coils obtained by GOSPEL and the original design that has been obtained with NESCOIL are compared in figure 11. The coils by GOSPEL have fewer bends than the original design, while both designs satisfy the engineering constraints on the minimum curvature radius and coil-to-coil distance. The dependence of convergence of f on the number of configurations to be excluded from the reproduction process, N ex , is investigated for the CFQS configuration. Figure 12 shows the evolution of the minimum value of f of each of the generations for different values of N ex . These calculations were done with N = 600 and N c = 8. The engineering constraints are imposed. The initial value of the minimum f is large (∼ 2.5) because of the engineering constraints. The minimum value of f rapidly decreases (2.5 → 0.25), indicating that a coil configuration that satisfies the engineering constraints is produced before the 10th generation for all cases. The acceleration of the convergence by using larger value of N ex can be seen  both before and after 10th generation. The minimum value of f 1 at the 100th generation is 0.08 for N ex = 0 and 0.02 for N ex = 500. This acceleration is almost canceled by the slowing down occurring after 100th generation, which may be due to the loss of the genetic diversity. The minimum values of f 1 at the 1000th generation reach almost the same level (∼ 0.015) for all cases.

Wendelstein 7-X
Next, we apply GOSPEL to Wendelstein 7-X. A boundary obtained from a magnetic surface in the vacuum magnetic field of the so-called standard-iota high-mirror configuration is used as a target. This configuration is characterized by a relative variation of the magnetic field of approximately 10% on the magnetic axis and the ι = 1 islands being used for divertor operation, where ι is the rotational transform.
The poloidal cross sections and the boundary magnetic surface are shown in figure 13. We can see the toroidal ripples of the field strength on the magnetic surface due to the discrete modular coils.
In the calculation, the number of modular coils is set to 50, which is the same as the actual device. The number of control points per coil is gradually increased from N c = 7 to 11. The number of individuals used are N/N ex = 960/160. Figure 14 shows the coil configuration obtained using N c = 11 without engineering constraints. The coils are drawn as those with square cross-section just to improve the readability of three-dimensional geometry of coils. The minimum distance between the coils is 22.7 cm and the minimum curvature radius is 15.6 cm, respectively. The minimum distance between coil and the target surface is 26.4 cm. Although we did not assume winding surface nor did we impose engineering constraints, interference between coils was not found. The normal field component is f 1 = 0.67 × 10 −2 . The Poincáre plot of the magnetic field generated by these coils are shown in figure 15 together with the equilibrium by VMEC. Although the tip at the outer side of the torus is not as sharp as the target surface at φ = π/5, the Poincáre plot is well reproduced by the obtained coil. The rotational transform profiles shown in figure 16 are also in a good agreement with that by VMEC.
These results show that the optimization scheme using three-dimensional B-spline curve and GA can properly work as a method of designing the external coils for stellarators.

Helical coil optimization
In this section, we explore helical coils that can reproduce the magnetic surface of the optimized stellarators, which are usually realized by modular coils.
Before targeting at the optimized stellarators, we firstly apply GOSPEL to a target magnetic surface that is produced by a known pattern of helical coils. We assume the coil configuration consisting of a pair of symmetric helical coils and two pairs of PF coils, which is similar to LHD without its inner-shaping PF coil pair. The major and minor radii of the helical coils are 3.9 m and 0.94 m, respectively. The major radius and height of the two PF coil pairs are (R, Z) =(5.5 m, 1.55 m) and (1.8 m, 0.8 m), respectively. The target magnetic surface was obtained by calculating the vacuum magnetic field assuming filamentary coils. We applied GOSPEL to this target surface assuming the same number of PF and helical coils. For the helical coils, paired (nonself-symmetric) helical coils are used. Engineering constraints were not imposed.
Calculation was done with N = 960, N ex = 160, and N c = 7. It was found that the coil configuration with helical coils having the same helical winding direction as the original coils were automatically selected as the best configuration at the early phase of the optimization. The helical coils of the best coil configuration at the 1000th generation and the original ones are compared in figure 17. The error field component of the optimized one is f 1 = 1.5 × 10 −2 . The shape of the original, self-symmetric helical coils are excellently reproduced. After several calculations with different random numbers, it was found that the reproduction of the shape of these helical coils is robust. On the other hand, the radius and the height of PF coils were different from the original coils and showed dependence on the random numbers. One of the reasons might be that these PF coils are located farther from the target surface than the helical coils and there may exist arbitrariness in the position and current of PF coils.  Next we use W7-X as the target. In the optimization of helical coils, the number of each type of coils must be specified. In this study, we assume a pair of PF coils and three helical coils: one self-symmetric helical coil and a pair of helical coils.
We found that different types of arrangement of helical coils show the best fitness in the early phase of the optimization. Typical patterns are summarized in figure 18. Note that PF coils are not drawn in the figure In the pattern A in the figure, the self-symmetric helical coil is located inside the magnetic axis at the vertically-elongated cross section. On the other hand, the self-symmetric coil of the pattern B has an opposite initial phase of helical winding, as it is located outside of the magnetic axis at the vertically-elongated cross section of the target surface. The paired helical coils of these patterns have the same direction of helical winding as the self-symmetric coil, and their initial phase is opposite to the self-symmetric coil. The coil currents in the helical coils of these patterns are comparable to each other and in the same direction. In the pattern C, the paired coils are not helical in the sense that they do not circulate around the magnetic axis. The coil current of these paired coils is an order smaller than that of the symmetric helical coil. In every pattern, the helical winding of the helical coil is in the same direction as the magnetic axis.
The convergence of minimum value of f up to the 500th generation for these patterns for N c = 7 are compared in figure 19. In these calculations, the number of individuals are set to N = 960. The large value for N ex (=720) is used to promote convergence to local minimum and avoid a transition to a different pattern. It is found that the pattern B tends to reproduce the target magnetic surface with the smallest error among these three patterns. It was also found that the convergence is much slower when using helical coils than modular coils if the same number of N ex is used. When we used a small number of N ex , 160, for example, transition between patterns and convergence to local minimum of pattern A or C occurred depending on the choice of the seed of the random numbers. This implies that the selection and/or the cross over method we introduced in section 2 may not fully utilize the inherent advantages of the GA for helical coil optimization.
To obtain an optimized helical coil configuration of the pattern B, we set N ex = 840. Engineering constraints were imposed on the length of helical coils (90 m, 66 m, 66 m), minimum value of R of helical coils (3.0 m), minimum distance between coils (40 cm), minimum distance between coils and target surface (30 cm), and the maximum height of PF coils (3.5 m). The number of control points per helical coil per toroidal period was gradually increased from N c = 7 to 15. After increasing the number of control points per helical coil to N c = 15, we obtained the helical coils of the pattern B with the error field component of f 1 = 1.7 × 10 −2 .
We also performed helical coil optimization for CFQS. The same tendency of coil patterns as in the case of W7-X was also seen. In a similar way, as in the case of the W7-X, we obtained optimized coil configuration of pattern B also for CFQS, with the error field component of f 1 = 1.9 × 10 −2 .
The optimized helical coil configurations for CFQS and W7-X are shown in figure 20 as volume coils with a square cross section, together with the vacuum magnetic surface that   is produced by filament coils. As expected for the helical coils, there is a wide distance between the coils. We can see a similar pattern in the shape of helical coils from the vertical views for these different magnetic configurations. This implies that the pattern B of helical coil configuration may be a generic pattern for the optimized stellarators.  The Poincáre plot of the vacuum magnetic field generated by these coils are shown in figures 21 and 22, together with those by the modular coils obtained in the previous subsections. In the figure, the field lines outside of the last-closed magnetic surface are also shown except for the modular coils of CFQS. These field lines are followed from near the lastclosed magnetic surface until they reach the virtual walls. These virtual walls are defined extending the last-closed magnetic surface not to interfere with the coils. For the modular coils of CFQS, because the last-closed magnetic surface was too close to the coils (∼5 cm), we did not perform the calculation. As we can see, despite the relatively large normal field components for helical coil configurations, the vacuum magnetic field generated by the modular and helical coils are in good agreement with each other in the core region for both magnetic configurations. On the other hand, the magnetic field near the coil shows quite different behavior depending on the type of the coil.
For CFQS, modular coils produce clean magnetic surfaces that extend very close to the coils, and there is almost no room between the last-closed magnetic surface and coils. On the other hand, helical coils produce a region of broken magnetic surfaces slightly away from the target surface that seem to construct divertor legs. For W7-X, the difference in the shape of the last-closed magnetic surface between the modular and the helical coils is not so significant. The modular coils also create broken magnetic surfaces outside of the last-closed magnetic surface and we can see structures that seem to be divertor legs. Note that in the W7-X operation the plasma edge is defined by the interaction of the magnetic field lines in the edge and the divertor plates. The divertor legs created by the helical coils, on the other hand, extend farther from the closed region because of the larger distance between the closed region and the virtual wall.
Finally, the rotational transform profiles are compared in figures 23 and 24. For W7-X, the errors from the VMEC result are almost at the same level for both modular and helical coils.

Conclusion
In this paper, we have developed a new numerical method for finding and optimizing the external coil of stellarators applying B-spline curve and the GA. We have introduced a cubic basis function for a closed B-spline curve to express a filamentary coil as a smooth interpolation of a set of discrete spatial control point. These control points, together with the coil position and coil current are optimized using the GA to minimize the objective function. The magnetic field component normal to the target magnetic surface is the main component of the objective function. The engineering constraints on the minimum curvature radius, minimum coil-to-coil distance, etc are also taken into account in the objective function.
We have developed a new numerical code GOSPEL based on the developed method. By introducing randomly-generated coil configuration for the initial generation, GOSPEL is designed to find the best coil arrangement without giving any specific initial guess. The implementation and applicability of the method has been verified applying GOSPEL to CFQS and W7-X. It has been shown that, as a general characteristic of GA, loss of the genetic diversity of the coil configurations tends to slow down the convergence of the objective function in the later phase of the optimization. We have successfully obtained modular coil configurations that well reproduce Poincáre plot and rotational transform profile for both of the target magnetic configurations. For CFQS, the error field produced by the obtained coils is at almost the same level as that by the original engineering design of CFQS obtained by NESCOIL, while the minimum curvature radius and coil-to-coil distance are increased. These results indicate not only the practicality but also an advantage of the three-dimensional B-spline curve in designing stellarator coil.
We have also explored helical coil configurations for these target magnetic configurations. Although LHD-type helical coils were robustly reproduced by GOSPEL using different seed value for the random number, it was found that GOSPEL tends to converge to three different types of coil arrangements depending on the seed value of the random number in the case of the optimized stellarators. Among these patterns, we found a specific pattern that can minimize the normal field component while being reasonably simple in shape. The optimized helical coils of this pattern generated magnetic configurations with divertor legs. In the region of the closed magnetic surfaces, the Poincáre plots and the rotational transform profiles were in a good agreements between the modular and helical coils.
In the optimization of helical coils targeting the optimized stellarators, a tendency of convergence to a local minimum was observed. Because of this tendency, we needed to manually adjust the calculation condition to obtain the coil configuration with small error field, despite using the GA. The cause of this problem is still not clear. However, it may be divided into three parts, namely, quality of initial generation, crossover law, and mutation law.
It is easily imaged that crossover between pattern A and pattern B in figure 18 may produce a child coil configuration with strange self-symmetric coil, because the phase of helical winding of the symmetric coils is opposite. A possible approach to avoid this problem is to explicitly introduce preset coil patterns and suppress crossover across each of the patterns. In this study, we fixed most of the parameters specific to the GA, such as α in the crossover, μ in the selection of parents, and the mutation rates. Dependency of convergence on these parameters, as well as the choice of crossover method and mutation laws, will be investigated in future.
From a physics and engineering point of view, it is of interest whether it is possible to design a feasible device based on helical coils obtained in this study. It should also be investigated if the divertor legs produced by these coils can be used as a robust guide of the heat and particle from the core plasma under the existence of plasma pressure and/or current. These points will be addressed in future works.