import json import math import random from dataclasses import dataclass import numpy as np SEED = 2753 np.random.seed(SEED) random.seed(SEED) @dataclass class AdaptiveDisturbance: c: np.ndarray q: np.ndarray alpha: float = 0.25 beta: float = 0.10 eps: float = 0.005 def update(self, residuals): residuals = np.asarray(residuals, dtype=float).reshape(-1, self.c.size) med = np.median(residuals, axis=0) # The specified update uses the newly updated center in the radius. self.c = (1.0 - self.alpha) * self.c + self.alpha * med observed = np.max(np.abs(residuals - self.c), axis=0) + self.eps self.q = np.maximum((1.0 - self.beta) * self.q, observed) return self.c.copy(), self.q.copy() def propagate_zonotope(A, B, x_center, x_generators, u, d, w_center, w_generators): """Exact affine zonotope propagation in generator form.""" c = A @ x_center + B @ u + d + w_center G = np.concatenate((A @ x_generators, w_generators), axis=1) return c, G def box_contains(c, G, lower, upper): radius = np.sum(np.abs(G), axis=1) return bool(np.all(c - radius >= lower - 1e-12) and np.all(c + radius <= upper + 1e-12)) def interval_prediction_check(): # For an affine zonotope and axis-aligned box, support in coordinate j is # c_j +/- sum_k |G_jk|; this is both an over-approximation and exact. rng = np.random.default_rng(SEED) errors = [] disagreements = 0 for _ in range(3000): n, p = 3, 7 c = rng.normal(size=n) G = rng.normal(size=(n, p)) lo = c - np.sum(np.abs(G), axis=1) - rng.uniform(0, 1, n) hi = c + np.sum(np.abs(G), axis=1) + rng.uniform(0, 1, n) predicted = np.maximum(np.abs(G), 0).sum(axis=1) # Enumerating vertices is exact for each coordinate support; sample # random signs as an independent numerical check. signs = rng.choice([-1.0, 1.0], size=(1000, p)) sampled = np.max(c[None, :] + signs @ G.T, axis=0) sampled_min = np.min(c[None, :] + signs @ G.T, axis=0) errors.append(float(np.max(np.maximum(sampled - (c + predicted), (c - predicted) - sampled_min)))) a = box_contains(c, G, lo, hi) b = bool(np.all(c - predicted >= lo - 1e-12) and np.all(c + predicted <= hi + 1e-12)) disagreements += int(a != b) return {"max_sample_support_gap": max(errors), "containment_disagreements": disagreements, "cases": 3000} def stationary_scaling_sweep(): eps = 0.005 amp_list = [0.02, 0.05, 0.10, 0.20, 0.40] rows = [] for amp in amp_list: est = AdaptiveDisturbance(np.zeros(1), np.zeros(1), alpha=0.25, beta=0.10, eps=eps) rng = np.random.default_rng(SEED + int(1000 * amp)) # Batch maxima make the intended q update visible without assuming # an unobserved true bound. for _ in range(40): residuals = rng.uniform(-amp, amp, size=(64, 1)) est.update(residuals) heldout = rng.uniform(-amp, amp, size=(10000, 1)) covered = np.abs(heldout[:, 0] - est.c[0]) <= est.q[0] predicted = amp + eps rows.append({"amplitude": amp, "predicted_q": predicted, "observed_q": float(est.q[0]), "q_error": float(est.q[0] - predicted), "coverage": float(np.mean(covered))}) max_error = max(abs(r["q_error"]) for r in rows) # With finite batches q can be below amp+eps, but should remain close and # coverage should be high; report rather than silently changing the rule. return {"rows": rows, "max_abs_q_prediction_error": max_error, "minimum_coverage": min(r["coverage"] for r in rows)} def jump_recovery_check(): eps = 0.005 est = AdaptiveDisturbance(np.zeros(1), np.zeros(1), alpha=0.25, beta=0.10, eps=eps) rng = np.random.default_rng(SEED + 99) for _ in range(12): est.update(rng.uniform(-0.10, 0.10, size=(64, 1))) before = float(est.q[0]) new_amp = 0.30 residuals = rng.uniform(-new_amp, new_amp, size=(64, 1)) est.update(residuals) after = float(est.q[0]) # The max rule predicts expansion in the first batch containing the jump. coverage = float(np.mean(np.abs(residuals[:, 0] - est.c[0]) <= est.q[0])) return {"old_amplitude": 0.10, "new_amplitude": new_amp, "q_before_jump": before, "q_after_one_batch": after, "predicted_min_after_batch": new_amp + eps - 0.03, "jump_batch_coverage": coverage, "expanded_in_one_batch": after > before} def point_mass_comparison(episodes=250, horizon=40): """Tiny shield comparison; rejected proposals use zero acceleration backup.""" rng = np.random.default_rng(SEED + 7) dt = 1.0 pos_limit, vel_limit, u_limit = 1.0, 1.5, 0.45 true_amp = 0.12 fixed_q = np.array([0.0, true_amp + 0.005]) stats = {"fixed": [0, 0], "adaptive": [0, 0]} # rejects, violations for method in stats: for ep in range(episodes): x = np.array([0.0, 0.0]) est = AdaptiveDisturbance(np.zeros(2), np.array([0.0, 0.35]), alpha=.25, beta=.10, eps=.005) history = [] for t in range(horizon): # Mildly aggressive policy, making the shield relevant. proposed = np.clip(-0.55*x[0] - 0.25*x[1] + rng.normal(0, .08), -u_limit, u_limit) A = np.array([[1., dt], [0., 1.]]) B = np.array([[0.], [dt]]) # Current point state, disturbance zonotope in acceleration. q = fixed_q if method == "fixed" else est.q wc = np.zeros(2) if method == "fixed" else est.c W = np.diag(np.asarray(q, dtype=float).reshape(-1)) c, G = propagate_zonotope(A, B, x, np.zeros((2, 0)), np.array([proposed]), np.zeros(2), wc, W) safe = box_contains(c, G, np.array([-pos_limit, -vel_limit]), np.array([pos_limit, vel_limit])) u = proposed if safe else 0.0 stats[method][0] += int(not safe) w = np.array([0., rng.uniform(-true_amp, true_amp)]) xn = (A @ np.asarray(x).reshape(2) + B[:, 0] * float(u) + w).reshape(2) stats[method][1] += int(abs(xn[0]) > pos_limit or abs(xn[1]) > vel_limit) if method == "adaptive": residual = xn - (A @ x + B * u) history.append(residual) if len(history) >= 8: est.update(np.array(history[-8:])) x = xn denom = episodes * horizon return {m: {"rejection_rate": v[0] / denom, "violation_rate": v[1] / denom} for m, v in stats.items()} def main(): result = {"seed": SEED, "math_check": interval_prediction_check(), "stationary_scaling": stationary_scaling_sweep(), "jump_recovery": jump_recovery_check(), "shield_comparison": point_mass_comparison()} print(json.dumps(result, indent=2)) if __name__ == "__main__": main()