CoronelBuendia
commented
1 year ago The previous attempt at solving this issue using optimisation was found to work well in some cases, but not be robust, also producing different results on different machines (possibly due to different RNG in COBYLA's internal random seed - but not sure).
I reproduce the key elements of that approach here for future reference:
def f_psi_norm(x):
return psi_norm_func(x) - psi_norm
search_range = mesh.hmax()
def find_extremum(
func: Callable[[np.ndarray], np.ndarray], x0: Iterable[float]
) -> np.ndarray:
"""
Extremum finding using constrained optimisation
"""
lower_bounds = x0 - search_range
upper_bounds = x0 + search_range
# NOTE: COBYLA appears to do a better job here, as it seems that the
# NLOpt implementation of SLSQP really requires a feasible starting
# solution, which is not so readily available with this tight equality
# constraint. The scipy SLSQP implementation apparently does not require
# such a good starting solution. Neither SLSQP nor COBYLA can guarantee
# convergence without a feasible starting point.
optimiser = Optimiser(
"COBYLA", 2, opt_conditions={"ftol_abs": 1e-10, "max_eval": 1000}
)
optimiser.set_objective_function(func)
optimiser.set_lower_bounds(lower_bounds)
optimiser.set_upper_bounds(upper_bounds)
f_constraint = OptimisationConstraint(
_f_constrain_psi_norm,
f_constraint_args={
"psi_norm_func": f_psi_norm,
"lower_bounds": lower_bounds,
"upper_bounds": upper_bounds,
},
constraint_type="equality",
)
optimiser.add_eq_constraints(f_constraint, tolerance=1e-10)
try:
x_star = optimiser.optimise(x0)
except ExternalOptError as e:
bluemira_warn(
f"calculate_plasma_shape_params::find_extremum failing at {x0}, defaulting to mesh value: {e}"
)
x_star = x0
return x_star
pi_opt = find_extremum(_f_min_radius, pi)
pl_opt = find_extremum(_f_min_vert, pl)
po_opt = find_extremum(_f_max_radius, po)
pu_opt = find_extremum(_f_max_vert, pu)with supporting functions:
def _f_max_radius(x, grad):
result = -x[0]
if grad.size > 0:
grad[0] = -1.0
grad[1] = 0.0
return result
def _f_min_radius(x, grad):
result = x[0]
if grad.size > 0:
grad[0] = 1.0
grad[1] = 0.0
return result
def _f_max_vert(x, grad):
result = -x[1]
if grad.size > 0:
grad[0] = 0.0
grad[1] = -1.0
return result
def _f_min_vert(x, grad):
result = x[1]
if grad.size > 0:
grad[0] = 0.0
grad[1] = 1.0
return result
def _f_constrain_psi_norm(
constraint: np.ndarray,
x: np.ndarray,
grad: np.ndarray,
psi_norm_func=None,
lower_bounds=None,
upper_bounds=None,
) -> np.ndarray:
"""
Constraint function for points on the psi_norm surface.
"""
result = psi_norm_func(x)
constraint[:] = result
if grad.size > 0:
grad[:] = approx_derivative(
psi_norm_func,
x,
f0=result,
bounds=[lower_bounds, upper_bounds],
method="3-point",
)
return np.array([result])
Description of issue / requirement to address
When finding the plasma shape parameters in
calculate_plasma_shape_paramswe presently just use the extrema of the flux surface as recovered from the tricontour operation on the mesh and the associated psi field.This is known to be inaccurate and sensitive to mesh size. There is also the very real possibility for oscillation, particularly in triangularity.
Proposed solution
In the past we have tried two different approaches here.
1) Refine the mesh and recalculate psi on the new mesh and return the extrema from the tricontours on that psi field 2) Use constrained optimisation to find the extrema. The latter approach was in fact used for some time before found wanting for a number of reasons (see implementation and reasons in the comment below).
Ideally, the constrained optimisation approach should be revisited to be more stable, probably not using COBYLA but SLSQP, as this is the most "correct" approach mathematically, and is not sensitive to mesh size (within reason). Failing that we could re-implement the refined mesh search as originally implemented.
Alternative solutions
Additional context
As part of a push towards including the procedure and associated module reshuffle in the next release, this issue has been separated out and will be treated in future releases.
@ivanmaione for info