freegsnke.inverse module

Implements the optimiser for the inverse Grad-Shafranov problem.

Copyright 2025 UKAEA, UKRI-STFC, and The Authors, as per the COPYRIGHT and README files.

This file is part of FreeGSNKE.

FreeGSNKE is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

FreeGSNKE is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

You should have received a copy of the GNU Lesser General Public License along with FreeGSNKE. If not, see <http://www.gnu.org/licenses/>.

class freegsnke.inverse.Inverse_optimizer(isoflux_set=None, null_points=None, psi_vals=None, coil_current_limits=None, psi_norm_limits=None, *, weight_isoflux=1.0, weight_nulls=1.0, weight_psi=1.0, mu_coils=100000.0, mu_psi_norm=1000000.0)[source]

Bases: object

This class implements a gradient based optimiser for the coil currents, used to perform (static) inverse GS solves.

__init__(isoflux_set=None, null_points=None, psi_vals=None, coil_current_limits=None, psi_norm_limits=None, *, weight_isoflux=1.0, weight_nulls=1.0, weight_psi=1.0, mu_coils=100000.0, mu_psi_norm=1000000.0)[source]

Initialise magnetic constraint definitions for inverse equilibrium optimisation.

This object stores magnetic configuration constraints that are enforced during inverse Grad-Shafranov optimisation.

Constraints supported include:

  • Isoflux constraints

  • Magnetic axis (O-point) and X-point constraints

  • Direct flux value constraints

  • Coil current bound constraints

  • Normalised flux inequality constraints

These constraints are used to construct objective and penalty terms during nonlinear current optimisation.

Parameters:

  • isoflux_set (list or ndarray, optional) –

    Collection of isoflux constraint objects.

    Each isoflux constraint is specified as:

    [Rcoords, Zcoords, weights]

    where:

    Rcoords : 1D array of radial coordinates Zcoords : 1D array of vertical coordinates weights : (optional, array of 1’s if not provided)1D array of weights to increase/decrease

    the influence of the constraint on the solution.

    All specified points within each set are required to share the same poloidal flux value.

  • null_points (list or ndarray, optional) –

    Magnetic null point constraints.

    Structure:

    [Rcoords, Zcoords]

    Specifies coordinates of:

    • X-points

    • O-points (magnetic axes)

  • psi_vals (list or ndarray, optional) –

    Direct flux value constraints.

    Structure:

    [Rcoords, Zcoords, psi_values]

    where:

    Rcoords, Zcoords, psi_values must have identical shapes.

    Used to enforce ψ(R,Z) = ψ_target at specified locations.

  • coil_current_limits (list, optional) –

    Hard inequality bounds on coil currents.

    Structure:

    [upper_limits, lower_limits]

    Each entry is a list with length equal to the number of controllable coils.

    Example:
    [

    [Imax1, Imax2, …], [Imin1, Imin2, …]

    ]

    Use None to indicate no bound.

  • psi_norm_limits (list or ndarray, optional) –

    Normalised flux inequality constraints.

    Structure:

    [Rcoord, Zcoord, normalised_psi_value, constraint_sign]

    Constraint form:

    If constraint_sign = 1:

    ψ_norm ≥ ψ_target

    If constraint_sign = -1:

    ψ_norm ≤ ψ_target

    Normalised flux is defined:

    ψ_norm = (ψ - ψ_axis) / (ψ_boundary - ψ_axis)

  • weight_isoflux (float) – The weight of the isoflux constraints in the least-squares optimisation problem (default = 1.0).

  • weight_nulls (float) – The weight of the null point (X-point) constraints in the least-squares optimisation problem (default = 1.0).

  • weight_psi (float) – The weight of the psi value constraints in the least-squares optimisation problem (default = 1.0).

  • mu_coils (float) – A penalty factor applied to violation of the coil current limits (default = 1e5).

  • mu_psi_norm (float) – A penalty factor applied to violation of the normalised psi limits (default = 1e5).

Notes

Increasing the weights/penalty factors causes the least-squares optimisation to proritise satisfying the higher-weighted/penaltied constraints.

build_control_coils(eq)[source]

Identify and cache coil systems available for active control.

This method constructs internal indexing and masking structures that separate:

  • Control coils (actively optimised)

  • Passive coils (fixed currents or structures)

These mappings are required for:

  • Current optimisation

  • Jacobian construction

  • Constraint enforcement

  • Vectorised current updates

