Finite-Candidate Neural Reference Shield / shield_experiment.py
Mechanism confirmed, baseline not beaten
1import json
2import itertools
3import time
4import numpy as np
5from scipy.optimize import minimize
6
7
8def constraints(r, R, B):
9 r = np.asarray(r)
10 return np.array([r[0]-B[0], -r[0]-B[0], r[1]-B[1], -r[1]-B[1], r@r-R*R])
11
12
13def shield(z, R, B, tol=1e-10):
14 """Exact Euclidean projection onto disk intersected with an axis-aligned box.
15 Candidates are the KKT active-set candidates for all subsets in 2D.
16 """
17 z = np.asarray(z, dtype=float)
18 candidates = [(z.copy(), ())]
19 # One box face: orthogonal projection onto the face.
20 faces = [(0, 1), (0, -1), (1, 1), (1, -1)]
21 for i, sign in faces:
22 rr = z.copy(); rr[i] = sign * B[i]
23 candidates.append((rr, (i, sign)))
24 # Disk alone (radial KKT candidate).
25 nz = np.linalg.norm(z)
26 if nz > R:
27 candidates.append((z * (R / nz), ('disk',)))
28 # Two box faces (corners).
29 for i, s1 in faces:
30 for j, s2 in faces:
31 if i < j:
32 rr = z.copy(); rr[i] = s1 * B[i]; rr[j] = s2 * B[j]
33 candidates.append((rr, (i, s1, j, s2)))
34 # Disk/box active intersections (the finite KKT roots).
35 for i, sign in faces:
36 c = sign * B[i]
37 rem = R*R-c*c
38 if rem >= -tol:
39 rem = max(0., rem)
40 for y in (np.sqrt(rem), -np.sqrt(rem)):
41 rr = np.array([0., 0.]); rr[i] = c; rr[1-i] = y
42 candidates.append((rr, (i, sign, 'disk')))
43 feasible = [(0.5*np.sum((rr-z)**2), rr, a)
44 for rr, a in candidates if np.max(constraints(rr, R, B)) <= tol]
45 if not feasible:
46 # Conservative fallback is the box center, which is feasible whenever the
47 # modeled actuator set is nonempty.
48 rr = np.zeros(2)
49 return rr, 0.5*np.sum((rr-z)**2), True
50 obj, rr, active = min(feasible, key=lambda t: t[0])
51 return rr, obj, False
52
53
54def pg_projection(z, R, B, steps=300, alpha=.08):
55 """A deliberately standard alternating projected-gradient feasibility repair."""
56 r = z.copy()
57 for _ in range(steps):
58 r = np.clip(r, -B, B)
59 n = np.linalg.norm(r)
60 if n > R: r *= R/n
61 return r
62
63
64def reference(z, R, B):
65 # SLSQP is only used for numerical verification, not the shield.
66 fun = lambda r: .5*np.sum((r-z)**2)
67 cons = ({'type':'ineq','fun':lambda r: R*R-r@r},
68 {'type':'ineq','fun':lambda r: B-r},
69 {'type':'ineq','fun':lambda r: B+r})
70 q = minimize(fun, np.clip(z, -B, B), constraints=cons, method='SLSQP',
71 options={'ftol':1e-12, 'maxiter':300})
72 return q.x, q.fun, q.success
73
74
75def run():
76 rng = np.random.default_rng(2906)
77 B = np.array([1.35, 1.10])
78 direction = np.array([.8, .6])
79 direction /= np.linalg.norm(direction)
80 # Prediction 1: every command inside the feasible set is unchanged.
81 interior_err = []
82 ref_err = []
83 max_violation = 0.
84 for _ in range(300):
85 R = rng.uniform(.65, 1.8)
86 z = rng.uniform(-1.5, 1.5, 2)
87 r, _, _ = shield(z, R, B)
88 if np.max(constraints(z, R, B)) <= 0:
89 interior_err.append(np.linalg.norm(r-z))
90 rr, _, ok = reference(z, R, B)
91 if ok: ref_err.append(np.linalg.norm(r-rr))
92 max_violation = max(max_violation, max(0., np.max(constraints(r, R, B))))
93 # Prediction 2: along a ray the transition is at min(box boundary, disk boundary).
94 R = 1.0
95 predicted_transition = min(B[0]/direction[0], B[1]/direction[1], R)
96 scales = np.linspace(.1, 1.7, 161)
97 errors = []
98 violations = []
99 for s in scales:
100 z = s*direction; r, _, _ = shield(z, R, B)
101 errors.append(np.linalg.norm(r-z)); violations.append(max(0., np.max(constraints(r,R,B))))
102 observed_transition = scales[np.argmax(np.asarray(errors) > 1e-8)]
103 # Prediction 3: changing disk radius shifts the transition linearly while disk is limiting.
104 radii = np.linspace(.45, 1.0, 12)
105 measured = []
106 for rad in radii:
107 es = []
108 for s in scales:
109 rr, _, _ = shield(s*direction, rad, B); es.append(np.linalg.norm(rr-s*direction))
110 measured.append(scales[np.argmax(np.asarray(es)>1e-8)])
111 slope = np.polyfit(radii[:8], measured[:8], 1)[0]
112 # Mini comparison on random commands: exact shield vs raw and iterative repair.
113 zs = rng.uniform(-2.0,2.0,(500,2)); Rs = rng.uniform(.65,1.25,500)
114 raw_v=[]; shield_v=[]; pg_v=[]; shield_dist=[]; pg_dist=[]
115 t0=time.perf_counter()
116 for z,rad in zip(zs,Rs):
117 rr,_,_=shield(z,rad,B); shield_v.append(max(0.,np.max(constraints(rr,rad,B))))
118 shield_dist.append(np.linalg.norm(rr-z)); raw_v.append(max(0.,np.max(constraints(z,rad,B))))
119 shield_ms=1000*(time.perf_counter()-t0)/len(zs)
120 t0=time.perf_counter()
121 for z,rad in zip(zs,Rs):
122 rr=pg_projection(z,rad,B); pg_v.append(max(0.,np.max(constraints(rr,rad,B)))); pg_dist.append(np.linalg.norm(rr-z))
123 pg_ms=1000*(time.perf_counter()-t0)/len(zs)
124 out={
125 'predictions': {
126 'interior_identity_max_error': float(max(interior_err)),
127 'predicted_identity_error': 0.0,
128 'transition_predicted': float(predicted_transition),
129 'transition_observed_grid': float(observed_transition),
130 'transition_abs_error': float(abs(observed_transition-predicted_transition)),
131 'radius_transition_slope_predicted': 1.0,
132 'radius_transition_slope_observed': float(slope),
133 },
134 'math_check': {'max_constraint_violation':float(max_violation), 'max_error_vs_SLSQP':float(max(ref_err)), 'n_reference':len(ref_err)},
135 'comparison': {
136 'raw_mean_violation':float(np.mean(raw_v)), 'shield_mean_violation':float(np.mean(shield_v)), 'pg_mean_violation':float(np.mean(pg_v)),
137 'shield_mean_distance':float(np.mean(shield_dist)), 'pg_mean_distance':float(np.mean(pg_dist)),
138 'shield_ms_per_command':shield_ms, 'pg_ms_per_command':pg_ms
139 }, 'seed':2906}
140 with open('results.json','w') as f: json.dump(out,f,indent=2)
141 print(json.dumps(out,indent=2))
142
143if __name__=='__main__': run()