[education/kstars] kstars/ekos/guide/offlinetrainer: offline_trainer: recommend base RA/DEC PID gains from step-response data
Jasem Mutlaq <[email protected]>
| Newsgroups | gmane.comp.kde.cvs |
|---|---|
| Message-ID | <[email protected]> |
Git commit 51f67c8336b6f9fadbc580942d57a2e76f979666 by Jasem Mutlaq.
Committed on 04/08/2026 at 02:55.
Pushed by mutlaqja into branch 'master'.
offline_trainer: recommend base RA/DEC PID gains from step-response data
M +19 -0 kstars/ekos/guide/offlinetrainer/README.md
M +21 -3 kstars/ekos/guide/offlinetrainer/train.py
M +208 -8 kstars/ekos/guide/offlinetrainer/train_harmonic.py
https://invent.kde.org/education/kstars/-/commit/51f67c8336b6f9fadbc580942d57a2e76f979666
diff --git a/kstars/ekos/guide/offlinetrainer/README.md b/kstars/ekos/guide/offlinetrainer/README.md
index 35c3297ad3..edb9061d98 100644
--- a/kstars/ekos/guide/offlinetrainer/README.md
+++ b/kstars/ekos/guide/offlinetrainer/README.md
@@ -32,6 +32,7 @@ python train.py --sysid-data ./sysid_data.json --output ./weights.json
* `--mount-type <type>`: Force a specific mount type (`WORM_GEAR`, `HARMONIC_DRIVE`, `DIRECT_DRIVE`). Usually not needed as it auto-detects from the sysid file.
* `--verbose`: Print detailed logs of the physics fitting and the MLP training epochs.
* `--gpu`: Use CUDA/MPS to accelerate training if available (not necessary for these small models).
+* `--pid-lambda-factor <float>`: HARMONIC_DRIVE only. Design parameter for the advisory PID auto-tune recommendation below (default `3.0`); larger values recommend a slower, more conservative gain.
## 4. Loading the Model into KStars
@@ -39,3 +40,21 @@ python train.py --sysid-data ./sysid_data.json --output ./weights.json
2. In KStars, open your Equipment Profile.
3. Under the **AI Guiding** section, set the **Weights File** path to point to your `weights.json`.
4. Start guiding. The AI will automatically engage!
+
+## 5. PID Auto-Tune Recommendation (HARMONIC_DRIVE, advisory only)
+
+For `HARMONIC_DRIVE` mounts, the trainer also derives a recommended base RA/DEC
+proportional (+ conservative integral) gain from the same `pulse_response`
+step-response sessions used to fit κ/τ, via a conservative SIMC/IMC-style
+step-response tuning rule. This is printed to the console and saved into
+`weights.json` as `recommended_ra_proportional_gain`, `recommended_ra_integral_gain`,
+`recommended_dec_proportional_gain`, `recommended_dec_integral_gain`, plus a
+`pid_autotune` block with the underlying fit diagnostics (process gain, τ, dead
+time, confidence).
+
+**This is advisory only and is never applied automatically.** Nothing in KStars
+reads these fields back — review the numbers (especially the `confidence` flag:
+`"unavailable"` means don't use them, `"low"` means the dead time could not be
+resolved below the guide-frame sampling floor) and, if you agree with them,
+manually update the RA/DEC Proportional Gain (and Integral Gain, if used) in
+KStars' Guide options yourself.
diff --git a/kstars/ekos/guide/offlinetrainer/train.py b/kstars/ekos/guide/offlinetrainer/train.py
index 7f5fdfd29c..9fa0662b9a 100644
--- a/kstars/ekos/guide/offlinetrainer/train.py
+++ b/kstars/ekos/guide/offlinetrainer/train.py
@@ -32,6 +32,11 @@ def parse_args():
help="Use GPU if available (optional — all models train fine on CPU)")
p.add_argument("--epochs", type=int, default=None,
help="Override default epoch count for neural models")
+ p.add_argument("--pid-lambda-factor", type=float, default=None,
+ help="HARMONIC_DRIVE only: SIMC design parameter for the advisory PID "
+ "auto-tune recommendation (lambda = max(tau, factor * dead_time)). "
+ "Larger = slower/more conservative gain recommendation. "
+ "Default: train_harmonic.SIMC_LAMBDA_L_FACTOR")
p.add_argument("--simulate", action="store_true",
help="Run closed-loop simulation after training completes")
p.add_argument("--plot", action="store_true",
@@ -83,9 +88,9 @@ def main():
epochs=args.epochs, verbose=args.verbose)
elif mount_type == "HARMONIC_DRIVE":
- from train_harmonic import train_harmonic
- weights = train_harmonic(sysid, gpu=args.gpu,
- epochs=args.epochs, verbose=args.verbose)
+ from train_harmonic import train_harmonic, SIMC_LAMBDA_L_FACTOR
+ weights = train_harmonic(sysid, gpu=args.gpu, epochs=args.epochs, verbose=args.verbose,
+ pid_lambda_factor=args.pid_lambda_factor or SIMC_LAMBDA_L_FACTOR)
else:
sys.exit(
@@ -101,6 +106,19 @@ def main():
print(f"[Ekos AI Trainer] ✓ Weights saved to {args.output}")
print(f"[Ekos AI Trainer] Load this file in KStars Equipment Profile → AI Guiding → Weights File")
+ pid = weights.get("pid_autotune")
+ if pid:
+ print(f"\n[Ekos AI Trainer] PID auto-tune recommendation (advisory -- not applied automatically):")
+ for axis in ("ra", "dec"):
+ r = pid[axis]
+ if r["confidence"] == "unavailable":
+ print(f" {axis.upper()}: unavailable ({r['reason']})")
+ else:
+ print(f" {axis.upper()}: proportional_gain={r['proportional_gain']:.3f} "
+ f"integral_gain={r['integral_gain']:.3f} confidence={r['confidence']}"
+ + (" [dead time unresolved below sampling floor]" if r["resolution_limited"] else ""))
+ print(f" Review before changing Options::rA/dECProportionalGain() in KStars -- see pid_autotune_plan.md.")
+
if args.simulate:
print(f"\n[Ekos AI Trainer] Launching simulation...")
import subprocess
diff --git a/kstars/ekos/guide/offlinetrainer/train_harmonic.py b/kstars/ekos/guide/offlinetrainer/train_harmonic.py
index 0fb96f0ebe..99a21d3017 100644
--- a/kstars/ekos/guide/offlinetrainer/train_harmonic.py
+++ b/kstars/ekos/guide/offlinetrainer/train_harmonic.py
@@ -28,6 +28,28 @@ except ImportError:
+# Default SIMC design parameter for the PID auto-tune recommendation (§ Phase 1b
+# below): lambda = max(tau, SIMC_LAMBDA_L_FACTOR * dead_time). Larger -> slower/
+# more robust closed loop; smaller -> faster/less margin. Kept as an explicit,
+# overridable constant per the plan's "expose lambda as a documented knob, not a
+# hidden constant" (pid_autotune_plan.md §3.2).
+SIMC_LAMBDA_L_FACTOR = 3.0
+
+# KStars' "integral gain" (Options::rAIntegralGain()/dECIntegralGain()) multiplies
+# a ~100-frame moving average of drift (gmath.cpp::processAxis(), drift_integral[k]),
+# not a classical accumulating integrator. SIMC's tau_I (a reset *time*) has no
+# principled mapping onto that EMA-style term, so rather than inventing one, the
+# integral recommendation is a conservative fixed fraction of the proportional
+# recommendation (pid_autotune_plan.md §1: "keep I conservative"). tau_I is still
+# reported in the output for reference.
+INTEGRAL_GAIN_CONSERVATIVE_FRACTION = 0.25
+
+# Below this many usable step-response fits, the recommendation is flagged "low"
+# confidence rather than withheld outright -- still informative, not to be applied
+# unattended (pid_autotune_plan.md §4).
+MIN_FITS_FOR_MEDIUM_CONFIDENCE = 6
+
+
def _effective_pixel_scale(sysid):
"""Pixel scale in arcsec/px. Older exports recorded it without the binning factor."""
eq = sysid.get("equipment", {})
@@ -47,11 +69,20 @@ def _effective_pixel_scale(sysid):
def train_harmonic(sysid: dict,
gpu: bool = False,
epochs: int = None,
- verbose: bool = False) -> dict:
+ verbose: bool = False,
+ pid_lambda_factor: float = SIMC_LAMBDA_L_FACTOR) -> dict:
"""
Fit κ/τ spring parameters, detect PE, fit drift, and train Q-net
for a harmonic drive mount.
Returns a weights dict compatible with HarmonicGuider::loadWeights().
+
+ Also computes a "pid_autotune" recommendation block (RA/DEC proportional +
+ conservative integral gain, derived from the same pulse_response sessions
+ via a SIMC/IMC-style step-response rule -- pid_autotune_plan.md §3-5). This
+ is advisory only: nothing in this module or HarmonicGuider ever reads it
+ back or applies it automatically. A human reviews the numbers and manually
+ updates Options::setRAProportionalGain()/etc. in KStars, the same
+ human-in-the-loop workflow already used for loading a new weights.json.
"""
eq = sysid["equipment"]
pixel_scale = _effective_pixel_scale(sysid)
@@ -67,6 +98,12 @@ def train_harmonic(sysid: dict,
if verbose:
print(f" κ_ra={kappa_ra:.3f} τ_ra={tau_ra:.2f}s")
print(f" κ_dec={kappa_dec:.3f} τ_dec={tau_dec:.2f}s")
+ print(f"\n--- Phase 1b: PID Auto-tune (recommendation only, not applied) ---")
+
+ # ── Step 1b: Recommend base P(+I) guide gains from the same step-response data ──
+ pid_autotune = _recommend_pid_gains(sysid, guide_exp, verbose, pid_lambda_factor)
+
+ if verbose:
print(f"\n--- Phase 2: PE Period Detection ---")
# ── Step 2: Detect PE period from free-drift data ──────────────────────
@@ -99,6 +136,9 @@ def train_harmonic(sysid: dict,
qnet_weights = _train_qnet(sysid, kappa_ra, tau_ra, kappa_dec, tau_dec,
pe_period, gpu, epochs, verbose)
+ def _recommended(axis_result, field):
+ return axis_result[field] if axis_result["confidence"] != "unavailable" else None
+
from train_direct_drive import _build_fingerprint
return {
"format_version": "1.0",
@@ -107,6 +147,13 @@ def train_harmonic(sysid: dict,
"mount_name": eq.get("mount_name", "unknown"),
"pixel_scale": pixel_scale,
"model_fingerprint": _build_fingerprint(sysid),
+ # Advisory PID auto-tune recommendation -- see train_harmonic()'s docstring.
+ # recommended_* fields are None when confidence is "unavailable" (don't apply).
+ "recommended_ra_proportional_gain": _recommended(pid_autotune["ra"], "proportional_gain"),
+ "recommended_ra_integral_gain": _recommended(pid_autotune["ra"], "integral_gain"),
+ "recommended_dec_proportional_gain": _recommended(pid_autotune["dec"], "proportional_gain"),
+ "recommended_dec_integral_gain": _recommended(pid_autotune["dec"], "integral_gain"),
+ "pid_autotune": pid_autotune,
"physical": {
"kappa_ra": float(kappa_ra),
"tau_ra": float(tau_ra),
@@ -131,17 +178,31 @@ def train_harmonic(sysid: dict,
# Phase 1: Spring parameter fitting from pulse_response sessions
# ═══════════════════════════════════════════════════════════════════════════════
-def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool):
+def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool,
+ return_fits: bool = False):
"""
Fit spring constant κ and time constant τ from pulse_response sessions.
Model: d(t) = P * (1 - κ * exp(-t/τ)) + v*t + c. Fits whose |P| is not
significantly above the residual noise are skipped.
- Returns: (kappa, tau_seconds)
+ Returns: (kappa, tau_seconds), or (kappa, tau_seconds, fit_info) if
+ return_fits is True. fit_info["fits"] is the list of per-fit records
+ (P_fit_px, tau_fit_s, residual_std_px, t_first_s, t_arr, pos_arr) this
+ function already computes and would otherwise discard — reused by
+ _recommend_pid_gains() as the step-response data for PID auto-tune, a
+ second, independent consumer of the same pulse_response sessions.
+ fit_info["sign_consistent"] mirrors the gate this function itself uses
+ to decide the fits are real mechanics rather than noise.
"""
# Unmeasured means unmodeled: the default kappa stays 0
DEFAULTS = (0.0, 1.5)
+
+ def _finish(kappa, tau, fit_records, sign_consistent):
+ if not return_fits:
+ return kappa, tau
+ return kappa, tau, {"fits": fit_records, "sign_consistent": sign_consistent}
+
pulse_sessions = [
s for s in sysid["sessions"]
if s.get("type") == "pulse_response" and s.get("pulse_axis", "").upper() == axis.upper()
@@ -150,7 +211,7 @@ def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool):
if not pulse_sessions:
if verbose:
print(f" [{axis}] No pulse_response sessions found. Using defaults (κ=0.2, τ=1.5s)")
- return DEFAULTS
+ return _finish(*DEFAULTS, [], None)
axis_key = "ra_raw_px" if axis.upper() == "RA" else "dec_raw_px"
@@ -208,6 +269,7 @@ def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool):
fit_signs = []
paired_signs = set()
skipped_noise = 0
+ fit_records = []
def accept_fit(kappa_fit, tau_fit, t_first):
# tau at the upper bound: exponential degenerate with the drift term
@@ -258,6 +320,11 @@ def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool):
continue
paired_signs.add(1.0 if P_fit > 0 else -1.0)
accept_fit(kappa_fit, tau_fit, t_arr[0])
+ fit_records.append({
+ "pulse_magnitude_ms": float(pulse_mag), "P_fit_px": float(P_fit),
+ "tau_fit_s": float(tau_fit), "residual_std_px": float(residual_std),
+ "t_first_s": float(t_arr[0]), "t_arr": t_arr, "pos_arr": diff,
+ })
if verbose:
print(f" [{axis}] Pulse {pulse_mag}ms paired {pos_dir}-{neg_dir}: "
f"κ={kappa_fit:.3f}, τ={tau_fit:.2f}s (P={P_fit:.2f}px, noise={residual_std:.2f}px)")
@@ -281,6 +348,11 @@ def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool):
continue
accept_fit(kappa_fit, tau_fit, t_arr[0])
fit_signs.append((s.get("pulse_direction", "?"), np.sign(P_fit)))
+ fit_records.append({
+ "pulse_magnitude_ms": float(pulse_mag), "P_fit_px": float(P_fit),
+ "tau_fit_s": float(tau_fit), "residual_std_px": float(residual_std),
+ "t_first_s": float(t_arr[0]), "t_arr": t_arr, "pos_arr": pos_arr,
+ })
if verbose:
print(f" [{axis}] Pulse {pulse_mag}ms {s.get('pulse_direction', '?')}: "
f"κ={kappa_fit:.3f}, τ={tau_fit:.2f}s (P={P_fit:.2f}px, noise={residual_std:.2f}px)")
@@ -290,7 +362,7 @@ def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool):
print(f" [{axis}] No pulse response measurable above noise "
f"({skipped_noise} skipped). Using defaults (κ=0.2, τ=1.5s). "
f"Consider larger protocol pulses.")
- return DEFAULTS
+ return _finish(*DEFAULTS, fit_records, None)
# Real responses have consistent signs per direction; paired diffs share one sign
by_dir = {}
@@ -304,7 +376,7 @@ def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool):
if verbose:
print(f" [{axis}] WARNING: response signs inconsistent across pulse directions "
f"— fits are noise, not mechanics. Using defaults (κ=0.2, τ=1.5s).")
- return DEFAULTS
+ return _finish(*DEFAULTS, fit_records, False)
kappa_result = float(np.median(kappas))
tau_result = float(np.median(taus)) if taus and kappa_result > 0.0 else DEFAULTS[1]
@@ -314,13 +386,141 @@ def _fit_spring_params(sysid: dict, axis: str, guide_exp: float, verbose: bool):
if verbose:
print(f" [{axis}] WARNING: fit pinned at bounds (κ={kappa_result:.3f}, "
f"τ={tau_result:.2f}s) — unphysical. Using defaults (κ=0.2, τ=1.5s).")
- return DEFAULTS
+ return _finish(*DEFAULTS, fit_records, consistent)
if verbose:
print(f" [{axis}] Final: κ={kappa_result:.3f} (from {len(kappas)} fits), "
f"τ={tau_result:.2f}s")
- return kappa_result, tau_result
+ return _finish(kappa_result, tau_result, fit_records, consistent)
+
+
+# ═══════════════════════════════════════════════════════════════════════════════
+# Phase 1b: PID auto-tune -- SIMC step-response gain recommendation
+#
+# Offline, calibration-time-only calculation (pid_autotune_plan.md). Reuses the
+# per-fit P_fit/tau_fit/residual_std/t_arr values _fit_spring_params() already
+# computes from pulse_response sessions -- no new data collection, just a second
+# consumer of the existing step-response fits.
+# ═══════════════════════════════════════════════════════════════════════════════
+
+def _estimate_dead_time_s(t_arr: np.ndarray, pos_arr: np.ndarray, residual_std: float) -> float:
+ """
+ First t at which |pos(t)| exceeds ~2.5x the fit's residual noise -- a proxy
+ for FOPDT dead time L (pid_autotune_plan.md §3.1). If even the first sample
+ already exceeds threshold (the common case at ~2-3s/frame cadence), this
+ returns t_arr[0]: L is only known to be <= the first sample, not resolved any
+ finer -- see _recommend_axis_pid_gain()'s "resolution_limited" flag.
+ """
+ threshold = max(2.5 * residual_std, 1e-6)
+ for t, p in zip(t_arr, pos_arr):
+ if abs(p) > threshold:
+ return float(t)
+ return float(t_arr[-1])
+
+
+def _calibration_ms_per_arcsec(sysid: dict, axis: str) -> float:
+ """
+ Median calibrated ms-per-arcsec pulse rate for this axis, from whichever
+ sessions recorded it (standard_guiding sessions carry ra_ms_per_arcsec/
+ dec_ms_per_arcsec). This is Calibration::ra/decPulseMillisecondsPerArcsecond()
+ at collection time -- the same normalization gmath.cpp::processAxis() applies
+ to proportional_gain, needed to convert a physical SIMC Kc (ms/arcsec) back
+ into KStars' dimensionless aggressiveness (pid_autotune_plan.md §3.3).
+ """
+ key = "ra_ms_per_arcsec" if axis.upper() == "RA" else "dec_ms_per_arcsec"
+ values = [s[key] for s in sysid["sessions"] if s.get(key, 0.0) and s[key] > 0.0]
+ return float(np.median(values)) if values else 0.0
+
+
+def _recommend_axis_pid_gain(sysid: dict, axis: str, guide_exp: float,
+ lambda_l_factor: float, verbose: bool) -> dict:
+ """
+ Derive a recommended base proportional/integral gain for one axis from its
+ pulse_response step-response fits, via a conservative SIMC/IMC-style PI rule
+ (pid_autotune_plan.md §3.2). Returns a dict; confidence "unavailable" means
+ the numbers (if present at all) should not be used.
+ """
+ _, _, fit_info = _fit_spring_params(sysid, axis, guide_exp, False, return_fits=True)
+ fits = fit_info["fits"]
+
+ if fit_info["sign_consistent"] is False:
+ return {"confidence": "unavailable",
+ "reason": "pulse-response signs inconsistent across directions -- fits are noise, not mechanics"}
+ if not fits:
+ return {"confidence": "unavailable", "reason": "no usable pulse-response step-response fits"}
+
+ cal_ms_per_arcsec = _calibration_ms_per_arcsec(sysid, axis)
+ if cal_ms_per_arcsec <= 0.0:
+ return {"confidence": "unavailable",
+ "reason": "no calibrated ms_per_arcsec recorded for this axis (need a standard_guiding session)"}
+
+ pixel_scale = _effective_pixel_scale(sysid)
+
+ K_samples, tau_samples, L_samples, t_first_samples = [], [], [], []
+ for f in fits:
+ if f["pulse_magnitude_ms"] <= 0.0:
+ continue
+ K_samples.append(abs(f["P_fit_px"]) * pixel_scale / f["pulse_magnitude_ms"])
+ t_first_samples.append(f["t_first_s"])
+ L_samples.append(_estimate_dead_time_s(f["t_arr"], f["pos_arr"], f["residual_std_px"]))
+ if f["tau_fit_s"] <= 9.8: # same "pinned at bound, degenerate with drift" guard as the spring fit
+ tau_samples.append(f["tau_fit_s"])
+
+ if not K_samples:
+ return {"confidence": "unavailable", "reason": "no fit had a usable pulse magnitude"}
+
+ K = float(np.median(K_samples)) # arcsec of steady-state response per ms of pulse
+ L = float(np.median(L_samples))
+ # tau can't be resolved below the sampling floor either: a fit whose tau_fit
+ # landed below t_first ("spring already released") is only known to be
+ # <= t_first, not physically ~0 -- floor it the same way L is floored, so a
+ # spuriously tiny fitted tau can't produce a dangerously aggressive Kc.
+ tau_raw = float(np.median(tau_samples)) if tau_samples else float(np.median(t_first_samples))
+ tau = max(tau_raw, L)
+ resolution_limited = bool(np.isclose(L, float(np.median(t_first_samples)), rtol=0.05))
+
+ lam = max(tau, lambda_l_factor * L)
+ Kc = (1.0 / K) * tau / (lam + L) # ms of pulse per arcsec of error
+ tau_i = min(tau, 4.0 * (lam + L)) # SIMC reset time, reported only -- see constant doc above
+
+ proportional_gain = Kc / cal_ms_per_arcsec
+ integral_gain = INTEGRAL_GAIN_CONSERVATIVE_FRACTION * proportional_gain
+ confidence = "low" if (resolution_limited or len(fits) < MIN_FITS_FOR_MEDIUM_CONFIDENCE) else "medium"
+
+ if verbose:
+ print(f" [{axis} PID] K={K:.5f} arcsec/ms tau={tau:.2f}s L={L:.2f}s lambda={lam:.2f}s "
+ f"(n={len(fits)} fits, cal={cal_ms_per_arcsec:.1f}ms/arcsec)")
+ print(f" [{axis} PID] Recommended proportional_gain={proportional_gain:.3f} "
+ f"integral_gain={integral_gain:.3f} confidence={confidence}")
+
+ return {
+ "confidence": confidence,
+ "proportional_gain": float(np.clip(proportional_gain, 0.0, 1.0)),
+ "integral_gain": float(np.clip(integral_gain, 0.0, 1.0)),
+ "process_gain_arcsec_per_ms": K,
+ "tau_s": tau,
+ "dead_time_s": L,
+ "lambda_s": lam,
+ "tau_i_s": tau_i,
+ "calibration_ms_per_arcsec": cal_ms_per_arcsec,
+ "n_fits": len(fits),
+ "resolution_limited": resolution_limited,
+ }
+
+
+def _recommend_pid_gains(sysid: dict, guide_exp: float, verbose: bool,
+ lambda_l_factor: float = SIMC_LAMBDA_L_FACTOR) -> dict:
+ """
+ Recommend base RA/DEC proportional (+ conservative integral) gains from the
+ pulse_response step-response data, via a conservative SIMC/IMC-style PI rule
+ (pid_autotune_plan.md §3). Calibration-time only, never auto-applied -- see
+ train_harmonic()'s docstring for how these surface in the weights JSON.
+ """
+ return {
+ "ra": _recommend_axis_pid_gain(sysid, "RA", guide_exp, lambda_l_factor, verbose),
+ "dec": _recommend_axis_pid_gain(sysid, "DEC", guide_exp, lambda_l_factor, verbose),
+ }
# ═══════════════════════════════════════════════════════════════════════════════