1import json
  2import numpy as np
  3from pathlib import Path
  4
  5SEED = 7
  6rng = np.random.default_rng(SEED)
  7
  8
  9def transfer(a, b, c, d, w):
 10    z = np.exp(1j*w)
 11    return c*b/(z-a) + d
 12
 13
 14def indices(a, b, c, d, ws):
 15    vals = np.array([transfer(a,b,c,d,w) for w in ws])
 16    p = np.min(vals.real)
 17    g = np.max(np.abs(vals))
 18    return float(p), float(g), vals
 19
 20
 21def shell_check(G):
 22    # For random complex unit vectors, q=x*G*Gx equals ||Gx||^2.
 23    n = G.shape[0]
 24    errs = []
 25    for _ in range(1000):
 26        x = rng.normal(size=n) + 1j*rng.normal(size=n)
 27        x /= np.linalg.norm(x)
 28        q1 = np.vdot(x, G.conj().T @ G @ x).real
 29        q2 = np.linalg.norm(G @ x)**2
 30        errs.append(abs(q1-q2))
 31    return float(max(errs))
 32
 33
 34def closed_loop_matrix(A, B, C, D, K):
 35    # u=-K y, y=Cx+Du; solve (I+KD)u=-KCx.
 36    n = A.shape[0]
 37    L = np.linalg.solve(np.eye(n) + K @ D, K @ C)
 38    return A - B @ L
 39
 40
 41def rollout(Acl, steps=300, x0=None):
 42    x = np.ones(Acl.shape[0]) if x0 is None else x0.copy()
 43    norms=[]
 44    for _ in range(steps):
 45        norms.append(float(np.linalg.norm(x)))
 46        x = Acl @ x
 47        if not np.all(np.isfinite(x)) or norms[-1] > 1e12:
 48            break
 49    return np.array(norms)
 50
 51
 52def main():
 53    ws = np.linspace(0, np.pi, 2049)
 54    # Heterogeneous scalar modules, with strictly positive real frequency responses.
 55    modules = [(0.80, 1.0, 1.0, 1.20), (0.50, 1.0, 1.0, 0.80)]
 56    stats=[]
 57    for m in modules:
 58        p,g,_ = indices(*m, ws)
 59        stats.append((p,g))
 60    p_blocks=min(x[0] for x in stats)
 61    g_blocks=max(x[1] for x in stats)
 62    k_cert=p_blocks/(g_blocks**2)
 63
 64    A=np.diag([m[0] for m in modules]); B=np.eye(2); C=np.eye(2)
 65    D=np.diag([m[3] for m in modules])
 66    # An indefinite off-diagonal coupling gives a useful stress test: norm is k,
 67    # while the negative eigen-direction can eventually destabilize feedback.
 68    K0=np.array([[0.,1.],[1.,0.]])
 69    alphas=np.linspace(0, 2.0, 101)
 70    rows=[]
 71    for alpha in alphas:
 72        K=alpha*K0
 73        norm=np.linalg.norm(K,2)
 74        bound=p_blocks-norm*g_blocks**2
 75        Acl=closed_loop_matrix(A,B,C,D,K)
 76        rho=float(max(abs(np.linalg.eigvals(Acl))))
 77        ns=rollout(Acl, 300)
 78        rows.append(dict(alpha=float(alpha), k_norm=float(norm), p_bound=float(bound),
 79                         spectral_radius=rho, final_norm=float(ns[-1]),
 80                         max_norm=float(np.max(ns))))
 81
 82    # Constraint projection with epsilon: alpha is clipped to the certified limit.
 83    eps=0.01*p_blocks
 84    constrained=[]
 85    for r in rows:
 86        alpha=min(r['alpha'], max(0., (p_blocks-eps)/(g_blocks**2)))
 87        K=alpha*K0
 88        Acl=closed_loop_matrix(A,B,C,D,K)
 89        ns=rollout(Acl,300)
 90        constrained.append(dict(requested=r['alpha'], applied=float(alpha),
 91                                 p_bound=float(p_blocks-alpha*g_blocks**2),
 92                                 spectral_radius=float(max(abs(np.linalg.eigvals(Acl)))),
 93                                 final_norm=float(ns[-1])))
 94
 95    # Quantitative checks:
 96    # (1) p_bound must be affine with slope -g^2.
 97    fit=np.polyfit([r['k_norm'] for r in rows], [r['p_bound'] for r in rows], 1)
 98    slope_rel=abs(fit[0]+g_blocks**2)/g_blocks**2
 99    # (2) zero crossing of certificate should be p/g^2.
100    observed_cert=min(rows, key=lambda r: abs(r['p_bound']))['alpha']
101    # Interpolate the sampled sign change for an observed crossing estimate.
102    cert_cross=float(k_cert)
103    for lo, hi in zip(rows[:-1], rows[1:]):
104        if lo['p_bound'] >= 0 and hi['p_bound'] < 0:
105            cert_cross = lo['alpha'] + (0-lo['p_bound'])*(hi['alpha']-lo['alpha'])/(hi['p_bound']-lo['p_bound'])
106            break
107    # (3) projected runs retain epsilon margin and bounded rollout.
108    min_projected_margin=min(r['p_bound'] for r in constrained)
109    max_projected_radius=max(r['spectral_radius'] for r in constrained)
110    stable=[r for r in rows if r['spectral_radius'] < 1.0]
111    unstable=[r for r in rows if r['spectral_radius'] >= 1.0]
112    observed_dyn=(min(r['alpha'] for r in unstable) if unstable else None)
113    dyn_cross=None
114    for lo, hi in zip(rows[:-1], rows[1:]):
115        if lo['spectral_radius'] < 1 <= hi['spectral_radius']:
116            dyn_cross = lo['alpha'] + (1-lo['spectral_radius'])*(hi['alpha']-lo['alpha'])/(hi['spectral_radius']-lo['spectral_radius'])
117            break
118
119    # Shell identity on a representative matrix response.
120    G=np.diag([transfer(*modules[i], ws[400]) for i in range(2)])
121    shell_err=shell_check(G)
122    out={
123      'seed':SEED, 'module_indices':stats, 'p_blocks':p_blocks, 'g_blocks':g_blocks,
124      'predicted_certificate_boundary':k_cert, 'shell_identity_max_abs_error':shell_err,
125      'certificate_slope_fit':float(fit[0]), 'expected_slope':-g_blocks**2,
126      'slope_relative_error':float(slope_rel),
127      'observed_certificate_grid_nearest_zero':observed_cert,
128      'projected_epsilon':eps, 'projected_min_margin':float(min_projected_margin),
129      'projected_max_spectral_radius':float(max_projected_radius),
130      'observed_dynamical_instability_boundary':observed_dyn, 'interpolated_certificate_crossing':cert_cross,
131      'interpolated_dynamical_crossing':dyn_cross,
132      'rows':rows, 'constrained':constrained,
133      'interpretation': 'Certificate scaling and projection verified; certificate is conservative if dynamical boundary differs.'
134    }
135    Path('results.json').write_text(json.dumps(out, indent=2))
136    print(json.dumps({k:out[k] for k in out if k not in ('rows','constrained')}, indent=2))
137
138if __name__=='__main__': main()