FPT Topology Control (pyna.control)

Note

pyna.control contains the generic FPT topology control framework, applicable to any area-preserving dynamical system. For toroidal-specific control helpers (gap response, island control, q-profile), use pyna.toroidal.control.

Magnetic topology state vector �observable quantities for control.

This module defines the complete set of observables that characterize

magnetic topology for both tokamaks (axisymmetric) and stellarators (3D).

Observable categories

---

Boundary : g_i (plasma-wall gap), X/O-point positions, DPm eigenvalues,

DPm eigenvectors, B_pol on LCFS

Core : q-profile, J-iota profile, O-point rotation transform

Derived : connection length L_c, flux expansion f_x, surface fate

class pyna.control.topology_state.OPointState(R, Z, A_matrix, DPm, DPm_eigenvalues, iota=0.0, q=inf)

Bases: object

State of a magnetic O-point (axis or island centre).

A_matrix: ndarray

DPm: ndarray

DPm_eigenvalues: ndarray

R: float

Z: float

iota: float = 0.0

property is_elliptic: bool

q: float = inf

class pyna.control.topology_state.SurfaceFate(value)

Bases: Enum

Classification of a flux surface under perturbation.

CHAOTIC = 4

DEFORMED = 2

INTACT = 1

ISLAND = 3

UNKNOWN = 5

class pyna.control.topology_state.TopologyState(xpoints=<factory>, opoints=<factory>, gap_gi=<factory>, Bpol_lcfs=None, theta_lcfs=None, q_samples=None, s_samples=None, J_iota_samples=None, surface_fate=<factory>, is_axisymmetric=True, phi_ref=0.0)

Bases: object

Complete magnetic topology state vector.

Contains all observables needed for multi-objective magnetic topology

control. Works for both axisymmetric (tokamak) and 3D (stellarator)

configurations.

Attributes

---

xpoints : list of XPointState

