Advanced Control — Phase 4 Modules (v0.16.0)¶

Demonstrates the ten new control modules added in scpn-control v0.16.0:

  1. Sliding-mode vertical stabilizer — super-twisting algorithm
  2. Gain-scheduled controller — regime-adaptive PID
  3. RWM feedback — resistive wall mode stabilization
  4. Mu-synthesis — structured uncertainty D-K iteration
  5. Fault-tolerant control — sensor/actuator fault detection + reconfiguration
  6. Shape controller — isoflux plasma boundary control

Dependencies: numpy, scipy, matplotlib (no GPU, no optional deps)

In [1]:
import sys
from pathlib import Path

sys.path.insert(0, str(Path(".").resolve().parent / "src"))

import numpy as np
import matplotlib.pyplot as plt

1. Sliding-Mode Vertical Stabilizer¶

The VerticalStabilizer wraps a super-twisting second-order sliding mode controller (Levant 1998) for tokamak vertical position control. Finite-time convergence and chattering-free continuous output.

In [2]:
from scpn_control.control.sliding_mode_vertical import (
    SuperTwistingSMC,
    VerticalStabilizer,
    lyapunov_certificate,
)

smc = SuperTwistingSMC(alpha=50.0, beta=100.0, c=0.5, u_max=1e5)
vs = VerticalStabilizer(
    n_index=1.5,
    Ip_MA=15.0,
    R0=6.2,
    m_eff=1e6,
    tau_wall=0.01,
    smc=smc,
)

print(f"Lyapunov certificate: {lyapunov_certificate(50.0, 100.0, 200.0)}")

dt = 1e-4
n_steps = 5000
Z = 0.05  # 5 cm initial displacement
dZ = 0.0
Z_ref = 0.0

z_hist, u_hist = [Z], []
for _ in range(n_steps):
    u = vs.step(Z, Z_ref, dZ, dt)
    # Simplified plant: M * ddZ = F_unstable + u
    gamma_growth = vs.n_index * 4e-7 * np.pi * vs.Ip**2 / (4 * np.pi**2 * vs.R0 * vs.m_eff)
    ddZ = gamma_growth * Z + u / vs.m_eff
    dZ += ddZ * dt
    Z += dZ * dt
    z_hist.append(Z)
    u_hist.append(u)

t = np.arange(len(z_hist)) * dt * 1e3  # ms

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6), sharex=True)
ax1.plot(t, np.array(z_hist) * 100, "b-", lw=1.5)
ax1.set_ylabel("Z displacement [cm]")
ax1.set_title("Super-twisting vertical stabilizer")
ax1.grid(True, alpha=0.3)

ax2.plot(t[:-1], np.array(u_hist) / 1e3, "r-", lw=1)
ax2.set_xlabel("Time [ms]")
ax2.set_ylabel("Control [kA-turns]")
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Lyapunov certificate: False
No description has been provided for this image

2. Gain-Scheduled Controller — Regime Detection¶

The GainScheduledController switches PID gains based on operating regime (ramp-up, L-mode, LH transition, H-mode, ramp-down). The RegimeDetector classifies plasma state automatically.

In [3]:
from scpn_control.control.gain_scheduled_controller import (
    RegimeDetector,
    OperatingRegime,
)

detector = RegimeDetector()

# Simulate a discharge: ramp-up -> L-mode -> H-mode -> ramp-down
regimes = []
n_phases = 4
steps_per_phase = 50

for phase in range(n_phases):
    for step in range(steps_per_phase):
        t_frac = (phase * steps_per_phase + step) / (n_phases * steps_per_phase)
        if phase == 0:  # ramp-up
            state = np.array([t_frac * 15, 0.5, 3.0, 0.8, 50.0, 5.0])
            dstate = np.array([15.0 / steps_per_phase, 0, 0, 0, 0, 0])
        elif phase == 1:  # L-mode flat
            state = np.array([15.0, 1.0, 3.0, 0.9, 100.0, 8.0])
            dstate = np.zeros(6)
        elif phase == 2:  # H-mode
            state = np.array([15.0, 2.5, 3.0, 0.85, 250.0, 10.0])
            dstate = np.zeros(6)
        else:  # ramp-down
            state = np.array([(1 - (step / steps_per_phase)) * 15, 1.0, 3.0, 0.9, 80.0, 6.0])
            dstate = np.array([-15.0 / steps_per_phase, 0, 0, 0, 0, 0])

        regime = detector.detect(state, dstate, tau_E=0.5, p_disrupt=0.01)
        regimes.append(regime)

