1import json, math
  2from pathlib import Path
  3import numpy as np
  4from scipy.linalg import eigvals
  5from scipy.optimize import minimize
  6from scipy.special import lambertw
  7
  8# Scalar delayed recurrent core:
  9# xdot(t) = -a*x(t) - b*x(t-h) + r(t), p=x, q=r.
 10# The stated LMI then certifies ||p||_2^2 <= gamma ||q||_2^2.
 11
 12def pade_system(a, b, h):
 13    if h == 0:
 14        return np.array([[-a-b]]), np.array([[1.]]), np.array([[1.]]), np.array([[0.]])
 15    # y_delay = -x + 2 z, z_dot = -2/h*z + 2/h*x
 16    A = np.array([[-a+b, -2*b], [2/h, -2/h]], float)
 17    B = np.array([[1.], [0.]])
 18    C = np.array([[1., 0.]])
 19    D = np.array([[0.]])
 20    return A, B, C, D
 21
 22def exact_roots(a,b,h,branches=range(-20,21)):
 23    if h == 0: return np.array([-a-b+0j])
 24    z = -b*h*np.exp(a*h)
 25    return np.array([-a + lambertw(z,k)/h for k in branches])
 26
 27def exact_stable(a,b,h):
 28    return np.max(exact_roots(a,b,h).real) < 0
 29
 30def critical_delay(a,b):
 31    # for b>a, xdot=-a x-b x(t-h), first imaginary crossing
 32    if b <= a: return float('inf')
 33    return math.acos(-a/b)/math.sqrt(b*b-a*a)
 34
 35def hinf(A,B,C,D, n=30000):
 36    if np.max(eigvals(A).real) >= -1e-8: return float('inf')
 37    poles = np.abs(eigvals(A))
 38    wmax = max(100., 100.*float(np.max(poles)))
 39    ws = np.concatenate(([0.], np.logspace(-5, math.log10(wmax), n-1)))
 40    best=0.
 41    I=np.eye(A.shape[0])
 42    for w in ws:
 43        H = (C @ np.linalg.solve(1j*w*I-A,B) + D)[0,0]
 44        best=max(best,abs(H))
 45    return float(best)
 46
 47def lmi_margin(A,B,C,D,P,gamma):
 48    """Standard bounded-real LMI: dV + ||p||^2 - gamma ||r||^2 <= 0.
 49    Here gamma is the squared gain bound and P is constrained positive definite.
 50    """
 51    n=A.shape[0]
 52    M=np.block([[A.T@P+P@A + C.T@C, P@B + C.T@D],
 53                [B.T@P + D.T@C, D.T@D - gamma*np.eye(D.shape[1])]])
 54    return float(np.max(np.linalg.eigvalsh(M)))
 55
 56def lmi_feasible(A,B,C,D,gamma):
 57    n=A.shape[0]
 58    # Cholesky-like parametrization P=L L' + eps I, optimized against largest eigenvalue.
 59    def unpack(v):
 60        L=np.zeros((n,n)); k=0
 61        for i in range(n):
 62            for j in range(i+1):
 63                L[i,j]=v[k]; k+=1
 64        return L@L.T + 1e-7*np.eye(n)
 65    x0=np.zeros(n*(n+1)//2)
 66    for i in range(n): x0[i*(i+1)//2+i]=1.
 67    fun=lambda v: max(0., lmi_margin(A,B,C,D,unpack(v),gamma))**2 + 1e-8*np.dot(v,v)
 68    res=minimize(fun,x0,method='Nelder-Mead',options={'maxiter':1000,'xatol':1e-9,'fatol':1e-12})
 69    P=unpack(res.x); m=lmi_margin(A,B,C,D,P,gamma)
 70    return bool(m <= -1e-6 and np.min(np.linalg.eigvalsh(P)) >= 1e-6), float(m), P
 71
 72def lmi_gamma(A,B,C,D):
 73    lo=max(1e-8,hinf(A,B,C,D)**2*0.95); hi=max(1.,hinf(A,B,C,D)**2*2+1e-3)
 74    if not np.isfinite(lo): return float('inf'), float('inf')
 75    # Feasibility is monotone in gamma for this bounded-real form.
 76    for _ in range(22):
 77        mid=(lo+hi)/2
 78        ok,_,_=lmi_feasible(A,B,C,D,mid)
 79        if ok: hi=mid
 80        else: lo=mid
 81    ok,m,_=lmi_feasible(A,B,C,D,hi*1.001)
 82    return float(hi),float(m)
 83
 84def run():
 85    np.random.seed(0); a=1.; b=1.4
 86    hc=critical_delay(a,b)
 87    # Prediction 1: direct characteristic roots cross at hc.
 88    delays=np.array([.5*hc,.8*hc,.95*hc,.99*hc,1.01*hc,1.2*hc,1.5*hc])
 89    root_rows=[{'h':float(h),'stable':bool(exact_stable(a,b,float(h))),'max_real':float(np.max(exact_roots(a,b,float(h)).real))} for h in delays]
 90    below=min(delays,key=lambda h: abs(h-.95*hc)); above=min(delays,key=lambda h: abs(h-1.05*hc))
 91    # Prediction 2: h=0 gain is exactly 1/(a+b), and certified gamma is its square.
 92    zero=[]
 93    for bb in [.2,.5,1.,2.,4.]:
 94        A,B,C,D=pade_system(a,bb,0.)
 95        g=hinf(A,B,C,D); zero.append({'b':bb,'measured_gain':g,'predicted_gain':1/(a+bb),'relative_error':abs(g-1/(a+bb))/(1/(a+bb))})
 96    # Prediction 3: gain increases toward the boundary and Padé remains stable below it.
 97    near=[]
 98    for h in [0.,.2,.4,.6,.8,.95*hc,1.02*hc]:
 99        A,B,C,D=pade_system(a,b,h); stable=bool(np.max(eigvals(A).real)<0)
100        g=hinf(A,B,C,D); gam= float('inf') if not np.isfinite(g) else g*g
101        near.append({'h':h,'exact_stable':exact_stable(a,b,h),'pade_stable':stable,'gain':g,'gamma_hinf':gam,'lmi_gamma':(None if not stable else lmi_gamma(A,B,C,D)[0])})
102    # Secondary same-system comparison: unconstrained near-boundary versus certified lower-delay block.
103    A1,B1,C1,D1=pade_system(a,b,.95*hc); A2,B2,C2,D2=pade_system(a,b,.35*hc)
104    comparison={'unconstrained_delay':.95*hc,'certified_delay':.35*hc,
105                'unconstrained_gain':hinf(A1,B1,C1,D1),'certified_gain':hinf(A2,B2,C2,D2),
106                'unconstrained_stable':bool(np.max(eigvals(A1).real)<0),'certified_stable':bool(np.max(eigvals(A2).real)<0)}
107    out={'parameters':{'a':a,'b':b,'predicted_critical_delay':hc},'prediction_1_boundary':{'rows':root_rows,'below':below,'above':above,'below_stable':exact_stable(a,b,below),'above_stable':exact_stable(a,b,above)},'prediction_2_zero_delay_scaling':zero,'prediction_3_delay_gain':near,'secondary_comparison':comparison}
108    Path('results.json').write_text(json.dumps(out,indent=2, default=lambda x: x.item() if hasattr(x, 'item') else str(x)))
109    print(json.dumps(out,indent=2, default=lambda x: x.item() if hasattr(x, 'item') else str(x)))
110
111if __name__=='__main__': run()