Parameters:

  • eq (freegsnke equilibrium object) – Provides tokamak coil system information.

  • attributes (Expected coil) –

    eq.tokamak.coils : list of coil objects Each coil must have:

    • coil.control : bool indicating controllability

Returns:

Updates solver internal coil control state.

Return type:

None

build_control_currents(eq)[source]

Extract coil current values for controllable coils from the equilibrium object.

This method builds a reduced current vector containing only the currents of actively optimised control coils.

Parameters:

eq (freegsnke equilibrium object) –

Equilibrium object providing current coil state via:

eq.tokamak.getCurrentsVec()

Returns:

Stores result in:

self.control_currents

Return type:

None

build_control_currents_Vec(full_currents_vec)[source]

Extract controllable coil currents from a full coil current vector.

This function performs projection:

I_control = M_control I_full

where M_control is a boolean masking operator selecting controllable coils.

Parameters:

full_currents_vec (ndarray) –

Full vector of coil currents.

Example:

eq.tokamak.getCurrentsVec()

Returns:

Stores result in:

self.control_currents

Return type:

None

build_full_current_vec(eq)[source]

Build full coil current vector from equilibrium object.

This retrieves all coil currents, including:

  • Control coils

  • Passive structure coils

Parameters:

eq (freegsnke equilibrium object) –

Provides coil current state via:

eq.tokamak.getCurrentsVec()

Returns:

Stores result in:

self.full_currents_vec

Return type:

None

build_greens(eq)[source]

Construct and cache magnetic Green’s function response operators.

This method precomputes magnetic response matrices used for:

  • Flux matching constraints

  • Null point field constraints

  • Isoflux surface optimisation

  • Flux inequality constraints

Green’s functions represent linear mappings between:

Coil current perturbations → Magnetic field / flux perturbations

These operators are used during:

  • Constraint optimisation

  • Jacobian approximation

  • Inverse equilibrium solving

Parameters:

eq (freegsnke equilibrium object) – Provides access to tokamak geometry and magnetic response calculation routines.

Returns:

Stores Green’s function operators internally.

Return type:

None

build_isoflux_lsq(full_currents_vec)[source]

Construct linear least-squares system for enforcing isoflux magnetic constraints.

This method assembles the linear optimisation problem:

A I_control ≈ b

where:

A = Green’s function response matrix I_control = vector of controllable coil currents b = plasma + vacuum flux vector

The least-squares system is derived from:

ψ_total = ψ_tokamak + ψ_plasma

and enforcing:

ψ_i − ψ_j ≈ constant (equal-flux constraints along magnetic surfaces)

Parameters:

full_currents_vec (ndarray) –

Full vector of all coil currents.

Example:

eq.tokamak.getCurrentsVec()

Returns:

  • A (list of ndarray) – List of arrays for each isoflux constraint set.

  • b (list of ndarray) – Right-hand side residual vectors.

  • loss (list of float) – Constraint violation magnitudes (L2 norms).

  • Mathematical formulation

  • ————————

  • For each isoflux constraint set – A_i = (dG_i)ᵀ P_control

    b_i = − ( Σ_k dG_i I_k + Δψ_plasma )

  • where – dG_i = pairwise Green’s function flux differences P_control = coil control projection operator

build_lsq(full_currents_vec)[source]

Assemble the global least-squares optimisation system combining all active magnetic and control constraints.

This method aggregates all available constraint types into a single linear optimisation problem of the form:

A I_control ≈ b

where:

A = stacked Jacobian / response matrices b = stacked residual constraint vectors

The constraint types that can be combined include:

  • Isoflux surface constraints

  • Magnetic null-point field constraints

  • Direct flux value constraints

  • Direct coil current constraints

Parameters:

full_currents_vec (ndarray) –

Full vector of coil currents.

Example:

eq.tokamak.getCurrentsVec()

Returns:

Stores assembled optimisation system internally:

self.A self.b self.loss

Return type:

None

Notes

This formulation solves the global constrained optimisation problem:

min_I || A I − b ||²

where each constraint type contributes a block to the system.

Constraint dimensional bookkeeping is stored as:

self.isoflux_dim self.nullp_dim self.psiv_dim self.curr_dim

build_null_points_lsq(full_currents_vec)[source]

Construct a least-squares system enforcing magnetic null-point constraints.

This method builds the linear optimisation system:

A I_control ≈ b

for magnetic field cancellation at prescribed null points.

Null point constraints enforce:

B_R(R_i, Z_i) ≈ 0 B_Z(R_i, Z_i) ≈ 0

where:

B_R = radial magnetic field component B_Z = vertical magnetic field component