All X-points (divertor, saddle, �.

opoints : list of OPointState

All O-points (magnetic axis + island O-points).

gap_gi : dict

Plasma-wall gaps g_i at named monitoring points {name: gap_m}.

Bpol_lcfs : ndarray or None

|B_pol| on LCFS at theta_lcfs grid points.

q_samples : ndarray or None

q(s) at s = s_samples.

J_iota_samples : ndarray or None

Current-density rotation transform at the same s grid.

surface_fate : dict

{psi_norm: SurfaceFate} for monitored flux surfaces.

is_axisymmetric : bool

If True use closed-form FPT; otherwise integrate along φ.

phi_ref : float

Reference toroidal angle for the Poincaré section (default 0).

Bpol_lcfs: ndarray | None = None

J_iota_samples: ndarray | None = None

gap_gi: Dict[str, float]

is_axisymmetric: bool = True

opoints: List[OPointState]

phi_ref: float = 0.0

q_samples: ndarray | None = None

s_samples: ndarray | None = None

surface_fate: Dict[float, SurfaceFate]

theta_lcfs: ndarray | None = None

to_vector(keys=None)

Flatten state into a 1-D numpy array for optimisation.

Returns

---

values : ndarray, shape (n_obs,)

labels : list of str, length n_obs

xpoints: List[XPointState]

class pyna.control.topology_state.XPointState(R, Z, A_matrix, DPm, DPm_eigenvalues, DPm_eigenvectors, D2Pm=None, D3Pm=None, Greene_residue=0.0)

Bases: object

State of a divertor X-point.

A_matrix: ndarray

D2Pm: ndarray | None = None

D3Pm: ndarray | None = None

DPm: ndarray

DPm_eigenvalues: ndarray

DPm_eigenvectors: ndarray

Greene_residue: float = 0.0

R: float

Z: float

property connection_length_scale: float

Inverse of max log-eigenvalue �proportional to L_c near X-point.

When eigenvalues are close to 1 the separatrix is nearly degenerate

and nearby field lines have a long connection length (good for

detachment). When they diverge from 1 the connection length drops.

property is_hyperbolic: bool

property stability_index: float

Half-trace of DPm.

pyna.control.topology_state.compute_topology_state(field_func, xpoint_guesses, opoint_guesses, is_axisymmetric=True, phi_ref=0.0, eps=0.0001)

Compute full topology state from a field function.

Parameters

---

field_func : callable

Magnetic field function: [R,Z,phi] �[dR/dl, dZ/dl, dphi/dl].

xpoint_guesses : list of (R, Z)

Initial guesses for X-point positions.

opoint_guesses : list of (R, Z)

Initial guesses for O-point positions.

is_axisymmetric : bool

Use closed-form exp(2πA) for DPm.

phi_ref : float

Reference toroidal angle.

eps : float

Finite-difference step for A-matrix.

Returns

---

TopologyState

Return type:

TopologyState

Response matrix ∂(topology_state)/∂(control_inputs).

R_ij = ∂(observable_i) / ∂(control_j)

For axisymmetric tokamaks:

∂(x_cyc)/∂(I_coil_k) = -A⁻¹ · ∂g(x_cyc)/∂(I_coil_k) ∂(DPm_eigval)/∂(I_coil_k) = from δDPm formula

For 3D/stellarator: same structure but needs φ-integration.

The response matrix enables the linear control problem:

δ_state ≈ R · δ_controls

hence:

min Σ_i w_i |state_i + R_ij δI_j − target_i|² s.t. constraints

pyna.control.response_matrix.build_full_response_matrix(base_field_func, coil_field_funcs, state, wall=None, field_func_key='default', observables=None, eps_current=1.0, equilibrium=None, plasma_response=False, R_grid=None, Z_grid=None)

Full response matrix combining all observable categories.

Combines: - X/O-point shifts and DPm eigenvalue changes (FPT closed-form) - Plasma-wall gap responses (FPT manifold shift, if wall is supplied) - q-profile (placeholder zeros, to be filled by pyna-qprofile-response)

Parameters:

• base_field_func (callable) – Base equilibrium field function.

• coil_field_funcs (list of callable) – Per-unit-current coil field functions.

• state (TopologyState) – Current topology state with at least one X-point.

• wall (WallGeometry or None) – First wall geometry. If None, gap rows are zero placeholders.

• field_func_key (str) – Hashable key identifying the equilibrium (for caching manifold growth).

• observables (list of str or None) – Reserved for future filtering.

• eps_current (float) – Current perturbation for per-unit normalisation.

• equilibrium (EquilibriumSolovev or compatible object, optional) – If provided and plasma_response=True, the plasma response δB_plasma is added to each coil’s δB_ext before computing topology/gap observables.

• plasma_response (bool) –

  Whether to include plasma response. - True (recommended for core quantities like q-profile):

  δB_total = δB_ext + solve_perturbed_gs(B0, J0, p0, δB_ext)

  • False (OK for edge quantities like gap_gi when coils are far):

    δB_total = δB_ext (vacuum approximation)

• R_grid (array-like or None) – Grid for plasma response computation. If None, a default grid around the equilibrium is constructed automatically.

• Z_grid (array-like or None) – Grid for plasma response computation. If None, a default grid around the equilibrium is constructed automatically.

Returns:

• R_full (ndarray, shape (n_obs_full, n_coils))

• labels_full (list of str)

pyna.control.response_matrix.build_response_matrix(base_field_func, coil_field_funcs, state, observables=None, eps_current=1.0, raw_coil_B_funcs=None, base_raw_B_func=None)

Build response matrix R[n_obs, n_coils] using FPT formulae.

For each coil k, computes δ_state when I_k increases by eps_current, using closed-form FPT for axisymmetric observables (X/O-point shifts and DPm changes) and zero placeholders for observables that require additional geometric data (gaps, q-profile).

Parameters:

• base_field_func (callable) – Base equilibrium field function: [R,Z,phi] → [BR/|B|, BZ/|B|, Bphi/(R|B|)].

• coil_field_funcs (list of callable) –

  Coil perturbation field functions (per unit current). Legacy form (only used when raw_coil_B_funcs is None):

  coil_field_funcs[k]([R,Z,phi]) → [δBR/|B0|, δBZ/|B0|, δBphi/(R·|B0|)]

  This form is first-order correct for |δB| << |B0| away from X-points. Near the poloidal-field null (X-point), |Bpol/B| ~ 1e-5, so normalized coil components are unreliable. Use raw_coil_B_funcs instead.

• state (TopologyState) – Current topology state (must be pre-computed).

• observables (list of str or None) – Reserved for future observable filtering; currently unused.

• eps_current (float) – Current perturbation amplitude used to define per-unit response. (The result is divided by eps_current so R is per-ampere.)

• raw_coil_B_funcs (list of callable or None) –

  If provided, preferred over coil_field_funcs for DPm eigenvalue computation. Each callable returns raw coil B components:

  raw_coil_B_funcs[k]([R,Z,phi]) → [δBR, δBZ, δBphi] (in Tesla)

  The base field (B0) must then also be recoverable. Supply base_raw_B_func alongside this argument.

• base_raw_B_func (callable or None) –

  Returns raw base B components:

  base_raw_B_func([R,Z,phi]) → [BR0, BZ0, Bphi0] (in Tesla)

  Required when raw_coil_B_funcs is provided.

Returns:

• R_mat (ndarray, shape (n_obs, n_coils))

• obs_labels (list of str, length n_obs)

Multi-objective magnetic topology controller.

Solves the weighted least-squares problem:

min Σ_i w_i · |state_i + R_ij · δI_j − target_i|² + λ · ||δI||² s.t. I_min_j ≤ δI_j ≤ I_max_j (coil current limits)

Uses scipy.optimize.minimize with the SLSQP method.

Design principle — “avoid whack-a-mole”

• All objectives are optimised simultaneously, not sequentially.

• Weights encode the physicist’s priorities.

• Hard bound constraints enforce safety limits.

• L2 regularisation on ΔI penalises large coil changes.

class pyna.control.optimizer.ControlConstraints(I_max=None, I_min=None, gap_min=<factory>, xpoint_R_bounds=None, xpoint_Z_bounds=None)

Bases: object

Hard constraints for the optimisation.

I_max: ndarray | None = None

I_min: ndarray | None = None

gap_min: Dict[str, float]

xpoint_R_bounds: Tuple[float, float] | None = None

xpoint_Z_bounds: Tuple[float, float] | None = None

class pyna.control.optimizer.ControlWeights(gap_gi=100.0, xpoint_position=10.0, DPm_eigenvalue=5.0, DPm_eigenvector=2.0, Bpol_lcfs=2.0, q_profile=1.0, opoint_position=3.0, iota_core=0.5, delta_I_regularization=0.1)

Bases: object

Weights for multi-objective topology optimisation.

Higher weight → higher priority. Safety-critical observables (gap_gi) should carry the highest weight.

Bpol_lcfs: float = 2.0

DPm_eigenvalue: float = 5.0

DPm_eigenvector: float = 2.0

delta_I_regularization: float = 0.1

gap_gi: float = 100.0

iota_core: float = 0.5

opoint_position: float = 3.0

q_profile: float = 1.0

xpoint_position: float = 10.0

class pyna.control.optimizer.TopologyController(n_coils)

Bases: object

Multi-objective magnetic topology controller.

Uses the linear response matrix R (∂state/∂controls) to find the optimal coil-current changes δI that move the topology toward a target state while respecting coil limits.

For nonlinear problems iterate:

compute state → build R → solve → apply δI → repeat.

Example

>>> ctrl = TopologyController(n_coils=6)
>>> weights = ControlWeights(gap_gi=100.0, DPm_eigenvalue=10.0)
>>> delta_I, result = ctrl.solve(
...     current_state, target_state,
...     response_matrix, obs_labels,
...     weights=weights,
... )

predict_response(current_state, delta_I, response_matrix, obs_labels)

Predict how each observable changes under δI.

Return type:

Dict[str, float]

Returns:

dict {label (predicted_change})

solve(current_state, target_state, response_matrix, obs_labels, weights=None, constraints=None)

Solve multi-objective optimisation for optimal δI.

Parameters:

• current_state (TopologyState)

• target_state (TopologyState)

• response_matrix (ndarray, shape (n_obs, n_coils))

• obs_labels (list of str)

• weights (ControlWeights or None)

• constraints (ControlConstraints or None)

Return type:

Tuple[ndarray, object]

Returns:

• delta_I (ndarray, shape (n_coils,))

• result (scipy OptimizeResult)

Flux surface fate classification under perturbation.

Classifies whether a flux surface (characterised by ψ or ι) will:

• Remain intact (KAM torus survives)

• Deform but stay closed

• Break into an island chain (resonance â†?residue crosses 0 or 1)

• Become chaotic (Chirikov overlap)

Criteria used:

• Greene’s residue R = (2 âˆ?Tr(DPm)) / 4

• Relative DPm perturbation magnitude vs ε_KAM / ε_chaos thresholds

• Higher-order tensors D²Pm, D³Pm (future extension)

pyna.control.surface_fate.Greene_residue(DPm)

Greene’s residue R = (2 âˆ?Tr(DPm)) / 4.

Interpretation

---

R < 0 : hyperbolic (X-point-like), field lines diverge

0 < R < 1 : elliptic (O-point-like), field lines rotate

R > 1 : inverse hyperbolic

R = 0 or 1 : parabolic (marginal / transition point)

Parameters

---

DPm : ndarray, shape (2,2)

Poincaré map Jacobian (one full period).

Returns

---

float

Return type:

float

pyna.control.surface_fate.classify_surface_fate(iota, delta_iota, DPm, delta_DPm, epsilon_KAM=0.05, epsilon_chaos=0.3)

Classify how a flux surface responds to a field perturbation.

Parameters

---

iota : float

Unperturbed rotation transform ι = 1/q.

delta_iota : float

Change in rotation transform under perturbation.

DPm : ndarray, shape (2,2)

Unperturbed Poincaré map Jacobian.

delta_DPm : ndarray, shape (2,2)

Change in DPm under perturbation.

epsilon_KAM : float

Relative DPm-change threshold below which the KAM torus survives.

epsilon_chaos : float

Relative DPm-change threshold above which chaos is expected.

Returns

---

SurfaceFate

Return type:

SurfaceFate

pyna.control.surface_fate.scan_surface_fates(field_func, delta_field_func, psi_values, q_values, field_period=1, epsilon_KAM=0.05, epsilon_chaos=0.3)

Scan flux surfaces from core to edge and classify each.

For each (psi, q) pair this function would compute the A-matrix and DPm

at the corresponding O-point (requiring the actual flux surface geometry).

The present implementation returns UNKNOWN for all surfaces as a stub â€? a full implementation needs the O-point (R,Z) positions for each ψ.

Parameters

---

field_func : callable

Base magnetic field function.

delta_field_func : callable

Perturbation field function.

psi_values : ndarray

Normalised poloidal flux values ψ_N âˆ?[0, 1].

q_values : ndarray

Safety factor q at each ψ.

field_period : int

Toroidal field period (1 for tokamak, N_fp for stellarator).

epsilon_KAM : float

KAM survival threshold.

epsilon_chaos : float

Chaos threshold.

Returns

---

dict {psi_norm: (SurfaceFate, Greene_residue_before, delta_residue)}

Return type:

Dict[float, Tuple[SurfaceFate, float, float]]