regime_names = [r.name for r in regimes]
regime_ints = [r.value for r in regimes]

plt.figure(figsize=(10, 3))
plt.plot(regime_ints, "k-", lw=1.5)
plt.yticks(
    range(len(OperatingRegime)),
    [r.name for r in OperatingRegime],
    fontsize=8,
)
plt.xlabel("Time step")
plt.title("Regime detection across discharge phases")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image

3. RWM Feedback — Resistive Wall Mode Stabilization¶

The RWMFeedbackController implements proportional-derivative feedback on radial magnetic field sensors to stabilize resistive wall modes. The RWMPhysics model computes growth rates from beta limits.

In [4]:
from scpn_control.control.rwm_feedback import RWMPhysics, RWMFeedbackController

# Scan beta_N and compute growth rates
beta_n_range = np.linspace(2.0, 4.5, 50)
beta_nowall = 2.8
beta_wall = 4.2
tau_wall = 0.005  # 5 ms

gamma_open = []
gamma_closed = []

ctrl_rwm = RWMFeedbackController(n_sensors=8, n_coils=6, G_p=5.0, G_d=0.1)

for bn in beta_n_range:
    rwm = RWMPhysics(bn, beta_nowall, beta_wall, tau_wall)
    gamma_open.append(rwm.growth_rate())
    gamma_closed.append(ctrl_rwm.effective_growth_rate(rwm))

plt.figure(figsize=(10, 5))
plt.plot(beta_n_range, gamma_open, "r-", lw=2, label="Open-loop")
plt.plot(beta_n_range, gamma_closed, "b-", lw=2, label="Closed-loop (G_p=5)")
plt.axhline(0, color="gray", ls="--", lw=0.5)
plt.axvline(beta_nowall, color="green", ls=":", label=f"beta_nowall={beta_nowall}")
plt.axvline(beta_wall, color="orange", ls=":", label=f"beta_wall={beta_wall}")
plt.xlabel(r"$\beta_N$")
plt.ylabel(r"Growth rate $\gamma$ [s$^{-1}$]")
plt.title("RWM feedback stabilization")
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(-500, 2000)
plt.show()

# Check stabilization at beta_N = 3.5
rwm_test = RWMPhysics(3.5, beta_nowall, beta_wall, tau_wall)
print(
    f"beta_N=3.5: unstable={rwm_test.is_unstable()}, "
    f"gamma_open={rwm_test.growth_rate():.1f} s^-1, "
    f"stabilized={ctrl_rwm.is_stabilized(rwm_test)}"
)
No description has been provided for this image 
beta_N=3.5: unstable=True, gamma_open=200.0 s^-1, stabilized=True

4. Mu-Synthesis — Structured Uncertainty Robust Control¶

The MuSynthesisController implements D-K iteration for controllers robust against structured parametric uncertainty. The compute_mu_upper_bound function computes the structured singular value via D-scaling optimization.

In [5]:
from scpn_control.control.mu_synthesis import (
    compute_mu_upper_bound,
    StructuredUncertainty,
    UncertaintyBlock,
)

# Define uncertainty blocks
uncertainty = StructuredUncertainty(
    [
        UncertaintyBlock("plasma_current", size=1, bound=0.1, block_type="real_scalar"),
        UncertaintyBlock("wall_position", size=1, bound=0.2, block_type="complex_scalar"),
    ]
)

# Compute mu upper bound for a sample transfer matrix
n = uncertainty.total_size()
rng = np.random.default_rng(42)

# Sweep frequency and compute mu
freqs = np.logspace(-1, 3, 50)
mu_vals = []
for omega in freqs:
    # Frequency-dependent transfer matrix (simplified)
    s = 1j * omega
    M = np.array(
        [
            [1.0 / (s + 10), 0.5 / (s + 20)],
            [0.3 / (s + 5), 1.0 / (s + 15)],
        ]
    )
    delta_struct = uncertainty.build_Delta_structure()
    mu = compute_mu_upper_bound(M, delta_struct)
    mu_vals.append(mu)

plt.figure(figsize=(8, 4))
plt.semilogx(freqs, mu_vals, "b-", lw=2)
plt.axhline(1.0, color="red", ls="--", label="Robustness boundary")
plt.xlabel("Frequency [rad/s]")
plt.ylabel(r"$\mu$ upper bound")
plt.title("Structured singular value vs frequency")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

print(f"Peak mu = {max(mu_vals):.4f} (< 1 means robust stability)")

