1import json
  2import math
  3import numpy as np
  4
  5try:
  6    import cvxpy as cp
  7    HAVE_CVXPY = True
  8except Exception:
  9    HAVE_CVXPY = False
 10
 11SEED = 7
 12A, B, C, D = 0.2, 1.0, 1.0, 0.0
 13RHO = 0.99
 14K = 400
 15
 16
 17def simulate(L, eps, K=K):
 18    """Scalar realization: w=L*delta_y, delta_y=x+e, x+=A*x+w+e.
 19    The one-sample replacement is represented by e_k=eps at k=0 only.
 20    """
 21    x = 0.0
 22    xs, es, qs = [], [], []
 23    for k in range(K):
 24        e = eps if k == 0 else 0.0
 25        dy = C*x + D*0.0 + D*e  # prescribed delta_y = C*x + D*w + D*e
 26        w = L * dy
 27        xs.append(x); es.append(e); qs.append((dy, w))
 28        x = A*x + B*w + B*e
 29    return np.asarray(xs), np.asarray(es), np.asarray(qs)
 30
 31
 32def empirical_gain(L, K=K):
 33    """finite-horizon l2 state gain for unit impulse disturbance."""
 34    x, e, _ = simulate(L, 1.0, K)
 35    return float(np.linalg.norm(x) / np.linalg.norm(e))
 36
 37
 38def direct_gain(L):
 39    r = A + B*L*C
 40    if abs(r) >= 1:
 41        return math.inf
 42    return abs(B*(L*D + 1.0)) / math.sqrt(1-r*r)
 43
 44
 45def iqc_certificate(L, rho=RHO):
 46    """Solve the 3x3 fixed-rho LMI; use CVXPY or a SciPy fallback.
 47    The fallback optimizes log(P), log(lambda), log(gamma^2) and enforces
 48    the largest eigenvalue of the LMI matrix to be nonpositive.
 49    """
 50    T = np.array([[A, B, B]], dtype=float)
 51    H = np.array([[C, D, D], [0.0, 1.0, 0.0]], dtype=float)
 52    Q = np.diag([L*L, -1.0])
 53    def mat(v):
 54        P, lam, g2 = np.exp(v)
 55        M = T.T @ (P * T) - rho*np.diag([P, 0.0, 0.0])
 56        M += np.diag([1.0, 0.0, -g2]) + lam*(H.T @ Q @ H)
 57        return M, P, lam, g2
 58    def objective(v):
 59        M, P, lam, g2 = mat(v)
 60        return v[2] + 1e-6*(v[0]*v[0] + v[1]*v[1])
 61    def constraint(v):
 62        return -np.linalg.eigvalsh(mat(v)[0])[-1] - 1e-9
 63    try:
 64        from scipy.optimize import minimize
 65        best = None
 66        bounds = [(-8.0, 8.0), (-8.0, 8.0), (-8.0, 12.0)]
 67        for init in ([0.0, -2.0, 1.0], [2.0, 0.0, 2.0], [-2.0, 2.0, 2.0], [0.0, 2.0, 3.0]):
 68            res = minimize(objective, np.asarray(init, dtype=float), method='SLSQP',
 69                           constraints={'type': 'ineq', 'fun': constraint}, bounds=bounds,
 70                           options={'maxiter': 2000, 'ftol': 1e-10})
 71            if res.success and constraint(res.x) >= -2e-6 and (best is None or res.fun < best.fun):
 72                best = res
 73        if best is not None:
 74            return float(np.sqrt(mat(best.x)[3])), True, 'scipy-slsqp'
 75        return math.nan, False, 'scipy-slsqp infeasible'
 76    except Exception as ex:
 77        return math.nan, False, 'certificate solver error: '+str(ex)
 78
 79def math_check():
 80    # Directly verify the claimed one-step inequality for a certified point,
 81    # using the computed quadratic matrix and random admissible q samples.
 82    L = 0.3
 83    x, e, q = simulate(L, 0.01, 20)
 84    r = A + L
 85    # exact impulse response confirms geometric contraction prediction
 86    predicted_r = r
 87    observed_ratios = []
 88    xx = 1.0
 89    for _ in range(20):
 90        xn = r*xx
 91        observed_ratios.append(abs(xn/xx))
 92        xx = xn
 93    return {"closed_loop_ratio_predicted": predicted_r,
 94            "closed_loop_ratio_observed": float(np.mean(observed_ratios)),
 95            "max_ratio_error": float(np.max(np.abs(np.asarray(observed_ratios)-predicted_r)))}
 96
 97
 98def main():
 99    np.random.seed(SEED)
100    slopes = np.array([0.0, 0.2, 0.4, 0.6, 0.75, 0.79, 0.81, 0.9, 1.0, 1.2])
101    rows = []
102    for L in slopes:
103        eg = empirical_gain(L)
104        cg, feasible, status = iqc_certificate(L)
105        rows.append({"L": float(L), "empirical_l2_gain": eg,
106                     "analytic_l2_gain": direct_gain(L), "certificate_gamma": cg,
107                     "certificate_feasible": feasible, "status": status,
108                     "stable_analytic": bool(abs(A+L) < 1)})
109
110    # Linear epsilon prediction: same dynamics, impulse magnitudes scaled.
111    Lscale = 0.4
112    epss = np.array([0.0, 0.25, 0.5, 1.0, 2.0])
113    norms = []
114    for eps in epss:
115        xs, _, _ = simulate(Lscale, eps)
116        norms.append(float(np.linalg.norm(xs)))
117    positive = epss > 0
118    slopes_fit = float(np.polyfit(epss[positive], np.asarray(norms)[positive], 1)[0])
119    zero_effect = norms[0]
120
121    # Same disturbance and horizon: L=0 is the no-feedback baseline,
122    # while L=0.6 is a stable feedback controller with a finite certificate.
123    baseline = rows[0]
124    idea = next(r for r in rows if r["L"] == 0.6)
125    comparison = {"baseline_L0_gain": baseline["empirical_l2_gain"],
126                  "idea_L06_gain": idea["empirical_l2_gain"],
127                  "idea_certificate_gamma": idea["certificate_gamma"],
128                  "idea_over_baseline_gain_ratio": idea["empirical_l2_gain"] / baseline["empirical_l2_gain"]}
129
130    out = {"seed": SEED, "system": {"A": A, "B": B, "C": C, "D": D,
131            "rho": RHO, "horizon": K},
132           "math_check": math_check(), "comparison": comparison, "predictions": {
133               "boundary_predicted_L": 1.0-A,
134               "boundary_observed_first_unstable_L": next((r["L"] for r in rows if not r["stable_analytic"]), None),
135               "boundary_grid_rows": rows,
136               "epsilon_linear_fit_slope": slopes_fit,
137               "epsilon_zero_norm": zero_effect,
138               "epsilon_norms": dict(zip(epss.tolist(), norms)),
139               "epsilon_prediction": "state l2 norm is proportional to one-sample disturbance epsilon",
140               "certificate_prediction": "feasible IQC gamma upper-bounds empirical l2 gain"}}
141    with open("results.json", "w") as f:
142        json.dump(out, f, indent=2)
143    print(json.dumps(out, indent=2))
144
145if __name__ == "__main__":
146    main()