simsopt.solve package

simsopt.solve.GPMO(pm_opt, algorithm='baseline', **kwargs)

GPMO is a greedy algorithm for the permanent magnet optimization problem.

GPMO is an alternative to to the relax-and-split algorithm. Full-strength magnets are placed one-by-one according to minimize the MSE (fB). Allows for a number of keyword arguments that facilitate some basic backtracking (error correction) and placing magnets together so no isolated magnets occur.

Parameters:

pm_opt – The grid of permanent magnets to optimize. PermanentMagnetGrid instance
algorithm –
The type of greedy algorithm to use. Options are

baseline:
the simple implementation of GPMO,

multi:
GPMO, but placing multiple magnets per iteration,

backtracking:
backtrack every few hundred iterations to improve the solution,

ArbVec:
the simple implementation of GPMO, but with arbitrary oriented polarization vectors,

ArbVec_backtracking:
same as above but w/ backtracking.

Easiest algorithm to use is ‘baseline’ but most effective algorithm is ‘ArbVec_backtracking’.
kwargs –
Keyword arguments for the GPMO algorithm and its variants. The following variables can be passed:

K: integer
Maximum number of GPMO iterations to run.

nhistory: integer
Every ‘nhistory’ iterations, the loss terms are recorded.

Nadjacent: integer
Number of neighbor cells to consider ‘adjacent’ to one another, for the purposes of placing multiple magnets or doing backtracking. Not to be used with ‘baseline’ and ‘ArbVec’.

dipole_grid_xyz: 2D numpy array, shape (ndipoles, 3).
XYZ coordinates of the permanent magnet locations. Needed for figuring out which permanent magnets are adjacent to one another. Not a keyword argument for ‘baseline’ and ‘ArbVec’.

max_nMagnets: integer.
Maximum number of magnets to place before algorithm quits. Only a keyword argument for ‘backtracking’ and ‘ArbVec_backtracking’ since without any backtracking, this is the same parameter as ‘K’.

backtracking: integer.
Every ‘backtracking’ iterations, a backtracking is performed to remove suboptimal magnets. Only a keyword argument for ‘backtracking’ and ‘ArbVec_backtracking’ algorithms.

thresh_angle: float.
If the angle between adjacent dipole moments > thresh_angle, these dipoles are considered suboptimal and liable to removed during a backtracking step. Only a keyword argument for ‘ArbVec_backtracking’ algorithm.

single_direction: int, must be = 0, 1, or 2.
Specify to only place magnets with orientations in a single direction, e.g. only magnets pointed in the +- x direction. Keyword argument only for ‘baseline’, ‘multi’, and ‘backtracking’ since the ‘ArbVec…’ algorithms have local coordinate systems and therefore can specify the same constraint and much more via the ‘pol_vectors’ argument.

reg_l2: float.
L2 regularization value, applied through the mmax argument in the GPMO algorithm. See the paper for how this works.

verbose: bool.
If True, print out the algorithm progress every ‘nhistory’ iterations. Also needed to record the algorithm history.

Returns:

Tuple of (errors, Bn_errors, m_history)

errors:: Total optimization loss values, recorded every ‘nhistory’ iterations.
Bn_errors:: \(|Bn|\) errors, recorded every ‘nhistory’ iterations.
m_history:: Solution for the permanent magnets, recorded after ‘nhistory’ iterations.

simsopt.solve.bnorm_obj_matrices(wframe, surf_plas, ext_field=None, area_weighted=True, bnorm_target=None, verbose=True)

Computes the normal field matrix and target field vector used for determining the squared-flux objective for wireframe current optimizations.

Parameters:

wframe (instance of the ToroidalWireframe class) – Wireframe whose segment currents are to be optimized
surf_plas (Surface class instance) – Surface of the target plasma, on which test points are placed for evaluation of the normal field.
ext_field (MagneticField class instance (optional)) – Constant external field assumed to be present in addition to the field produced by the wireframe.
bnorm_target (1d double array (optional)) – Target value of the normal field on the plasma boundary to be produced by the combination of the wireframe and ext_field. Zero by default. Test points on the plasma boundary corresponding to the elements of bnorm_target must agree with the test points of surf_plas.
area_weighted (boolean (optional)) – Determines whether the normal field matrix and target normal field vector elements are weighted by the square root of the area ascribed to each test point on the plasma boundary. If true, the weightings will be applied and therefore the optimization objective will be the square integral of the normal field on the plasma boundary. If false, the weightings will NOT be applied, and the optimization objective will be the sum of squares of the normal field at each test point.
verbose (boolean (optional)) – If true, will print progress to screen with durations of certain steps

