Ω + SE44 implementation mapping (discrete-time control demo)
import numpy as np
import matplotlib.pyplot as plt
# -----------------------------
# Ω + SE44 implementation mapping (discrete-time control demo)
# -----------------------------
rng = np.random.default_rng(7)
# Plant: x_{k+1} = A x_k + B u_k + w_k
A = np.array([[1.02, 0.05],
[0.00, 0.98]])
B = np.array([[0.10],
[0.08]])
# LQR-ish stabilizing gain (picked for demo; can be derived formally)
K = np.array([[2.6, 1.4]]) # u = -K x
# Ω structure:
# Ω_k = (state_k + bias_k) * alpha_k
# state_k := y_k (a measurable scalar derived from x_k)
# bias_k := disturbance estimate (simple EWMA of residual)
# alpha_k := adaptive gain (bounded), updated from coherence/entropy
#
# u_k = sat( -K x_k + Ω_k )
#
# SE44 gate computed from measurable signals:
# C_k (coherence) = exp(-||e_k|| / c_scale) where e_k is tracking error to 0
# S_k (entropy) = min(1, ||Δx_k|| / s_scale)
# RMS_drift_k = RMS of recent Δx norms over a window
# Gate:
# ACCEPT if C>=Cmin and S<=Smax and RMS<=RMSmax
# else REBIND: u_k := u_{k-1} (safe hold)
Cmin = 0.985
Smax = 0.01
RMSmax = 0.001
c_scale = 0.06 # sets how strict coherence is vs state magnitude
s_scale = 0.50 # sets entropy normalization
alpha_min, alpha_max = 0.0, 0.8 # bounded amplification
alpha = 0.4
# Bias estimator (EWMA of residuals on y)
beta = 0.05
bias = 0.0
# Saturation on control
u_max = 2.0
# Drift window
W = 25
dx_hist = []
def sat(u, umax):
return np.clip(u, -umax, umax)
def coherence(x):
# "coherence" as semantic-shape alignment → here: small state magnitude
return float(np.exp(-np.linalg.norm(x) / c_scale))
def entropy(dx):
# "entropy" as permissible novelty → here: normalized step change magnitude
return float(min(1.0, np.linalg.norm(dx) / s_scale))
def rms_drift(dx_hist):
if len(dx_hist) == 0:
return 0.0
arr = np.array(dx_hist)
return float(np.sqrt(np.mean(arr**2)))
# Simulation
T = 800
x = np.array([[0.15],
[-0.10]])
u_prev = 0.0
# logs
Xs = np.zeros((T, 2))
Us = np.zeros(T)
Omegas = np.zeros(T)
Cs = np.zeros(T)
Ss = np.zeros(T)
RMSs = np.zeros(T)
accept = np.zeros(T, dtype=int)
for k in range(T):
# measurable scalar state proxy y (engineers: choose a measurement model; here y = x1)
y = float(x[0, 0])
# residual / innovation for bias estimator (pretend desired y* = 0)
resid = y
bias = (1 - beta) * bias + beta * resid
# omega signal
omega = (y + bias) * alpha
# candidate control (baseline stabilizing feedback + omega)
u_cand = float(-K @ x + omega)
u_cand = float(sat(u_cand, u_max))
# Predict dx for gating (using candidate u_cand)
w = rng.normal(0.0, 0.0002, size=(2, 1)) # small bounded disturbance
x_next_pred = A @ x + B * u_cand + w
dx = x_next_pred - x
# Update drift history (use predicted dx norm for gate)
dx_norm = float(np.linalg.norm(dx))
dx_hist.append(dx_norm)
if len(dx_hist) > W:
dx_hist.pop(0)
Ck = coherence(x)
Sk = entropy(dx)
RMSk = rms_drift(dx_hist)
# SE44 gate
if (Ck >= Cmin) and (Sk <= Smax) and (RMSk <= RMSmax):
uk = u_cand
accept[k] = 1
else:
uk = u_prev # REBIND: hold last safe control
accept[k] = 0
# Apply to real plant with disturbance
w = rng.normal(0.0, 0.0002, size=(2, 1))
x = A @ x + B * uk + w
u_prev = uk
# adaptive alpha update (bounded): increase only when gate passes; decay otherwise
if accept[k] == 1:
alpha = min(alpha_max, alpha + 0.001)
else:
alpha = max(alpha_min, alpha - 0.004)
# log
Xs[k] = x.ravel()
Us[k] = uk
Omegas[k] = omega
Cs[k] = Ck
Ss[k] = Sk
RMSs[k] = RMSk
# Plot results (no explicit colors)
t = np.arange(T)
plt.figure()
plt.plot(t, Xs[:, 0], label="x1")
plt.plot(t, Xs[:, 1], label="x2")
plt.legend()
plt.title("State trajectory under Ω + SE44-gated control")
plt.xlabel("k")
plt.ylabel("state")
plt.figure()
plt.plot(t, Us, label="u")
plt.legend()
plt.title("Control input (with saturation)")
plt.xlabel("k")
plt.ylabel("u")
plt.figure()
plt.plot(t, Omegas, label="Ω contribution")
plt.legend()
plt.title("Ω signal injected into control law")
plt.xlabel("k")
plt.ylabel("Ω")
plt.figure()
plt.plot(t, Cs, label="C (coherence)")
plt.plot(t, Ss, label="S (entropy)")
plt.plot(t, RMSs, label="RMS drift")
plt.legend()
plt.title("SE44 gate signals (measurable)")
plt.xlabel("k")
plt.ylabel("value")
plt.figure()
plt.plot(t, accept.astype(float), label="SE44 ACCEPT=1 / REBIND=0")
plt.ylim(-0.1, 1.1)
plt.legend()
plt.title("Gate decisions over time")
plt.xlabel("k")
plt.ylabel("accept")
plt.show()
# Return a few numeric summaries for engineering inspection
summary = {
"final_state": Xs[-1].tolist(),
"state_norm_max": float(np.max(np.linalg.norm(Xs, axis=1))),
"control_abs_max": float(np.max(np.abs(Us))),
"accept_rate": float(np.mean(accept)),
"alpha_bounds": [alpha_min, alpha_max],
}
summary
STDOUT/STDERR
Matplotlib is building the font cache; this may take a moment.
Result
{'final_state': [171657.31534174748, 0.0006874320312023781],
'state_norm_max': 171657.31534174748,
'control_abs_max': 0.0,
'accept_rate': 0.0,
'alpha_bounds': [0.0, 0.8]}
Comments
Post a Comment