Prescribed-Performance Hidden-State Observer / ppo_observer_experiment.py
Mechanism failed
1"""Prescribed-performance hidden-state observer MVP.
2
3Run: /home/maxwelhelp/main/bin/python3 ppo_observer_experiment.py
4"""
5import json, math, random
6from pathlib import Path
7import numpy as np
8
9SEED = 7
10np.random.seed(SEED); random.seed(SEED)
11
12
13def envelope(k, r0=1.0, rinf=.08, lo=.035, dt=.02):
14 return (r0-rinf)*np.exp(-lo*k*dt)+rinf
15
16
17def ppo_update(hat, e, r, dt=.02, gamma=(2.0,), rinf=.08, eps=1e-6):
18 """Euler PPO update, with safe transformed innovation clipping."""
19 n = len(hat)
20 xi = float(np.clip((e-hat[0])/r, -1+eps, 1-eps))
21 T = math.log((1+xi)/(1-xi))
22 old = hat.copy()
23 for i in range(n-1):
24 hat[i] += dt*(old[i+1] + gamma[i]*r*T/(rinf**(i+1)))
25 hat[-1] += dt*gamma[-1]*r*T/(rinf**n)
26 return hat, xi, T
27
28
29def fixed_update(hat, e, gain=1.0, dt=.02):
30 # A simple conventional fixed-gain first-order residual observer.
31 hat[0] += dt*gain*(e-hat[0])
32 return hat
33
34
35def math_check():
36 # Constant error: PPO should rapidly reduce normalized innovation while
37 # retaining a shrinking prescribed envelope; compare a deliberately modest
38 # fixed gain, which cannot track the shrinking envelope as well.
39 steps=1800; e=.18
40 hp=np.zeros(1); hf=np.zeros(1); pxis=[]; fxis=[]; violations=0; fviol=0
41 for k in range(steps):
42 r=envelope(k)
43 hp, xi, _ = ppo_update(hp,e,r,gamma=(1.4,),rinf=.08)
44 hf=fixed_update(hf,e,gain=1.0)
45 pxis.append(abs((e-hp[0])/r)); fxis.append(abs((e-hf[0])/r))
46 violations += pxis[-1] >= 1
47 fviol += fxis[-1] >= 1
48 return {
49 'ppo_max_normalized_innovation': float(max(pxis)),
50 'fixed_max_normalized_innovation': float(max(fxis)),
51 'ppo_final_normalized_innovation': float(pxis[-1]),
52 'fixed_final_normalized_innovation': float(fxis[-1]),
53 'ppo_violations': int(violations), 'fixed_violations': int(fviol),
54 'ppo_improves_over_fixed': float(pxis[-1]) < float(fxis[-1]),
55 }
56
57
58def gain_sweep():
59 """Probe the Euler gain boundary claimed by the construction."""
60 out = {}
61 for g in (0.2, 1.4, 5.0, 20.0, 100.0, 500.0):
62 h = np.zeros(1); max_xi = 0.0
63 for k in range(600):
64 h, xi, _ = ppo_update(h, .18, envelope(k), dt=.02,
65 gamma=(g,), rinf=.08)
66 max_xi = max(max_xi, abs(xi))
67 out[str(g)] = {'max_abs_xi': float(max_xi),
68 'final_hat_e1': float(h[0])}
69 return out
70
71
72def make_data(n=2600):
73 # Observable noisy sinusoid with slowly varying frequency/amplitude.
74 t=np.arange(n+1,dtype=np.float64)
75 clean=np.sin(.075*t)+.25*np.sin(.013*t+1.1)
76 noisy=clean + .10*np.random.randn(n+1)
77 return noisy.astype(np.float32), clean.astype(np.float32)
78
79
80def train_gru(train, epochs=5, hidden=20):
81 import torch
82 from torch import nn
83 torch.manual_seed(SEED)
84 device='cuda' if torch.cuda.is_available() else 'cpu'
85 try:
86 model=nn.GRU(1,hidden,batch_first=True).to(device)
87 head=nn.Linear(hidden,1).to(device)
88 opt=torch.optim.Adam(list(model.parameters())+list(head.parameters()),lr=3e-3)
89 # Short windows make the comparison cheap and reproducible.
90 x=[]; y=[]; W=24
91 for i in range(0,len(train)-W-1,2):
92 x.append(train[i:i+W]); y.append(train[i+W])
93 x=torch.tensor(np.asarray(x),device=device).unsqueeze(-1)
94 y=torch.tensor(np.asarray(y),device=device).view(-1,1)
95 for _ in range(epochs):
96 perm=torch.randperm(len(x),device=device)
97 for j in range(0,len(x),128):
98 ix=perm[j:j+128]; out,_=model(x[ix]); pred=head(out[:,-1])
99 loss=((pred-y[ix])**2).mean(); opt.zero_grad(); loss.backward(); opt.step()
100 return model,head,device
101 except Exception as ex:
102 # CUDA OOM or another backend issue: retry on CPU as required.
103 torch.cuda.empty_cache() if torch.cuda.is_available() else None
104 device='cpu'; model=nn.GRU(1,hidden,batch_first=True).to(device); head=nn.Linear(hidden,1).to(device)
105 opt=torch.optim.Adam(list(model.parameters())+list(head.parameters()),lr=3e-3)
106 W=24; x=[]; y=[]
107 for i in range(0,len(train)-W-1,2): x.append(train[i:i+W]); y.append(train[i+W])
108 x=torch.tensor(np.asarray(x)).unsqueeze(-1); y=torch.tensor(np.asarray(y)).view(-1,1)
109 for _ in range(epochs):
110 for j in range(0,len(x),128):
111 out,_=model(x[j:j+128]); pred=head(out[:,-1]); loss=((pred-y[j:j+128])**2).mean(); opt.zero_grad(); loss.backward(); opt.step()
112 return model,head,device
113
114
115def predict_compare(data, clean, model, head, device):
116 import torch
117 W=24; hp=np.zeros(1); hf=np.zeros(1); raw=[]; ppo=[]; fixed=[]; truth=[]; xis=[]
118 with torch.no_grad():
119 for k in range(W,len(data)-1):
120 seq=torch.tensor(data[k-W:k],device=device).view(1,W,1)
121 pred=float(head(model(seq)[0][:,-1]).cpu().item())
122 # e is the currently observable residual; observer estimates its
123 # slowly varying component and corrects the next output.
124 e=float(data[k]-pred); r=envelope(k-W,r0=.55,rinf=.08,lo=.08)
125 hp,xi,_=ppo_update(hp,e,r,dt=.05,gamma=(1.0,),rinf=.08)
126 hf=fixed_update(hf,e,gain=1.2,dt=.05)
127 raw.append(pred); ppo.append(pred+hp[0]); fixed.append(pred+hf[0]); truth.append(float(data[k+1])); xis.append(abs(xi))
128 def rmse(a,b): return float(np.sqrt(np.mean((np.asarray(a)-np.asarray(b))**2)))
129 return {'raw_gru_rmse':rmse(raw,truth),'ppo_corrected_rmse':rmse(ppo,truth),'fixed_observer_rmse':rmse(fixed,truth),
130 'ppo_mean_abs_xi':float(np.mean(xis)),'ppo_max_abs_xi':float(np.max(xis)),
131 'ppo_fraction_envelope_violations':float(np.mean(np.asarray(xis)>=1.0))}
132
133
134def main():
135 check=math_check()
136 train,clean=make_data(); split=1900
137 model,head,device=train_gru(train[:split],epochs=5)
138 pred=predict_compare(train[split-24:],clean[split-24:],model,head,device)
139 result={'seed':SEED,'device':device,'math_check':check,'gain_sweep':gain_sweep(),'prediction':pred}
140 Path('results.json').write_text(json.dumps(result,indent=2))
141 print(json.dumps(result,indent=2))
142
143if __name__=='__main__': main()