Returns:

Amat (2d double array) – Inductance matrix relating normal field at test points on the plasma boundary to currents in each segment.
bvec (double array (column vector)) – Vector giving the target values of the normal field at each test point on the plasma boundary.

simsopt.solve.constrained_mpi_solve(prob: ConstrainedProblem, mpi: MpiPartition, grad: bool = False, abs_step: float = 1e-07, rel_step: float = 0.0, diff_method: str = 'forward', opt_method: str = 'SLSQP', options: dict = None)

Solve a constrained minimization problem using MPI. All MPI processes (including group leaders and workers) should call this function.

Parameters:

prob – ConstrainedProblem object defining the objective function, parameter space, and constraints.
mpi – A MpiPartition object, storing the information about how the pool of MPI processes is divided into worker groups.
grad – Whether to use a gradient-based optimization algorithm, as opposed to a gradient-free algorithm. If unspecified, a a gradient-free algorithm will be used by default. If you set grad=True finite-difference gradients will be used.
abs_step – Absolute step size for finite difference jac evaluation
rel_step – Relative step size for finite difference jac evaluation
diff_method – Differentiation strategy. Options are "centered" and "forward". If "centered", centered finite differences will be used. If "forward", one-sided finite differences will be used. For other values, an error is raised.
opt_method – Constrained solver to use: One of "SLSQP", "trust-constr", or "COBYLA". Use "COBYLA" for derivative-free optimization. See scipy.optimize.minimize for a description of the methods.
options –
dict, options keyword which is passed to scipy.optimize.minimize.

simsopt.solve.constrained_serial_solve(prob: ConstrainedProblem, grad: bool = None, abs_step: float = 1e-07, rel_step: float = 0.0, diff_method: str = 'forward', opt_method: str = 'SLSQP', options: dict = None)

Solve a constrained minimization problem using scipy.optimize, and without using any parallelization.

Parameters:

prob – ConstrainedProblem object defining the objective function, parameter space, and constraints.
grad – Whether to use a gradient-based optimization algorithm, as opposed to a gradient-free algorithm. If unspecified, a a gradient-free algorithm will be used by default. If you set grad=True for a problem, finite-difference gradients will be used.
abs_step – Absolute step size for finite difference jac evaluation
rel_step – Relative step size for finite difference jac evaluation
diff_method – Differentiation strategy. Options are "centered" and "forward". If "centered", centered finite differences will be used. If "forward", one-sided finite differences will be used. For other settings, an error is raised.
opt_method –
Constrained solver to use: One of "SLSQP", "trust-constr", or "COBYLA". Use "COBYLA" for derivative-free optimization. See scipy.optimize.minimize for a description of the methods.
options –
dict, options keyword which is passed to scipy.optimize.minimize.

simsopt.solve.get_gsco_iteration(iteration, res, wframe)

Returns the intermediate solution obtained by GSCO at a given iteration.

Parameters:

iteration (integer) – The iteration at which data are requested (0 corresponds to the initialization)
res (dictionary) – Dictionary returned by optimize_wireframe
wframe (wireframe class instance) – Wireframe whose segment currents were optimized

Returns:

x_iter – Solution at the requested iteration

Return type:

double array (1d column vector)

simsopt.solve.gsco_wireframe(wframe, A, c, lambda_S, no_crossing, match_current, default_current, max_current, max_iter, print_interval, no_new_coils=False, max_loop_count=0, x_init=None, loop_count_init=None, verbose=True)

Runs the Greedy Stellarator Coil Optimization algorithm to optimize the currents in a wireframe.

Parameters:

