PIPO-PITO bounded recurrent gain / pipo_pito_experiment.py
Failed on benchmark
1"""Numerical MVP for the PIPO-PITO bounded recurrent gain idea.
2
3The scalar positive plant is x' = (-d + w)x + p, y=x. A fixed recurrent
4input gain uses w=gamma. The proposed controller uses w'=b-a*y*w.
5This is intentionally a small mechanism test, not a neural benchmark.
6
7Run:
8 /home/maxwelhelp/main/bin/python3 pipo_pito_experiment.py
9"""
10import math
11import json
12import numpy as np
13
14
15def controller_constant_y(a, b, V, w0, dt=1e-4):
16 """Test w'=b-a*V*w against the PITO upper-bound prediction."""
17 eps = 2.0*b/(a*V)
18 predicted_T = max(0.0, math.log(2*w0/eps)/(a*V)) if w0 > eps else 0.0
19 horizon = max(2*predicted_T, 4/(a*V), 0.1)
20 n = int(math.ceil(horizon/dt))
21 w = float(w0)
22 crossing = 0.0 if w <= eps else None
23 trace = []
24 for i in range(1, n+1):
25 w += dt*(b-a*V*w)
26 if crossing is None and w <= eps:
27 crossing = i*dt
28 if i > 10 and w > eps/2 + 1e-8:
29 trace.append((i*dt, math.log(w-eps/2)))
30 slope = float('nan')
31 if len(trace) > 3:
32 xy = np.asarray(trace)
33 slope = -float(np.polyfit(xy[:, 0], xy[:, 1], 1)[0])
34 # Exact crossing is useful to distinguish the conservative theorem bound.
35 exact = (0.0 if w0 <= eps else
36 max(0., math.log((w0-eps/2)/(eps/2))/(a*V)))
37 return {"epsilon": eps, "V": V, "V_times_epsilon": V*eps,
38 "predicted_T_bound": predicted_T, "observed_T": crossing,
39 "exact_T": exact, "rate_predicted": a*V,
40 "rate_observed": slope, "rate_relative_error": abs(slope-a*V)/(a*V)}
41
42
43def closed_loop(gamma, regulated, seed=0, dt=0.002, horizon=80.0):
44 rng = np.random.default_rng(seed)
45 d, p, a, b = 0.4, 0.02, 1.0, 0.2
46 x = 0.1 + 0.01*rng.random()
47 w = gamma if not regulated else 2.0
48 max_x = x
49 for i in range(int(horizon/dt)):
50 if regulated:
51 w = max(0.0, w + dt*(b-a*x*w))
52 gain = w if regulated else gamma
53 x = max(0.0, x + dt*((-d+gain)*x+p))
54 max_x = max(max_x, x)
55 if x > 1e100:
56 break
57 # For p>0, gamma<d converges, gamma=d grows linearly, gamma>d grows exponentially.
58 unbounded = (not regulated and gamma >= d and x > 10.0)
59 return {"final_x": float(x), "max_x": float(max_x), "final_w": float(w),
60 "unbounded_by_80": bool(unbounded),
61 "predicted_regime": ("bounded" if gamma < d else "unbounded (marginal at equality)")}
62
63
64def main():
65 a, b, w0 = 1.3, 0.7, 3.0
66 epsilons = [0.2, 0.4, 0.8, 1.2]
67 pito = [controller_constant_y(a, b, 2*b/(a*eps), w0) for eps in epsilons]
68 # Fixed plant threshold: coefficient -d+gamma changes sign at gamma=d.
69 gammas = [0.2, 0.39, 0.4, 0.6, 0.8]
70 fixed = {str(g): closed_loop(g, False) for g in gammas}
71 regulated = closed_loop(0.8, True)
72 max_rate_error = max(r["rate_relative_error"] for r in pito)
73 max_bound_slack = max((r["observed_T"] or 0)-r["predicted_T_bound"] for r in pito)
74 result = {
75 "parameters": {"a": a, "b": b, "d": 0.4, "p": 0.02,
76 "w0": w0, "dt": 0.002, "horizon": 80.0},
77 "pito_sweep": pito,
78 "fixed_gain_sweep": fixed,
79 "regulated_at_nominal_gamma_0.8": regulated,
80 "quantitative_checks": {
81 "PITO_rate_prediction": {"predicted": "a*V", "max_relative_error": max_rate_error},
82 "PITO_threshold_prediction": {"predicted": "V*epsilon=2*b/a",
83 "observed_values": [r["V_times_epsilon"] for r in pito],
84 "target": 2*b/a},
85 "crossing_bound": {"predicted": "observed_T <= predicted_T_bound",
86 "maximum_observed_minus_predicted": max_bound_slack},
87 "plant_boundary": {"predicted": "gamma=d=0.4; gamma<d bounded, gamma>=d unbounded",
88 "observed": {g: fixed[str(g)]["unbounded_by_80"] for g in gammas}},
89 "regulated_boundedness": {"predicted": "feedback reduces gain and bounds x",
90 "observed_max_x": regulated["max_x"],
91 "observed_final_w": regulated["final_w"]}
92 }
93 }
94 print(json.dumps(result, indent=2))
95
96
97if __name__ == "__main__":
98 main()