In silico optimization of critical currents in superconductors

TL;DR

This paper compares various optimization strategies to identify the best defect structures in superconductors that maximize their critical current capacity, using simulations of different defect geometries.

Contribution

It introduces a systematic comparison of local and global optimization methods for designing defect structures in superconductors to enhance critical currents.

Findings

Global optimization outperforms local methods in finding optimal defect configurations.

Optimal defect structures vary significantly with defect shape and placement.

The methods successfully predict defect arrangements that could improve superconductor performance.

Abstract

For many technological applications of superconductors the performance of a material is determined by the highest current it can carry losslessly - the critical current. In turn, the critical current can be controlled by adding non-superconducting defects in the superconductor matrix. Here we report on systematic comparison of different local and global optimization strategies to predict optimal structures of pinning centers leading to the highest possible critical currents. We demonstrate performance of these methods for a superconductor with randomly placed spherical, elliptical, and columnar defects.

Tables

Table 1. Table 1: Optimal PSO parameters 𝐪 opt = { S \mathbf{q}_{\mathrm{opt}}=\{S , ω 𝜔 \omega , ϕ p subscript italic-ϕ p \phi_{\mathrm{p}} , ϕ g } \phi_{\mathrm{g}}\} and F f , α ​ ( 𝐪 opt ) subscript 𝐹 𝑓 𝛼 subscript 𝐪 opt F_{f,\alpha}(\mathbf{q}_{\mathrm{opt}}) for the sphere function.

$n$	$\{S$	$\omega$	${\varphi }_p$	${\varphi }_g\}$	$F_{f,\alpha }$
2	5	0.22	0.93	1.93	84.7
3	5	0.36	1.35	1.68	131.3
4	8	0.23	0.80	1.96	183.7
5	10	0.18	0.99	1.96	235.5
6	9	0.36	1.55	1.55	290.5
7	12	0.27	1.18	1.75	348.4

Table 2. Table 2: Optimal adaptive pattern search parameters 𝐪 opt = { μ \mathbf{q}_{\mathrm{opt}}=\{\mu , σ 𝜎 \sigma , k s subscript 𝑘 s k_{\mathrm{s}} , k u } k_{\mathrm{u}}\} and F f , α ​ ( 𝐪 opt ) subscript 𝐹 𝑓 𝛼 subscript 𝐪 opt F_{f,\alpha}(\mathbf{q}_{\mathrm{opt}}) for the sphere function.

$n$	$\{\mu$	$\sigma$	$k_s$	$k_u\}$	$F_{f,\alpha }$
2	1	0.21	1.00	0.24	44.3
3	2	0.22	1.00	0.33	64.0
4	4	0.12	1.00	0.29	82.0
5	2	0.36	1.00	0.29	99.5
6	3	0.21	1.00	0.33	114.9
7	13	0.25	1.00	0.33	128.7

Table 3. Table 3: Optimal PSO parameters 𝐪 opt = { S \mathbf{q}_{\mathrm{opt}}=\{S , ω 𝜔 \omega , ϕ p subscript italic-ϕ p \phi_{\mathrm{p}} , ϕ g } \phi_{\mathrm{g}}\} , F f , α ​ ( 𝐪 opt ) subscript 𝐹 𝑓 𝛼 subscript 𝐪 opt F_{f,\alpha}(\mathbf{q}_{\mathrm{opt}}) , r f ​ ( 𝐪 opt ) subscript 𝑟 𝑓 subscript 𝐪 opt r_{f}(\mathbf{q}_{\mathrm{opt}}) , and N f , α ​ ( 𝐪 opt ) subscript 𝑁 𝑓 𝛼 subscript 𝐪 opt N_{f,\alpha}(\mathbf{q}_{\mathrm{opt}}) for the Rosenbrock function.

$n$	$\{S$	$\omega$	${\varphi }_p$	${\varphi }_g\}$	$F_{f,\alpha }$	$r_f$	$N_{f,\alpha }$
2	28	0.25	$-0.21$	1.58	823	0.99	1
3	34	0.29	$-0.15$	1.73	2942	0.95	2
4	40	0.35	$-0.17$	1.59	9352	0.73	4
5	49	0.28	$-0.21$	1.67	15650	0.61	5
6	61	0.30	$-0.20$	1.70	22339	0.56	6
7	73	0.27	$-0.19$	1.68	39277	0.53	7

Table 4. Table 4: Optimal adaptive pattern search parameters 𝐪 opt = { μ \mathbf{q}_{\mathrm{opt}}=\{\mu , σ 𝜎 \sigma , k s subscript 𝑘 s k_{\mathrm{s}} , k u } k_{\mathrm{u}}\} , F f , α ​ ( 𝐪 opt ) subscript 𝐹 𝑓 𝛼 subscript 𝐪 opt F_{f,\alpha}(\mathbf{q}_{\mathrm{opt}}) , r f ​ ( 𝐪 opt ) subscript 𝑟 𝑓 subscript 𝐪 opt r_{f}(\mathbf{q}_{\mathrm{opt}}) , and N f , α ​ ( 𝐪 opt ) subscript 𝑁 𝑓 𝛼 subscript 𝐪 opt N_{f,\alpha}(\mathbf{q}_{\mathrm{opt}}) for the Rosenbrock function.

$n$	$\{\mu$	$\sigma$	$k_s$	$k_u\}$	$F_{f,\alpha }$	$r_f$	$N_{f,\alpha }$
2	2	0.28	1.67	0.43	227	1.00	1
3	3	0.15	1.48	0.46	392	1.00	1
4	4	0.22	1.81	0.38	781	0.94	2
5	5	0.14	1.71	0.45	1384	0.92	2
6	6	0.37	1.96	0.36	1775	0.93	2
7	7	0.10	2.77	0.34	2993	0.92	2

Table 5. Table 5: Optimal PSO parameters 𝐪 opt = { S \mathbf{q}_{\mathrm{opt}}=\{S , ω 𝜔 \omega , ϕ p subscript italic-ϕ p \phi_{\mathrm{p}} , ϕ g } \phi_{\mathrm{g}}\} , F f , α ​ ( 𝐪 opt ) subscript 𝐹 𝑓 𝛼 subscript 𝐪 opt F_{f,\alpha}(\mathbf{q}_{\mathrm{opt}}) , r f ​ ( 𝐪 opt ) subscript 𝑟 𝑓 subscript 𝐪 opt r_{f}(\mathbf{q}_{\mathrm{opt}}) , and N f , α ​ ( 𝐪 opt ) subscript 𝑁 𝑓 𝛼 subscript 𝐪 opt N_{f,\alpha}(\mathbf{q}_{\mathrm{opt}}) for Rastrigin function.