wframe (instance of the ToroidalWireframe class) – Wireframe whose segment currents are to be optimized. This function will update the currents instance variable to the solution. The optimizer will obey the constraints set in wframe.
Amat (contiguous 2d double array) – Inductance matrix relating normal field at test points on the plasma boundary to currents in each segment.
bvec (contiguous double array (column vector)) – Vector giving the target values of the normal field at each test point on the plasma boundary.
lambda_S (double) – Weighting factor for the objective function proportional to the number of active segments
no_crossing (boolean) – If true, the solution will forbid currents from crossing within the wireframe
match_current (boolean) – If true, added loops of current will match the current of the loop(s) adjacent to where they are added
default_current (double) – Loop current to add during each iteration to empty loops or to any loop if not matching existing current
max_current (double) – Maximum magnitude for the current in each segment
max_iter (integer) – Maximum number of iterations to perform
print_interval (integer) – Number of iterations between subsequent progress prints to screen
no_new_coils (boolean (optional)) – If true, the solution will forbid the addition of currents to loops where no segments already carry current (although it is still possible for an existing coil to be split into two coils); False by default
max_loop_count (integer (optional)) – If nonzero, sets the maximum number of current increments to add to a given loop in the wireframe (default is zero)
x_init (contiguous double array (optional)) – Initial current values to impose on the wireframe segments; will overwrite wframe.currents if used
loop_count_init (contiguous integer array (optional)) – Signed number of loops of current added to each loop in the wireframe prior to the optimization (optimization will add to these numbers); zero by default
verbose (boolean (optional)) – If true, will print progress to screen with durations of certain steps

Returns:

x (double array (1d column vector)) – Solution (currents in each segment of the wirframe)
loop_count (integer array (1d columnn vector)) – Signed number of current loops added to each loop in the wireframe
iter_hist (integer array) – Array with the iteration numbers of the data recorded in the history arrays. The first index is 0, corresponding to the initial guess x_init; the last is the final iteration.
curr_hist (double array) – Array with the signed loop current added at each iteration, taken to be zero for the initial guess (iteration zero)
loop_hist (integer array) – Array with the index of the loop to which current was added at each iteration, taken to be zero for the initial guess (iteration zero)
f_B_hist (1d double array (column vector)) – Array with values of the f_B objective function at each iteration
f_S_hist (1d double array (column vector)) – Array with values of the f_S objective function at each iteration
f_hist (1d double array (column vector)) – Array with values of the f_S objective function at each iteration
x_init (1d double array (column vector)) – Copy of the initial guess provided to the optimizer

simsopt.solve.least_squares_mpi_solve(prob: LeastSquaresProblem, mpi: MpiPartition, grad: bool = False, abs_step: float = 1e-07, rel_step: float = 0.0, diff_method: str = 'forward', save_residuals: bool = False, **kwargs)

Solve a nonlinear-least-squares minimization problem using MPI. All MPI processes (including group leaders and workers) should call this function.

This function solves min f(x) subject to lb <= x <= ub where the bounds are taken from the prob.bounds attribute.

Parameters:

prob – Optimizable object defining the objective function(s) and parameter space.
mpi – A MpiPartition object, storing the information about how the pool of MPI processes is divided into worker groups.
grad – Whether to use a gradient-based optimization algorithm, as opposed to a gradient-free algorithm. If unspecified, a a gradient-free algorithm will be used by default. If you set grad=True finite-difference gradients will be used.
abs_step – Absolute step size for finite difference jac evaluation
rel_step – Relative step size for finite difference jac evaluation
diff_method – Differentiation strategy. Options are “centered”, and “forward”. If centered, centered finite differences will be used. If forward, one-sided finite differences will be used. Else, error is raised.
save_residuals – Whether to save the residuals at each iteration. This may be useful for debugging, although the file can become large.
kwargs – Any arguments to pass to scipy.optimize.least_squares. For instance, you can supply max_nfev=100 to set the maximum number of function evaluations (not counting finite-difference gradient evaluations) to 100. Or, you can supply method to choose the optimization algorithm.

simsopt.solve.least_squares_serial_solve(prob: LeastSquaresProblem, grad: bool = None, abs_step: float = 1e-07, rel_step: float = 0.0, diff_method: str = 'forward', save_residuals: bool = False, **kwargs)

Solve a nonlinear-least-squares minimization problem using scipy.optimize, and without using any parallelization.

Parameters:

