1import json, math, os, random
 2import numpy as np
 3import torch
 4
 5SEED=7
 6np.random.seed(SEED); random.seed(SEED); torch.manual_seed(SEED)
 7torch.set_num_threads(4)
 8DT=0.12
 9P=np.diag([1.0, 0.35]).astype(np.float32)
10C2=0.025
11
12# True plant has an unmodeled bounded force; nominal model omits it.
13def true_step(x,u, amp=0.22):
14    x=np.asarray(x); u=np.asarray(u)
15    return np.stack([x[...,0]+DT*x[...,1], x[...,1]+DT*(u+0.18*np.sin(1.7*x[...,0])+amp*np.sin(2.3*x[...,0]))],-1)
16def nominal_step(x,u):
17    x=np.asarray(x); u=np.asarray(u)
18    return np.stack([x[...,0]+DT*x[...,1], x[...,1]+DT*(u+0.18*np.sin(1.7*x[...,0]))],-1)
19def V(x): return 0.5*np.einsum('...i,ij,...j->...',x,P,x)
20def gradV(x): return np.asarray(x)@P.T
21
22def analytic_lipschitz(z,delta):
23    # Valid on the uncertainty ball for quadratic V: ||Pz|| + ||P||_2 delta.
24    return np.linalg.norm(gradV(z),axis=-1)+np.max(np.linalg.eigvalsh(P))*delta
25
26def math_check():
27    rng=np.random.default_rng(3); rows=[]
28    for d in [0.,.03,.08,.15,.25]:
29        z=rng.normal(size=(4000,2)); e=rng.normal(size=(4000,2)); e=e/np.linalg.norm(e,axis=1,keepdims=True)*d*rng.random((4000,1))
30        lhs=np.abs(V(z+e)-V(z)); rhs=analytic_lipschitz(z,d)*np.linalg.norm(e,axis=1)
31        rows.append((d,float(np.max(lhs-rhs)),float(np.mean(lhs/rhs.clip(1e-9)))))
32    return rows
33
34def policy_torch(K,x):
35    return -(x@K.T).clamp(-3,3)
36def train(robust, delta, n=1000, steps=180):
37    rng=np.random.default_rng(SEED+n+int(delta*1000)+(1 if robust else 0))
38    x=rng.uniform(-2,2,size=(n,2)).astype(np.float32)
39    # held-out-like nominal transition samples, using the nominal model as learned dynamics
40    xt=torch.tensor(x); K=torch.tensor([[1.0,1.0]],requires_grad=True)
41    opt=torch.optim.Adam([K],lr=.045)
42    for _ in range(steps):
43        u=policy_torch(K,xt)
44        # differentiable nominal transition for policy gradients
45        z= torch.stack([xt[:,0]+DT*xt[:,1], xt[:,1]+DT*(u[:,0]+.18*torch.sin(1.7*xt[:,0]))],1)
46        dv=0.5*torch.sum((z@torch.tensor(P))*z,1)-0.5*torch.sum((xt@torch.tensor(P))*xt,1)
47        lv=torch.sqrt(torch.sum((z@torch.tensor(P))**2,1)+1e-8)+delta
48        target=C2*torch.sum(xt*xt,1)
49        margin=dv+(lv*delta if robust else 0)+target
50        loss=torch.mean(u[:,0]**2)*.012 + torch.mean(torch.nn.functional.softplus(8*margin))/8 + .0005*torch.sum(K*K)
51        opt.zero_grad(); loss.backward(); opt.step()
52    return K.detach().numpy()
53
54def evaluate(K, delta_eval, amp=.22):
55    rng=np.random.default_rng(900+int(delta_eval*1000)); x=rng.uniform(-2,2,size=(1800,2)).astype(np.float32)
56    u=-(x@K.T).clip(-3,3); z=nominal_step(x,u[:,0]); dv=V(z)-V(x); lv=analytic_lipschitz(z,delta_eval)
57    robust_margin=dv+lv*delta_eval+C2*np.sum(x*x,1)
58    nominal_margin=dv+C2*np.sum(x*x,1)
59    # random bounded errors test the actual next state and certificate
60    e=true_step(x,u[:,0],amp=amp)-z
61    actual_margin=V(z+e)-V(x)+C2*np.sum(x*x,1)
62    conv=0
63    for j in range(400):
64        xx=rng.uniform(-2,2,size=2)
65        for t in range(60):
66            uu=float(np.clip(-(K@xx)[0],-3,3)); xx=true_step(xx,uu,amp=amp)
67            if np.linalg.norm(xx)<.12: conv+=1; break
68    return dict(min_robust=float(np.min(-robust_margin)), median_robust=float(np.median(-robust_margin)), min_nominal=float(np.min(-nominal_margin)), actual_violation=float(np.mean(actual_margin>0)), convergence=conv/1500, K=K.tolist())
69
70def mechanism_sweep():
71    # Prediction 1: for a fixed policy, m(delta)=m(0)-L*delta-(lambda_max P)delta^2;
72    # observed fit should have positive linear coefficient L and small quadratic term.
73    K=np.array([[1.25,1.55]],dtype=np.float32); rng=np.random.default_rng(11); x=rng.uniform(-1.8,1.8,(12000,2)); u=-(x@K.T).clip(-3,3)[:,0]; z=nominal_step(x,u)
74    base=V(z)-V(x)+C2*np.sum(x*x,1); L=np.linalg.norm(gradV(z),axis=1); lam=np.max(np.linalg.eigvalsh(P))
75    ds=np.array([0,.02,.05,.10,.18,.28]); observed=[]; predicted=[]
76    for d in ds:
77        observed.append(float(np.min(-(base+(L+d)*d))))
78        predicted.append(float(np.min(-base-L*d-lam*d*d)))
79    # Prediction 2: actual bounded error is certified whenever predicted robust margin >= 0.
80    # Sweep scalar c2 to locate the sign transition and compare actual violations.
81    c2s=np.array([0,.01,.02,.03,.05,.08]); boundary=[]; actual=[]
82    e=true_step(x,u,amp=.22)-z; ae=V(z+e)-V(x)
83    for c in c2s:
84        mm=ae+c*np.sum(x*x,1); boundary.append(float(np.min(-mm))); actual.append(float(np.mean(mm>0)))
85    return {'delta_sweep':{'delta':ds.tolist(),'observed_min_certificate':observed,'predicted_min_certificate':predicted},'c2_sweep':{'c2':c2s.tolist(),'observed_min_actual_certificate':boundary,'actual_violation_fraction':actual}}
86
87def main():
88    check=math_check(); mech=mechanism_sweep()
89    # Same data-size proxy: uncertainty radius grows with scarcity, calibrated conservatively.
90    rounds=[]
91    for frac,delta in [(1.0,.055),(.5,.11),(.2,.18)]:
92        kb=train(False,delta,n=int(500*frac)); kr=train(True,delta,n=int(500*frac))
93        rounds.append({'fraction':frac,'delta':delta,'baseline':evaluate(kb,delta),'robust':evaluate(kr,delta)})
94    out={'math_check':check,'mechanism_predictions':mech,'training_comparison':rounds}
95    with open('results.json','w') as f: json.dump(out,f,indent=2)
96    print(json.dumps(out,indent=2))
97if __name__=='__main__': main()