1import json
  2import numpy as np
  3from scipy.linalg import eigvalsh
  4
  5# Scalar latent plant: z_dot = A z, y = C z + w.
  6# The observation disturbance w enters the observer error through -K w.
  7A, C = -0.4, 1.0
  8alpha, gamma = 0.2, 0.9
  9
 10
 11def lmi_matrix(K, P):
 12    Ac = A - K * C
 13    Bc = -K                 # Bw=0, Dw=1
 14    Ce, De = 1.0, 0.0       # performance output q=e
 15    return np.array([[2 * Ac * P + Ce**2 + 2 * alpha * P,
 16                     P * Bc + Ce * De],
 17                     [P * Bc + Ce * De,
 18                      De**2 - gamma**2]])
 19
 20
 21def best_storage(K):
 22    # Search the one-dimensional positive storage metric, as a transparent
 23    # numerical analogue of optimizing P=LL^T+p_min I.
 24    ps = np.logspace(-4, 3, 600)
 25    vals = np.array([eigvalsh(lmi_matrix(K, p))[-1] for p in ps])
 26    i = int(vals.argmin())
 27    return float(ps[i]), float(vals[i])
 28
 29
 30def exact_error(K, e0, times, w):
 31    # Piecewise-constant forcing, exact propagation over each small interval.
 32    e = np.empty(len(times)); e[0] = e0
 33    dt = times[1] - times[0]
 34    Ac = A - K
 35    for j in range(len(times) - 1):
 36        # e_dot = Ac e - K w[j]
 37        e[j+1] = np.exp(Ac*dt)*e[j] + (-K*w[j]) * np.expm1(Ac*dt)/Ac
 38    return e
 39
 40
 41def dissipation_check(K, seed=7):
 42    rng = np.random.default_rng(seed)
 43    dt, T = 0.001, 8.0
 44    t = np.arange(0, T + dt, dt)
 45    # Smooth, piecewise disturbance, including an impulsive-like short pulse.
 46    w = 0.35 * rng.standard_normal(len(t)-1)
 47    w = np.convolve(w, np.ones(20)/20, mode='same')
 48    w[1000:1030] += 5.0
 49    e0 = 1.3
 50    e = exact_error(K, e0, t, w)
 51    P, lmax = best_storage(K)
 52    V = P * e**2
 53    lhs = (V[-1] - V[0] + np.sum((2*alpha*V[:-1] + e[:-1]**2 - gamma**2*w**2)*dt))
 54    # Direct induced-energy ratio (excluding initial-condition energy).
 55    zero = exact_error(K, 0.0, t, w)
 56    ratio = np.sum(zero[:-1]**2)*dt / (np.sum(w**2)*dt)
 57    return dict(P=P, lmi_max_eigenvalue=lmax, dissipation_residual=float(lhs),
 58                induced_energy_ratio=float(ratio), initial_V=float(V[0]),
 59                final_V=float(V[-1]))
 60
 61
 62def exact_hinf_gain(K):
 63    # For e_dot=-(A_abs+K)e-Kw, q=e, the stable transfer is -K/(s+0.4+K).
 64    # Its H-infinity norm is attained at zero frequency.
 65    return float(K / (K - A))
 66
 67
 68def zero_disturbance_decay(K):
 69    dt, T, e0 = 0.002, 8.0, 1.0
 70    t = np.arange(0, T+dt, dt)
 71    e = exact_error(K, e0, t, np.zeros(len(t)-1))
 72    # Compare observed decay exponent with required V decay exponent 2 alpha.
 73    P, _ = best_storage(K)
 74    V = P*e**2
 75    observed_rate = -np.log(V[-1]/V[0]) / T
 76    return dict(observed_V_decay_rate=float(observed_rate), required_rate=2*alpha,
 77                V_initial=float(V[0]), V_final=float(V[-1]))
 78
 79
 80def main():
 81    # Baseline: unconstrained high-gain innovation correction.
 82    baseline_K = 10.0
 83    # Idea: choose the largest gain on a coarse grid that is LMI feasible,
 84    # then refine by the storage search. This is a tiny proxy for training K,P.
 85    candidates = np.linspace(0.05, 10.0, 300)
 86    feasible = [(K, best_storage(K)) for K in candidates if best_storage(K)[1] <= 1e-8]
 87    idea_K, (idea_P, idea_lmi) = feasible[-1]
 88    result = {
 89        'config': {'A': A, 'C': C, 'alpha': alpha, 'gamma': gamma},
 90        'baseline': {'K': baseline_K, 'lmi': best_storage(baseline_K),
 91                     'exact_hinf_gain': exact_hinf_gain(baseline_K),
 92                     'zero_disturbance': zero_disturbance_decay(baseline_K),
 93                     'disturbance': dissipation_check(baseline_K)},
 94        'idea': {'K': float(idea_K), 'lmi': (idea_P, idea_lmi),
 95                 'exact_hinf_gain': exact_hinf_gain(idea_K),
 96                 'zero_disturbance': zero_disturbance_decay(idea_K),
 97                 'disturbance': dissipation_check(idea_K)},
 98        'claim_check': {
 99            'idea_lmi_feasible': bool(idea_lmi <= 1e-8),
100            'baseline_lmi_feasible': bool(best_storage(baseline_K)[1] <= 1e-8),
101            'idea_dissipation_holds': bool(dissipation_check(idea_K)['dissipation_residual'] <= 1e-5),
102            'idea_exact_hinf_below_gamma': bool(exact_hinf_gain(idea_K) <= gamma),
103            'baseline_exact_hinf_below_gamma': bool(exact_hinf_gain(baseline_K) <= gamma)
104        }
105    }
106    with open('results.json', 'w') as f:
107        json.dump(result, f, indent=2)
108    print(json.dumps(result, indent=2))
109
110if __name__ == '__main__':
111    main()