1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2907
  6
  7
  8def aug_matrix(A, B, KC, KI, rho):
  9    return np.array([[A, -B * KC], [KI, rho]], dtype=float)
 10
 11
 12def spectral_radius(M):
 13    return float(np.max(np.abs(np.linalg.eigvals(M))))
 14
 15
 16def unsaturated_decay(A, B, KC, KI, rho, n=100):
 17    M = aug_matrix(A, B, KC, KI, rho)
 18    v = np.array([1.0, 0.0])
 19    norms = []
 20    for _ in range(n):
 21        norms.append(float(np.linalg.norm(v)))
 22        v = M @ v
 23    # Dominant asymptotic factor estimated from log(norm) after transients.
 24    fit = np.polyfit(np.arange(25, n), np.log(np.maximum(np.asarray(norms[25:]), 1e-15)), 1)
 25    return spectral_radius(M), float(np.exp(fit[0]))
 26
 27
 28def stability_sweep():
 29    A, B, KC, rho = 0.82, 1.0, 0.32, 0.68
 30    kis = np.linspace(0.0, 1.8, 181)
 31    rows = []
 32    for ki in kis:
 33        sr = spectral_radius(aug_matrix(A, B, KC, ki, rho))
 34        M = aug_matrix(A, B, KC, ki, rho)
 35        v = np.array([1.0, 0.0])
 36        for _ in range(500):
 37            v = M @ v
 38        # A deliberately strict finite-horizon divergence test, not the theory.
 39        empirical_stable = bool(np.linalg.norm(v) < 1e-3)
 40        rows.append((float(ki), sr, empirical_stable))
 41    stable = [x for x in rows if x[2]]
 42    empirical_boundary = max(x[0] for x in stable) if stable else float('nan')
 43    crossing = [rows[i] for i in range(1, len(rows))
 44                if rows[i-1][1] < 1.0 <= rows[i][1]]
 45    predicted_boundary = crossing[0][0] if crossing else float('nan')
 46    return {
 47        'A': A, 'B': B, 'KC': KC, 'rho': rho,
 48        'predicted_KI_boundary_rho_eq_1': predicted_boundary,
 49        'observed_KI_boundary_500_step_state_decay_below_1e-3': empirical_boundary,
 50        'boundary_abs_difference': abs(predicted_boundary - empirical_boundary),
 51        'spectral_stability_classification_agreement': float(np.mean([x[2] == (x[1] < 1.0) for x in rows])),
 52    }
 53
 54
 55def saturation_episode(ki, ka, sat=0.25, n=180, hold_end=95, rho=0.92, kp=0.28, kc=0.55):
 56    y, c = 0.0, 0.0
 57    ys, rs, cs, es = [], [], [], []
 58    for t in range(n):
 59        yd = 1.0 if 10 <= t < hold_end else 0.0
 60        e = yd - y
 61        tilde = kp * e + kc * c  # nominal neural reference z_t is zero
 62        r = float(np.clip(tilde, -sat, sat))
 63        d = r - tilde
 64        y = 0.92 * y + 0.08 * r
 65        c = rho * c + ki * e + ka * d
 66        ys.append(y); rs.append(r); cs.append(c); es.append(e)
 67    es, cs, rs = np.asarray(es), np.asarray(cs), np.asarray(rs)
 68    recovery = None
 69    for t in range(hold_end, n):
 70        if np.max(np.abs(es[t:min(t + 8, n)])) < 0.1:
 71            recovery = t - hold_end
 72            break
 73    return {
 74        'rms_error': float(np.sqrt(np.mean(es ** 2))),
 75        'post_release_rms_error': float(np.sqrt(np.mean(es[hold_end:] ** 2))),
 76        'max_abs_c': float(np.max(np.abs(cs))),
 77        'release_abs_c': float(abs(cs[hold_end])),
 78        'recovery_steps': recovery if recovery is not None else n - hold_end,
 79        'max_constraint_violation': float(np.max(np.maximum(np.abs(rs) - sat, 0.0))),
 80        'terminal_abs_error': float(abs(es[-1])),
 81    }
 82
 83
 84def saturation_sweep():
 85    durations = np.array([10, 25, 40, 55, 70])
 86    ordinary = [saturation_episode(.075, 0.0, hold_end=10 + int(d)) for d in durations]
 87    antiwindup = [saturation_episode(.075, .9, hold_end=10 + int(d)) for d in durations]
 88    ordinary_c = np.array([x['release_abs_c'] for x in ordinary])
 89    anti_c = np.array([x['release_abs_c'] for x in antiwindup])
 90    # Prediction: residual feedback lowers the slope of windup versus saturation duration.
 91    ordinary_slope = float(np.polyfit(durations, ordinary_c, 1)[0])
 92    anti_slope = float(np.polyfit(durations, anti_c, 1)[0])
 93    return {
 94        'durations': durations.tolist(), 'ordinary_integral': ordinary,
 95        'residual_antiwindup': antiwindup,
 96        'predicted_ordering': 'antiwindup slope < ordinary slope',
 97        'observed_release_c_slope_ordinary': ordinary_slope,
 98        'observed_release_c_slope_antiwindup': anti_slope,
 99        'observed_slope_ratio_anti_over_ordinary': anti_slope / ordinary_slope,
100    }
101
102
103def unsaturated_sweep():
104    A, B, KC, rho = .80, 1.0, .30, .70
105    kis = [.02, .08, .16, .24]
106    rows = []
107    for ki in kis:
108        predicted, observed = unsaturated_decay(A, B, KC, ki, rho)
109        rows.append({'KI': ki, 'predicted_spectral_radius': predicted,
110                     'observed_asymptotic_decay_factor': observed,
111                     'relative_error': abs(observed-predicted)/predicted})
112    return rows
113
114
115def main():
116    sat = saturation_sweep()
117    result = {
118        'seed': SEED,
119        'math_prediction_1_unsaturated_decay': unsaturated_sweep(),
120        'math_prediction_2_stability_boundary': stability_sweep(),
121        'math_prediction_3_saturation_windup_scaling': sat,
122        'mini_experiment_baseline_no_compensation': saturation_episode(0.0, 0.0),
123        'mini_experiment_ordinary_integral': saturation_episode(.075, 0.0),
124        'mini_experiment_idea_residual_antiwindup': saturation_episode(.075, .9),
125    }
126    Path('results.json').write_text(json.dumps(result, indent=2))
127    print(json.dumps(result, indent=2))
128
129if __name__ == '__main__':
130    main()