/tmp/ipykernel_3092279/3814760742.py:31: DeprecationWarning: StructuredUncertainty.build_Delta_structure is deprecated and will be removed in 0.25.0; use scpn_control.control.static_mu_analysis.StructuredUncertainty.build_delta_structure.
  delta_struct = uncertainty.build_Delta_structure()
/tmp/ipykernel_3092279/3814760742.py:32: DeprecationWarning: compute_mu_upper_bound is deprecated and will be removed in 0.25.0; use scpn_control.control.static_mu_analysis.compute_static_mu_upper_bound.
  mu = compute_mu_upper_bound(M, delta_struct)
No description has been provided for this image 
Peak mu = 0.1272 (< 1 means robust stability)

5. Fault-Tolerant Control — Detection & Reconfiguration¶

The FDIMonitor (Fault Detection & Isolation) detects sensor dropouts, drifts, and stuck actuators from residual analysis. The ReconfigurableController nulls out faulty channels.

In [6]:
from scpn_control.control.fault_tolerant_control import FDIMonitor

fdi = FDIMonitor(n_sensors=4, n_actuators=3, threshold_sigma=3.0, n_alert=3)

rng = np.random.default_rng(42)
n_steps = 100
fault_step = 40

faults_detected = []
residuals_s2 = []

for i in range(n_steps):
    y_pred = np.array([1.0, 2.0, 3.0, 4.0])
    noise = rng.normal(0, 0.1, 4)

    if i >= fault_step:
        # Inject sensor 2 dropout (stuck at last good value)
        y_meas = y_pred + noise
        y_meas[2] = 3.0  # stuck
    else:
        y_meas = y_pred + noise

    reports = fdi.update(y_meas, y_pred, t=float(i))
    faults_detected.append(len(reports))
    residuals_s2.append(abs(y_meas[2] - y_pred[2]))

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
ax1.plot(residuals_s2, "b-", lw=1)
ax1.axvline(fault_step, color="red", ls="--", label="Fault injected")
ax1.set_ylabel("|Residual| sensor 2")
ax1.set_title("Sensor fault detection")
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2.bar(range(n_steps), faults_detected, width=1.0, color="orange")
ax2.set_xlabel("Time step")
ax2.set_ylabel("Faults detected")
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

first_detection = next((i for i, f in enumerate(faults_detected) if f > 0), None)
print(f"Fault injected at step {fault_step}, detected at step {first_detection}")
No description has been provided for this image 
Fault injected at step 40, detected at step None

6. Shape Controller — Isoflux Plasma Boundary¶

The PlasmaShapeController computes coil current corrections to maintain the desired plasma boundary shape using a Tikhonov-regularized pseudoinverse of the coil-to-shape Jacobian.

In [7]:
from scpn_control.control.shape_controller import (
    PlasmaShapeController,
    ShapeTarget,
    CoilSet,
)

target = ShapeTarget(
    isoflux_points=[(8.2, 0.0), (7.8, 1.5), (6.5, 3.0), (5.0, 2.5), (4.5, 0.0), (5.0, -2.5), (6.5, -3.0), (7.8, -1.5)],
    gap_points=[(8.5, 0.0, 1.0, 0.0), (4.2, 0.0, -1.0, 0.0)],
    gap_targets=[0.1, 0.08],
)

coils = CoilSet(n_coils=6)
# kernel=None works — ShapeJacobian uses a mock Jacobian for the demo
shape_ctrl = PlasmaShapeController(target=target, coil_set=coils, kernel=None)

# Simulate shape correction over several iterations
psi = np.ones((33, 33)) * 0.5  # mock flux grid
coil_currents = np.zeros(6)

corrections = []
for i in range(10):
    delta_I = shape_ctrl.step(psi, coil_currents)
    coil_currents = np.clip(coil_currents + delta_I, -coils.max_currents, coils.max_currents)
    corrections.append(np.linalg.norm(delta_I))

plt.figure(figsize=(8, 4))
plt.semilogy(corrections, "bo-", lw=2)
plt.xlabel("Iteration")
plt.ylabel("|delta_I| [A]")
plt.title("Shape controller convergence")
plt.grid(True, alpha=0.3)
plt.show()

print(f"Final coil currents: {np.round(coil_currents, 1)}")
No description has been provided for this image 
Final coil currents: [-33.5  38.9  -3.5 -86.   72.6 -38.9]

Summary: scpn-control v0.16.0 adds production-grade controllers for vertical stability, shape control, RWM stabilization, fault tolerance, and robust synthesis under structured uncertainty. These compose with the existing H-infinity, MPC, and SNN controllers for multi-layered plasma control architectures.

See docs/api.md for the complete API reference.