Structured-μ Robust Optimizer / structured_robust_optimizer.py
Unverified
1"""MVP verification of a structured-mu feedback optimizer on a scalar quadratic.
2
3The local error dynamics are x[t+1] = A x[t] + w[t], z[t] = r x[t],
4where A=a-c is the nominal closed-loop pole, c is the optimizer feedback
5(g_t=a*x_t, update u=-c*x), and w=delta*z, |delta|<=1. Thus the uncertain
6pole is A+delta*r. The transfer from w to z is M(z)=r/(z-A), and the
7scalar structured singular value is |M|. Its frequency maximum is
8r/(1-|A|), so the certificate is mu<1 iff |A|+r<1.
9"""
10import json
11import math
12import numpy as np
13from dataclasses import dataclass
14
15@dataclass
16class StructuredMuController:
17 """Low-order (zero-state) controller with a conservative sampled mu guard."""
18 nominal_curvature: float
19 uncertainty_radius: float
20 target_mu: float = 0.90
21
22 def choose_gain(self):
23 # c=a makes A=0 and is the minimizer of the exact scalar certificate.
24 # If the uncertainty itself is too large, no static controller can certify it.
25 c = self.nominal_curvature
26 achieved = self.mu_estimate(c)
27 if achieved >= self.target_mu:
28 # Keep the nominal pole at zero and expose an honest failure.
29 return c, achieved, False
30 return c, achieved, True
31
32 def mu_estimate(self, c, nfreq=4097):
33 A = self.nominal_curvature - c
34 omega = np.linspace(0, 2*np.pi, nfreq)
35 M = self.uncertainty_radius / (np.exp(1j*omega) - A)
36 return float(np.max(np.abs(M)))
37
38def exact_mu(a, r, c):
39 A = a-c
40 return r/(1-abs(A)) if abs(A) < 1 else float("inf")
41
42def simulate(a, r, c, delta, steps=120, x0=1.0):
43 x = float(x0)
44 traj = [x]
45 for _ in range(steps):
46 x = (a-c+delta*r)*x
47 traj.append(x)
48 if not np.isfinite(x) or abs(x) > 1e12:
49 break
50 return np.asarray(traj)
51
52def numerical_core_check():
53 # Prediction 1: sampled frequency max equals r/(1-|A|), including its worst case.
54 rows=[]
55 for a,c,r in [(0.8,0.35,0.2),(0.8,0.8,0.65),(0.4,0.1,0.3)]:
56 A=a-c
57 w=np.linspace(0,2*np.pi,200001)
58 sampled=float(np.max(np.abs(r/(np.exp(1j*w)-A))))
59 analytic=exact_mu(a,r,c)
60 rows.append({"a":a,"c":c,"r":r,"sampled_mu":sampled,"analytic_mu":analytic,
61 "relative_error":abs(sampled-analytic)/analytic})
62 return rows
63
64def boundary_sweep():
65 # Prediction 2: for fixed A, transition is r*=1-|A|; measure actual trajectories.
66 a,c=0.8,0.35
67 A=a-c
68 predicted=1-abs(A)
69 rs=np.linspace(predicted-0.18,predicted+0.18,19)
70 records=[]
71 for r in rs:
72 # delta=+1 is the worst structured perturbation for positive A.
73 x=simulate(a,r,c,delta=1.0,steps=100)
74 stable=bool(np.max(np.abs(x)) < 1e6 and abs(x[-1]) < abs(x[0]))
75 records.append((float(r),stable,float(abs(x[-1]))))
76 stable_rs=[r for r,s,_ in records if s]
77 unstable_rs=[r for r,s,_ in records if not s]
78 observed=(max(stable_rs)+min(unstable_rs))/2 if stable_rs and unstable_rs else float("nan")
79 return {"A":A,"predicted_r_boundary":predicted,"observed_r_boundary":observed,
80 "records":records,"boundary_error":abs(observed-predicted)}
81
82def gain_sweep():
83 # Prediction 3: c=a (zero nominal pole) minimizes worst-case contraction and mu.
84 a,r=0.8,0.45
85 gains=np.linspace(0.0,1.6,161)
86 mus=np.array([exact_mu(a,r,c) for c in gains])
87 contractions=np.array([abs(a-c)+r for c in gains])
88 best_mu_gain=float(gains[np.argmin(mus)])
89 best_contraction_gain=float(gains[np.argmin(contractions)])
90 robust=StructuredMuController(a,r)
91 c,mu,cert=robust.choose_gain()
92 # baseline is a common under-corrected fixed step c=0.35
93 base_c=0.35
94 base_worst=abs(a-base_c)+r
95 idea_worst=abs(a-c)+r
96 return {"a":a,"r":r,"predicted_optimal_gain":a,
97 "observed_mu_minimizer":best_mu_gain,
98 "observed_contraction_minimizer":best_contraction_gain,
99 "controller_gain":c,"controller_mu":mu,"certificate":cert,
100 "baseline_gain":base_c,"baseline_worst_pole":base_worst,
101 "idea_worst_pole":idea_worst,
102 "mu_at_zero_uncertainty":exact_mu(a,0.0,a)}
103
104def divergence_probability_sweep(seed=11):
105 # Prediction: random bounded structured perturbations become unstable when
106 # the worst-case certificate mu=r/(1-|A|) exceeds one.
107 rng=np.random.default_rng(seed)
108 a,c=0.8,0.35
109 A=a-c
110 rs=np.linspace(0.40,0.70,16)
111 rows=[]
112 for r in rs:
113 divergent=0
114 trials=1000
115 for _ in range(trials):
116 delta=float(rng.uniform(-1,1))
117 pole=abs(A+delta*r)
118 # Long-horizon growth from x0=1, with a conservative numerical cutoff.
119 x=abs(pole)**80
120 divergent += int(x>10)
121 rows.append({"r":float(r),"mu":float(exact_mu(a,r,c)),
122 "divergence_probability":divergent/trials})
123 return {"predicted_mu_boundary":1.0,"rows":rows}
124
125def quadratic_mini_experiment(seed=7):
126 # Same tiny setup, adversarial bounded curvature/noise represented by delta.
127 rng=np.random.default_rng(seed)
128 a=0.8; r=0.45; horizon=80; trials=200
129 baseline_c=0.35
130 idea=StructuredMuController(a,r); idea_c,_,_=idea.choose_gain()
131 outcomes={"baseline":[],"idea":[]}
132 for _ in range(trials):
133 # Fixed but unknown structured perturbation plus small gradient noise.
134 delta=float(rng.uniform(-1,1))
135 noise=float(rng.normal(0,0.01))
136 for name,c in [("baseline",baseline_c),("idea",idea_c)]:
137 x=1.0; maxabs=1.0
138 for _ in range(horizon):
139 x=(a-c+delta*r)*x + noise
140 maxabs=max(maxabs,abs(x))
141 outcomes[name].append((abs(x),maxabs))
142 summary={}
143 for k,v in outcomes.items():
144 v=np.asarray(v)
145 summary[k]={"median_final_abs":float(np.median(v[:,0])),
146 "median_max_abs":float(np.median(v[:,1])),
147 "fraction_max_over_10":float(np.mean(v[:,1]>10))}
148 summary["controller_gain"]=idea_c
149 return summary
150
151def main():
152 out={"core_math":numerical_core_check(),"boundary":boundary_sweep(),
153 "gain_prediction":gain_sweep(),"divergence_sweep":divergence_probability_sweep(),"mini_experiment":quadratic_mini_experiment()}
154 print(json.dumps(out,indent=2))
155
156if __name__ == "__main__": main()