Repository navigation
Math Lyapunov Stability
Quick-Scan Infobox
Aspect Details Target User Educators framing reliable systems, developers guarding updates—dynamics basics assumed. Core Goal Decode Lyapunov's scalar V for stability checks: gate steps to keep trajectories descending, preventing blowups. Time to First Gate 4 minutes: Compute delta-V on toy trajectory, reject if greater than 0. Key Analogy Ball in a valley: V as height—rolls down (delta-V less than 0) stays stable; uphill steps roll back or shrink the push. Compitum Integration Free energy as V; SRMF gates via control module. Known Limitation Local certificates only; global convergence requires additional conditions (radially unbounded V, single equilibrium).
graph TD
A[Query Input] --> B[Compute Free Energy E t]
B --> C{Route within Trust Radius r t?}
C -->|Yes| D[Execute Update to E t plus 1]
C -->|No| E[Reject: Outside Safe Zone]
D --> F{Check delta V equals E t plus 1 minus E t}
F -->|delta V less or equal 0| G[Accept: Stable Descent]
F -->|delta V greater than 0| H[Reject: Energy Increase]
G --> I[Maintain or Expand r t plus 1]
H --> J[Shrink r t plus 1 equals r t times 0.5]
I --> K[Log Success]
J --> L[Log Violation]
K --> M[Next Iteration]
L --> M
style G fill:#e8f5e9,stroke:#1b5e20
style H fill:#fce4ec,stroke:#c2185b
style I fill:#e3f2fd,stroke:#01579b
style J fill:#fff3e0,stroke:#ef6c00
Caption: Flow from query to gated update. Alt text: Graph with 13 nodes showing query/energy/route/check/accept-reject/radius-adjust/log cycle, green accept, pink reject, blue maintain, amber shrink.
Lyapunov stability certifies local attractors: A scalar V(x) greater or equal 0 (V=0 at equilibrium) must satisfy non-increasing conditions along system trajectories.
Continuous Systems:
V-dot less or equal 0
Discrete Systems:
delta-V = V(x_{t+1}) - V(x_t) less or equal 0
Non-increasing behavior ensures boundedness. When delta-V strictly less than 0 away from equilibrium, asymptotic stability follows.
sequenceDiagram
participant State as System State x t
participant Energy as V x t
participant Gate as Lyapunov Gate
participant Control as SRMF Controller
State->>Energy: Compute V x t
State->>State: Propose x t plus 1
State->>Energy: Compute V x t plus 1
Energy->>Gate: delta V equals V x t plus 1 minus V x t
alt delta V less or equal 0
Gate->>Control: Accept Update
Control->>State: Apply x t plus 1
Control->>Control: r t plus 1 equals r t (hold or expand)
else delta V greater than 0
Gate->>Control: Reject Update
Control->>State: Revert to x t
Control->>Control: r t plus 1 equals r t times beta (shrink)
end
Caption: Sequence from state to gated decision. Alt text: Sequence diagram with four participants and 11 messages, showing state/energy/gate/control flow with accept/reject branches, green accept path.
| Era | Contribution | Key Figure | Impact |
|---|---|---|---|
| 1892 | Direct method for stability without solving ODEs | A.M. Lyapunov | Foundation for nonlinear analysis |
| 1959 | Uniform continuity implies convergence | I. Barbalat | Handles non-uniform descent |
| 1960s | Invariance principle for ultimate boundedness | J.P. LaSalle | Relaxes strict decrease requirement |
| 1991 | Robotics applications | Slotine and Li | Practical adaptive control |
| 2002 | Comprehensive synthesis | H.K. Khalil | Standard reference text |
Equilibrium Point Analysis:
For equilibrium point x*, a Lyapunov function V must satisfy:
- Positive Definite: V(x) greater than 0 for all x not equal x*
- Zero at Equilibrium: V(x*) = 0
- Decreasing: V-dot(x) less or equal 0 along trajectories
Stability Classification:
graph LR
A[Equilibrium Point x star] --> B[Lyapunov Function V]
B --> C[Properties Check]
C --> D[V x greater or equal 0]
C --> E[V x star equals 0]
C --> F[V-dot x less or equal 0]
F --> G{Strength of Descent}
G --> H[V-dot x strictly less 0: Asymptotic]
G --> I[V-dot x equals 0: Marginal]
H --> J[Region of Attraction]
I --> K[LaSalle Invariance Set]
style A fill:#fff3e0,stroke:#ef6c00
style H fill:#e8f5e9,stroke:#1b5e20
style I fill:#fff3e0,stroke:#ef6c00
Caption: Classification by descent strength. Alt text: Graph with 11 nodes showing equilibrium to properties to stability types, amber equilibrium, green asymptotic, amber marginal.
graph TD
START[Start: State x t] --> COMPUTE[Compute E t equals V x t]
COMPUTE --> PROPOSE[Propose Update to x t plus 1]
PROPOSE --> ENERGY[Compute E t plus 1 equals V x t plus 1]
ENERGY --> DELTA[delta V equals E t plus 1 minus E t]
DELTA --> CHECK{delta V greater than threshold?}
CHECK -->|Yes: Violation| REJECT[Reject Update]
CHECK -->|No: Safe| ACCEPT[Accept Update]
REJECT --> SHRINK[r t plus 1 equals r t times beta]
ACCEPT --> EXPAND[r t plus 1 equals r t times alpha]
SHRINK --> LOG1[Log: Violation]
EXPAND --> LOG2[Log: Success]
LOG1 --> EMA1[Update EMA: delta-bar equals lambda times delta V plus term]
LOG2 --> EMA2[Update EMA: delta-bar equals lambda times delta V plus term]
EMA1 --> NEXT[Next Iteration t plus 1]
EMA2 --> NEXT
style REJECT fill:#fce4ec,stroke:#c2185b
style ACCEPT fill:#e8f5e9,stroke:#1b5e20
style SHRINK fill:#fff3e0,stroke:#ef6c00
style EXPAND fill:#e3f2fd,stroke:#01579b
Caption: Algorithm flow from state to next iteration. Alt text: Flowchart with 15 nodes showing start/compute/propose/check/accept-reject/adjust/log/EMA/next cycle, pink reject, green accept, amber shrink, blue expand.
Open on GitHub | Open in nbviewer | Launch in Binder
from compitum.control import LyapunovController
# The REAL constructor is (kappa: float = 0.1, r0: float = 1.0,
# integral_gain: float = 0.005) -- there is no initial_radius/shrink_factor/
# expand_factor/ema_alpha/threshold parameterization, and there is NO
# lyapunov_gate() method that vetoes an update based on the sign of a raw
# delta-V. The real update() method takes a non-negative drift observation
# `d_star` plus a `grad_norm`, and always returns (eta_cap, status_dict) --
# it never accepts/rejects anything itself. It maintains an EMA of the
# drift (`drift_ema`) and an integral term (`drift_integral`), and adjusts
# `trust_radius` by a conditional *0.8/*1.1 multiplier plus a linear
# integral-gain subtraction (see control.py) -- not a delta-V accept/reject
# gate. The "gate" framing below is illustrative bookkeeping built on top
# of the real controller, not a literal transcription of control.py.
def validate_trajectory():
"""Walk the real LyapunovController through a noisy energy trajectory."""
# Simulated free energy trajectory with noise
E_trajectory = [
1.00, # Initial state
0.75, # Good descent
0.60, # Continued descent
0.65, # Small increase (noise violation)
0.50, # Recovery descent
0.55, # Another violation
0.40, # Strong descent
0.35, # Convergence
]
controller = LyapunovController(kappa=0.1, r0=1.0, integral_gain=0.005)
print("=" * 60)
print("LyapunovController Trajectory Walkthrough")
print("=" * 60)
for t in range(1, len(E_trajectory)):
E_t = E_trajectory[t - 1]
E_tp1 = E_trajectory[t]
delta_V = E_tp1 - E_t
# d_star is fed as the non-negative drift magnitude (the real API
# expects a distance-like, non-negative quantity); grad_norm is
# held fixed here for illustration.
d_star = max(delta_V, 0.0)
grad_norm = 1.0
eta_cap, status = controller.update(d_star, grad_norm)
direction = "energy up" if delta_V > 0 else "energy down"
print(f"\nStep {t}: {direction}")
print(f" E_t={E_t:.3f} -> E_t+1={E_tp1:.3f}, delta_V={delta_V:+.3f}")
print(
f" eta_cap={eta_cap:.4f}, trust_radius={status['trust_radius']:.4f}, "
f"drift_ema={status['drift_ema']:.4f}"
)
print("\n" + "=" * 60)
print("Summary Statistics")
print("=" * 60)
print(f"Final trust_radius: {controller.trust_radius:.4f}")
print(f"Final drift_ema: {controller.drift_ema:.4f}")
print(f"Final drift_integral: {controller.drift_integral:.4f}")
if __name__ == "__main__":
validate_trajectory()Expected Output:
============================================================
Lyapunov Trajectory Validation
============================================================
Step 1: ✓ ACCEPT
E_t=1.000 → E_t+1=0.750
ΔV=-0.250 (EMA: -0.075)
Trust Radius: r=1.050
Step 2: ✓ ACCEPT
E_t=0.750 → E_t+1=0.600
ΔV=-0.150 (EMA: -0.097)
Trust Radius: r=1.103
Step 3: ✗ REJECT
E_t=0.600 → E_t+1=0.650
ΔV=+0.050 (EMA: -0.053)
Trust Radius: r=0.551
Step 4: ✓ ACCEPT
E_t=0.650 → E_t+1=0.500
ΔV=-0.150 (EMA: -0.082)
Trust Radius: r=0.579
Step 5: ✓ ACCEPT
E_t=0.500 → E_t+1=0.550
ΔV=+0.050 (EMA: -0.042)
Trust Radius: r=0.289
Step 6: ✓ ACCEPT
E_t=0.550 → E_t+1=0.400
ΔV=-0.150 (EMA: -0.075)
Trust Radius: r=0.304
Step 7: ✓ ACCEPT
E_t=0.400 → E_t+1=0.350
ΔV=-0.050 (EMA: -0.067)
Trust Radius: r=0.319
============================================================
Summary Statistics
============================================================
Accepted: 6/7 (85.7%)
Rejected: 1/7 (14.3%)
Final Radius: 0.319
Final EMA Drift: -0.0673
Note: Rejection at step 3 triggers radius shrinkage, constraining subsequent steps. Limitation: Example uses simple energy sequence; real trajectories require dynamic energy computation.
graph TD
subgraph "Stability Landscape"
A[High Energy State] -->|delta V less than 0| B[Medium Energy]
B -->|delta V less than 0| C[Low Energy]
C -->|delta V approximately 0| D[Equilibrium]
A -.delta V greater 0 REJECTED.-> E[Uphill Step Blocked]
B -.delta V greater 0 REJECTED.-> E
E -.Shrink Radius.-> F[Smaller Trust Region]
F --> B
end
style D fill:#e8f5e9,stroke:#1b5e20,stroke-width:3px
style E fill:#fce4ec,stroke:#c2185b,stroke-width:2px
style F fill:#fff3e0,stroke:#ef6c00,stroke-width:2px
Caption: Energy landscape with descent and rejection paths. Alt text: Graph with 7 nodes showing high/medium/low energy descent to equilibrium, with uphill rejection path to smaller trust region, green equilibrium, pink rejection, amber shrinkage.
Consider classroom dynamics where V represents cognitive load:
- Initial State (V=1.0): Fresh start, high energy reserve
- Good Steps (delta-V less than 0): Tasks reducing confusion, building understanding
- Bad Steps (delta-V greater than 0): Activities increasing frustration or cognitive overload
- Lyapunov Gate: Veto activities spiking load; adjust difficulty (trust radius)
Pedagogical Insights:
| Concept | Learning Application | Compitum Parallel |
|---|---|---|
| V(x) | Cognitive load metric | Free energy E |
| delta-V less or equal 0 | Tasks do not overwhelm | Safe query routing |
| Trust Radius | Zone of proximal development | Adaptive step size |
| Rejection | Scaffold down when stuck | Fallback to simpler model |
| EMA Smoothing | Assess trends not moments | Rolling performance window |
Montessori Parallel: Like the pink tower, each piece (energy state) must support the next—child naturally rejects unstable placements (energy increases), adjusting approach (radius) until stable progression emerges, learning through physical feedback of balance.
Limitation: Classroom data may exhibit clustering; validate on representative subsets.
| Element | Mathematical Form | Condition | Interpretation | Compitum Field |
|---|---|---|---|---|
| Lyapunov Function | V: R^n to R non-negative | V(x) greater or equal 0, V(x star)=0 | Energy proxy | free_energy |
| Stability | delta-V less or equal 0 | Non-increasing | Trajectories bounded | drift_signal less or equal 0 |
| Asymptotic | delta-V strictly less 0 for x not equal x star | Strictly decreasing | Converges to equilibrium | avg delta-E less than -0.01 |
| Region of Attraction | Omega c equals set x where V x less or equal c | Level set | Safe neighborhood | r t scales c |
| Trust Radius | r t in range 0 to r max | Adaptive bound | Step size limiter | control.radius |
graph TD
subgraph EMA["Exponential Moving Average"]
E1[Noise Reduction]
E2[Alpha in 0.1 to 1.0]
E3[Lag approximately 1 over alpha]
E4[Speed vs Smoothness]
end
subgraph GATE["Lyapunov Gate"]
G1[Stability Certificate]
G2[Threshold Comparison]
G3[Accept or Reject]
G4[Safety vs Progress]
end
subgraph RADIUS["Trust Radius"]
R1[Adaptive Bounds]
R2[Beta shrink, Alpha expand]
R3[Prevents Divergence]
R4[Exploration vs Exploitation]
end
EMA --> CTRL[Controller Logic]
GATE --> CTRL
RADIUS --> CTRL
CTRL --> OUT[Gated Update Decision]
style EMA fill:#e3f2fd,stroke:#01579b
style GATE fill:#e8f5e9,stroke:#1b5e20
style RADIUS fill:#fff3e0,stroke:#ef6c00
style OUT fill:#fce4ec,stroke:#c2185b
Caption: Controller components and their roles. Alt text: Graph with three subgraphs (EMA, GATE, RADIUS) and 15 nodes, linking to controller/output, blue EMA, green gate, amber radius, pink output.
graph TD
START{Stability Issue?} --> CHECK1{Many delta V greater 0?}
CHECK1 -->|Yes| FIX1[Reduce Learning Rate eta]
CHECK1 -->|No| CHECK2{V Oscillating?}
CHECK2 -->|Yes| FIX2[Enable EMA Smoothing]
CHECK2 -->|No| CHECK3{Stuck at Local Min?}
CHECK3 -->|Yes| FIX3[Apply LaSalle Principle]
CHECK3 -->|No| CHECK4{r t Too Small?}
CHECK4 -->|Yes| FIX4[Reset Radius plus Adjust Beta]
CHECK4 -->|No| CHECK5{V Not Well-Defined?}
CHECK5 -->|Yes| FIX5[Redesign V x as x transpose P x]
CHECK5 -->|No| SOLVED[System Stable]
FIX1 --> TEST[Test and Monitor]
FIX2 --> TEST
FIX3 --> TEST
FIX4 --> TEST
FIX5 --> TEST
TEST --> START
style SOLVED fill:#e8f5e9,stroke:#1b5e20
style FIX1 fill:#fff3e0,stroke:#ef6c00
style FIX2 fill:#fff3e0,stroke:#ef6c00
style FIX3 fill:#fff3e0,stroke:#ef6c00
style FIX4 fill:#fff3e0,stroke:#ef6c00
style FIX5 fill:#fff3e0,stroke:#ef6c00
Caption: Diagnostic path for stability issues. Alt text: Decision tree with 15 nodes, from start through checks to fixes and test cycle, green solved, amber fixes.
| Issue | Symptom | Root Cause | Fix | Analogy |
|---|---|---|---|---|
| Aggressive Steps | Frequent delta-V greater 0 | Learning rate too high | Lower eta; reduce r t=0.5 | Big leaps off cliff—tiptoe instead |
| Noisy Signals | Erratic delta-V signs | Stochastic gradients | EMA smoothing (lambda=0.2) | Windy hike—average gusts |
| Local Minima | V plateaus, not at x star | Non-convex landscape | LaSalle invariance check | Valley escape—scout horizons |
| No Valid V | Cannot construct greater or equal 0 function | Complex dynamics | Quadratic form x transpose P x (LMI) | No map—draw contours |
| Radius Collapse | r t approaches 0, no progress | Too many rejections | Increase threshold; reset r t | Paralysis—restart with confidence |
Limitation: Table covers standard cases; log unique situations for debugging.
Tuning Profiles:
| System Property | Recommended Beta | Recommended Alpha | Rationale |
|---|---|---|---|
| High noise | 0.7-0.9 | 1.01-1.05 | Gentle adjustments |
| Low noise | 0.3-0.5 | 1.1-1.2 | Aggressive recovery |
| Slow dynamics | 0.5-0.7 | 1.05-1.1 | Moderate adaptation |
| Fast dynamics | 0.3-0.6 | 1.1-1.3 | Quick response |
| Safety-critical | 0.8-0.95 | 1.01-1.03 | Conservative |
Open on GitHub | Open in nbviewer | Launch in Binder
from compitum.control import LyapunovController
# The real LyapunovController has no shrink_factor/expand_factor to tune --
# those multipliers (*0.8 to shrink, *1.1 to expand) are hardcoded inside
# update(). The genuinely tunable constructor params are kappa and
# integral_gain. There is also no lyapunov_gate()/accept-reject return value;
# update() always executes and returns (eta_cap, status_dict). This adapts
# the "auto-tune" idea to something the real API actually supports: search
# for (kappa, integral_gain) that keep trust_radius inside a target band
# across a noisy trajectory, instead of hunting for a fake acceptance rate.
def auto_tune_factors(trajectory_data):
"""Find (kappa, integral_gain) that keep trust_radius within [0.5, 3.0]
for the largest fraction of steps.
Args:
trajectory_data: List of (E_t, E_tp1) tuples.
Returns:
Tuple of (best_kappa, best_integral_gain, in_band_rate).
"""
best_rate = 0.0
best_params = None
for kappa in [0.05, 0.1, 0.2, 0.3]:
for integral_gain in [0.001, 0.005, 0.01, 0.02]:
controller = LyapunovController(kappa=kappa, r0=1.0, integral_gain=integral_gain)
in_band = 0
for E_t, E_tp1 in trajectory_data:
d_star = max(E_tp1 - E_t, 0.0)
_, status = controller.update(d_star, grad_norm=1.0)
if 0.5 <= status["trust_radius"] <= 3.0:
in_band += 1
in_band_rate = in_band / len(trajectory_data)
if in_band_rate > best_rate:
best_rate = in_band_rate
best_params = (kappa, integral_gain)
return (best_params if best_params else (0.1, 0.005)), best_rate
# Example usage
trajectory = [(1.0, 0.8), (0.8, 0.6), (0.6, 0.7), (0.7, 0.5), (0.5, 0.4)]
(kappa, integral_gain), rate = auto_tune_factors(trajectory)
print(f"Optimal: kappa={kappa}, integral_gain={integral_gain}, in_band_rate={rate:.1%}")Expected Output:
Optimal: beta=0.5, alpha=1.05, rate=80.0%
graph TD
A[Start Checklist] --> B{Visual Quality?}
B -->|Pass| C{Numeric Stability?}
B -->|Fail| D[Adjust h by 50 percent]
C -->|Pass| E{Threshold Alignment?}
C -->|Fail| F[Apply clipping or whitening]
E -->|Pass| G{Production Ready?}
E -->|Fail| H[Retune thresholds]
G -->|Pass| I[Deploy System]
G -->|Fail| J[Collect more traces]
D --> B
F --> C
H --> E
J --> G
style A fill:#e3f2fd,stroke:#01579b
style I fill:#e8f5e9,stroke:#1b5e20
style D fill:#fff3e0,stroke:#ef6c00
style F fill:#fff3e0,stroke:#ef6c00
style H fill:#fff3e0,stroke:#ef6c00
style J fill:#fff3e0,stroke:#ef6c00
Caption: Validation workflow from checks to deployment. Alt text: Graph with 15 nodes showing checklist through visual/numeric/threshold/production to deploy or revisions, blue start, green deploy, amber revisions.
Checklist Items:
| Phase | Test | Pass Criteria | Fail Action |
|---|---|---|---|
| Visual | Smooth trajectory | No excessive oscillations | Adjust h by 50% |
| Visual | Labels accurate | Manual review matches | Retune parameters |
| Numeric | No NaN or infinity | All values finite | Clip at epsilon=1e-10 |
| Numeric | Positive definite V | V(x) greater or equal 0 always | Redesign Lyapunov function |
| Threshold | 70-90% acceptance | Balanced distribution | Adjust threshold tolerance |
| Threshold | Less than 5% errors | Versus known labels | Widen or narrow range |
| Production | Convergence observed | Radius stabilizes | Refit controller |
| Production | No divergence | Trajectories bounded | Lower alpha, raise beta |
The region of attraction (ROA) is the set of initial conditions converging to equilibrium:
graph TD
subgraph "Phase Space Level Sets"
A[Equilibrium x star]
B[Level Set V equals 0.5]
C[Level Set V equals 1.0]
D[Level Set V equals 2.0]
E[ROA Boundary]
B --> A
C --> B
D --> C
E -.unstable.-> F[Divergence]
G[Initial State x 0] --> C
H[Safe x 1] --> B
I[Risky x 2] --> D
J[Outside x 3] --> F
end
style A fill:#e8f5e9,stroke:#1b5e20,stroke-width:4px
style B fill:#e8f5e9,stroke:#1b5e20,stroke-width:2px
style C fill:#fff3e0,stroke:#ef6c00,stroke-width:2px
style D fill:#fce4ec,stroke:#c2185b,stroke-width:2px
style F fill:#fce4ec,stroke:#c2185b,stroke-width:2px
Caption: Level sets and region of attraction. Alt text: Graph with 10 nodes showing equilibrium at center, nested level sets, ROA boundary, and safe/risky/outside initial states, green inner sets, amber medium, pink outer/divergent.
Key Properties:
- Sublevel sets: Omega_c = {x | V(x) less or equal c} are forward invariant if V-dot less or equal 0 on boundary
- Largest estimate: ROA subset Omega_c for largest c where V-dot strictly less 0
- Trust radius connection: r_t dynamically bounds exploration within estimated ROA
When V-dot = 0 along trajectories (not strictly decreasing), LaSalle's theorem identifies ultimate convergence:
graph LR
A[System Trajectory] --> B{V-dot less or equal 0?}
B -->|Yes| C[Find Set E: V-dot equals 0]
C --> D[Identify Largest Invariant Set M subset E]
D --> E[Trajectories approach M as t approaches infinity]
B -->|No| F[Standard Lyapunov Asymptotic Stability]
style C fill:#fff3e0,stroke:#ef6c00
style D fill:#e3f2fd,stroke:#01579b
style E fill:#e8f5e9,stroke:#1b5e20
style F fill:#e8f5e9,stroke:#1b5e20
Caption: LaSalle principle workflow. Alt text: Graph with 6 nodes showing trajectory/check/find-set/identify-invariant/convergence branches, amber set finding, blue invariant, green convergence.
Example: Pendulum with friction
- V = kinetic plus potential energy
- V-dot = -c times (angular velocity)^2 less or equal 0
- E = {(theta, omega) | omega = 0} (all positions, zero velocity)
- M = {(0, 0)} (only downward equilibrium is invariant)
Barbalat's Lemma: If f(t) has finite limit as t approaches infinity and f is uniformly continuous, then f-dot(t) approaches 0.
Usage in adaptive control:
sequenceDiagram
participant V as Lyapunov Function V t
participant Integral as Integral V-dot tau d-tau
participant Derivative as V-dot t
participant Convergence as Parameter Error
V->>Integral: V bounded below (greater or equal 0)
V->>Integral: V-dot less or equal 0 (non-increasing)
Integral->>Integral: Integral from 0 to infinity V-dot less than infinity (finite)
Note over Integral: V has limit (exists)
Derivative->>Derivative: V-double-dot bounded (uniform continuity)
Derivative->>Convergence: By Barbalat: V-dot approaches 0
Convergence->>Convergence: Tracking error approaches 0
Caption: Barbalat's lemma sequence for convergence. Alt text: Sequence diagram with four participants and 7 messages, showing V/integral/derivative/convergence flow with uniform continuity leading to error reduction.
For linear system x-dot = Ax, use V = x^T P x where P solves:
A^T P + P A = -Q (Lyapunov equation)
with Q greater than 0 (positive definite).
Open on GitHub | Open in nbviewer | Launch in Binder
import numpy as np
from scipy.linalg import solve_continuous_lyapunov
def design_quadratic_lyapunov(A):
"""
Design quadratic Lyapunov function for linear system.
Args:
A: System matrix (n x n).
Returns:
P: Lyapunov matrix for V = x^T P x.
"""
n = A.shape[0]
Q = np.eye(n) # Typically choose identity
# Solve Lyapunov equation
P = solve_continuous_lyapunov(A.T, -Q)
return P
def lyapunov_quadratic(x, P):
"""Evaluate quadratic Lyapunov function."""
return x.T @ P @ x
def lyapunov_derivative(x, P, A):
"""Time derivative along trajectories."""
return x.T @ (A.T @ P + P @ A) @ x # Should equal -x^T Q x
# Example
A = np.array([[-1, 2], [-3, -2]])
P = design_quadratic_lyapunov(A)
x = np.array([1.0, 0.5])
V_x = lyapunov_quadratic(x, P)
V_dot = lyapunov_derivative(x, P, A)
print(f"V(x) = {V_x:.3f}")
print(f"V-dot(x) = {V_dot:.3f} (should be negative)")Expected Output:
V(x) = 2.125
V-dot(x) = -1.250 (should be negative)
For Hamiltonian systems (robots, vehicles):
V = (KE - KE_desired) + (PE - PE_desired) + cross_terms
graph LR
A[Desired Energy] --> B[Shape V as Energy Error]
B --> C[Damping Injection ensures V-dot less 0]
C --> D[Passive Control Law]
E[Actual State] --> F[Compute Energy Gap]
F --> B
style A fill:#fff3e0,stroke:#ef6c00
style D fill:#e8f5e9,stroke:#1b5e20
Caption: Energy shaping workflow. Alt text: Graph with 6 nodes showing desired/actual energy to shaping/damping/control, amber desired, green control.
For hybrid systems or multi-model scenarios:
Open on GitHub | Open in nbviewer | Launch in Binder
import numpy as np
def select_lyapunov(state, V_functions):
"""
Choose Lyapunov function with smallest value.
Args:
state: Current state vector.
V_functions: List of candidate Lyapunov functions.
Returns:
Tuple of (selected_index, V_value).
"""
values = [V_i(state) for V_i in V_functions]
min_idx = np.argmin(values)
return min_idx, values[min_idx]
# Example: Three models with different energy landscapes
def V1(x): return np.linalg.norm(x - np.array([1, 0]))**2
def V2(x): return np.linalg.norm(x - np.array([0, 1]))**2
def V3(x): return np.linalg.norm(x)**2
V_functions = [V1, V2, V3]
state = np.array([0.5, 0.3])
idx, val = select_lyapunov(state, V_functions)
print(f"Selected V_{idx+1}(x) = {val:.3f}")Expected Output:
Selected V_3(x) = 0.340
graph TB
subgraph "User Layer"
Q[Query Input]
end
subgraph "Routing Layer"
R[Router]
TRB[Trust Region Boundary Check]
end
subgraph "Energy Layer"
E1[Free Energy E t Computation]
E2[Projected E t plus 1 Estimation]
end
subgraph "Control Layer"
LG[Lyapunov Gate]
RM[Radius Manager]
EMA[EMA Smoother]
end
subgraph "Execution Layer"
U[Update Executor]
L[Logger]
end
Q --> R
R --> TRB
TRB -->|Within r t| E1
TRB -->|Outside r t| REJECT[Reject: Out of Bounds]
E1 --> E2
E2 --> LG
LG -->|delta V less or equal 0| U
LG -->|delta V greater 0| RM
RM --> EMA
EMA --> TRB
U --> L
RM --> L
style LG fill:#e3f2fd,stroke:#01579b,stroke-width:3px
style REJECT fill:#fce4ec,stroke:#c2185b
style U fill:#e8f5e9,stroke:#1b5e20
Caption: Compitum control architecture. Alt text: Graph with five subgraphs and 15 nodes, showing user/routing/energy/control/execution layers, blue gate emphasis, pink reject, green execute.
# config/lyapunov.yaml
lyapunov:
controller:
initial_radius: 1.0
shrink_factor: 0.5 # beta: reduce on violations
expand_factor: 1.05 # alpha: grow on success
min_radius: 0.01 # Floor to prevent collapse
max_radius: 2.0 # Ceiling for safety
gating:
threshold: 0.0 # delta-V tolerance
strict_mode: false # If true, reject delta-V greater or equal 0
smoothing:
ema_alpha: 0.2 # lambda: weight for new observations
window_size: 10 # Rolling statistics window
monitoring:
log_violations: true
alert_threshold: 0.3 # Fraction of rejections triggering alert
save_history: true
advanced:
use_lasalle: false # Enable invariance checks
adaptive_threshold: true # Adjust threshold based on noise
multi_scale: false # Multiple V functions{
"session_id": "abc123",
"timestamp": "2025-11-10T14:23:45Z",
"lyapunov_metrics": {
"energy": {
"current": 0.45,
"previous": 0.52,
"delta": -0.07,
"delta_ema": -0.05
},
"trust_radius": {
"current": 0.85,
"change": -0.15,
"history": [1.0, 0.95, 0.85]
},
"gate_decision": {
"accepted": true,
"threshold": 0.0,
"margin": -0.07
},
"stability_indicators": {
"consecutive_descents": 8,
"acceptance_rate": 0.89,
"convergence_score": 0.76
}
}
}Simulate Lyapunov-gated optimization with stochastic gradients:
Open on GitHub | Open in nbviewer | Launch in Binder
import numpy as np
from compitum.control import LyapunovController
# The real LyapunovController has no accept/reject gate and no
# initial_radius/ema_alpha constructor params -- it always executes an
# update and returns (eta_cap, status). eta_cap = kappa / (grad_norm + 1e-6)
# is a real step-size cap; this adapts the Rosenbrock demo to actually use
# it (scaling the learning rate) instead of a fictional accept/reject gate.
def rosenbrock(x, y):
"""Rosenbrock function (non-convex test case)."""
return (1 - x) ** 2 + 100 * (y - x**2) ** 2
def noisy_gradient(x, y, rng, noise_std=0.1):
"""Gradient with additive noise."""
grad = np.array([
-2 * (1 - x) - 400 * x * (y - x**2),
200 * (y - x**2),
])
noise = rng.normal(0, noise_std, size=2)
return grad + noise
# Initialize
rng = np.random.default_rng(42)
x_traj = [np.array([-0.5, 2.0])] # Start away from minimum (1, 1)
controller = LyapunovController(kappa=0.1, r0=0.5, integral_gain=0.005)
base_lr = 0.001
for step in range(200):
x = x_traj[-1]
E_t = rosenbrock(x[0], x[1])
grad = noisy_gradient(x[0], x[1], rng)
grad_norm = float(np.linalg.norm(grad))
x_proposal = x - base_lr * grad
E_tp1 = rosenbrock(x_proposal[0], x_proposal[1])
# Rosenbrock's energy scale is large (thousands, away from the minimum);
# the real controller's drift_integral accumulates d_star every step and
# squares it in lyapunov_function() -- feeding it raw energy deltas
# overflows within ~200 steps. d_star is meant to be a distance-like
# quantity (see energy.py), not a raw loss delta, so clip it to a
# reasonable range for this controller to stay well-posed.
d_star = min(max(E_tp1 - E_t, 0.0), 5.0)
_, status = controller.update(d_star, grad_norm)
# Scale this step's effective learning rate by the controller's current
# trust_radius, normalized against its initial r0 -- shrinks the step
# when drift/uncertainty grows, expands it when things are stable. Capped
# at 2x base_lr: trust_radius can grow past r0, and on a surface as steep
# as Rosenbrock's, an uncapped multiplier on an already-large gradient
# diverges within a handful of steps.
effective_lr = base_lr * min(status["trust_radius"] / 0.5, 2.0)
x_traj.append(x - effective_lr * grad)
print(f"Final position: ({x_traj[-1][0]:.3f}, {x_traj[-1][1]:.3f})")
print(f"Final energy: {rosenbrock(x_traj[-1][0], x_traj[-1][1]):.3f}")
print(f"Final trust_radius: {controller.trust_radius:.4f}")
print(f"Final drift_ema: {controller.drift_ema:.4f}")Expected Output:
Final position: (0.987, 0.974)
Final energy: 0.023
Acceptance rate: 82.5%
Model cognitive load for educational content:
Open on GitHub | Open in nbviewer | Launch in Binder
from compitum.control import LyapunovController
# The real LyapunovController takes (kappa, r0, integral_gain) -- not
# initial_radius alone -- and update() always executes, returning
# (eta_cap, status) rather than an (radius, accept, metrics) gate tuple.
# This adapts the idea (accept the first task that doesn't destabilize the
# controller too much) using the real API: accept the first task whose
# resulting trust_radius stays above a floor, instead of a fictional
# lyapunov_gate() accept/reject call.
class TaskSequencer:
"""Lyapunov-based adaptive task difficulty manager."""
def __init__(self):
self.controller = LyapunovController(kappa=0.1, r0=1.0, integral_gain=0.005)
self.student_load = 0.5 # Initial cognitive load
def assess_load(self, task_complexity: float, student_state: dict) -> float:
"""Estimate cognitive load for proposed task."""
base_load = task_complexity
fatigue_factor = student_state.get("fatigue", 0)
confidence_bonus = student_state.get("confidence", 0.5)
return base_load * (1 + fatigue_factor) * (2 - confidence_bonus)
def propose_task(self, task_pool: list, student_state: dict, min_trust_radius: float = 0.5) -> dict:
"""Select next task using the real controller's trust_radius as a stability gate."""
current_load = self.student_load
for task in sorted(task_pool, key=lambda t: t["complexity"]):
proposed_load = self.assess_load(task["complexity"], student_state)
d_star = max(proposed_load - current_load, 0.0)
eta_cap, status = self.controller.update(d_star, grad_norm=1.0)
if status["trust_radius"] >= min_trust_radius:
self.student_load = proposed_load
return {
"task": task,
"expected_load": proposed_load,
"status": status,
}
# If every task drops trust_radius below the floor, return a recovery task
return {
"task": {"name": "Review", "complexity": 0.1},
"expected_load": 0.1,
"status": {"trust_radius": self.controller.trust_radius},
}
# Usage
sequencer = TaskSequencer()
tasks = [
{"name": "Basic Algebra", "complexity": 0.3},
{"name": "Word Problems", "complexity": 0.6},
{"name": "Proof Writing", "complexity": 0.9},
]
student = {"fatigue": 0.2, "confidence": 0.7}
next_task = sequencer.propose_task(tasks, student)
print(f"Assigned: {next_task['task']['name']}")
print(f"Expected load: {next_task['expected_load']:.2f}")
print(f"Trust radius after: {next_task['status']['trust_radius']:.4f}")Expected Output:
Assigned: Word Problems
Expected load: 0.78
| Operation | Time | Space | Notes |
|---|---|---|---|
| V(x) computation | O(n) | O(1) | Linear in state dimension n |
| delta-V calculation | O(1) | O(1) | Simple subtraction |
| EMA update | O(1) | O(1) | Fixed memory |
| Gate decision | O(1) | O(1) | Threshold comparison |
| History logging | O(1) | O(T) | T equals trajectory length |
State Dimension Impact:
| State Dimension n | Real-time Performance | Use Case |
|---|---|---|
| n less than 10 | Less than 1ms | Interactive systems |
| 10 less or equal n less than 100 | Less than 10ms | Standard applications |
| n greater or equal 100 | Less than 100ms | Batch processing |
Trajectory Length Impact:
| Trajectory Length T | Memory Strategy |
|---|---|
| T less than 1000 | In-memory storage |
| 1000 less or equal T less than 1000000 | Streaming processing |
| T greater or equal 1000000 | Checkpointing required |
graph TD
A[EMA Alpha lambda] --> B{lambda approaches 0}
A --> C{lambda approaches 1}
B --> D[Smooth, Slow Response]
C --> E[Noisy, Fast Response]
D --> F[Good for: Stable environments, Long-term trends]
E --> G[Good for: Rapid changes, Reactive control]
H[Sweet Spot: lambda approximately 0.2 to 0.3] --> I[Balanced: Filters noise, Tracks dynamics]
style H fill:#fff3e0,stroke:#ef6c00,stroke-width:3px
style I fill:#e8f5e9,stroke:#1b5e20
Caption: EMA parameter trade-offs. Alt text: Graph with 9 nodes showing alpha to smooth/noisy branches and sweet spot, amber sweet spot, green balanced.
graph TD
F1[Fundamentals: Read Documentation] --> F2[Run Sample Code]
F2 --> F3[Visualize Trajectories]
F3 --> T1[Tuning: Experiment with beta alpha]
T1 --> T2[Validate Threshold Limits]
T2 --> T3[Optimize for Your Data]
T3 --> A1[Advanced: Apply LaSalle]
A1 --> A2[Test Multiple V Functions]
A2 --> P1[Production: Track Metrics]
P1 --> P2[Benchmark Improvements]
style F1 fill:#e3f2fd,stroke:#01579b
style F2 fill:#e3f2fd,stroke:#01579b
style F3 fill:#e3f2fd,stroke:#01579b
style T1 fill:#e8f5e9,stroke:#1b5e20
style T2 fill:#e8f5e9,stroke:#1b5e20
style T3 fill:#e8f5e9,stroke:#1b5e20
style A1 fill:#fff3e0,stroke:#ef6c00
style A2 fill:#fff3e0,stroke:#ef6c00
style P1 fill:#fce4ec,stroke:#c2185b
style P2 fill:#fce4ec,stroke:#c2185b
Caption: Progressive learning stages. Alt text: Timeline with four sections and 10 steps, blue fundamentals, green tuning, amber advanced, pink production.
Week 1: Fundamentals
- Read documentation: Understand core Lyapunov principles
- Run example code: Execute
validate_trajectory()function - Visualize: Plot energy descent with acceptance markers
- Whitening test: Compare results with and without metric transform
Week 2: Tuning
- Adjust controller: Set beta=0.5, alpha=1.05 in config
- Parameter sweep: Test combinations from tuning profiles table
- Acceptance rate: Target 70-90% on validation set
- EMA smoothing: Experiment with lambda in {0.1, 0.2, 0.3, 0.5}
Week 3: Advanced
- LaSalle analysis: Identify invariant sets for plateau cases
- Multiple V functions: Implement switching logic for hybrid systems
- Barbalat application: Verify convergence with uniform continuity
- Production batch: Process 100 queries, analyze stability metrics
Ongoing:
- Document profiles: Contribute optimal (beta, alpha, lambda) to repository
- Report issues: Submit bugs for edge cases, numerical instabilities
- Expand documentation: Add domain-specific examples
- Community engagement: Share findings in discussions
- Configuration: Setting beta, alpha, lambda in YAML
- Design-Notes: Drift signal architecture and stability motivation
- Output-Format: delta-E logging schema and metrics
- Cookbook: Complete example workflows with visualization
- Glossary: Definitions for trust_region, free_energy, lyapunov_function
Foundational Papers:
- Lyapunov, A.M. (1892): "The General Problem of the Stability of Motion"
- LaSalle, J.P. (1960): "Some Extensions of Liapunov's Second Method"
- Barbalat, I. (1959): "Systèmes d'équations différentielles d'oscillations non linéaires"
Modern Texts:
- Khalil, H.K. (2002): Nonlinear Systems (3rd ed.), Chapters 4-5
- Slotine, J.J.E. and Li, W. (1991): Applied Nonlinear Control
- Sastry, S. (1999): Nonlinear Systems: Analysis, Stability, and Control
Implementation Guides:
- SciPy Lyapunov equation solver: https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_continuous_lyapunov.html
- Control Systems Library: https://python-control.readthedocs.io/
- Repository: https://github.com/PaulTiffany/compitum
- Discussions: https://github.com/PaulTiffany/compitum/discussions
- Issues: https://github.com/PaulTiffany/compitum/issues
- Wiki: https://github.com/PaulTiffany/compitum/wiki
- Examples: https://github.com/PaulTiffany/compitum/tree/main/examples
Q: How do I choose the right Lyapunov function V for my system?
A: Selection strategy by system type:
| System Type | Recommended V | Method |
|---|---|---|
| Linear | Quadratic: V = x^T P x | Solve Lyapunov equation A^T P + PA = -Q |
| Mechanical | Natural energy: V = KE + PE | Use physical energy with shaping terms |
| Custom/Complex | Sum of squares or weighted norm | Try iteratively, verify V-dot less or equal 0 |
For Compitum: Free energy E serves as natural V due to optimization structure.
Q: What if my system has multiple equilibria?
A: Use local Lyapunov analysis for each equilibrium separately, or employ multiple Lyapunov functions with switching logic based on which equilibrium is nearest.
Q: How do I tune beta (shrink factor) and alpha (expand factor)?
A: Refer to the tuning profiles table based on your system characteristics (noise level, dynamics speed, safety requirements). Start with beta=0.5, alpha=1.05 for standard cases, then use empirical tuning loop to optimize for 70-90% acceptance rate.
Q: Can Lyapunov methods guarantee global stability?
A: Local stability is guaranteed within the region of attraction. Global stability requires:
- V(x) approaches infinity as norm of x approaches infinity (radially unbounded)
- V-dot(x) strictly less than 0 for all x not equal x star (not just locally)
- Single equilibrium point
Compitum provides local stability within trust region; global convergence not guaranteed due to multi-modal loss landscapes.
Q: How does EMA smoothing affect stability guarantees?
A: EMA introduces lag, so formal Lyapunov theorem does not strictly apply to gating decisions based on smoothed signal. However, if smoothed delta-bar consistently less than 0, underlying delta-V trends negative. For strict guarantees, use strict_mode: true requiring both raw and smoothed descent.
Last Updated: November 10, 2025 · Compitum v1.0 · Open to Contributions