1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2915
  6rng = np.random.default_rng(SEED)
  7
  8
  9def certificate(a, gain, gamma, eps=None, kappa=None):
 10    # Scalar A=-a, C=1, P=1, phi_gamma(x)=gamma*tanh(x).
 11    # OSL rho=gamma and QIB alpha=gamma^2 are globally valid.
 12    if eps is None:
 13        eps = gamma if gamma > 0 else 1.0
 14    M = -2.0 * (a + gain) + 2.0 * gamma + gamma * gamma / eps + eps
 15    # largest certifiable margin for this fixed P and selected epsilon
 16    return float(M), float(-M)
 17
 18
 19def dynamics_error(e, a, gain, gamma):
 20    # true x and estimate differ by e; nonlinear incremental term is exact.
 21    return -(a + gain) * e + gamma * np.tanh(e)
 22
 23
 24def rk4_step(e, h, a, gain, gamma):
 25    f = lambda z: dynamics_error(z, a, gain, gamma)
 26    k1 = f(e); k2 = f(e + .5*h*k1); k3 = f(e + .5*h*k2); k4 = f(e + h*k3)
 27    return e + h*(k1 + 2*k2 + 2*k3 + k4)/6
 28
 29
 30def estimate_rate(times, vals, start=0.2):
 31    mask = (times >= start) & (vals > 1e-14)
 32    if mask.sum() < 3: return float('nan')
 33    return float(np.polyfit(times[mask], np.log(vals[mask]), 1)[0])
 34
 35
 36def verify_bounds():
 37    # Numerical check of OSL/QIB over a broad pair grid, including unequal signs.
 38    gamma = 1.7
 39    xs = np.linspace(-8, 8, 401)
 40    worst_os = -np.inf
 41    worst_q = -np.inf
 42    for x in xs[::4]:
 43        for xhat in xs[::4]:
 44            e = xhat-x
 45            if abs(e) < 1e-9: continue
 46            dp = gamma*(np.tanh(xhat)-np.tanh(x))
 47            worst_os = max(worst_os, e*dp/e**2)
 48            worst_q = max(worst_q, dp**2/e**2)
 49    return {'rho_theory': gamma, 'rho_sample_max': float(worst_os),
 50            'alpha_theory': gamma**2, 'alpha_sample_max': float(worst_q),
 51            'bounds_hold': bool(worst_os <= gamma+1e-10 and worst_q <= gamma**2+1e-10)}
 52
 53
 54def boundary_sweep(a=1.0, gain=0.5):
 55    # With epsilon=gamma, M=-2(a+L)+4 gamma; sufficient certificate boundary is (a+L)/2.
 56    rows=[]
 57    for gamma in np.linspace(0, 2.2, 23):
 58        M, kap = certificate(a, gain, gamma)
 59        e=2.0; h=.002; n=int(8/h); vals=[]
 60        for _ in range(n):
 61            vals.append(abs(e)); e=rk4_step(e,h,a,gain,gamma)
 62        # observed contraction if final error is below initial by a meaningful amount
 63        observed = bool(vals[-1] < vals[0]*1e-3)
 64        rows.append((gamma, M, observed, vals[-1]))
 65    # Certificate sign transition (not finite-time threshold): M first becomes nonnegative.
 66    cert_boundary = next((float(g) for g,m,o,v in rows if m >= 0), float('nan'))
 67    # Exact local tanh stability boundary is gamma=a+L because tanh'(0)=1.
 68    exact_boundary = a + gain
 69    return rows, cert_boundary, exact_boundary
 70
 71
 72def rate_sweep(a=1.0, gain=0.5):
 73    # Prediction: log V slope is bounded above by M (the certificate matrix scalar).
 74    rows=[]
 75    for gamma in [0.0, .2, .5, .8, 1.0, 1.2, 1.4]:
 76        M, kap = certificate(a,gain,gamma)
 77        h=.001; times=np.arange(0, 3+h/2, h); e=0.02; vals=[]
 78        for t in times:
 79            vals.append(e*e); e=rk4_step(e,h,a,gain,gamma)
 80        obs=estimate_rate(times,np.asarray(vals),.3)
 81        rows.append({'gamma':gamma,'certificate_upper_bound_logV_slope':M,'observed_logV_slope':obs,'bound_holds':bool(obs <= M + 1e-3)})
 82    return rows
 83
 84
 85def gain_sweep(a=1.0, gamma=1.2, target_kappa=0.8):
 86    # From the optimized Young certificate: kappa=2(a+L)-4 gamma, so L_required=2 gamma-a+kappa/2.
 87    predicted = 2*gamma-a+target_kappa/2
 88    rows=[]
 89    for gain in [0,.1,.2,.3,.4,.5,.7,1.0,1.5,1.8,2.0]:
 90        M,kap=certificate(a,gain,gamma)
 91        rows.append({'gain':gain,'certificate_kappa':kap,'meets_target':bool(kap>=target_kappa-1e-9)})
 92    observed=min((r['gain'] for r in rows if r['meets_target']), default=float('nan'))
 93    return predicted, rows, observed
 94
 95
 96def noisy_comparison(a=1.0, gamma=0.5):
 97    # Same Euler observer setup; baseline has no correction, idea uses certified L.
 98    # State is driven by a known bounded input, while only noisy state output is observed.
 99    h=.01; steps=1000; noise=.08; trials=40
