1import json
  2import numpy as np
  3from scipy.optimize import minimize
  4
  5SEED = 1436
  6rng = np.random.default_rng(SEED)
  7
  8
  9def sector_matrix(alpha=0.0, beta=1.0):
 10    return np.array([[-alpha * beta, 0.5 * (alpha + beta)],
 11                     [0.5 * (alpha + beta), -1.0]])
 12
 13
 14def solve_certificate(A, B1, B2, C1, D11, D12, C2, D21, D22,
 15                      alpha=-1.0, beta=1.0, eps=1e-7):
 16    """Fast scalar certificate search: exact LMI eigenvalues on a log grid."""
 17    X=np.array([[float(np.asarray(A).ravel()[0]),float(B1),float(B2)]])
 18    V=np.array([[float(C1),float(D11),float(D12)]])
 19    E=np.array([[float(C2),float(D21),float(D22)]])
 20    W=np.array([[0.,1.,0.]])
 21    Q=np.vstack([V,W]).T @ sector_matrix(alpha,beta) @ np.vstack([V,W])
 22    best=None
 23    # Homogeneity permits P=1 in the scalar state dimension; scan multiplier and gain.
 24    for P in np.logspace(-3,3,9):
 25     for lam in np.r_[0., np.logspace(-2,2,8)]:
 26      for g2 in np.r_[np.logspace(-2,4,15)]:
 27        L=X.T@(P*X)-np.diag([P,0.,0.])+E.T@E-g2*np.diag([0.,0.,1.])+lam*Q
 28        ev=float(np.linalg.eigvalsh(L).max())
 29        if ev < -eps and (best is None or g2 < best[0]): best=(g2,lam,ev,P)
 30    if best is None: return False,np.inf,np.nan,np.nan,np.inf
 31    return True,float(np.sqrt(best[0])),float(best[3]),float(best[1]),float(best[2])
 32
 33
 34def collect_and_identify(a, b, n=5000, noise=0.0):
 35    """Exciting trajectory and least-squares reconstruction of [A B1 B2], [C1 D11 D12]."""
 36    u = rng.normal(0.0, 0.35, n)
 37    x = np.zeros(n + 1); w = np.zeros(n); v = np.zeros(n)
 38    for k in range(n):
 39        v[k] = x[k]
 40        w[k] = np.tanh(v[k])
 41        x[k + 1] = a*x[k] + b*w[k] + u[k] + noise*rng.normal()
 42    # exact hidden state is available in this controlled verification; identify linear lifted model
 43    Z = np.column_stack([x[:-1], w, u])
 44    Xhat = np.linalg.lstsq(Z, x[1:], rcond=None)[0]
 45    Vhat = np.linalg.lstsq(Z, v, rcond=None)[0]
 46    return Xhat, Vhat, (x, u, v, w)
 47
 48
 49def nonlinear_metrics(a, b, n=250, delta=1e-5):
 50    """Local perturbation decay and finite-horizon impulse gain."""
 51    x1 = delta; x2 = -delta
 52    diffs = []
 53    for _ in range(n):
 54        diffs.append(abs(x1-x2))
 55        x1 = a*x1 + b*np.tanh(x1)
 56        x2 = a*x2 + b*np.tanh(x2)
 57    ratio = (diffs[-1] / diffs[0]) ** (1.0 / max(n-1, 1))
 58    # disturbance impulse, zero initial state; output is x
 59    x = 0.0; ss = 0.0
 60    for k in range(n):
 61        d = 1.0 if k == 0 else 0.0
 62        x = a*x + b*np.tanh(x) + d
 63        ss += x*x
 64    return ratio, np.sqrt(ss)
 65
 66
 67def main():
 68    # b is fixed and a is swept. Tanh has sector [-1,1]; robust sector boundary is |a|+b=1.
 69    b = 0.45
 70    a_values = np.array([-1.45,-1.1,-0.8,-0.5,-0.2,0.0,0.2,0.4,0.5,0.6,0.75])
 71    rows = []
 72    for a in a_values:
 73        cert = solve_certificate(a, b, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0)
 74        identX, identV, _ = collect_and_identify(a, b, n=500)
 75        data_cert = solve_certificate(identX[0], identX[1], identX[2],
 76                                      identV[0], identV[1], identV[2],
 77                                      1.0, 0.0, 0.0)
 78        rho, impulse = nonlinear_metrics(a, b)
 79        rows.append({"a": float(a), "local_rho_pred": abs(a+b),
 80                     "measured_rho": float(rho), "impulse_l2": float(impulse),
 81                     "cert": bool(cert[0]), "gamma": float(cert[1]),
 82                     "data_cert": bool(data_cert[0]), "data_gamma": float(data_cert[1]),
 83                     "lmi_maxeig": float(cert[4])})
 84
 85    # A finer boundary search reports the first certified stable point and compares it
 86    # with the analytical local and sector worst-case boundaries.
 87    fine = np.linspace(-1.2, 0.7, 10)
 88    fine_cert = [solve_certificate(a, b, 1.0, 1.0, 0, 0, 1, 0, 0)[0] for a in fine]
 89    cert_stable = [a for a, ok in zip(fine, fine_cert) if ok]
 90    observed_boundary = max([a for a in cert_stable if a >= 0.0], default=np.nan)
 91    # For tanh sector [-1,1], the positive-a robust boundary is a+b=1.
 92    predicted_boundary = 1.0 - b
 93    # Fit gain scaling in stable region against exact linearized Hinf prediction 1/(1-rho).
 94    stable_rows = [r for r in rows if r["local_rho_pred"] < 0.92 and r["cert"]]
 95    for r in stable_rows:
 96        r["linear_gain_pred"] = 1.0 / (1.0 - r["local_rho_pred"])
 97    out = {"seed": SEED, "b": b, "rows": rows,
 98           "prediction_boundary_a": predicted_boundary,
 99           "observed_certificate_boundary_a": float(observed_boundary),
100           "boundary_error": float(observed_boundary-predicted_boundary),
101           "positive_boundary_definition": "largest nonnegative certified a on the tested grid",
102           "stable_gain_rows": stable_rows,
103           "notes": [
104             "Prediction 1: local perturbations contract for |a+b|<1 and grow for |a+b|>1.",
105             "Prediction 2: tanh sector certificate boundary is a+b=1 on the positive-a side.",
106             "Prediction 3: disturbance gain grows as 1/(1-(a+b)) near the boundary."
107           ]}
108    with open("results.json", "w") as f: json.dump(out, f, indent=2)
109    print(json.dumps({"predicted_boundary": predicted_boundary,
110                      "observed_certificate_boundary": observed_boundary,
111                      "boundary_error": observed_boundary-predicted_boundary,
112                      "rows": rows}, indent=2))
113
114if __name__ == "__main__":
115    main()