$n$	$\{S$	$\omega$	${\varphi }_p$	${\varphi }_g\}$	$F_{f,\alpha }$	$r_f$	$N_{f,\alpha }$
2	50	0.25	2.00	1.00	934	1.00	1
3	145	0.58	2.08	0.89	4514	1.00	1
4	265	0.58	2.21	0.75	10953	1.00	1
5	365	0.48	2.22	0.77	20410	1.00	1
6	650	0.5	2.25	0.67	43939	0.998	1
7	945	0.51	2.29	0.65	81779	0.993	1

Table 6. Table 6: Optimal adaptive search parameters 𝐪 opt = { S \mathbf{q}_{\mathrm{opt}}=\{S , ω 𝜔 \omega , ϕ p subscript italic-ϕ p \phi_{\mathrm{p}} , ϕ g } \phi_{\mathrm{g}}\} , F f , α ​ ( 𝐪 opt ) subscript 𝐹 𝑓 𝛼 subscript 𝐪 opt F_{f,\alpha}(\mathbf{q}_{\mathrm{opt}}) , r f ​ ( 𝐪 opt ) subscript 𝑟 𝑓 subscript 𝐪 opt r_{f}(\mathbf{q}_{\mathrm{opt}}) , and N f , α ​ ( 𝐪 opt ) subscript 𝑁 𝑓 𝛼 subscript 𝐪 opt N_{f,\alpha}(\mathbf{q}_{\mathrm{opt}}) for Rastrigin function.

$n$	$\{\mu$	$\sigma$	$k_s$	$k_u\}$	$F_{f,\alpha }$	$r_f$	$N_{f,\alpha }$
2	3	0.25	1.79	0.33	1051	0.22	19
3	1	0.37	1.94	0.31	7972	0.04	113
4	1	0.39	1.69	0.29	27017	0.01	327
5	2	0.62	1.00	0.45	132134	0.004	1149
6	6	0.3	6.59	0.29	197052	0.003	1532
7	7	0.41	5.93	0.11	584564	0.001	4603

Table 7. Table 7: Optimal parameters of the pinscape for maximal critical current in the superconductor.

$n$	$𝐱=\{{\rho }_s$	$D_s\}$	$J_c$
2	0.22	$3.5$	0.0235

Equations

$${\mathbf{x}}_{\mathrm{opt}}=\arg \underset{\mathbf{x}\in \Omega }{\min }\{f(\mathbf{x})\},f_{\mathrm{opt}}=\underset{\mathbf{x}\in \Omega }{\min }\{f(\mathbf{x})\},$$

