import json, math, random from pathlib import Path import numpy as np # Feasible-action mapping safety layer for a 1D point mass. # x=[position, velocity], u=p is held constant over H predictions. class ConstantAccelerationSafety: def __init__(self, dt=0.2, H=5, u_min=-1.0, u_max=1.0, pos_min=-2.0, pos_max=2.0, vel_min=-2.0, vel_max=2.0, terminal_pos_min=-0.4, terminal_pos_max=0.4, terminal_vel_min=-0.5, terminal_vel_max=0.5): self.dt, self.H = dt, H self.u_min, self.u_max = u_min, u_max self.pos_min, self.pos_max = pos_min, pos_max self.vel_min, self.vel_max = vel_min, vel_max self.tpos_min, self.tpos_max = terminal_pos_min, terminal_pos_max self.tvel_min, self.tvel_max = terminal_vel_min, terminal_vel_max def interval(self, state): """Return exact feasible interval [lo,hi] for constant p.""" x0, v0 = map(float, state) lo, hi = self.u_min, self.u_max # x_k = x0 + k dt v0 + .5 (k dt)^2 p; v_k=v0+k dt p for k in range(1, self.H + 1): A = 0.5 * (k*self.dt)**2; B = x0 + k*self.dt*v0 for bound, is_lower in [(self.pos_min, True), (self.pos_max, False)]: if is_lower: if A > 0: lo = max(lo, (bound-B)/A) else: if A > 0: hi = min(hi, (bound-B)/A) A = k*self.dt; B = v0 for bound, is_lower in [(self.vel_min, True), (self.vel_max, False)]: if is_lower: if A > 0: lo = max(lo, (bound-B)/A) else: if A > 0: hi = min(hi, (bound-B)/A) # Explicit terminal set (at H; duplicate state constraints are harmless). k=self.H; A=0.5*(k*self.dt)**2; B=x0+k*self.dt*v0 lo=max(lo,(self.tpos_min-B)/A); hi=min(hi,(self.tpos_max-B)/A) A=k*self.dt; B=v0 lo=max(lo,(self.tvel_min-B)/A); hi=min(hi,(self.tvel_max-B)/A) return lo, hi def project(self, state, z): lo, hi = self.interval(state) if lo > hi: return None, float('inf'), (lo,hi) p = min(max(float(z), lo), hi) return p, abs(p-float(z)), (lo,hi) def rollout(self, state, p): x,v=map(float,state); traj=[] for _ in range(self.H): x=x+self.dt*v+0.5*self.dt*self.dt*p; v=v+self.dt*p traj.append((x,v,p)) return np.asarray(traj) def certified(self,state,p): if p is None: return False q=self.rollout(state,p) return bool(np.all((q[:,0]>=self.pos_min-1e-9)&(q[:,0]<=self.pos_max+1e-9)& (q[:,1]>=self.vel_min-1e-9)&(q[:,1]<=self.vel_max+1e-9)& (q[:,2]>=self.u_min-1e-9)&(q[:,2]<=self.u_max+1e-9)& (q[-1,0]>=self.tpos_min-1e-9)&(q[-1,0]<=self.tpos_max+1e-9)& (q[-1,1]>=self.tvel_min-1e-9)&(q[-1,1]<=self.tvel_max+1e-9))) def run(): np.random.seed(1341); random.seed(1341) sl=ConstantAccelerationSafety() out={"predictions":{},"simulation":{}} # Prediction 1: projection is identity throughout feasible interval. state=np.array([0.0,0.0]); lo,hi=sl.interval(state) zs=np.linspace(lo,hi,101); ds=np.array([sl.project(state,z)[1] for z in zs]) out["predictions"]["inside_identity"]={"predicted":"d=0 for z in [lo,hi]","observed_max_d":float(ds.max()),"interval":[lo,hi]} # Prediction 2: outside distance has slope one (weighted scalar projection). zleft=np.linspace(lo-3,lo-.05,40); zright=np.linspace(hi+.05,hi+3,40) dl=np.array([sl.project(state,z)[1] for z in zleft]); dr=np.array([sl.project(state,z)[1] for z in zright]) slope_l=np.polyfit(zleft,dl,1)[0]; slope_r=np.polyfit(zright,dr,1)[0] out["predictions"]["outside_linear_distance"]={"predicted":"d=|z-boundary|, slope=1","observed_left_slope":float(slope_l),"observed_right_slope":float(slope_r),"max_abs_slope_error":float(max(abs(slope_l+1),abs(slope_r-1)))} # Prediction 3: the active terminal-velocity boundary shifts linearly with current velocity. states=np.array([[0.,0.],[0.,0.1],[0.,-0.1],[0.,0.3],[0.,-0.3]]) bounds=np.array([sl.interval(s) for s in states]) # v_H=v0+H*dt*p and v_H <= terminal_vel_max, so hi shifts by -dv/(H dt). predicted_shift=-(states[:,1]-states[0,1])/(sl.H*sl.dt) observed_shift=bounds[:,1]-bounds[0,1] residual=observed_shift-predicted_shift out["predictions"]["state_dependent_boundary"]={"predicted":"upper boundary shift=-delta_v/(H dt)","predicted_shifts":predicted_shift.tolist(),"observed_hi_shifts":observed_shift.tolist(),"max_abs_residual":float(np.max(np.abs(residual)))} # Certified-regime random test: exact projection must satisfy every constraint. cert=0; feasible=0 for _ in range(10000): s=np.random.uniform([-0.35,-0.35],[0.35,0.35]); z=np.random.uniform(-4,4) p,d,b=sl.project(s,z) if p is not None: feasible+=1; cert += int(sl.certified(s,p)) out["predictions"]["zero_violation"]={"predicted":"0 violations whenever solver feasible","feasible_cases":feasible,"certified_fraction":cert/max(feasible,1)} # One-step stress test: compare coordinate clipping with the predictive map. clip_bad = safe_bad = 0 for _ in range(10000): s=np.random.uniform([-1.8,-1.8],[1.8,1.8]); z=np.random.uniform(-4,4) pc=float(np.clip(z,sl.u_min,sl.u_max)); q=sl.rollout(s,pc)[0] clip_bad += int(not (-2<=q[0]<=2 and -2<=q[1]<=2)) ps,_,_=sl.project(s,z) if ps is not None: q=sl.rollout(s,ps)[0] safe_bad += int(not (-2<=q[0]<=2 and -2<=q[1]<=2)) out["predictions"]["clipping_comparison"]={"predicted":"predictive projection has no certified one-step violations","clipping_violations":clip_bad,"projection_violations":safe_bad} # Same tiny control task: clipping only limits action; safety projection enforces predicted constraints. def target_policy(s): return 2.2*(0.65-s[0])-0.8*s[1] def episode(method, n=80): s=np.array([-0.8,0.0]); violations=0; reward=0.; projections=[] for t in range(n): z=target_policy(s)+np.random.default_rng(100+t).normal(0,.18) if method=='clip': p=float(np.clip(z,sl.u_min,sl.u_max)); d=0. else: p,d,b=sl.project(s,z) if p is None: p=0.; d=abs(z-p) # execute one step, then assess actual constraints x,v=s; ns=np.array([x+sl.dt*v+.5*sl.dt**2*p,v+sl.dt*p]) bad=not (-2<=ns[0]<=2 and -2<=ns[1]<=2 and -1<=p<=1) violations += int(bad) reward -= float((ns[0])**2+0.1*ns[1]**2+0.01*p*p) projections.append(d); s=ns return reward,violations,float(np.mean(projections)) reps=30 base=np.array([episode('clip') for _ in range(reps)]); idea=np.array([episode('safe') for _ in range(reps)]) out["simulation"]={"episodes":reps,"baseline_clip":{"mean_reward":float(base[:,0].mean()),"violations_per_episode":float(base[:,1].mean()),"mean_projection":float(base[:,2].mean())},"idea_projection":{"mean_reward":float(idea[:,0].mean()),"violations_per_episode":float(idea[:,1].mean()),"mean_projection":float(idea[:,2].mean())}} Path('results.json').write_text(json.dumps(out,indent=2)) print(json.dumps(out,indent=2)) assert ds.max()<1e-10 and abs(slope_l+1)<1e-10 and abs(slope_r-1)<1e-10 assert cert==feasible if __name__=='__main__': run()