These fields are expressed using Green’s function response matrices:

B(R,Z) = Σ_k G_k(R,Z) I_k + B_plasma

Parameters:

full_currents_vec (ndarray) –

Full vector of coil currents.

Example:

eq.tokamak.getCurrentsVec()

Returns:

  • A (ndarray) – Combined Jacobian matrix for null-point field constraints.

  • b (ndarray) – Residual field mismatch vector (negated for optimisation solving).

  • loss (list of float) –

    Constraint violation magnitudes:

    [ ||B_R||₂ , ||B_Z||₂ ]

Notes

This formulation is used to stabilise:

  • Magnetic axis positioning

  • X-point control

  • Divertor topology shaping

Mathematical formulation

The optimisation problem solved is:

min_I || G I + B_plasma ||²

where:

G = magnetic Green’s response operator I = coil current vector

build_plasma_isoflux_lsq(full_currents_vec, trial_plasma_psi)[source]

Assemble the least-squares system for plasma-only isoflux optimisation.

This method constructs a linearised least-squares problem of the form:

A_plasma x ≈ b_plasma

where x ∈ ℝ² represents plasma transformation parameters (e.g. normalisation and nonlinear shaping terms).

The residual enforces pairwise isoflux constraints:

Δψ_total = Δψ_coils + Δψ_plasma ≈ 0

using both raw flux differences and nonlinear transformed (ψ̂ log ψ̂) contributions.

Parameters:

  • full_currents_vec (ndarray) – Full vector of coil currents.

  • trial_plasma_psi (ndarray) – Plasma flux array defined on the equilibrium grid.

Notes

The constructed system solves for two plasma parameters:

x[0] → linear normalisation scaling x[1] → nonlinear shaping contribution

Results are stored internally:

self.A_plasma self.b_plasma self.loss_plasma

build_plasma_vals(trial_plasma_psi)[source]

Compute and cache plasma-dependent magnetic quantities from a candidate flux solution.

This method evaluates the plasma flux field and derives quantities required for constraint optimisation and nonlinear solve diagnostics, including:

  • Magnetic field components at null points

  • Isoflux surface flux differences

  • Pointwise flux constraint values

The plasma flux field is interpolated using bicubic spline interpolation.

Parameters:

trial_plasma_psi (ndarray) –

Candidate plasma poloidal flux solution.

Must have shape compatible with:

len(self.eqR) × len(self.eqZ)

Returns:

Updates internal solver state variables:

self.brp self.bzp self.d_psi_plasma_vals_iso self.psi_plasma_vals

Return type:

None

build_psi_vals_lsq(full_currents_vec)[source]

Construct a least-squares system enforcing direct flux value constraints.

This method builds the optimisation system:

A I_control ≈ b

where:

ψ_model(R_i, Z_i) = ψ_target(R_i, Z_i)

is enforced by matching magnetic flux values at specified locations.

The optimisation residual is defined as:

b = ψ_tokamak(I) + ψ_plasma

  • ψ_target

  • ⟨b⟩

Mean flux removal is applied to remove arbitrary vertical flux offsets, since the Grad–Shafranov equation is invariant under constant flux shifts.

Parameters:

full_currents_vec (ndarray) –

Full coil current vector.

Example:

eq.tokamak.getCurrentsVec()

Returns:

  • A (ndarray) – Jacobian matrix mapping coil current perturbations → flux changes.

  • b (ndarray) – Flux mismatch residual vector.

  • normalised_loss (list of float) – Normalised constraint violation magnitude.

Notes

This constraint formulation is commonly used for:

  • Magnetic axis pinning

  • Boundary flux matching

  • Profile shape control

Mathematical formulation

Solve:

min_I || G I + ψ_plasma − ψ_target ||²

optimize_currents(eq, profiles, full_currents_vec, trial_plasma_psi, l2_reg)[source]

Solve the constrained least-squares optimisation problem for coil current updates.

This method computes optimal coil current corrections by solving:

min_I || A I − b ||² + λ || I ||²

where:

A = combined constraint Jacobian matrix b = combined constraint residual vector λ = Tikhonov (L2) regularisation parameter

The optimisation accounts for:

  • Magnetic topology constraints

  • Flux surface matching

  • Null point field cancellation

  • Hardware current constraints