100    out={}
101    for name,gain in [('baseline_no_observer',0.0),('certified_observer',1.0)]:
102        errs=[]
103        for tr in range(trials):
104            x=1.5; xhat=0.0; se=0.
105            for k in range(steps):
106                u=.7*np.sin(.025*k)
107                # nonlinear plant: stable linear part plus bounded nonlinear residual and input
108                x += h*(-a*x + gamma*np.tanh(x) + u)
109                y=x + rng.normal(0,noise)
110                xhat += h*(-a*xhat + gamma*np.tanh(xhat) + u + gain*(y-xhat))
111                if k >= steps//2: se += (xhat-x)**2
112            errs.append(se/(steps//2))
113        out[name]={'mean_mse':float(np.mean(errs)),'std_mse':float(np.std(errs))}
114    return out
115
116
117def main():
118    a=1.0; gain=.5
119    bounds=verify_bounds()
120    boundary_rows, observed_boundary, exact_boundary=boundary_sweep(a,gain)
121    rates=rate_sweep(a,gain)
122    pred_gain,gain_rows,observed_gain=gain_sweep()
123    noisy=noisy_comparison()
124    result={
125      'seed':SEED,
126      'model':'scalar stable nonlinear plant with phi_gamma(x)=gamma*tanh(x), P=1, C=1',
127      'math':{'a':a,'observer_gain':gain,'rho': 'gamma','alpha':'gamma^2',
128              'certificate_at_eps_gamma':'M=-2*(a+L)+4*gamma',
129              'predicted_certificate_boundary_gamma':(a+gain)/2,
130              'predicted_exact_local_boundary_gamma':a+gain},
131      'bound_check':bounds,
132      'boundary_sweep':{'predicted_certificate_boundary':(a+gain)/2,'observed_certificate_sign_grid_gamma':observed_boundary,
133                        'predicted_exact_local_boundary':a+gain,'observed_exact_local_boundary':exact_boundary,
134                        'rows':[{'gamma':g,'M':m,'contracted':o,'final_abs_error':v} for g,m,o,v in boundary_rows]},
135      'rate_sweep':rates,
136      'gain_sweep':{'predicted_min_gain_for_kappa_0.8':pred_gain,'observed_grid_min_gain':observed_gain,'rows':gain_rows},
137      'noisy_comparison':noisy}
138    Path('results.json').write_text(json.dumps(result,indent=2))
139    print(json.dumps({'bound_check':bounds,'predicted_certificate_boundary':(a+gain)/2,
140      'observed_certificate_sign_grid':observed_boundary,
141      'predicted_exact_local_boundary':a+gain,'rate_sweep':rates,
142      'predicted_gain':pred_gain,'observed_gain_grid':observed_gain,'noisy':noisy},indent=2))
143
144if __name__=='__main__': main()