prob – LeastSquaresProblem object defining the objective function(s) and parameter space.
grad – Whether to use a gradient-based optimization algorithm, as opposed to a gradient-free algorithm. If unspecified, a a gradient-free algorithm will be used by default. If you set grad=True for a problem, finite-difference gradients will be used.
abs_step – Absolute step size for finite difference jac evaluation
rel_step – Relative step size for finite difference jac evaluation
diff_method – Differentiation strategy. Options are "centered", and "forward". If "centered", centered finite differences will be used. If "forward", one-sided finite differences will be used. Else, error is raised.
save_residuals – Whether to save the residuals at each iteration. This can be useful for debugging, although the file can become large.
kwargs –
Any arguments to pass to scipy.optimize.least_squares. For instance, you can supply max_nfev=100 to set the maximum number of function evaluations (not counting finite-difference gradient evaluations) to 100. Or, you can supply method to choose the optimization algorithm.

simsopt.solve.optimize_wireframe(wframe, algorithm, params, surf_plas=None, ext_field=None, area_weighted=True, bnorm_target=None, Amat=None, bvec=None, verbose=True)

Optimizes the segment currents in a wireframe class instance.

The currents field of the wireframe class instance will be updated with the solution.

There are two different modes of operation depending on which optional input parameters are set:

The user supplies a target plasma boundary and, if needed, a constant external magnetic field; in this case, the function calculates the normal field matrix and required normal field at the test points in preparation for performing the least-squares solve. In this mode, the parameter surf_plas must be supplied and, if relevant, ext_field.

The user supplies a pre-computed normal field matrix and vector with the required normal field on the plasma boundary, in which case it is not necessary to perform a field calculation before the least- squares solve. In this mode, the parameters Amat and bvec must be supplied.

IMPORTANT: for Regularized Constrained Least Squares (‘rcls’) optimizations, the parameter assume_no_crossings MUST be set to True if the wireframe is constrained to allow no crossing currents; otherwise, the optimization will not work properly!

Parameters:

wframe (instance of the ToroidalWireframe class) – Wireframe whose segment currents are to be optimized
algorithm (string) –
Optimization algorithm to use. Options are
- "rcls": Regularized Constrained Least Squares
- "gsco": Greedy Stellarator Coil Optimization
params (dictionary) – As specified in the lists below under Parameters for RCLS optimizations or Parameters for GSCO optimizations
surf_plas (Surface class instance (optional)) – Surface of the target plasma, on which test points are placed for evaluation of the normal field. If supplied, a magnetic field calculation will be performed to determine the normal field matrix and target normal field vector prior to performing the least- squares solve as described in mode (1) above.
ext_field (MagneticField class instance (optional)) – Constant external field assumed to be present in addition to the field produced by the wireframe. Used to calculate the target normal field vector in mode (1) as described above.
bnorm_target (double array (optional)) – Target value of the normal field on the plasma boundary to be produced by the combination of the wireframe and ext_field. Zero by default. Test points on the plasma boundary corresponding to the elements of bnorm_target must agree with the test points of surf_plas.
area_weighted (boolean (optional)) – Determines whether the normal field matrix and target normal field vector elements are weighted by the square root of the area ascribed to each test point on the plasma boundary. If true, the weightings will be applied and therefore the optimization will minimize the square integral of the normal field on the plasma boundary. If false, the weightings will NOT be applied, and the optimization will minimize the sum of squares of the normal field at each test point.
Amat (2d double array (optional)) – Inductance matrix relating normal field at test points on the plasma boundary to currents in each segment. This can be supplied along with bvec to skip the field calculation as describe in mode (2) above. Must have dimensions (n_test_points, wFrame.n_segments), where n_test_points is the number of test points on the plasma boundary.
bvec (double array (optional)) – Vector giving the target values of the normal field at each test point on the plasma boundary. This can be supplied along with Amat to skip the field calculation. Must have dimensions (n_test_points, 1).
verbose (boolean (optional)) – If true, will print progress to screen with durations of certain steps

Returns:

results – As specified below under Results dictionary

Return type:

dictionary

Parameters for RCLS optimizations:

reg_W: (scalar, 1d array, or 2d array) - Scalar, array, or matrix for regularization. If a 1d array, must have the same number of elements as wframe.n_segments. If a 2d array, both dimensions must be equal to wframe.n_segments.
assume_no_crossings: (boolean (optional)) - If true, will assume that the wireframe is constrained such that its free segments form single-track loops with no forks or crossings. Default is False.