Parameters:

  • eq (FreeGSNKE equilibrium object) – Equilibrium providing geometry and tokamak configuration.

  • profiles (FreeGSNKE profile object) – Plasma profile properties used for constraint evaluation.

  • full_currents_vec (ndarray) –

    Full vector of all coil currents.

    Example:

    eq.tokamak.getCurrentsVec()

  • trial_plasma_psi (ndarray) –

    Candidate plasma flux solution.

    Must have same shape as:

    eq.R

  • l2_reg (float or ndarray) –

    Tikhonov regularisation parameter.

    If float:

    λ I² penalty is applied uniformly.

    If array:

    Allows coil-wise regularisation weighting.

Returns:

  • delta_current (ndarray) – Optimal coil current update vector.

  • loss (float) – Residual loss norm.

optimize_currents_grad(full_currents_vec, trial_plasma_psi, isoflux_weight=1.0, null_points_weight=1.0, psi_vals_weight=1.0)[source]

Compute the gradient of the magnetic least-squares objective with respect to control coil currents.

This method evaluates the gradient of the unconstrained objective:

J(ΔI) = ½ || A ΔI − b ||²

evaluated at ΔI = 0, i.e.

∇J = Aᵀ b

Optional weighting factors allow different constraint classes (isoflux, null points, direct flux values) to be scaled prior to gradient evaluation.

Parameters:

  • full_currents_vec (ndarray) –

    Full vector of all coil current values.

    Example:

    eq.tokamak.getCurrentsVec()

  • trial_plasma_psi (ndarray) – Plasma flux contribution. Must have same shape as eq.R.

  • isoflux_weight (float, optional) – Weight applied to isoflux constraint residuals.

  • null_points_weight (float, optional) – Weight applied to null-point field constraint residuals.

  • psi_vals_weight (float, optional) – Weight applied to direct flux value constraint residuals.

Returns:

  • grad (ndarray) – Gradient of the weighted least-squares objective with respect to control coil current updates.

  • loss (float) – Combined magnetic constraint loss (unweighted).

optimize_currents_quadratic(eq, profiles, full_currents_vec, reg_matrix, *, mu_coils=None, mu_psi_norm=None, A=None, b=None)[source]

Solve the regularised constrained least-squares problem using convex optimisation.

This method computes coil current updates by solving the quadratic program:

minimise_ΔI

|| A ΔI − b ||²

  • ΔIᵀ R ΔI

  • penalty(slack variables)

subject to optional inequality constraints:

  • Coil current upper/lower bounds

  • Normalised flux (ψ_norm) inequality constraints

Slack variables are introduced to allow temporary constraint violation, penalised in the objective function. This improves optimisation robustness and prevents infeasibility during nonlinear solve iterations.

Parameters:

  • eq (FreeGSNKE equilibrium object) – Equilibrium object used to evaluate plasma quantities and ψ_norm.

  • profiles (FreeGSNKE profile object) – Provides boundary flux (psi_bndry) and related quantities.

  • full_currents_vec (ndarray) – Full vector of coil currents.

  • reg_matrix (ndarray) – Regularisation matrix (n_control_coils × n_control_coils), typically diagonal, encoding Tikhonov L2 regularisation.

  • mu_coils (float | None) – Scaling factor for coil current limit slack penalties. If None, defaults to self.mu_coils (default ~1e5).

  • mu_psi_norm (float | None) – Scaling factor for ψ_norm slack penalties. If None, defaults to self.mu_psi_norm.

  • A (ndarray | None) – Sensitivity matrix. If None, uses self.A.

  • b (ndarray | None) – Target residual vector. If None, uses self.b.

Returns:

  • delta (ndarray) – Optimal coil current update vector ΔI.

  • loss (float) –

    Combined loss including:

    • Magnetic constraint residual norm

    • Slack variable penalties

Notes

The optimisation problem solved is a convex quadratic program (QP).

Base objective:

min_ΔI ||A ΔI − b||² + ΔIᵀ R ΔI

With optional inequality constraints:

I_min ≤ I_current + ΔI ≤ I_max ψ_norm constraints (≥ or ≤ type)

Slack variables s ≥ 0 allow:

constraint_violation ≤ s

with large quadratic penalties to discourage violation.

optimize_plasma_psi(full_currents_vec, trial_plasma_psi, l2_reg)[source]

Solve the regularised least-squares problem for plasma parameters.

This solves:

min_x || A_plasma x − b_plasma ||² + xᵀ R x

where:

x ∈ ℝ² R = Tikhonov regularisation matrix

Parameters:

  • full_currents_vec (ndarray) – Full vector of coil currents.

  • trial_plasma_psi (ndarray) – Plasma flux array defined on the equilibrium grid.

  • l2_reg (float or ndarray (length 2)) –

    Tikhonov regularisation strength.

    If float:

    Uniform regularisation applied.

    If array:

    Diagonal regularisation weights.

