1import json
  2import numpy as np
  3from pathlib import Path
  4
  5
  6def transport_queue_check(seed=0):
  7    rng = np.random.default_rng(seed)
  8    d = 7
  9    u = rng.normal(size=40)
 10    # Discrete upwind transport: q_{k+1}[s]=q_k[s-1], q_k[0]=u_k.
 11    q = np.zeros(d)
 12    observed = []
 13    for k, uk in enumerate(u):
 14        q[1:] = q[:-1].copy()
 15        q[0] = uk
 16        observed.append(q[-1])
 17    expected = [0.0] * (d - 1) + list(u[:-(d - 1)])
 18    err = float(np.max(np.abs(np.asarray(observed) - np.asarray(expected))))
 19    return {"delay": d, "max_abs_error": err}
 20
 21
 22def closed_form_predictor(A, B, x, queued):
 23    d = len(queued)
 24    out = np.linalg.matrix_power(A, d) @ x
 25    for j, uj in enumerate(queued):
 26        out += np.linalg.matrix_power(A, d - 1 - j) @ B @ uj
 27    return out
 28
 29
 30def predictor(A, B, x, queued):
 31    """Compute A^d x + sum A^(d-1-j) B queued[j]."""
 32    p = np.asarray(x, dtype=float).copy()
 33    for uj in queued:  # queued is chronological: oldest input first
 34        p = A @ p + B @ np.asarray(uj, dtype=float)
 35    return p
 36
 37
 38def simulate(A, B, K, d, steps=250, seed=123, compensated=True, noise_std=0.0):
 39    rng = np.random.default_rng(seed)
 40    n, m = B.shape
 41    x = np.array([1.0, -0.7], dtype=float)
 42    # pipeline entries are controls which will be applied over the next d steps
 43    queue = [np.zeros(m) for _ in range(d)]
 44    norms = []
 45    controls = []
 46    for k in range(steps):
 47        queued = [q.copy() for q in queue]
 48        applied = queue.pop(0)
 49        if compensated:
 50            target = predictor(A, B, x, queued)
 51        else:
 52            target = x
 53        u = K @ target
 54        queue.append(np.asarray(u).copy())
 55        x = A @ x + B @ applied + noise_std * rng.normal(size=n)
 56        norms.append(float(np.linalg.norm(x)))
 57        controls.append(float(np.linalg.norm(u)))
 58    return {
 59        "rms_last_50": float(np.sqrt(np.mean(np.square(norms[-50:])))),
 60        "max_norm": float(np.max(norms)),
 61        "final_norm": float(norms[-1]),
 62        "norms": norms,
 63        "controls": controls,
 64    }
 65
 66
 67def lyapunov_margin(F):
 68    # P from the discrete Lyapunov series; this is a constructive numerical check.
 69    n = F.shape[0]
 70    P = np.zeros((n, n))
 71    term = np.eye(n)
 72    for _ in range(10000):
 73        P += term
 74        term = F.T @ term @ F
 75        if np.linalg.norm(term) < 1e-12:
 76            break
 77    eig = np.linalg.eigvalsh(P)
 78    residual = F.T @ P @ F - P
 79    # For Q=I, residual should be -I up to truncation.
 80    return {"rho": float(np.max(np.abs(np.linalg.eigvals(F)))),
 81            "P_min_eig": float(np.min(eig)),
 82            "lyapunov_max_eig": float(np.max(np.linalg.eigvalsh(residual))),
 83            "series_terms": _}
 84
 85
 86def main():
 87    # A mildly unstable plant and a gain with a stable no-delay target.
 88    A = np.array([[1.08, 0.12], [0.0, 1.03]])
 89    B = np.array([[0.0], [1.0]])
 90    K = np.array([[-0.72, -0.48]])
 91    F = A + B @ K
 92    rng = np.random.default_rng(9)
 93    Aq = np.array([[0.91, 0.13], [-0.04, 0.86]])
 94    Bq = np.array([[0.2], [0.7]])
 95    xq = rng.normal(size=2)
 96    uq = [rng.normal(size=1) for _ in range(5)]
 97    algebraic_error = float(np.max(np.abs(predictor(Aq, Bq, xq, uq) - closed_form_predictor(Aq, Bq, xq, uq))))
 98    out = {"transport_check": transport_queue_check(), "predictor_formula_check": {"max_abs_error": algebraic_error}, "target": lyapunov_margin(F), "runs": []}
 99    # Equal setup, fixed seeds; baseline is stale-state feedback.
100    for d in [1, 2, 4, 8, 16, 32]:
101        base = simulate(A, B, K, d, compensated=False)
102        pred = simulate(A, B, K, d, compensated=True)
103        out["runs"].append({"d": d, "baseline": {k: v for k, v in base.items() if k not in ("norms", "controls")},
104                            "predictor": {k: v for k, v in pred.items() if k not in ("norms", "controls")}})
105    # Noise sensitivity is reported at a representative long delay.
106    out["noise"] = []
107    for sigma in [0.0, 0.005, 0.02]:
108        b = simulate(A, B, K, 16, compensated=False, noise_std=sigma)
109        p = simulate(A, B, K, 16, compensated=True, noise_std=sigma)
110        out["noise"].append({"sigma": sigma, "baseline_rms_last_50": b["rms_last_50"],
111                             "predictor_rms_last_50": p["rms_last_50"]})
112    Path("results.json").write_text(json.dumps(out, indent=2))
113    print(json.dumps(out, indent=2))
114
115
116if __name__ == "__main__":
117    main()