Parameters for GSCO optimizations:

lambda_S: (double) - Weighting factor for the objective function proportional to the number of active segments
max_iter: (integer) - Maximum number of iterations to perform
print_inteval: (integer) - Number of iterations between subsequent prints of progress to screen
no_crossing: (boolean (optional)) - If true, the solution will forbid currents from crossing within the wireframe; default is false
match_current: (boolean (optional)) - If true, added loops of current will match the current of the loop(s) adjacent to where they are added; default is false
default_current: (double (optional)) - Loop current to add during each iteration to empty loops or to any loop if not matching existing current; default is 0
max_current: (double (optional)) - Maximum magnitude for the current in each segment; default is infinity
max_loop_count: (integer (optional)) - If nonzero, sets the maximum number of current increments to add to a given loop in the wireframe (default is zero)
x_init: (double array (optional)) - Initial current values to impose on the wireframe segments; will overwrite wframe.currents if used
loop_count_init: (integer array (optional)) - Signed number of loops of current added to each loop in the wireframe prior to the optimization (optimization will add to these numbers); zero by default

Results dictionary:

x: (1d double array) - Array of currents in each segment according to the solution. The elements of wframe.currents will be set to these values.
Amat: (2d double array) - Inductance matrix used for optimization, area weighted if requested
bvec: (1d double array (column vector)) - Target values of the normal field on the plasma boundary, area weighted if requested
wframe_field: (WireframeField class instance) - Magnetic field produced by the wireframe
f_B: (double) - Values of the sub-objective function f_B
f: (double) - Values of the total objective function f

For RCLS optimizations only:

f_R: (double) - Value of the sub-objective function f_R

For GSCO optimizations only:

loop_count: (integer array (1d columnn vector)) - Signed number of current loops added to each loop in the wireframe
iter_hist: (integer array) - Array with the iteration numbers of the data recorded in the history arrays. The first index is 0, corresponding to the initial guess x_init; the last is the final iteration.
curr_hist: (double array) - Array with the signed loop current added at each iteration, taken to be zero for the initial guess (iteration zero)
loop_hist: (integer array) - Array with the index of the loop to which current was added at each iteration, taken to be zero for the initial guess (iteration zero)
f_B_hist: (1d double array (column vector)) - Array with values of the f_B objective function at each iteration
f_S_hist: (1d double array (column vector)) - Array with values of the f_S objective function at each iteration
f_hist: (1d double array (column vector)) - Array with values of the f_S objective function at each iteration
x_init: (1d double array (column vector)) - Copy of the initial guess provided to the optimizer
f_S: (double) - Values of the sub-objective function f_S

simsopt.solve.rcls_wireframe(wframe, Amat, bvec, reg_W, assume_no_crossings, verbose)

Performs a Regularized Constrained Least Squares optimization for the segments in a wireframe.

Parameters:

wframe (instance of the ToroidalWireframe class) – Wireframe whose segment currents are to be optimized. This function will update the currents instance variable to the solution. The optimizer will obey the constraints set in wframe.
Amat (2d double array) – Inductance matrix relating normal field at test points on the plasma boundary to currents in each segment.
bvec (double array (column vector)) – Vector giving the target values of the normal field at each test point on the plasma boundary.
reg_W (scalar, 1d array, or 2d array) – Scalar, array, or matrix for regularization. If a 1d array, must have the same number of elements as wframe.n_segments. If a 2d array, both dimensions must be equal to wframe.n_segments.
assume_no_crossings (boolean (optional)) – If true, will assume that the wireframe is constrained such that its free segments form single-track loops with no forks or crossings. Default is False.
verbose (boolean (optional)) – If true, will print progress to screen with durations of certain steps

Returns:

x (double array (1d column vector)) – Solution (currents in each segment of the wirframe)
f_B, f_R, f (doubles) – Values of the partial objective functions (f_B for field accuracy, f_R for regularization) and the total objective function (f = f_B + f_R)

simsopt.solve.regularized_constrained_least_squares(A, b, W, C, d)

Solves a linear least squares problem with Tikhonov regularization subject to linear equality constraints on the variables.

