feat: Initialization of Cyber-Physical Systems library , with Lyapuno…

Cslib.lean

@@ -2,6 +2,8 @@ module  -- shake: keep-all
 
 public import Cslib.Algorithms.Lean.MergeSort.MergeSort
 public import Cslib.Algorithms.Lean.TimeM
+public import Cslib.CPS.DiscreteLinearTime.AsymptoticStability
+public import Cslib.CPS.DiscreteLinearTime.Basic
 public import Cslib.Computability.Automata.Acceptors.Acceptor
 public import Cslib.Computability.Automata.Acceptors.OmegaAcceptor
 public import Cslib.Computability.Automata.DA.Basic

Cslib/CPS/DiscreteLinearTime/AsymptoticStability.lean

@@ -0,0 +1,255 @@
+/-
+Copyright (c) 2026 Bashar Hamade. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Bashar Hamade
+-/
+
+module
+
+public import Cslib.Init
+
+
+public import Mathlib.Analysis.Normed.Group.Basic
+public import Mathlib.Analysis.Normed.Module.RCLike.Real
+public import Mathlib.Analysis.Normed.Operator.ContinuousLinearMap
+public import Mathlib.Analysis.Complex.Order
+public import Mathlib.Analysis.Normed.Algebra.Spectrum
+public import Mathlib.Order.Basic
+public import Mathlib.Analysis.Normed.Algebra.GelfandFormula
+public import Cslib.CPS.DiscreteLinearTime.Basic
+
+@[expose] public section
+
+/-!
+# Asymptotic Stability Analysis
+
+This module analyzes the asymptotic stability of discrete-time linear systems.
+It uses the Gelfand spectral radius formula to show that if the
+spectral radius of the system matrix is less than one,
+the state converges to zero under zero input.
+
+## Main Theorems
+* `DiscreteLinearSystemState.asymptotic_stability_discrete`: The main stability theorem.
+* `DiscreteLinearSystemState.gelfand_le_one_when_spectral_radius_le_one`:
+Relates spectral radius to Gelfand limit.
+
+## References
+https://en.wikipedia.org/wiki/Lyapunov_stability
+https://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-public_24Jul2020.pdf
+-/
+
+
+
+
+open scoped ComplexOrder
+
+variable {σ ι : Type*}
+variable [NormedAddCommGroup σ] [NormedSpace ℂ σ]
+variable [NormedAddCommGroup ι] [NormedSpace ℂ ι]
+
+namespace DiscreteLinearSystemState
+
+/-- System evolution function. -/
+noncomputable def systemEvolutionAlt (sys : DiscreteLinearSystemState σ ι) : ℕ → σ
+  | 0 => sys.x₀
+  | k + 1 => sys.A (systemEvolutionAlt sys k) + sys.B (sys.u k)
+
+/-- Discrete State space representation property. -/
+def stateSystemEquation (sys : DiscreteLinearSystemState σ ι) : Prop :=
+  ∀ k : ℕ, sys.x (k + 1) = sys.A (sys.x k) + sys.B (sys.u k)
+
+/-- Power definition: A^(k+1) = A^k * A. -/
+lemma system_power_multiplication (a : σ →L[ℂ] σ) (k : ℕ) :
+    a ^ (k + 1) = (a ^ k).comp a := by
+  induction k with
+  | zero =>
+    simp only [zero_add, pow_one, pow_zero]
+    exact ContinuousLinearMap.id_comp a
+  | succ k ih =>
+    rw [pow_succ]
+    rw [ih]
+    rfl
+
+/-- Power definition commutation: A^(k+1) = A * A^k. -/
+lemma system_power_multiplication_flopped (a : σ →L[ℂ] σ) (k : ℕ) :
+    a ^ (k + 1) = a.comp (a^k) := by
+  induction k with
+  | zero =>
+    simp only [zero_add, pow_one, pow_zero]
+    exact ContinuousLinearMap.id_comp a
+  | succ k ih =>
+    rw [pow_succ]
+    rw [ih]
+    simp only [←ContinuousLinearMap.mul_def]
+    rw [mul_assoc]
+    congr 1
+
+
+
+
+
+/-- Lemma: State evolution under zero input. -/
+lemma state_evolution_zero_input (sys : DiscreteLinearSystemState σ ι)
+ (h_init : sys.x 0 = sys.x₀)
+ (h_state : stateSystemEquation sys)
+ (h_zero_input : sys.u = zeroInput) :
+ ∀ k, sys.x k = (sys.A ^ k) sys.x₀ := by
+  intro k
+  induction k with
+ | zero =>
+   simp [pow_zero, h_init]
+ | succ k ih =>
+   have h1 : sys.x (k + 1) = sys.A (sys.x k) + sys.B (sys.u k) := h_state k
+   rw [ih] at h1
+   rw [h_zero_input] at h1
+   unfold zeroInput at h1
+   simp only [Function.const, ContinuousLinearMap.map_zero, add_zero] at h1
+   rw [h1]
+   rw [←ContinuousLinearMap.comp_apply]
+   congr 1
+   rw [system_power_multiplication_flopped]
+
+
+
+
+
+
+/-- Bound on state norm. -/
+theorem bound_x_norm
+    (sys : DiscreteLinearSystemState σ ι)
+    (hx : ∀ k, sys.x k = (sys.A ^ k) sys.x₀) :
+    ∀ k, ‖sys.x k‖ ≤ ‖sys.A ^ k‖ * ‖sys.x₀‖ := by
+  intro k
+  rw [hx k]
+  exact ContinuousLinearMap.le_opNorm (sys.A ^ k) sys.x₀
+
+/-- Spectral radius strictly less than one. -/
+def spectralRadiusLTOne (a : σ →L[ℂ] σ) : Prop :=
+  spectralRadius ℂ a < 1
+
+/-- If the spectral radius of `a` is less than one, then the
+Gelfand formula limits to less than one. -/
+theorem gelfand_le_one_when_spectral_radius_le_one
+    [CompleteSpace σ] (a : σ →L[ℂ] σ) :
+    spectralRadiusLTOne a →
+    Filter.limsup (fun (n : ℕ) => (‖a ^ n‖₊ : ENNReal) ^ (1 / n : ℝ)) Filter.atTop < 1 := by
+  intro h_spectral
+  unfold spectralRadiusLTOne at h_spectral
+  have h_gelfand := spectrum.limsup_pow_nnnorm_pow_one_div_le_spectralRadius a
+  convert lt_of_le_of_lt h_gelfand h_spectral
+
+
+
+/-- If the Gelfand limit is < 1, then eventually the root-term is bounded by some `r < 1`. -/
+theorem gelfand_eventually_bounded
+    (a : σ →L[ℂ] σ) (h : Filter.limsup (fun n :
+    ℕ ↦ (‖a ^ n‖₊ : ENNReal) ^ (1 / n : ℝ)) Filter.atTop < 1) :
+    ∃ (r : ENNReal) (N : ℕ), 0 < r ∧ r < 1
+    ∧ ∀ (k : ℕ), N < k → (‖a ^ k‖₊ : ENNReal) ^ (1 / k : ℝ) < r :=
+by
+  obtain ⟨r, h_limsup_lt_r, h_r_lt_one⟩ := exists_between h
+  have r_pos : 0 < r := lt_of_le_of_lt (zero_le _) h_limsup_lt_r
+  have eventually_lt := Filter.eventually_lt_of_limsup_lt h_limsup_lt_r
+  rw [Filter.eventually_atTop] at eventually_lt
+  obtain ⟨N, hN⟩ := eventually_lt
+  use r, N
+  refine ⟨r_pos, h_r_lt_one, fun k hk => hN k (Nat.le_of_lt hk)⟩
+
+
+/-- If terms are eventually bounded by `r^k` with `r < 1`, the sequence tends to zero. -/
+theorem power_to_zero (sys : DiscreteLinearSystemState σ ι)
+  (r : ENNReal) (N : ℕ)
+  (r_pos : 0 < r)
+  (r_lt_one : r < 1)
+  (h_power : ∀ (k : ℕ), N < k → (‖sys.A ^ k‖₊ : ENNReal) < r ^ k) :
+  Filter.Tendsto (fun k => ‖sys.A ^ k‖) Filter.atTop (nhds 0) := by
+  have r_real_zero : Filter.Tendsto (fun n => (r ^ n).toReal) Filter.atTop (nhds 0) := by
+    simp only [ENNReal.toReal_pow, tendsto_pow_atTop_nhds_zero_iff, ENNReal.abs_toReal]
+    cases r with
+    | top => simp
+    | coe x =>
+      simp only [ENNReal.coe_toReal, NNReal.coe_lt_one]
+      exact ENNReal.coe_lt_one_iff.mp r_lt_one
+  rw [Metric.tendsto_atTop]
+  intro ε ε_pos
+  obtain ⟨N₁, hN₁⟩ := Metric.tendsto_atTop.mp r_real_zero ε ε_pos
+  use max N N₁ + 1
+  intro k hk
+  have hkN : N < k := lt_of_le_of_lt (le_max_left N N₁) (Nat.lt_of_succ_le hk)
+  have hkN₁ : N₁ ≤ k := le_trans (le_max_right N N₁) (Nat.le_of_succ_le hk)
+  have h_bound := h_power k hkN
+  have h_r_small := hN₁ k hkN₁
+  simp only [ gt_iff_lt]
+  simp only [ENNReal.toReal_pow] at h_r_small
+  have h_r_ennreal : (r ^ k).toReal < ε := by
+    rw [ENNReal.toReal_pow]
+    rw [Real.dist_eq] at h_r_small
+    simp only [sub_zero, abs_pow, ENNReal.abs_toReal] at h_r_small
+    exact h_r_small
+  have h_finite : (r ^ k) ≠ ⊤ := by
+    simp only [ne_eq, ENNReal.pow_eq_top_iff]
+    intro h
+    cases h with
+    | intro h_r_top h_k_ne_zero =>
+      have : r < ⊤ := r_lt_one.trans_le le_top
+      exact ne_of_lt this h_r_top
+  rw [Real.dist_eq]
+  simp only [sub_zero, abs_norm, gt_iff_lt]
+  calc ‖sys.A ^ k‖
+  _ = (↑‖sys.A ^ k‖₊ : ENNReal).toReal := by simp
+  _ < (r ^ k).toReal := (ENNReal.toReal_lt_toReal ENNReal.coe_ne_top h_finite).mpr h_bound
+  _ < ε := h_r_ennreal
+
+
+
+
+
+
+
+/-- Main theorem: Asymptotic Stability. If the spectral radius is < 1,
+the zero-input response converges to zero. -/
+theorem asymptotic_stability_discrete [CompleteSpace σ] (sys : DiscreteLinearSystemState σ ι)
+  (h_init : sys.x 0 = sys.x₀)
+  (h_state : stateSystemEquation sys)
+  (h_zero_input : sys.u = zeroInput)
+  (h_spectral : spectralRadiusLTOne sys.A) :
+  Filter.Tendsto (fun k => ‖sys.x k‖) Filter.atTop (nhds 0) := by
+  have h_gelfand := spectrum.limsup_pow_nnnorm_pow_one_div_le_spectralRadius sys.A
+  have h_gelfand_le_one : Filter.limsup (fun (n : ℕ) =>
+  (‖sys.A ^ n‖₊ : ENNReal) ^ (1 / n : ℝ)) Filter.atTop < 1 := by
+      unfold spectralRadiusLTOne at h_spectral
+      refine lt_of_le_of_lt ?_ h_spectral
+      exact h_gelfand
+  have eventually_bounded := gelfand_eventually_bounded sys.A h_gelfand_le_one
+  obtain ⟨r, N, r_pos, r_lt_one, h_bound⟩ := eventually_bounded
+  have h_power : ∀ (k : ℕ), N < k → ↑‖sys.A ^ k‖₊ < r ^ k := by
+      intros k' hk'
+      specialize h_bound k' hk'
+      have h_k'_pos : 0 < k' := Nat.zero_lt_of_lt hk'
+      have h_inv_k' : (k' : ℝ)⁻¹ * k' = 1 := by
+        field_simp
+      have h_pow : (↑‖sys.A ^ k'‖₊ ^ (1 / k' : ℝ)) ^ (k' : ℝ) < r ^ k' := by
+        rw [← ENNReal.rpow_natCast r k']
+        exact ENNReal.rpow_lt_rpow h_bound (Nat.cast_pos.mpr h_k'_pos)
+      rw [← ENNReal.rpow_natCast, ← ENNReal.rpow_mul] at h_pow
+      simp only [one_div, ENNReal.rpow_natCast] at h_pow
+      rw [h_inv_k'] at h_pow
+      simp only [ENNReal.rpow_one] at h_pow
+      exact h_pow
+  have hx : ∀ k, sys.x k = (sys.A ^ k) sys.x₀ :=
+      state_evolution_zero_input sys h_init h_state h_zero_input
+  have pow_zero : Filter.Tendsto (fun k => ‖sys.A ^ k‖) Filter.atTop (nhds 0) :=
+      power_to_zero sys r N r_pos r_lt_one h_power
+  simp only [hx]
+  have h_norm_eq : ∀ k, ‖(sys.A ^ k) sys.x₀‖ ≤ ‖sys.A ^ k‖ * ‖sys.x₀‖ :=
+    fun k => ContinuousLinearMap.le_opNorm (sys.A ^ k) sys.x₀
+  have h_prod_zero : Filter.Tendsto (fun k => ‖sys.A ^ k‖ * ‖sys.x₀‖) Filter.atTop (nhds 0) := by
+    rw [← zero_mul ‖sys.x₀‖]
+    exact Filter.Tendsto.mul_const ‖sys.x₀‖ pow_zero
+  exact tendsto_of_tendsto_of_tendsto_of_le_of_le
+    (tendsto_const_nhds)
+    h_prod_zero
+    (fun k => norm_nonneg _)
+    h_norm_eq
+
+end DiscreteLinearSystemState

Cslib/CPS/DiscreteLinearTime/Basic.lean

@@ -0,0 +1,105 @@
+ /-
+Copyright (c) 2026 Bashar Hamade. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Bashar Hamade
+-/
+
+module
+
+public import Cslib.Init
+public import Mathlib.Analysis.Normed.Module.Basic
+public import Mathlib.Analysis.Normed.Operator.ContinuousLinearMap
+public import Mathlib.LinearAlgebra.Span.Basic
+public import Mathlib.LinearAlgebra.FiniteDimensional.Defs
+public import Mathlib.Data.Finset.Basic
+public import Mathlib.Data.Complex.Basic
+public import Mathlib.Order.CompleteLattice.Basic
+public import Mathlib.Analysis.Complex.Order
+public import Mathlib.Analysis.Normed.Field.Lemmas
+public import Mathlib.Analysis.Normed.Field.Basic
+public import Mathlib.Analysis.Complex.Exponential
+public import Mathlib.Analysis.Normed.Group.Basic
+public import Mathlib.Analysis.Normed.Module.RCLike.Real
+public import Mathlib.Analysis.Normed.Algebra.Spectrum
+public import Mathlib.Order.Basic
+public import Mathlib.LinearAlgebra.Charpoly.Basic
+public import Mathlib.LinearAlgebra.Matrix.Charpoly.LinearMap
+
+@[expose] public section
+
+open scoped ComplexOrder
+
+
+set_option linter.style.emptyLine false
+set_option linter.deprecated.module false
+
+universe u v
+
+/-!
+# Basic definitions for Discrete Linear Time Systems
+
+This module defines the state space representation of a discrete-time linear dynamical system.
+It includes the definition of the system state, the evolution function,
+and the property of satisfying the state equation.
+
+## Main Definitions
+* `DiscreteLinearSystemState`: Structure representing the system matrices (A and B),
+the current state, input, and initial state.
+* `DiscreteLinearSystemState.systemEvolution`: Function computing the state at time `k`
+given an input sequence.
+* `DiscreteLinearSystemState.satisfies_state_equation`: Proposition stating that
+the sequence `x` satisfies the linear difference equation `x(k+1) = A x(k) + B u(k)`.
+
+
+## References
+
+https://en.wikipedia.org/wiki/State-space_representation
+https://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-public_24Jul2020.pdf
+
+-/
+
+
+
+variable {σ : Type u} {ι : Type v}
+variable [TopologicalSpace σ] [NormedAddCommGroup σ] [NormedSpace ℂ σ]
+variable [TopologicalSpace ι] [NormedAddCommGroup ι] [NormedSpace ℂ ι]
+
+/-- Discrete-time linear dynamical system with state equation x(k+1) = A·x(k) + B·u(k). -/
+structure DiscreteLinearSystemState (σ : Type u) (ι : Type v)
+    [TopologicalSpace σ] [NormedAddCommGroup σ] [NormedSpace ℂ σ]
+    [TopologicalSpace ι] [NormedAddCommGroup ι] [NormedSpace ℂ ι] where
+  /-- State transition matrix (A), mapping the current state to the next state component (n×n). -/
+  A : σ →L[ℂ] σ
+  /-- Input matrix (B), mapping the current input to the next state component (n×p). -/
+  B : ι →L[ℂ] σ
+  /-- Initial state -/
+  x₀ : σ
+  /-- State sequence -/
+  x : ℕ → σ
+  /-- Input sequence -/
+  u : ℕ → ι
+
+namespace DiscreteLinearSystemState
+
+variable {sys : DiscreteLinearSystemState σ ι}
+
+/-- System evolution function from initial state -/
+noncomputable def systemEvolution (u : ℕ → ι) : ℕ → σ
+  | 0 => sys.x₀
+  | k + 1 => sys.A (systemEvolution u k) + sys.B (u k)
+
+/-- Discrete state space representation property -/
+def satisfies_state_equation : Prop :=
+  ∀ k : ℕ, sys.x (k + 1) = sys.A (sys.x k) + sys.B (sys.u k)
+
+
+/-- Evolution from zero initial state with given input -/
+noncomputable def evolveFromZero
+   (u : ℕ → ι) (sys : DiscreteLinearSystemState σ ι) : ℕ → σ
+  | 0 => 0
+  | k + 1 => sys.A (evolveFromZero u sys k) + sys.B (u k)
+
+/-- Zero input sequence -/
+def zeroInput : ℕ → ι := Function.const _ 0
+
+end DiscreteLinearSystemState

references.bib

@@ -254,3 +254,14 @@ @article{         ShepherdsonSturgis1963
   publisher    = {Association for Computing Machinery},
   address      = {New York, NY, USA}
 }
+
+@book{            AstromMurray2008,
+  author       = {Karl J. {\AA}str{\"o}m and Richard M. Murray},
+  title        = {Feedback Systems: An Introduction for Scientists and Engineers},
+  year         = {2008},
+  publisher    = {Princeton University Press},
+  address      = {Princeton, NJ},
+  url          = {https://www.cds.caltech.edu/~murray/books/AM08/pdf/fbs-public_24Jul2020.pdf},
+  urldate      = {2020-07-24},
+  isbn         = {978-0691135762}
+}