$$\mathcal{V}({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{n+1})=\frac{|\text{det}(\mathbf{V})|}{n!{\Delta }^n},$$

$${\bar{E}}_f(\mathbf{q})=\frac{1}{M}\sum _{i=1}^M{E}_{f,i}(\mathbf{q}),$$

$${\mathbf{q}}_{\mathrm{opt}}$$

$$F_{f,\alpha }(\mathbf{q})$$

$$N_{f,\alpha }(\mathbf{q})$$

$$f(\mathbf{x})=\sum _{i=1}^n{x}_i^2.$$

$$f(\mathbf{x})=\sum_{i=1}^{n-1}\bigl{[}100(x_{i+1}-x_{i}^{2})^{2}+(1-x_{i})^{2}\bigr{]}.$$

$$f(x,y)=100(y-x^2{)}^2+(1-x{)}^2.$$

$$f(\mathbf{x})=10n+\sum_{i=1}^{n}\bigl{[}x_{i}^{2}-10\cos(2\pi x_{i})\bigr{]}.$$

$$({\partial }_t+i\mu )\psi =ϵ(\mathbf{r})\psi -|\psi {|}^2\psi +(\nabla -i\mathbf{A}{)}^2\psi +\zeta (\mathbf{r},t),$$

$$\mathbf{J}=\mathrm{Im}\bigl{[}\psi^{*}(\nabla-i\mathbf{A})\psi\bigr{]}-\nabla\mu.$$

$$f(\mathbf{x})=-J_{\mathrm{c}}(\mathbf{x})$$

$$\mathbf{x}=(\rho_{\mathrm{s}},D_{\mathrm{s}}),\qquad N_{\mathrm{s}}=\biggl{\lfloor}\frac{6L_{x}L_{y}L_{z}\rho_{\mathrm{s}}}{\pi D_{\mathrm{s}}^{3}}\biggr{\rfloor},$$

$$\mathbf{x}=(\rho_{\mathrm{e}},D_{xy},D_{z}),\qquad N_{\mathrm{e}}=\biggl{\lfloor}\frac{4L_{x}L_{y}L_{z}\rho_{\mathrm{e}}}{\pi\kappa D_{xy}^{2}\tilde{D}_{z}}\biggr{\rfloor},$$

$$\mathbf{x}=(\rho_{\mathrm{s}},D_{\mathrm{s}},\rho_{\mathrm{c}},D_{\mathrm{c}}),\qquad N_{\mathrm{c}}=\biggl{\lfloor}\frac{4L_{x}L_{y}\rho_{\mathrm{s}}}{\pi D_{xy}^{2}}\biggr{\rfloor},$$

Full text

I Introduction

The most important feature of superconductors for high-performance applications Kwok et al. (2016) is its ability to carry large currents with almost no dissipative loss. At the same time, recent advances have made it possible to manufacture type-II-superconductor-based cables on industrial scales making large-scale application possible.Holesinger et al. (2008); Malozemoff (2012) The main dissipative mechanism in these type-II superconductors stems from the motion of magnetic vortices within the superconducting matrix. These are normal elastic filaments, carrying quantized magnetic fluxes,Abrikosov (1957) which appear in these superconductors in magnetic fields larger than the first critical field — in contrast to type-I superconductors, which lose their superconducting properties once the magnetic field enters the material. The highest amount of current that can be passed through a superconductor without dissipation is known as the critical current, which strongly depends on the vortex dynamics and its interaction with non-superconducting defects. Strictly speaking, the critical current in the presence of vortices is always zero at finite temperatures due to thermal creep, but it is conventionally defined as the current at which the dissipation/voltage reaches some threshold value. At sufficiently large currents the dissipation and associated heating of the superconductor will ultimately lead to the loss of superconductivity.

Despite all technological advances, typical critical currents are still way below the theoretical limit and from an economic point of view, too low to be competitive with conventional cables for large-distance energy transport. An obvious way to improve the critical current in existing superconductors is to impede the vortex motion more effectively. This can be accomplished by placing defects, which can “pin” magnetic vortices and thus prevent their motion,Blatter et al. (1994); Blatter and Geshkenbein (2003) in a sophisticated way. The entirety of all defects is also called the pinning landscape (or pinscape), which in conjunction with intrinsic material inhomogeneities and the sample geometry define the critical current of a given sample.

This defines the task, which we consider in this work: the optimization of the pinscape for highest possible critical currents. More precisely, we concentrate on defects or pinning centers, which can be controlled during the fabrication of a sample, like self-assembled inclusionsOrtalan et al. (2009); Sadovskyy et al. (2016a) or irradiation defects.Sadovskyy et al. (2016b); Leroux et al. (2015); Eley et al. (2017) The efficiency of the pinscape strongly depends on the shape, size, and arrangement of individual pinning centers. Indeed, bigger defects ensure a larger pinning force, but, at the same time, reduce the effective cross-section of the superconductor needed for current flow. The optimal pinscape also depends on the intended application, particularly on the type of superconductor and on the value and the direction of external magnetic fields. The question we address in this paper is, how one can find the best pinscape most effectively. This systematic prediction of optimal pinscapes for a given application aims at replacing the traditional trial-and-error approach.Sadovskyy et al. (2016b)

Here, we test several optimization strategies needed for a systematic improvement of the critical current in superconductors. We compare the efficiency of a global method (particle swarm optimization) and three local methods (Nelder-Mead method, pattern search, and adaptive pattern search) for a typical critical current optimization problem. The critical current for a given pinscape is calculated using a GPU-based iterative solver for the time-dependent Ginzburg-Landau (TDGL) equation describing type-II superconductors.Sadovskyy et al. (2015) This model correctly captures the vortex dynamicsBlatter et al. (1994); Blatter and Geshkenbein (2003); Vodolazov (2013) in superconductors in the vicinity of the critical temperature and is capable of reproducing experimental critical currents for a given pinscape.Berdiyorov et al. (2006a, b); Sadovskyy et al. (2016a, b, 2017) In addition, we provide a detailed analysis of these methods applied to several benchmark functions for comparison.

The article is organized as follows. In Sec. II we formulate the general optimization problem and describe the optimization methods studied here. In Sec. III we present a detailed comparison of the efficiency of the chosen optimization strategies on benchmark functions. In Sec. IV we briefly describe the TDGL model for superconductors and define three physically relevant optimization problems in Sec. V. The results of the optimization for these problems are discussed in Sec. VI. Finally, we summarize our results in Sec. VII.

II Optimization methods and problem formulation

Optimization methods are divided into two classes: local and global methods. Examples of global search methods include particle swarm optimization (PSO), cuckoo search, simulated annealing, etc. Here we focus on PSO because of its straight forward parallelizability. We compare this global PSO method to three local methods: pattern search (aka coordinate descent), adaptive pattern search, and the Nelder-Mead method (aka downhill simplex method). The pattern search and Nelder-Mead methods are standard local methods and their analysis, convergence properties, and pitfalls have been widely studied.Hooke and Jeeves (1961); Torczon (1997); Audet and Dennis Jr (2002); Dolan et al. (2003); Lagarias et al. (1998); McKinnon (1998); Rios and Sahinidis (2013) The adaptive pattern search method is a recent improvement on the traditional pattern search.Loshchilov et al. (2011) These methods can be used in conjunction with some more sophisticated methods, or routines, such as multi-level single linkage, which begins with many different initial starting points, and collects them into multiple sets, which depend on whether they are sufficiently close to a previously found local optimum. If they are not close to any previously found optimum, a local search is started.Rinnooy Kan and Timmer (1987a, b) The primary use behind this type of routine is to terminate local searches, which are falling into the basin of attraction of an optimum point already found, which reduces the computation time and number of function evaluations.

The general optimization problem of a function $f$ can be formulated as follows

[TABLE]

where $\mathbf{x}\in\Omega$ is a set of parameters in parameter space $\Omega$. $\mathbf{x}$ defines all variable properties in the system, which in turn, defines the objective function $f(\mathbf{x})$. In the case of a superconductor, the parameter set $\mathbf{x}$ can, e.g., determine the pinscape of the system and the corresponding objective function is $f(\mathbf{x})=-J_{\mathrm{c}}(\mathbf{x})$, where $J_{\mathrm{c}}(\mathbf{x})$ is the critical current of the superconductor with pinscape $\mathbf{x}$. Note, that the minimization of this $f(\mathbf{x})$ is equivalent to the maximization of $J_{\mathrm{c}}(\mathbf{x})$.

II.1 Particle swarm optimization

The PSO algorithm is a meta-heuristic global optimization algorithm.Kennedy (2011); Poli et al. (2007) Its convergence properties have been studied in a simplified form, where a single particle was used and the randomness in the algorithm was replaced by its averaged value.Trelea (2003) It performs well on all test problems, but is typically outshone by local methods when there was only a single minimum. The utility of this method is best seen on functions with multiple local optima, see Fig. LABEL:sub@fig:Rastrigin_function. This test function is relevant for problems having multiple local extrema, such as the problem of periodically arranged pinning centers in superconductors. In that case, many local maxima of $J_{\mathrm{c}}$ exist at integer values of the ratio of number of pinning centers to number of vortices.Sadovskyy et al. (2017)

The PSO has four main control parameters given by $\mathbf{q}=\{S$, $\omega$, $\phi_{\mathrm{p}}$, $\phi_{\mathrm{g}}\}$, where $S$ is the swarm size, $\omega$ the inertia of the individual particle (its tendency to move in its current direction), $\phi_{\mathrm{p}}$ and $\phi_{\mathrm{g}}$ are the weights for the particle to move towards its personal and global best point in parameter space, respectively. The pseudocode for PSO is presented in Listing 1.

The objective function is updated independently for each particle, which makes PSO parallelizable. In this way, a large architecture can make this method highly efficient for multimodal surfaces, where a large swarm size can be way more efficient in converging towards a global solution then local methods.

The biggest challenge is in determining an appropriate exit criterion for the routine. We use

1. (i)

the change, $|f(\mathbf{g}_{\mathrm{best}})-\langle f(\mathbf{g})\rangle_{\scriptscriptstyle M}|$, in the best found objective function value, $f(\mathbf{g}_{\mathrm{best}})$, with $f(\mathbf{g})$ averaged over the last $M$ iterations of PSO, $\langle f(\mathbf{g})\rangle_{\scriptscriptstyle M}$; or 2. (ii)

the average distance of the particle from the swarm global best position $K_{\mathrm{exit}}^{-1}\sum_{i=1}^{K_{\mathrm{exit}}}\|(\mathbf{x}_{i}-\mathbf{g})/(\mathbf{U}-\mathbf{L})\|_{2}$, where the sum is taken over $K_{\mathrm{exit}}$ particles with the best objective values, see lines 21, 24, and 25 in Listing 1. A small variant would be to use $\mathbf{p}_{i}$ instead $\mathbf{x}_{i}$ of each respective particle, which would be less sensitive to the exploratory nature of PSO (through the inertial term $\omega\mathbf{v}_{i}$).

The first exit criterion avoids taking an unnecessary amount of evaluations if there is no sufficient improvement (the method could then be terminated, switched to a local method, or restarted). The second was done to avoid cases where there may be many local optimum and so a few particles get trapped among them. This avoid premature stopping if many particles are spread out in parameter space and the improvement in $\mathbf{g}$ is slow. We stop when a certain proportion of the swarm has been attracted to the best found point. A higher proportion means more evaluations but also a higher probability that we converged to the correct optimum. Therefore, there is a balance between speed and accuracy. In the testing on the benchmark functions, we found no significant difference in accuracy for $K_{\mathrm{exit}}>S/2$, but a difference in speed. In the simulations presented below we chose $K_{\mathrm{exit}}=0.7S$.

II.2 Pattern search

Pattern search is a straightforward method, which starting from a random point, evaluates $2n+1$ points (including the initial point) where $n$ is the dimension of the parameter space, by moving a distance along each dimension in the search space. It then moves to the point which improves the function the most. The method then evaluates $2n-1$ new points (it does not need to re-evaluate the point its on or the point it came from). If no improvement is made, the step size is reduced. Once the step size is below some threshold, it exits out of the loop. Along with its simplicity comes its ability to converge to non-stationary problems on some relatively simple problems.Powell (1973) It has a particularly difficult time with functions with coordinate systems which are highly correlated, e.g., roughly speaking, the function gradient is not along any of the main coordinate axes, such as the Rosenbrock function shown in Fig. LABEL:sub@fig:Rosenbrock_function.

II.3 Adaptive pattern search

The adaptive pattern search method is a recent modification of pattern search. The algorithm works similar to pattern search, however, while searching the parameter space, it adapts the coordinate system to achieve faster convergence, see Listing 2 adapted from Ref. Loshchilov et al., 2011. The presented adaptive encoding is similar to principal component analysis (PCA). The covariance matrix $\mathbf{C}$ and the transformation matrix $\mathbf{B}$ are updated using the most successful $\mu$ points. The transformation matrix modifies the coordinates to make them as uncorrelated as possible. Unlike PCA which typically looks to reduce dimensionality by retaining only eigenvectors corresponding to the largest eigenvalues, adaptive encoding retains all components.

Apart from the adaptive encoding, adaptive pattern search is very similar to pattern search. Suppose, we have an $n$-dimensional optimization problem, it initially searches in each coordinate direction, independently recording the new optimal as it progresses (in contrast to the generic pattern search we employed, which is not done sequentially). Then it keeps the best $\mu$ points, $\mu<2n+1$, and Listing 2 is called. This transformation is updated after each sweep of all $n$ dimensions ($2n-1$ points) and applied in the next iteration. There are four parameters $\mathbf{q}=\{\mu$, $\sigma$, $k_{\mathrm{s}}$, $k_{\mathrm{u}}\}$ which control the success and efficiency of this method: $\mu$ is the number of points used in the adaptive encoding function call, $\sigma$ is the initial step size, and $k_{\mathrm{s}}$ ($k_{\mathrm{u}}$) are the increase (decrease) of the step size upon successful (unsuccessful) improvement of the function value. The convergence rate was observed to be most sensitive to the parameter $\mu$. The utility of this method is best seen when applied to the Rosenbrock function, where pattern search needs over 20,000 iterations to converge in a two-dimensional parameter space, while the adaptive pattern search method typically converged in less than 1,000 iterations.

The key for the method’s improved performance on Rosenbrock-type functions, is the adaptive encoding part of the algorithm. At this point we remark that the adaptive pattern search is typically a more efficient/performing method compared to the generic pattern search, but comparison to the latter offers a useful way to measure the overall shape of the optimal solution. We can logically deduce from the ratio of performance between these two methods, that if adaptive pattern search is orders of magnitude better than pattern search, then this implies that the parameters of interest are highly correlated.

II.4 Nelder-Mead method

The Nelder-Mead method (or downhill simplex) was chosen for its relative simple structure and its independence on the choice of the chosen coordinate system (it does not move along each dimension sequentially). This is most easily seen in comparing the method against the Rosenbrock function. This method utilizes simplexes, which are polytopes with $n+1$ points (or vertices) in $n$ dimensions (i.e., triangles for $n=2$, or tetrahedrons in $n=3$).

Over time, many variants of Nelder-Mead have been conceived.Price et al. (2002); Nazareth and Tseng (2002); Conn et al. (2009) Following Ref. Conn et al., 2009, Listing 3 describes the general Nelder-Mead method. Here we used standard values for parameters $\alpha=1$, $\gamma=2$, $\rho=1/2$, and $\sigma=1/2$. The exit criteria in line 25 is satisfied when the points are within a certain distance and is defined by the “volume”

[TABLE]

where $\Delta=\min\|\mathbf{x}_{i}-\mathbf{x}_{j}\|$ and matrix $\mathbf{V}=(\mathbf{\dot{x}}_{1}$, $\ldots$, $\mathbf{\dot{x}}_{n})$ consists of $n$ vectors $\mathbf{\dot{x}}_{i}=\mathbf{x}_{i+1}-\mathbf{x}_{1}$. This choice of the exit criteria was justified in Ref. Conn et al., 2009.

III Testing on benchmark functions

The above mentioned optimization methods are tested on the three benchmark functions shown in Fig. 1: the sphere function, the Rosenbrock function, and the Rastrigin function. These were chosen for their similarity to previously obtained critical-current surfaces in superconductors for low-dimensional sets of parameters. We expect scenarios where either a single optimum or multiple local with one global extrema exists. As a particular and practical example, the behavior of the solvers on the Rosenbrock function is related to the optimization of a pinscape consisting of spherical defects with two parameters being the number and diameter of the defects. In this case, it turns out that the Rosenbrock-type $J_{\mathrm{c}}$ surface can be removed by replacing the number of defects by the volume fraction. In general, however, one cannot be sure that an appropriate transformation exists or can even be found.

The study was broken up as follows: optimal parameters for the PSO and adaptive pattern search are obtained using the three benchmark functions. These were found by overlaying each respective optimization routine with PSO. For example, adaptive pattern search has tuning parameters $\mathbf{q}=\{\mu$, $k_{\mathrm{s}}$, $k_{\mathrm{u}}$, $\sigma\}$, see Sec. II.3. The nested PSO algorithm then searches through this parameter space in an attempt to find the optimal parameters for the algorithm. We initially considered the function

[TABLE]

where $E_{f,i}$ is the number of function evaluations required to find the global optimum for a function $f$, $\mathbf{q}$ are the parameters used for the optimization routine, and $M$ are the total number of simulations which successfully found the global optimum. However, this is not the most useful measure as it does not take into account the rate at which the algorithm successfully finds the global optimum. Indeed, defining $r_{f}(\mathbf{q})$ as the rate at which a correct solution is found (to within a specified tolerance) for a set of the method parameters $\mathbf{q}$. Then, it may turn out that $\bar{E}_{f}(\mathbf{q}_{1})<\bar{E}_{f}(\mathbf{q}_{2})$, but $r_{f}(\mathbf{q}_{1})\ll r_{f}(\mathbf{q}_{2})$ for some certain $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$. It then may happen that we would require many more runs for $\mathbf{q}_{1}$ so that $\mathbf{q}_{2}$ was actually the better choice. Therefore, we make a refinement to Eq. (2), and obtain the optimal parameters $\mathbf{q}_{\mathrm{opt}}$ for each optimization method by solving the following auxiliary optimization problem:

[TABLE]

where $N_{f,\alpha}$ is the number of iterations needed to be at least $\alpha$ sure that we have found the global solution and ${\bar{E}}_{f}(\mathbf{q})$ was defined by Eq. (2). In this work, we used $\alpha=0.99$. The dimensionality in the function was absorbed into the optimization problem, and an analysis of the problem dimensionality and number of iterations was tested.

We sampled $10^{3}$ different starting configurations $\mathbf{x}$, where $x_{i}\in[-10,10]$ for the algorithms and then ran the nested PSO algorithm 10 times for each dimension and each benchmark function. The best parameters were recorded. Once these were obtained, we tested all the algorithms mentioned using the same starting points (when possible) and compared the performance. The algorithms each had a maximum iteration number of $10^{3}n^{2}$ where $n=|\Omega|$ (the exception being PSO, where we used $\max(10^{3}n^{2},100S)$ to ensure at least 100 possible iterations of PSO) is the dimension of the original optimization problem (1), and would exit out of the loop with a tolerance of $10^{-3}n^{2}$.

III.1 Sphere function

The sphere function shown in Fig. LABEL:sub@fig:sphere_function is defined by

[TABLE]

Tables 1 and 2 show a tabulated view of the effectiveness of the chosen methods to find the extremum of this function in the given number of iterations. The function is very simple and the coordinates are uncorrelated. Thus, a wide range of parameters actually turned out to be similarly effective. A comparison of the method’s performances is presented in Fig. 2 by using 100 random initial starting configurations and employing Eq. (3).

III.2 Rosenbrock function

The Rosenbrock function shown in Fig. LABEL:sub@fig:Rosenbrock_function is a standard test for optimization methods and is given by

[TABLE]

Optimization of this function demonstrates the utility of coordinate independent local search methods (like Nelder-Mead or adaptive pattern search). This is due to the fact that the minimum is contained inside a parabolic valley which requires constant shrinking of the step size for pattern search to make progress. Consider the 2D case

[TABLE]

A change of variables to the $(u,v)$ plane given by $u=1-x$ and $v=y-x^{2}$ leads to elliptical level sets $g(u,v)=u^{2}+100v^{2}$, which is much more favorable to coordinate-dependent methods. In practice it is usually difficult or impossible to find the appropriate transform converting to elliptical level sets. Nevertheless, this gives a useful test of the morphology of the surface of $J_{\mathrm{c}}$ by comparing iterations between coordinate-dependent and coordinate-independent methods.

Tables 3 and 4 show the effectiveness of the PSO and adaptive pattern search methods to minimize the Rosenbrock function in a given number of iterations. As we can see from the table, the PSO method is much slower than Nelder-Mead and adaptive pattern search. The optimization of parameters for PSO have revealed that the dimensionality and swarm size are (perhaps not surprisingly) correlated. As the dimensionality increases, the optimal swarm size (holding other parameters fixed) increases. To verify this, we sampled 100 random starting points for swarm sizes between 10–200; the results are shown in Fig. 4.

III.3 Rastrigin function

The third test function we consider here, is the Rastrigin function [see Fig. LABEL:sub@fig:Rastrigin_function], which is the sum of a sphere function and periodic oscillations, having multiple minima. This type of objective function is similar to the objective function of superconductors having periodic pinscapes. It is defined by

[TABLE]

This is where PSO performs particularly well compared to the local search methods. Table 5 shows its effectiveness to solve the Rastrigin function optimization problem in the given number of iterations. The local methods do a comparatively poor job as can be seen in Table 6 and Fig. 5. The optimal swarm size also grows rapidly with increasing dimensionality and we sampled 100 random starting points for swarm sizes between 50 and 1,000; the results are shown in Fig. 6.

IV Model for the critical current in superconductors

Here we define the objective function for the optimization problem in superconductors mentioned before. For the description of the vortex dynamics in strong type-II superconductors, we use the TDGL equation for the superconducting order parameter $\psi=\psi(\mathbf{r},t)$,

[TABLE]

which is solved numerically. Here $\mu=\mu(\mathbf{r},t)$ is the chemical potential, $\mathbf{A}$ is the vector potential associated with an external magnetic field $\mathbf{B}$ as $\mathbf{B}=\nabla\times\mathbf{A}$, and $\zeta(\mathbf{r},t)$ is the temperature-dependent $\delta$-correlated Langevin term. The unit of length is given by the superconducting coherence length $\xi$ and the unit of magnetic field is given by the upper critical field $H_{\mathrm{c2}}$. See Ref. Sadovskyy et al., 2015 for the details of TDGL model implementation and physical units. The current density is given by the expression

[TABLE]

To determine the critical current value — the maximal current, which can flow through the superconductor without dissipation, — we use a finite-electrical-field criterion. Specifically, we choose a certain small external electric field, $E_{\mathrm{c}}=10^{-4}$, which measures the dissipation and adjust the applied external current, $J$, to reach this electrical-field/dissipation level on average during the simulation. The time-averaged value of external current over a steady state gives the critical current, $J_{\mathrm{c}}=\langle J\rangle$.

The critical current in the presence of an external magnetic field is mostly defined by the pattern of non-superconducting defects, including their sizes, shapes, and spatial distribution, which prevents vortices from moving under the influence of the Lorentz force, $\mathbf{f}_{\mathrm{\scriptscriptstyle L}}=\mathbf{J}\times\mathbf{B}$.

The pinscape is characterized by a set of parameters $\mathbf{x}$, which corresponds to the objective function

[TABLE]

used in optimization problem (1). Each element $\mathbf{x}$ of the parameter space $\Omega$ describes the pinscape in the superconductor, e.g. the shape of each defect and their spatial distribution. The optimal configuration of the defects $\mathbf{x}_{\mathrm{opt}}$ corresponds to the minimization of the objective function, $f_{\mathrm{opt}}=f(\mathbf{x}_{\mathrm{opt}})$.

Knowledge of the shape and behavior of the function $J_{\mathrm{c}}(\mathbf{x})$ is not known a priori. In addition if we consider, for example, a random placement of defects in the domain, each realization can yield slightly different values for $J_{\mathrm{c}}$ — for the same $\mathbf{x}$ (if it does not describe the positions of all defects explicitly), such that averaging is required. In that case one can expect, as the number of random simulations tends to infinity, $J_{\mathrm{c}}$ approaches a “true” value due to self-averaging. For a finite number of trials it leads, however, to a noisy $J_{\mathrm{c}}$ surface. This can create difficulties for local methods to converge to the global solution. A modification of the local methods to multi-level starting points, helps to overcome the noise in these types of problems.

To obtain the critical current $J_{\mathrm{c}}(\mathbf{x})$, we solve the TDGL equation in the domain of interest with a specified parameter set $\mathbf{x}$. Each evaluation of $J_{\mathrm{c}}$ is relatively expensive (typically, a few GPU-hours), but can be performed independently for different $\mathbf{x}$. In order to reduce noise in $J_{\mathrm{c}}$, one can (i) average it over several realizations of random positions of defects and/or (ii) increase the system size. Both techniques naturally increase the computation time of $J_{\mathrm{c}}$. The final value defines the objective function for the optimization problem.

V Optimization of the critical current

The objective function in Eq. (9) is then used to solve the general optimization problem for the pinscape defined in Eq. (1). Typical pinscapes with like defects are described by parameter spaces $\Omega$ having $n=2$ to $8$ dimensions. Here we consider three particular and important cases for $n=2$, $3$, $4$ to analyze the described optimization strategies. For practical applications also the robustness of an optimum is important, since e.g. the size and shape of defects cannot be controlled exactly in an experiment. Therefore, the found optimal critical current should be relatively insensitive to small perturbations in the optimal pinscape $\mathbf{x}_{\mathrm{opt}}$. In the examples studied here, this condition is fulfilled.

For the simulations needed to obtain the critical current for a given pinscape, we use a three-dimensional superconducting sample with dimensions $64\times 64\times 64$ in units of coherence length $\xi$ with (quasi-)periodic boundary conditions. The external current $J$ is applied in the $x$-direction, perpendicular to the magnetic field $B$ in the $z$-direction. All defects are modeled by local variation of the critical temperature, such that defects are regions being in the normal state. To this end the system is discretized by a regular cubic grid with grid resolution of half a coherence length. As demonstrated in earlier works,Sadovskyy et al. (2016a, b) this resolution is sufficient to capture the involved physical processes correctly. A pattern generator then creates the pinscape according to the defining parameter set $\mathbf{x}$, by assigning each grid point either to the superconductor or to the defect, distinguished by a high or low local $T_{\mathrm{c}}$ value. Defect regions with low $T_{\mathrm{c}}$ values typically consist of many connected grid point.

First, we consider randomly placed spherical defects (typically self-assembled inclusions in real sample) in a three-dimensional system. All inclusions have the same diameter, such that the parameter space for this problem is defined by the defect density $\rho_{\mathrm{s}}$ within the sample and their diameter $D_{\mathrm{s}}$. As mentioned before, this set of parameters is not unique, such that the same pinscape could be described by a different $\mathbf{x}$ (e.g., the diameter and total number of defects). This two-parameter pinscape problem was considered recently in Ref. Koshelev et al., 2016 by sampling of the parameter space, which was possible in that case.

Another example of a two-parameter problem would be a system with two kinds of defects having fixed shape, e.g., a superconductor having columnar defects of fixed diameter imprinted (typically introduced by high-energy heavy-ion irradiation in a real sample), which already had intrinsic (usually chemically grown) nano-rod pinning centers.Sadovskyy et al. (2016b) In this case the parameter space will be defined by the concentration $\rho_{\mathrm{c}}$ of columnar defects and the concentration $\rho_{\mathrm{n}}$ of the intrinsic nano-rod defects, i.e., $\mathbf{x}=(\rho_{\mathrm{c}},\rho_{\mathrm{n}})$.

For the first problem we consider here, the pinscape is defined by

[TABLE]

where $N_{\mathrm{s}}$ was defined in Eq. (10a). In the following we discuss the results for these three optimization problems as well as the benchmark function optimizations.

VI Results and discussions

Discussion of benchmark function optimizations.

The optimization of parameters for adaptive pattern search and PSO has led to some interesting results. First, the nested PSO algorithm finds many different “good” choices for the parameters needed. This is to be expected, since each problem is unique, and it would be surprising if multiple sets of parameters did not yield similar results. What was surprising, is the sensitivity of the parameters to small deviations from optima that were obtained on Rastrigin’s function. A good example would be in a particular optimization run of adaptive pattern search’s parameters for $n=6$. We found that the set $\{\mu$, $\sigma$, $k_{\mathrm{s}}$, $k_{\mathrm{u}}\}=\{6$, $0.3$, $6.59$, $0.29\}$ was the best for that particular run found with $2\cdot 10^{5}$ function evaluations needed. Just changing one of these parameters slightly to $\{6$, $0.3$, $6.59$, $0.39\}$, needs $4.25\cdot 10^{5}$ function evaluations.

This sensitivity was not observed in the other two benchmark functions, meaning that the sensitivity of these parameters for the optimization routines are highly dependent on the type of function, being particularly sensitive to functions with multiple optima. At the other end of this spectrum, it was observed that many different sets of parameter values led to a minimal amount of evaluations. Essentially, there was not a large difference between the best observed set, which is reported in Tables 1 and 2, and other choices. From this we can conclude that surfaces with uncorrelated dimensional space and with small number of optima are more likely to be insensitive to the choice of parameters and hence do not need to be tuned for the particular problem.

The Rastrigin function is particularly useful in analyzing the effect multiple optima have on local methods. One can find the number of extrema as a function of dimension size $n$. The derivative of the Rastrigin function in each respective coordinate is set to 0 and is given by $x_{i}+10\pi\sin(2\pi x_{i})=0$. We recall that the bounds for our benchmark functions are $x_{i}\in[-10,10]$. It is easy to see that due to symmetry, we only need to consider $x_{i}>0$. We further note that solutions are not possible for $x_{i}>10\pi$ and for each oscillation period, one obtains two extrema. Since the period is 1, this implies 10 minima (20 extremum which by symmetry has 10 maxima and 10 minima) within the above bounds. Furthermore, symmetry implies 10 more in the negative direction plus the obvious $x_{i}=0$ (global) minimum, which results in 21 minima total. It then follows that the number of minima in our problem for $n$ dimensions is $20^{n}+1\approx 20^{n}$.

Another aspect of the optimization is how high the probability for a found optimum being the global optimum is. For this we can consider any local method which is run $R$ times, and suppose the method has an equal probability to converge to any of the $m$ minima. We also assume $m\gg 1$, and we want to be $\alpha=0.99$ sure we have found the global solution. This would require $R\approx 4.6m$ [or $R=O(m)$]. Thus the number of local runs should scale approximately linearly with increasing minima. With just $n=2$ for the Rastrigin function for example would imply that a local method would need to be run $R\approx 2\cdot 10^{3}$ times to be 99% sure we have converged to the global solution. Again, it is important to note that this number assumed each extremum was equally likely, that is, each of their respective basins of attractions were the same size. This assumption in general is not true, but it provides an insight into how large the number of runs of a local method would be needed. For example, suppose a local method converges in 50 evaluations on average to a local extremum. In two dimensions this would require on average close to $4\cdot 10^{5}$ total evaluations of the objective function. Compare this with the number needed for PSO. In reality this turns out to be an overestimation, the total number of evaluations for local methods are far lower, but this analysis provides some insight as to why they suffer at this level and how the number of evaluations will scale with increasing number of optimums.

In contrast, PSO and other global methods are more robust to the number of local optima as can be seen by comparing Figs. 3 and 5, the number of function evaluations is on the same order for both Rosenbrock and Rastrigin functions. A possible reason why PSO behaves similarly for Rosenbrock and Rastrigin is because there is only a ”1D path” towards the minimum for Rosenbrock, while the sphere function has a 50% chance to lower the objective in each direction. Therefore, particles are as likely to become stuck in the Rosenbrock valley as they are in a Rastrigin local minimum.

The PSO optimal parameters also provide some insight into the nature of the function. It helps to first interpret what each parameter does. Increasing $S$ obviously leads to a greater chance of finding the best solution, the cost is an increase in function evaluations. The inertial weight $\omega<1$ is not a coincidence. As the particles approach the global best, $\mathbf{v}_{i}\leftarrow\omega\mathbf{v}_{i}$ and this has a fixed point at $\mathbf{v}_{i}=0$. The $\phi_{\mathrm{p}}$ and $\phi_{\mathrm{g}}$ are quantities which affect the exploratory nature of the algorithm and the weight associated with the global information that the swarm provides to each individual particle. As $\phi_{\mathrm{g}}\to 0$ each of the particles behave more independently and their movement becomes more individualized. As $\phi_{\mathrm{p}}\to 0$, we are essentially telling each particle to only trust the group on the correct direction to go.

Armed with this interpretation, we can understand some of the results from the tables in Sec. III. The sphere function and Rosenbrock both have few minima (sphere has only one, and the Rosenbrock has no more than two for $n\leqslant 7$). Therefore we expect $\phi_{\mathrm{g}}\geqslant\phi_{\mathrm{p}}$. However, the results are different in Table 5, which shows $\phi_{\mathrm{p}}$ increasing and $\phi_{\mathrm{g}}$ decreasing as $n$ (number of optima) increases. This makes sense, with multiple optima, we want the swarm to behave more independently and not get caught among the many local optima. However we don’t want $\phi_{\mathrm{g}}=0$ or else we would take very long to converge.

Discussion of critical current optimizations.

Figures 8, 9, and 10 show the optimization “paths” for all four optimization methods for the 2D [Eq. (10a) with parameter space dimension $n=2$], 3D [Eq. (10b) with $n=3$], and 4D [Eq. (10c) with $n=4$] problems, respectively. These plots show the values of the currently best critical current and the associated parameter set $\mathbf{x}$ as function of number of evaluations of the objective function (solid lines). If a new optimal value is found it is marked by circles, while all other evaluations (with objective values above the current minimum) are marked by gray dots. The optimal values are shown by dashed horizontal lines. Note, that the pattern search algorithm is not sequential, meaning that a minimum is only calculated after a certain number of evaluations (up to $2n+1$), which defines an iteration of the method. The iterations are marked by vertical dotted lines in the plots. If the number of evaluations per iteration is less that $2n+1$, the remaining objective values were calculated before and are just looked up in a database. Iterations are shaded according to the step sizes in parameter space. The optimization results are summarized in Table 7. In the 2D case for randomly placed spherical inclusions, the optimal values reproduce the ones of Ref. Koshelev et al., 2016 at magnetic field $B=0.1H_{\mathrm{c2}}$. In the 3D case for spheroidal inclusions, the optimal pinscape corresponds to an infinite size in $z$ direction (larger then the system size), i.e., when the spheroids have “evolved” to columnar inclusion, and the optimal value of the critical current is increased by $\sim 50\%$. In the 4D case, for a combination of the columnar and spherical defects, the 3D result is confirmed, since the volume fraction for the spheres vanish and the best pinscape turns out to be consisting of columnar defects in $z$ direction only. The diameter and volume fraction of the best columnar defects in the 3D and 4D cases is accordingly the same.

As mentioned before, in order to obtain a good critical current value for randomly placed inclusions, either a large number of random realizations or much larger systems are needed. While the latter is easier to run, we chose the former for performance reasons, such that a single objective function evaluation requires averaging the critical current over many pinscapes, described by the same $\mathbf{x}$. However, even with a large number of realizations, the critical current as a function of $\mathbf{x}$ can still be quite noisy. In the worst case this leads to local methods failing to reach the optimal solution. In contrast to that, global methods like PSO are much more resistant to noise in the objective function.

First, let us analyze the two parameter optimization more closely (see Fig. 8). All optimization methods converge to the same optimal pinscape. However, the optimization “paths” are quite different as one could expect: although all methods reach a critical current value close to the optimum after about 10–20 evaluations — indicating a robust, relatively flat peak in parameter space, the actual optimal configuration takes quite a while to be reached: It takes about twice as many evaluations for the local methods and an order of magnitude more for the global PSO method. Note, that defining an appropriate exit criterion for the optimization can be challenging, since the improvement in the critical current can be marginal, while the configuration can still change noticeably, as is clearly seen in the Nelder-Mead path.

Increasing the parameter space to three dimensions, (Fig. 9), shows that all methods need more evaluations as expected (scaling roughly linear with the dimension) and again the local methods reach a close to optimal critical current significantly faster than PSO. However, one can see that the local methods already start to have convergence problems, which is most likely a result of the noise in the objective function. The figure shows that the adaptive pattern search is a bit more efficient dealing with this noise than simple coordinate decent. However, in this case all methods still reach a configuration close to the optimal one.

An obvious observation is that increasing the parameter space, cannot lower the optimal value of the objective function. Indeed, the chosen examples are in a sense supersets of each other with higher dimensionality of $\Omega$, such that the best $J_{\mathrm{c}}$ value either increases or remains the same. When going over from the 3D to 4D optimization problem, we kept the result for the optimal $z$-diameter of the defects, i.e., that the optimal defects are cylinders, but added back spherical defects to the system. In the 4D case one could then imagine three different scenarios for the global optimal configuration: coexistence of both defect types, only cylindrical defects remain, or only spherical are optimal. Due to the results of the 3D problem, the latter is rather unlikely, since vortices and columnar defects are well aligned and the pinning potential of the defects are therefore “optimal.”

As seen from Fig. 10, the local methods converge to different optimal solutions. The coexistence scenario appears to be a good pinscape. However, the largest $J_{\mathrm{c}}$ value is still below the cylinder-only solution. One can expect that the optimal configuration is dependent on the angle between applied magnetic field and the main axis of the cylindrical defects. If they deviate from the parallel setup we studied here, the coexistence of spherical and columnar defects might be optimal, but it is clear that at a 90∘ angle between field and columns only spherical defects would “survive” in an optimization problem (for isotropic superconductors).

When we compare our results for the $J_{\mathrm{c}}$ optimization to the studied test functions before, we can conclude, that the $J_{\mathrm{c}}$ surface for the chosen parameters is roughly uncorrelated, since the number of evaluations required for adaptive pattern search and pattern search to converge are nearly the same.

Overall, we see that the PSO algorithm works reliable in all cases, at the cost of about 10 times more evaluations needed to reach the optimum, where often many evaluations, on the order $10^{2}$–$10^{3}$, are needed for a marginal gain in the objective function ($<1\%$). The local methods start to fail in higher dimensional parameter spaces, most likely due to noise in the objective function. As mentioned before, combining the local searches with a global method can potentially reduce the overall number of evaluations. Another way to mitigate this problem is to use more realizations or larger system sizes, which should smooth out the error from the random placement of defects.

Another benefit of using a global method such as PSO, is that it consistently finds better solution than the local methods. In higher dimensional spaces $n\gtrsim 4$, taking a larger swarm size $S$ becomes necessary to ensure convergence. The local methods, however, converge faster. But in any case there is no guarantee that the found configuration is the optimum. One could only estimate a probability that the found result is close to the optimal one.

Finally, we remark that larger oscillations in the optimization paths for the local methods in higher dimensions suggest a flatter global maximum of $J_{\mathrm{c}}$ or a more robust critical current when one deviates from the optimal configuring. This is particularly important for practical applications.

VII Conclusions

We have tested various methods for optimization problems and determined that particle swarm optimization and adaptive pattern search performed the best on the benchmark functions as well as in the physical optimization problem presented. Particle swarm optimization in particular does a relatively good job at handling the noisy surface of the objective function, but takes a long time to converge, where a significant amount of time is used for marginal improvement. This is a result of the difficulty to define an exit criterion. On the other side, adaptive pattern search can get caught in local extrema, which are either physical or a result of the noise. However, the number of evaluations are often much less (on the order of 7–10 times). As mentioned before, a multi-level local method could be used to mitigate this, which should make adaptive pattern search competitive with particle swarm optimization. Although we presented only some results for the pinscape optimization in superconductors, our studies of different optimization strategies can be important for a variety of different physical systems where the evaluation of the objective function is expensive, e.g., requiring first principle calculations or molecular dynamics simulations.

Acknowledgements.

This work was supported by the Scientific Discovery through Advanced Computing (SciDAC) program funded by U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Science, Division of Materials Science and Engineering. For the computational work we used resources of the Argonne Leadership Computing Facility (DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357) and the GPU cluster GAEA at Northern Illinois University.