In other words, minimizes:

0.5 * ((A * x - b)**2 + (W * x)**2

such that:

C * x = d

Here, A is the design matrix, b is the target vector, W is the regularization matrix (normally diagonal), and C and d contain the coefficients and constants of the constraint equations, respectively.

Parameters:

A (array with dimensions m*n) – Design matrix
b (array with dimension m) – Target vector
W (scalar, array with dimension n, or array with dimension n*n) – Regularization matrix
C (array with dimension p*n) – Coefficients of the solution vector elements in each of the p constraint equations. Must be full rank.
d (array with dimension p) – Constants appearing on the right-hand side of the constraint equations, i.e. B*x = d.

Returns:

x – Solution to the least-squares problem

Return type:

array with dimension n

simsopt.solve.relax_and_split(pm_opt, m0=None, **kwargs)

Uses a relax-and-split algorithm for solving the permanent magnet optimization problem, which solves a convex and nonconvex part separately.

Defaults to the MwPGP convex step and no nonconvex step. If a nonconvexity is specified, the associated prox function must be defined in this file. Relax-and-split allows for speedy algorithms for both steps and the imposition of convex equality and inequality constraints (including the required constraint on the strengths of the dipole moments).

Parameters:

pm_opt – The grid of permanent magnets to optimize.
m0 – Initial guess for the permanent magnet dipole moments. Defaults to a starting guess of all zeros. This vector must lie in the hypersurface spanned by the L2 ball constraints. Note that if algorithm is being used properly, the end result should be independent of the choice of initial condition.
kwargs –
Keyword arguments to pass to the algorithm. The following arguments can be passed to the MwPGP algorithm:

epsilon:
Error tolerance for the convex part of the algorithm (MwPGP).

nu:
Hyperparameter used for the relax-and-split least-squares. Set nu >> 1 to reduce the importance of nonconvexity in the problem.

reg_l0:
Regularization value for the L0 nonconvex term in the optimization. This value is automatically scaled based on the max dipole moment values, so that reg_l0 = 1 corresponds to reg_l0 = np.max(m_maxima). It follows that users should choose reg_l0 in [0, 1].

reg_l1:
Regularization value for the L1 nonsmooth term in the optimization,

reg_l2:
Regularization value for any convex regularizers in the optimization problem, such as the often-used L2 norm.

max_iter_MwPGP:
Maximum iterations to perform during a run of the convex part of the relax-and-split algorithm (MwPGP).

max_iter_RS:
Maximum iterations to perform of the overall relax-and-split algorithm. Therefore, also the number of times that MwPGP is called, and the number of times a prox is computed.

verbose:
Prints out all the loss term errors separately.

Returns:

A tuple of optimization loss, solution at each step, and sparse solution.

The tuple contains

errors:: Total optimization loss after each convex sub-problem is solved.
m_history:: Solution for the permanent magnets after each convex sub-problem is solved.
m_proxy_history:: Sparse solution for the permanent magnets after each convex sub-problem is solved.

simsopt.solve.serial_solve(prob: Optimizable | Callable, grad: bool = None, abs_step: float = 1e-07, rel_step: float = 0.0, diff_method: str = 'centered', **kwargs)

Solve a general minimization problem (i.e. one that need not be of least-squares form) using scipy.optimize.minimize, and without using any parallelization.

Parameters:

prob – Optimizable object defining the objective function(s) and parameter space.
grad – Whether to use a gradient-based optimization algorithm, as opposed to a gradient-free algorithm. If unspecified, a gradient-based algorithm will be used if prob has gradient information available, otherwise a gradient-free algorithm will be used by default. If you set grad=True in which gradient information is not available, finite-difference gradients will be used.
abs_step – Absolute step size for finite difference jac evaluation
rel_step – Relative step size for finite difference jac evaluation
diff_method – Differentiation strategy. Options are "centered", and "forward". If "centered", centered finite differences will be used. If "forward", one-sided finite differences will be used. Else, error is raised.
kwargs –
Any arguments to pass to scipy.optimize.least_squares. For instance, you can supply max_nfev=100 to set the maximum number of function evaluations (not counting finite-difference gradient evaluations) to 100. Or, you can supply method to choose the optimization algorithm.