Returns:

  • delta_current (ndarray) – Optimal update vector for plasma parameters (length 2).

  • loss_plasma (float) – Isoflux constraint violation norm (pre-update).

Notes

The solution is computed via normal equations:

x = (AᵀA + R)⁻¹ Aᵀ b

Since the system dimension is small (2×2), direct inversion is computationally inexpensive.

plot(axis=None, show=True)[source]

Visualise the active coil control constraints.

This method provides a graphical representation of the magnetic and operational constraints currently configured in the object.

Parameters:

  • axis (matplotlib.axes.Axes | None, optional) – Axis object on which to draw the plot. If None, a new figure and axis are created internally.

  • show (bool, optional) – If True, calls matplotlib.pyplot.show() before returning.

Returns:

Axis containing the generated plot.

Return type:

matplotlib.axes.Axes

Notes

This is a thin wrapper around:

freegs4e.plotting.plotIOConstraints

and exists primarily for convenience and API consistency.

prepare_for_plasma_optimization(eq)[source]

Prepare geometry- and Green’s-function-dependent quantities required for plasma optimisation.

This method ensures that:

  • Source-domain geometric properties are updated

  • Magnetic Green’s operators are constructed

Parameters:

eq (FreeGSNKE equilibrium object) – Provides geometry, coil configuration, and grid data.

Notes

This is typically called prior to plasma-only optimisation steps to ensure response operators and geometric mappings are consistent with the current equilibrium configuration.

prepare_for_solve(eq)[source]

Prepare constraint solver quantities after object instantiation.

This method must be called before performing optimisation or nonlinear solve steps. It builds:

  • Control coil selection masks

  • Green’s function response operators

These quantities are required for computing:

  • Constraint residuals

  • Gradient/Jacobian approximations

  • Loss function evaluations during optimisation

Parameters:

eq (freegsnke equilibrium object) –

Equilibrium object providing:

  • Coil geometry

  • Coil current state

  • Magnetic Green’s function operators

Returns:

Updates internal constraint solver state.

Return type:

None

prepare_plasma_psi(trial_plasma_psi)[source]

Preprocess plasma flux values for normalisation and constraint evaluation.

This method computes shifted plasma flux extrema used for normalised flux calculations and ψ_norm-based constraints.

Parameters:

trial_plasma_psi (ndarray) – Plasma flux array.

prepare_plasma_vals_for_plasma(trial_plasma_psi)[source]

Precompute plasma-dependent quantities required for plasma optimisation.

This method evaluates the plasma flux at constraint locations and constructs pairwise flux-difference terms used in isoflux constraints.

In addition to raw flux differences, a nonlinear transformed version of the normalised flux is computed:

ψ̂ = (ψ − ψ_min) / ψ0 f(ψ̂) = ψ̂ log(ψ̂)

which is then used to build weighted pairwise differences.

Parameters:

trial_plasma_psi (ndarray) – Plasma flux array defined on the equilibrium grid (eqR, eqZ).

Notes

For each isoflux constraint set:

  1. Interpolate ψ at constraint coordinates.

  2. Construct pairwise differences:

    ψ_i − ψ_j

  3. Construct transformed differences:

    ψ0 [ f(ψ̂_i) − f(ψ̂_j) ]

These quantities are cached for use in plasma optimisation routines.

The logarithmic transform enhances sensitivity near ψ̂ → 0 and improves conditioning when enforcing plasma-based constraints.

rebuild_full_current_vec(control_currents, filling=0)[source]

Reconstruct a full coil current vector from control coil values.

This performs the inverse mapping:

I_full = P I_control + I_background

where:
P = projection operator embedding control currents

into full coil system ordering.

Non-controlled coils are assigned the value specified by:

filling

Parameters:

  • control_currents (ndarray) – Current values for actively controlled coils.

  • filling (float, optional) – Default value used to fill non-controlled coil entries.

Returns:

full_current_vec – Complete coil current vector in system ordering.

Return type:

ndarray

source_domain_properties(eq)[source]

Source and cache computational domain geometry from the equilibrium object.

This method extracts and stores the spatial grid coordinates used for solving the Grad–Shafranov equation.

Parameters:

eq (FreeGSNKE equilibrium object) –

Equilibrium object providing computational domain information:

  • eq.R : 2D radial grid coordinates

  • eq.Z : 2D vertical grid coordinates

Returns:

Stores domain geometry internally.

Return type:

None