state-space model

state-space model

訊號學家自創一個詞彙state-space model。
嚴格來說是三件事情，但是整包稱作state-space model：
一、系統追加內部狀態（internal state）。
二、多項式推廣為矩陣（matrix representation）。
三、特殊系統改寫成矩陣（state-space representation /
　　　　　　　　　　　　phase-variable representation）。
本章介紹一。下章介紹二三。

state原義是指特定時刻，各個變數的數值，所組成的數組。
state space原義是指所有時刻，全部的數組。
上述三件事情，僅第三件事情的一部分跟state space有關。
訊號學家沒有考慮清楚，胡亂造詞。成為歷史共業。

system with internal state

system with internal state (no input signal)

system:            system with internal state:

y = f(θ)           ⎰ x +⃡ 1 = f(x)
                   ⎱ y     = g(x)

                   ╭╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╮
┌───┐              ╎  x ┌───┐     ┌───┐ ╎
│ f │──→ y         ╎ ┌─→│ f │──┬──│ g │──→ y
└───┘              ╎ │  └───┘  │  └───┘ ╎
                   ╎ └─────────┘        ╎
                   ╰╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╯
                              f'

                   x: internal state
y: output signal   y: output signal

系統g的輸入訊號是遞迴數列x。
系統g的輸出訊號是數列y。

y是輸出訊號。g是系統（開迴路）。
x是內部狀態。f是另一個系統（閉迴路）。
整體可以視作一個生成器f'。

數學當中，函數必須擁有輸入變數。
訊號學當中，系統不一定有輸入訊號x，但是通常有參數θ。
嚴謹起見，生成器的數學式子必須明確寫出參數θ，當作輸入變數。

system with internal state

system:             system with internal state:

y = f(x)            ⎰ x +⃡ 1 = f(x, u)
                    ⎱ y     = g(x, u)

                       ╭╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╮
     ┌───┐          u ────┬─────────┐        ╎
x ──→│ f │──→ y        ╎  └─→┌───┐  └─→┌───┐ ╎
     └───┘             ╎  ┌─→│ f │──┬──│ g │──→ y
                       ╎ x│  └───┘  │  └───┘ ╎
                       ╎  └─────────┘        ╎
                       ╰╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╯
                                   f'

                    x: internal state
x: input signal     u: input signal
y: output signal    y: output signal

u是輸入訊號。y是輸出訊號。g是系統（開迴路）。
x是內部狀態。f是另一個系統（閉迴路）。

系統g追加第二道輸入訊號x。
複製一份u，套用另一個系統f，當作第二道輸入訊號x。

符號意義被更動，u和x地位對調，f和g地位對調。
阿就古人當初沒有考慮清楚。成為歷史共業。

LTI system with internal state (no input signal)

LTI system:         LTI system with internal state:

y = a ∗ θ           ⎰ x +⃡ 1 = a ∗ x
                    ⎱ y     = c ∗ x

LTI system with internal state

LTI system:         LTI system with internal state:

y = a ∗ x           ⎰ x +⃡ 1 = a ∗ x + b ∗ u
                    ⎱ y     = c ∗ x + d ∗ u

範例

sensor (a = d = zero)

   ╭╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╮
 u ╎ ┌───────────┐  x  ┌────────────┐ ╎ y
────→│ f ≜ plant │────→│ g ≜ sensor │────→
   ╎ └───────────┘     └────────────┘ ╎
   ╰╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╯

使用感測器測量系統。
真實輸入訊號u，真實輸出訊號x，測量而得的輸出訊號y。
x無法得知，視作內部狀態。名副其實。

random number generator (b = d = zero and c = identity)

╭╌╌╌╌╌╌╌╌╌╌╌╌╌╮
╎  x ┌───┐    ╎
╎ ┌─→│ f │──┬──→ y
╎ │  └───┘  │ ╎
╎ └─────────┘ ╎
╰╌╌╌╌╌╌╌╌╌╌╌╌╌╯

密碼學當中，隨機數生成器經常採用AR system。
如果還想實施線性變換，那就接上MA system。不過一般來說沒有必要。
如果還想利用密鑰影響隨機數，那就接上輸入訊號。

inertial motion measurement (b = d = zero)

╭╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╮
╎  x⃗ ≜ ⎡x⎤                 ╎
╎      ⎣v⎦ ┌───┐     ┌───┐ ╎
╎ ┌───────→│ f │──┬──│ g │──→ y ≜ v
╎ │        └───┘  │  └───┘ ╎
╎ └───────────────┘        ╎
╰╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╯

⎰ x[n+1] = v[n] Δt + x[n]     牛頓運動定律（慣性運動）
⎱ y[n] = v[n]                 感測器（測速槍）

觀察自然現象x⃗，使用測量儀器g，收集數據y⃗。
自然現象源自自迴歸模型f，合乎科學定律。
測量儀器是移動平均模型g，合乎電路原理。

x是物體位置，v是物體速度。
x和v都是內部狀態，兩道訊號拼成向量x⃗，形成MIMO system。
y是用測速槍觀察到的速度。y是輸出訊號。

accelerated motion measurement (d = zero)

       ╭╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╮
u ≜ a ────┐                  ╎
       ╎  └─→┌───┐     ┌───┐ ╎
       ╎  ┌─→│ f │──┬──│ g │──→ y ≜ v
       ╎ x⃗│  └───┘  │  └───┘ ╎
       ╎  └─────────┘        ╎
       ╰╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╯

⎰ x[n+1] = v[n] Δt + x[n] + a[n]    牛頓運動定律（施力）
⎱ y[n] = v[n]                       感測器（測速槍）

a是加速度，來自外力。a = f/m。
a是輸入訊號。

system model

discrete-time system / continuous-time system

discrete-time system:
⎰ x[n+1] = fₙ(x[n], ..., x[n-p], u[n], ..., u[n-q])
⎱ y[n]   = gₙ(x[n], ..., x[n-r], u[n], ..., u[n-s])

continuous-time system:
⎰ x′(t) = f(t, x(t), u(t))
⎱ y(t)  = g(t, x(t), u(t))

causal system / time-invariant system / linear system

time-variant causal system:
⎰ x[n+1] = fₙ(x[n], ..., x[n-p], u[n], ..., u[n-q])
⎱ y[n]   = gₙ(x[n], ..., x[n-r], u[n], ..., u[n-s])

time-invariant causal system:
⎰ x[n+1] = f(x[n], ..., x[n-p], u[n], ..., u[n-q])
⎱ y[n]   = g(x[n], ..., x[n-r], u[n], ..., u[n-s])

linear time-variant causal system:
⎧ x[n+1] = aₙ₀ x[n] + aₙ₁ x[n-1] + ... + aₙₚ x[n-p]
⎨        + bₙ₀ u[n] + bₙ₁ u[n-1] + ... + bₙ₉ u[n-q]
⎪ y[n]   = cₙ₀ x[n] + cₙ₁ x[n-1] + ... + cₙᵣ x[n-r]
⎩        + dₙ₀ u[n] + dₙ₁ u[n-1] + ... + dₙₛ u[n-s]

linear time-invariant causal system:
⎧ x[n+1] = a₀ x[n] + a₁ x[n-1] + ... + aₚ x[n-p]
⎨        + b₀ u[n] + b₁ u[n-1] + ... + b₉ u[n-q]
⎪ y[n]   = c₀ x[n] + c₁ x[n-1] + ... + cᵣ x[n-r]
⎩        + d₀ u[n] + d₁ u[n-1] + ... + dₛ u[n-s]

LTI system

system                    LTI system
with internal state:      with internal state:

  ╭╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╮       ╭╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╮
  ╎     ┌───┐     ╎       ╎             ┌───┐             ╎
  ╎ ┌──→│ f │     ╎       ╎           ┌─│ a │─┐           ╎
  ╎ │ ┌→└───┘─┐   ╎       ╎           │ └───┘ │           ╎
  ╎ │ ├── x ←─┘   ╎       ╎    ┌───┐  ↓+      │  ┌───┐    ╎
  ╰╌│╌│╌╌╌╌╌╌╌╌─╌╌╯       ╎ ┌─→│ b │─→⊕─→ x ──┴─→│ c │──┐ ╎
    │ └→┌───┐             ╎ │  └───┘ +           └───┘  │ ╎
u ──┴──→│ g │─────→ y     ╰╌│╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌╌│╌╯
        └───┘               │           ┌───┐           ↓+
                        u ──┴──────────→│ d │──────────→⊕─→ y
                                        └───┘          +

系統模型畫成方塊圖，形成新觀點。
一、內部狀態x，自迴歸a。不斷改變自身，不受外界影響。
二、接著讓x受到外界影響。
　　引入輸入訊號。進來之前先過b，出去之前先過c。
三、線性非時變系統，只有三種運算，位移、倍率、加法。
　　因此可以直接把兩道輸入訊號加在一起。
　　因此可以直接將輸入訊號改成輸出訊號。
　　出去之前先過c再過d，兩者複合當作新c。
　　阿就前饋看起來比較工整啊。

維基百科的方塊圖，畫成其他造型。看你喜歡哪種都行。
https://commons.wikimedia.org/wiki/File:Typical_State_Space_model.svg

system representation

polynomial representation / matrix representation

polynomial representation:
⎧ x[n+1] = a₀ x[n] + a₁ x[n-1] + ... + aₚ x[n-p]
⎨        + b₀ u[n] + b₁ u[n-1] + ... + b₉ u[n-q]
⎪ y[n]   = c₀ x[n] + c₁ x[n-1] + ... + cᵣ x[n-r]
⎩        + d₀ u[n] + d₁ u[n-1] + ... + dₛ u[n-s]

matrix representation:
⎰ x⃗[n+1] = A x⃗[n] + B u⃗[n]
⎱ y⃗[n]   = C x⃗[n] + D u⃗[n]

四串系統參數可以改寫成四個常對角矩陣ABCD。

state-space representation / phase-variable representation

state-space model也有這兩種表示法。不再贅述。

companion matrix realization

模仿ARMA system，湊出矩陣ABCD。

https://people.duke.edu/~hpgavin/SystemID/CourseNotes/LTI.pdf

expansion

專著《Filtering and System Identification: A Least Squares Approach》。

reachability matrix / observability matrix

可到達性矩陣：x[n+1]遞迴函數，改寫成矩陣。
可觀察性矩陣：y[n]遞迴函數，改寫成矩陣。

訊號從數值u[n] x[n] y[n]推廣為向量u⃗[n] x⃗[n] y⃗[n]。
簡單起見，這個小節的u x y不添加向量符號、也不使用粗體字。

k是向量長度。
rank(A) = k，那麼A的直條的線性組合，得以形成每一種向量。
注意到，數列索引值n／數列長度n、向量長度k，兩者意義不同。

expansion of x (polynomial representation):

⎧ x[1] = A x[0] + B u[0]
⎨ x[2] = A² x[0] + AB u[0] + B u[1]
⎪   :
⎩ x[k] = Aᵏ x[0] + Aᵏ⁻¹B u[0] + ⋯ + AB u[k-2] + B u[k-1]

formula of x (matrix representation):

                                      ⎡ u[0]   ⎤
x[k] = Aᵏ x[0] + [ Aᵏ⁻¹B ⋯ A²B AB B ] ⎢   :    ⎥
                                      ⎣ u[k-1] ⎦
                           Ɔ

reachability matrix (controllability matrix):

𝐶 = [ A⁰B ⋯ Aᵏ⁻¹B ]

reachability:

    any y can be generated by x[0] and u
<=> any x[k] can be generated by x[0] and u[0⋯k-1]
<=> any (x[k] - Aᵏ x[0]) can be generated by u[0⋯k-1]
<=> rank(𝐶) = k

expansion of y (polynomial representation):

⎧ y[0] = C x[0] + D u[0]
⎨ y[1] = CA x[0] + CB u[0] + D u[1]
⎪   :
⎩ y[k-1] = CAᵏ⁻¹ x[0] + CAᵏ⁻²B u[0] + ... + CAB u[k-1] + D u[k-1]

expansion of y (matrix representation):

⎡ y[0]   ⎤   ⎡ C     ⎤        ⎡ D    0   ⋯ ⎤ ⎡ u[0]   ⎤
⎢   :    ⎥ = ⎢ CA    ⎥ x[0] + ⎢ CB   D   ⋯ ⎥ ⎢   :    ⎥
⎢   :    ⎥   ⎢  :    ⎥        ⎢ CAB  CB  ⋯ ⎥ ⎢   :    ⎥
⎣ y[k-1] ⎦   ⎣ CAᵏ⁻¹ ⎦        ⎣ :    :     ⎦ ⎣ u[k-1] ⎦
                 𝑂                   𝐺

observability matrix:

    ⎡ CA⁰   ⎤
𝑂 = ⎢  :    ⎥
    ⎣ CAᵏ⁻¹ ⎦

observability:

    unique x can be reconstructed by y and u
<=> unique x[0] can be reconstructed by y[0⋯k-1] and u[0⋯k-1]
<=> unique x[0] can be reconstructed by (y[0⋯k-1] - 𝐺 u[0⋯k-1])
<=> rank(𝑂) = k

Krylov subspace：湊足k種連續次方，得以構成k維空間。
一、可到達性矩陣、可觀察性矩陣至少需要列出0次方到k-1次方。
二、訊號長度至少是k個數字，得以檢查可到達性、可觀察性。

stability

專著《Optimal State Estimation: Kalman, H_∞ and Nonlinear Approaches》。

Hurwitz matrix

如果A穩定，那麼整個系統穩定。
https://en.wikipedia.org/wiki/Hurwitz-stable_matrix

stabilizability

output controllablility

output controllable if and only if the output controllability
Gramian is positive definite.

https://arxiv.org/pdf/2306.08523
http://maecourses.ucsd.edu/~mdeolive/mae280a/lecture12.pdf

dual system

專著《Linear systems》。

rank / nullity

矩陣邊長是n，而實際維度不足n。
套用此矩陣、實施變換，有些維度完全消失。
例如投影矩陣，有些維度完整保留，有些維度完全消失。
消失的維度，特徵值是零，特徵向量未定義。
（數學家規定必須湊滿n個特徵值。但是不必湊滿n個特徵向量。）

stability / instability

繼續深入。
實際維度是k，其中包含了縮小、不變、放大。
矩陣邊長是n，其中包含了消失、縮小、不變、放大。
特徵值，其中包含了等於0，絕對值小於1、絕對值等於1、絕對值大於1。
不斷套用此矩陣、不斷實施變換，
有些維度完全消失、趨近消失、保持不變、變成正負無限大。

消失、縮小、不變，最終形成穩態。
放大，最終不會形成穩態。

reachability / controllability / stabilizability

可到達性：總是存在控制訊號，各種初始狀態都可以得到各種狀態。
可控制性：總是存在控制訊號，可以得到各種狀態。
可穩定性：總是存在控制訊號，讓系統穩定。

可到達性：每種狀態皆可抵達。
　　　　　所有初始條件，皆可抵達特定一種狀態。
可控制性：每種狀態皆可抵達。
　　　　　存在初始條件，可以抵達特定一種狀態。
可穩定性：每種狀態不見得皆可抵達。
　　　　　所有初始條件，輸出訊號皆可抵達穩態。

可到達性：已知初始狀態x[0]、當前狀態x[n]，
　　　　　總是存在控制訊號u，
　　　　　可以從初始狀態x[0]得到當前狀態x[n]。
　　　　（可以得到狀態數列x、輸出數列y。畢竟已知系統公式。）
可控制性：已知當前狀態x[n]。
　　　　　總是存在控制訊號u、初始狀態x[0]，
　　　　　可以得到當前狀態x[n]。
　　　　（可以得到狀態數列x。畢竟已知系統公式。）
可穩定性：已知初始狀態x[0]，
　　　　　總是存在控制訊號u，
　　　　　使得每種輸出訊號y到達穩態。
　　　　（即便不滿足可控制性。例如ABCD矩陣全是零。）

上述定義經過推導得到下述結論。
可到達性：初始狀態x[0] = 0可以到達所有狀態。
可控制性：存在初始狀態x[0]可以到達狀態x[n] = 0。

連續時間系統，可到達性<=>可控制性。
離散時間系統，可到達性=>可控制性。
離散時間系統，A有反矩陣，可到達性<=>可控制性。

observability / constructability / detectability

可觀察性：已知控制訊號u、輸出訊號y，
　　　　　總是存在唯一一種初始狀態x[0]。
　　　　　可以用控制訊號u得到輸出訊號y。
　　　　（可以得到唯一一種狀態數列x。畢竟已知系統公式。）
可建設性：已知控制訊號u、輸出訊號y，
　　　　　總是存在一種初始狀態x[0]。
　　　　　可以用控制訊號u得到輸出訊號y。
　　　　（可以得到一種狀態數列x。畢竟已知系統公式。）
可偵測性：已知控制訊號u、輸出訊號y，
　　　　　可以讓輸出訊號y到達穩態。
　　　　（即便不滿足可觀察性。例如不知道狀態但知道會穩定。）

連續時間系統，可觀察性<=>可建設性。
離散時間系統，可觀察性=>可建設性。
離散時間系統，C有反矩陣，可觀察性<=>可建設性。

dual system

四個矩陣ABCD通通轉置，仍然是LTI state-space model。
四個矩陣ABCD通通轉置，reachability與observability恰好對調。
原理是線性代數基本定理：kernel(A) = image(Aᵀ)⟂
就這樣。

┌──────────────────┐ dual ┌──────────────────┐
│   reachability   │←────→│  observability   │
└──────────────────┘      └──────────────────┘
         ⇓                          ⇓
┌──────────────────┐ dual ┌──────────────────┐
│ controllability  │←────→│ constructability │
└──────────────────┘      └──────────────────┘

controllable canonical form / observable canonical form

以矩陣A為主角，
實施companion matrix realization，
得到controllable canonical form。

四個矩陣ABCD通通轉置，得到observable canonical form。
即是dual system。

http://www.dii.unimo.it/~zanasi/didattica/Teoria_dei_Sistemi/Luc_TDS_ING_2016_Reachability_and_Controllability.pdf
http://www.dii.unimo.it/~zanasi/didattica/Teoria_dei_Sistemi/Luc_TDS_ING_2016_Observability_and_Constructability.pdf

minimal realization

similarity transformation

A' = T⁻¹AT
B' = T⁻¹B
C' = CT
D' = D
x' = T⁻¹x
u' = u
y' = y

minimal realization

minimal realization is both controllable and observable.

eigensystem realization

eigensystem realization (Ho–Kálmán realization)

教科書將卷積核從f改成g。但是我想要用f。【尚待確認】

1. get (g[1], ⋯, g[n])
   by impulse response

2. get 𝑂,𝐶
   by compact SVD of Hankel((g[1], ⋯, g[n]))

3. get A,B,C,D
   A: extract from Hankel((g[2], ⋯, g[k+1]))
   B: first column of 𝐶
   C: first row of 𝑂
   D: orignal D (by similarity transformation)

使用impulse response。
輸入是脈衝函數，輸出稱作脈衝響應。
LTI系統當中，脈衝響應恰是卷積核。
最後，卷積核進一步改寫成矩陣。

https://math.stackexchange.com/questions/3275985/
https://par.nsf.gov/servlets/purl/10312038

Markov parameters

卷積核（脈衝響應）改寫成矩陣。然後省略g[0]。

Markov parameters:

(CA⁰B, CA¹B, ...)

theorem:

Hankel((g[1], ⋯, g[k])) = Hankel((CA⁰B, ⋯, CAᵏ⁻¹B)) = 𝑂𝐶

⎡ g[1] g[2]   ⋯ g[k]    ⎤   ⎡ CA⁰B   CA¹B ⋯ CAᵏ⁻¹B  ⎤
⎢ g[2] g[3]   ⋯ g[k+1]  ⎥   ⎢ CA¹B   CA²B ⋯ CAᵏB    ⎥
⎢   :    :        :     ⎥ = ⎢  :      :      :      ⎥ = 𝑂𝐶
⎢   :    :        :     ⎥   ⎢  :      :      :      ⎥
⎣ g[k] g[k+1] ⋯ g[2k-1] ⎦   ⎣ CAᵏ⁻¹B CAᵏB ⋯ CA²ᵏ⁻²B ⎦

Toeplitz matrix = constant diagonal matrix
Hankel matrix = constant skew-diagonal matrix

可到達性矩陣故意被左右顛倒，就是為了產生常反對角矩陣。
如果我沒搞錯，這與卷積互相呼應。
卷積的其中一個數列也是需要左右顛倒。

Hankel singular values

Hankel((g[1], ⋯, g[n]))
做共軛分解conjugate decomposition，得到𝑂與𝐶。
從中間剖開一人分一半。

共軛分解採用compact SVD，順便降維，以便形成minimal realization。
去掉奇異值是零的多餘維度。
方便起見，新維度還是標記成k，新矩陣則追加下標k。
A = UΣVᵀ = U√Σ√ΣVᵀ  compact singular value decomposition
Aₖ = Uₖ√Σₖ√ΣₖVₖᵀ    minimal realization (rank k)
𝑂 = Uₖ√Σₖ
𝐶 = √ΣₖVₖᵀ

A = BᵀB。
對稱半正定矩陣A，分解成矩陣內積。
沒有正式學術名稱。
少數文獻稱作共軛分解conjugate decomposition。
概念宛如將一個平方值，分解成共軛複數相乘。
分解方式有無限多種。其中有兩種知名方式：
一、Cholesky decomposition，分解成上三角矩陣。
    A = LLᵀ and B = Lᵀ
二、eigendecomposition。特徵分解。
    A = EΛE⁻¹ = EΛEᵀ = E√Λ√ΛEᵀ and B = √ΛEᵀ
對稱半正定矩陣的情況下，特徵分解EVD等同奇異值分解SVD！
不是對稱半正定矩陣的情況下，只好改用SVD。
Hankel((g[1], ⋯, g[k]))是對稱矩陣，但是通常不是半正定矩陣。

system matrices

算A：移位一個時刻，湊出A。然後𝑂和𝐶移項即得。
⎡ g[2]   g[3]   ⋯ g[k+1] ⎤   ⎡ CA¹B CA²B   ⋯ CAᵏB    ⎤
⎢ g[3]   g[4]   ⋯ g[k+2] ⎥   ⎢ CA²B CA²B   ⋯ CAᵏ⁺¹B  ⎥
⎢   :      :        :    ⎥ = ⎢  :      :      :      ⎥ = 𝑂A𝐶
⎢   :      :        :    ⎥   ⎢  :      :      :      ⎥
⎣ g[k+1] g[k+2] ⋯ g[2k]  ⎦   ⎣ CAᵏB CAᵏ⁺¹B ⋯ CA²ᵏ⁻¹B ⎦

let Hₖ = Hankel((g[1], ⋯, g[k]))
Aₖ = 𝑂ₖ⁺(Hₖ+⃡1)𝐶ₖ⁺

算B：𝐶ₖ第零個直條就是Bₖ。
算C：𝑂ₖ第零個橫條就是Cₖ。
算D：原本的D降維。【尚待確認】

balanced realization

balanced realization

四個矩陣ABCD的方程組，對角化比較複雜。
四個矩陣不能一起對角化。
只能以𝐶為主、或者以𝑂為主。

有人發現取平方（矩陣內積和外積）再開根號，𝐶和𝑂就可以一起對角化。
宛如SVD的背後原理。SVD即是EVD取平方（矩陣內積和外積）再開根號。

1. get P,Q by Lyapunov equation
   AP + PAᵀ + BBᵀ = 0
   AᵀQ + QA + CᵀC = 0
2. get 𝑂 by Cholesky decomposion
   Q = 𝑂ᵀ𝑂
3. get U,Σ by eigendecomposition
   𝑂P𝑂ᵀ = 𝑂𝐶𝐶ᵀ𝑂ᵀ= HHᵀ = UΣ²Uᵀ
4. get T by conjugate decomposition 
   T   = 𝑂⁻¹U√Σ
   T⁻¹ = √ΣUᵀ𝑂

其實步驟1和2是多餘的。
H = Hankel((g[1], ⋯, g[k]))
使用矩陣內積HᵀH、矩陣外積HHᵀ，即可完成平衡實現。
步驟1和2主要是為了介紹Lyapunov equation。

https://ocw.mit.edu/courses/6-241j-dynamic-systems-and-control-spring-2011/7f6754029d6945bd3bdc745080890bf8_MIT6_241JS11_chap26.pdf

image

線性代數當中，
image：一個矩陣，嘗試各種輸入向量，所能得到的各種輸出向量。輸出向量的集合。

theorem:
image(AAᵀ) = image(A)

直觀解讀：
一、Aᵀx是輸入向量x實施變換，獲得一部分的輸入向量。比原本少。
　　再經過A變換，獲得一部分的輸出向量。比原本少。
　　因此image(AAᵀ) ⊆ image(A)。
二、Aᵀ是垂直投影。Aᵀ消滅的維度，A消滅的維度，基本相等。
　　換句話說，rank(Aᵀ) = rank(A)。
　　即便額外套用Aᵀ，剩下的維度還是一樣多。
　　因此image(AAᵀ) = image(A)。

嚴謹證明，分成三個階段：
(1) AᵀAx = 0 <=> Ax = 0
(2) kernel(AᵀA) = kernel(A)
(3) image(AᵀA) = image(Aᵀ)

(1)(⟹) 等號兩邊同乘xᵀ
   AᵀAx = 0
=> xᵀAᵀAx = 0
=> ‖Ax‖² = 0
=> Ax = 0

(1)(⟸) 等號兩邊同乘Aᵀ。
   Ax = 0
=> AᵀAx = 0

(2) 因為兩個方程式等價，所以兩個解集合也等價。

(3) 線性代數基本定理：kernel(A) = image(Aᵀ)⟂
   kernel(AᵀA) = kernel(A)
=> image((AᵀA)ᵀ)⟂ = image(Aᵀ)⟂
=> image(AᵀA)⟂ = image(Aᵀ)⟂
=> image(AᵀA) = image(Aᵀ)

https://math.stackexchange.com/questions/2411508

reachability

圖論當中，
reachable：從一個起點可以到達一個終點。
connected：從各種起點可以到達各種終點。

訊號學當中，
reachable set：從各種起點所能到達的各種終點。終點的集合。

使用矩陣來描述相鄰關係，
訊號學的reachable set＝線性代數的image。

先備知識請見本站文件「adjacency matrix」。
圖論當中，A是圖。Aᵀ是反向圖（顛倒所有邊）。
套用Aᵀ：順走1步。
套用A：逆走1步。
套用AAᵀ：先順走1步、再逆走1步（矩陣乘法由右往左讀）。

考慮所有可能的起點，以及所有可以到達的終點。
套用A、套用AAᵀ，所有可以到達的終點一樣多。

套用R = (A⁰+A¹+⋯+Aᵏ)：
逆走0步到k步，所到之處通通聯集。
套用RRᵀ = (A⁰+A¹+⋯+Aᵏ)(A⁰+A¹+⋯+Aᵏ)ᵀ：
先順走再逆走（矩陣乘法由右往左讀）。
展開式子，變成兩兩相乘，得到各種情況。例如順1逆2。
套用P = A⁰(A⁰)ᵀ+A¹(A¹)ᵀ+⋯+Aᵏ(Aᵏ)ᵀ：
順0逆0、順1逆1、……、順k逆k，所到之處通通聯集。

考慮所有可能的起點，以及所有可以到達的終點。
套用R、套用P，所有可以到達的終點一樣多。

reachability Gramian

Gram matrix            MᵀM

k-step                 Aᵏ⁻¹B
outer product          (Aᵏ⁻¹B)(Aᵏ⁻¹B)ᵀ = Aᵏ⁻¹BBᵀ(Aᵏ⁻¹)ᵀ
image                  image(Aᵏ⁻¹B) = image((Aᵏ⁻¹B)(Aᵏ⁻¹B)ᵀ)

k-reachability matrix         𝐶ₖ = [ A⁰B ⋯ Aᵏ⁻¹B ]
k-reachability Gramian matrix [ A⁰BBᵀ(A⁰)ᵀ ⋯ Aᵏ⁻¹BBᵀ(Aᵏ⁻¹)ᵀ ]

k-reachability         R[k] = sum Aᵏ⁻¹B
                             k=0⋯∞
k-reachability Gramian P[k] = sum Aᵏ⁻¹BBᵀ(Aᵏ⁻¹)ᵀ = sum 𝐶ₖ 𝐶ₖᵀ
                             k=0⋯∞                k=0⋯∞
k-reachable set        image(R[k]) = image(P[k])

https://control.asu.edu/Classes/MAE507/507Lecture06.pdf
https://control.asu.edu/Classes/MAE507/507Lecture17.pdf

Lyapunov equation

k-reachability Gramian P[k+1] = A P[k] Aᵀ + BBᵀ
Lyapunov stability     image(P[k+1]) = image(P[k]) when k→∞
Lyapunov equation      X = AXAᵀ + Q
                       where X = P[k] and k→∞
                             Q = BBᵀ
Lyapunov theorem       Lyapunov stability <=> Lyapunov equation

Hankel singular value

reachability Gramian   P = 𝐶𝐶ᵀ  =>  AP + PAᵀ + BBᵀ = 0
observability Gramian  Q = 𝑂ᵀ𝑂  =>  AᵀQ + QA + CᵀC = 0
Hankel singular value  Σ = sqrt(eigenvalue(PQ))
balanced realization   P' = Q' = Σ

similarity transformation

P' = T⁻¹P(Tᵀ)⁻¹
Q' = TᵀQT

表格

專著《Linear Optimal Control: H₂ and H_∞ Methods》。

                   LTI state-space model
╭────────────────────────────┬────────────────────────────╮
│ polynomial representation  │ matrix representation      │
╞════════════════════════════╧════════════════════════════╡
│ time domain                                             │
╞════════════════════════════╤════════════════════════════╡
│ sequence                   │ companion matrix           │
│ a,b,c,d                    │ A,B,C,D                    │
│ a = (a₀, a₁, ...)          │                            │
│ b = (b₀, b₁, ...)          │                            │
│ c = (c₀, c₁, ...)          │                            │
│ d = (d₀, d₁, ...)          │                            │
├────────────────────────────┤────────────────────────────┤
│ linear recurrence          │ linear recurrence          │
│ ⎧ x[n+1] = a₀ x[n] + ⋯     │ ⎰ x⃗[n+1] = A x⃗[n] + B u⃗[n] │
│ ⎨        + b₀ u[n] + ⋯     │ ⎱ y⃗[n]   = C x⃗[n] + D u⃗[n] │
│ ⎪ y[n]   = c₀ x[n] + ⋯     │                            │
│ ⎩        + d₀ u[n] + ⋯     │                            │
├────────────────────────────┤────────────────────────────┤
│ expansion                  │ expansion                  │
│ x[n+1] = ......            | x[n+1] = Aⁿ⁺¹ x[0] + Ɔ u⃗[n]|
│ y[n]   = ......            │ y⃗[n]   = 𝑂 x[0] + 𝐺 u⃗[n]   │
├────────────────────────────┤────────────────────────────┤
│ convolution                │ Hankel matrix              │
│ ⎰ x+⃡1 = a ∗ x + b ∗ u      │ 𝑂𝐶                         │
│ ⎱ y   = c ∗ x + d ∗ u      │                            │
├────────────────────────────┤────────────────────────────┤
│ impulse response           │ Markov parameters          │
│ g[n] = ......              │ g[n] = CAⁿ⁻¹B + D δ[n]     │
│                            │      = ⎰ D      , if n = 0 │
│                            │        ⎱ CAⁿ⁻¹B , if n > 0 │
├────────────────────────────┤────────────────────────────┤
│ expansion                  │ expansion                  │
│ y[n] = ......              │ y[n] = C Aⁿ x[0]           │
│                            │      +  sum  CAⁿ⁻¹⁻ⁱB u[i] │
│                            │       i=0⋯n-1              │
│                            │      + D u[0]              │
╞════════════════════════════╧════════════════════════════╡
│ z-domain                                                │
╞════════════════════════════╤════════════════════════════╡
│ generating function        │ resolvent                  │
│ A(z) = a₀z⁰ + a₁z⁻¹ + ⋯    │ G(z) = C(zI-A)⁻¹B + D      │
│ B(z) = b₀z⁰ + b₁z⁻¹ + ⋯    │          adj(zI-A)         │
│ C(z) = c₀z⁰ + c₁z⁻¹ + ⋯    │      = C ————————— B + D   │
│ D(z) = d₀z⁰ + d₁z⁻¹ + ⋯    │          det(zI-A)         │
├────────────────────────────┤────────────────────────────┤
│ BIBO stability             │ BIBO stability             │
│ A(z) = 0                   │ det(zI-A) = 0              │
╰────────────────────────────┴────────────────────────────╯

system identification

專著《Numerical Methods for Linear Control Systems》。

  min   sum ║⎡x[n+1]⎤ _ ⎡ A B ⎤ ⎡x[n]⎤║²
A,B,C,D  n  ║⎣y[n]  ⎦   ⎣ C D ⎦ ⎣u[n]⎦║ꜰ

linear least squares有三種解法：
1. normal equation
2. QR decomposition
3. singular value decomposition

大家習慣使用QR分解，時間複雜度最低。

訊號學家自創一個詞彙subspace identification method。
其實就是這三種解法。

知名演算法N4SID。
N4SID將訊號切成三段：過去／現在／未來。
我不知道有何效果。請自行參考專著。

system functionality【尚無正式名稱】

verb     | action noun | agent noun
---------|-------------|-----------
observe  | observation | observer
control  | control     | controller

觀測器：系統輸出之後，追加一個系統，用來觀測系統內部狀態。
　　　　輸入訊號需要前饋。
控制器：系統輸入之前，追加一個系統，用來調整輸入訊號暨輸出訊號。
　　　　輸出訊號需要回饋。

       ┌──────────┐
u ──┬─→│  system  │──┬──────────────────→ y
    │  └──────────┘  └─→┌──────────┐
    ╰──────────────────→│ observer │────→ x̂
                        └──────────┘

       ┌──────────┐  u  ┌──────────┐
r ────→│controller│────→│  system  │──┬─→ y
    ╭─→└──────────┘     └──────────┘  │
    ╰─────────────────────────────────╯

linear quadratic observer（Kalman filter）

linear least squares

專著《Optimal Estimation of Dynamic Systems》。

linear least squares              線性系統（公式解）
weighted linear least squares     追加權重（能量範數）
constrained linear least squares  追加約束條件（拉格朗日乘數法）
nonlinear least squares           非線性系統（最佳化演算法）

https://books.google.com.tw/books?id=ITKKkBFxgNsC&pg=PA53

linear least squares:

given y = Ax
find x = argmin ‖y - Ax‖²
where A is overdetermined system (k ≥ n)

⎡ y[0]   ⎤   ⎡ A[0][0]   ... A[0][n-1]   ⎤ ⎡ x[0]   ⎤
⎢   :    ⎥   ⎢     :             :       ⎥ ⎢   :    ⎥
⎢   :    ⎥ = ⎢     :             :       ⎥ ⎢   :    ⎥
⎢   :    ⎥   ⎢     :             :       ⎥ ⎢   :    ⎥
⎣ y[k-1] ⎦   ⎣ A[k-1][0] ... A[k-1][n-1] ⎦ ⎣ x[n-1] ⎦
    y                     A                    x

error:

e = y - Ax

objective function:

J = ‖y - Ax‖²
  = e₀² + e₁² + ... + eₙ₋₁²
  = eᵀe
  = (y - Ax)ᵀ(y - Ax)
  = yᵀy - xᵀAᵀy - yᵀAx + xᵀAᵀAx

solution:

dJ/dx = - 2Aᵀy + 2AᵀAx = 0   極值位於一次微分等於零的地方
x = (AᵀA)⁻¹Aᵀy               normal equation
K = (AᵀA)⁻¹Aᵀ                pseudoinverse
x is least-squares solution if A has full column rank.

weighted linear least squares:

given y = Ax
find x = argmin ‖y - Ax‖²

objective function:

J = w₀e₀² + w₁e₁² + ... + wₙ₋₁eₙ₋₁²
  = eᵀWe
  = (y - Ax)ᵀW(y - Ax)

weight:

    ⎡ w₀        ⎤
W = ⎢    ⋱      ⎥
    ⎣      wₙ₋₁ ⎦

solution:

dJ/dx = - 2AᵀWy + 2AᵀWAx = 0
x = (AᵀWA)⁻¹AᵀWy
K = (AᵀWA)⁻¹AᵀW

linear least squares (with noise):

given y = Ax + v
find x = argmin ‖y - Ax - v‖²
where v is zero-mean white Gaussian noise

theorem:

normal regression = linear regression

solution (the same):

x = (AᵀA)⁻¹Aᵀy

weighted linear least squares (with noise):

given y = Ax + v
find x = argmin ‖y - Ax‖²
where E[vᵢ²] = σᵢ² and E[vᵢvⱼ] = 0

assumptions:

1. noises are zero-mean: correlation = covariance
2. noises has its variance
3. noises are independent: independent => uncorrelated

error covariance matrix:

             ⎡ σ₀²         ⎤
R = E[vvᵀ] = ⎢     ⋱       ⎥
             ⎣       σₙ₋₁² ⎦

objective function:

J = e₀²/σ₀² + e₁²/σ₁² + ... + eₙ₋₁²/σₙ₋₁²
  = eᵀR⁻¹e
  = (y - Ax)ᵀR⁻¹(y - Ax)

誤差變異數通通調整成一樣，
以便讓normal regression = linear regression。

solution:

dJ/dx = - 2AᵀR⁻¹y + 2AᵀR⁻¹Ax = 0
x = (AᵀR⁻¹A)⁻¹AᵀR⁻¹y

sequential linear least squares (recursive least squares)

專著《Adaptive Filter Theory》。

原始名稱recursive least squares。
雖然是online algorithm、iterative method，
但是原始作者取名recursive。
畢竟當時online這個詞彙尚未流行。
畢竟數值分析領域所有演算法都是iterative method，只好換個詞彙。
另外，數值分析領域習慣將online algorithm冠上前綴sequential。

linear least squares            -> offline algorithm
sequential linear least squares -> online algorithm

凡是線性迴歸，都可以改成線上演算法。舉例來說：
AR system求得系統參數，改成線上演算法，稱作Wiener filter。
state-space model求得內部狀態，改成線上演算法，稱作Kalman filter。

一、時間有限。改成online algorithm，即時計算當前最佳解。
　　目標函數是以前到當前所有誤差總和。
二、記憶體有限。隨時記錄最新k個訊號、形成最新k×n矩陣。
　　理論上，向量與矩陣越長，答案越精確。
　　方便起見，每回合向量x、矩陣A，保持相同長度k、而且盡量長。

sequential linear least squares:

yₙ = Aₙ xₙ

objective function:

Jₙ = sum ‖yᵢ - Aᵢ xᵢ‖²
    i=0⋯n

initial value:

x₀ = (A₀ᵀA₀)⁻¹A₀ᵀy₀
P₀ = (A₀ᵀA₀)⁻¹

iteration:

Kₙ = Pₙ₋₁ Aₙᵀ (I + AₙPₙ₋₁Aₙᵀ)⁻¹
xₙ = xₙ₋₁ + Kₙ (yₙ - Aₙxₙ₋₁)
Pₙ = (I - KₙAₙ) Pₙ₋₁

K: pseudoinverse of A              K = (AᵀA)⁻¹Aᵀ
P: inverse of Gram matrix of A     P = (AᵀA)⁻¹

改成online algorithm，必須找到遞迴公式。

xₙ的遞迴公式：xₙ = xₙ₋₁ + Kₙ (yₙ - Aₙ xₙ₋₁)。
利用上回合的輸入xₙ₋₁，
求得這回合的輸入xₙ。
因為我們只知輸出yₙ，所以利用輸出反推輸入。
Kₙ是Aₙ的虛擬反矩陣。

解釋一下遞迴公式的意義。
觀察(yₙ - Aₙ xₙ₋₁)。
  yₙ - Aₙ xₙ₋₁
= Aₙ xₙ - Aₙ xₙ₋₁
= Aₙ (xₙ - xₙ₋₁)
輸入變化xₙ - xₙ₋₁。
輸出變化Aₙ (xₙ - xₙ₋₁)。
利用輸出變化反推輸入變化。

Kₙ的遞迴公式：Pₙ₋₁ Aₙᵀ (I + AₙPₙ₋₁Aₙᵀ)⁻¹。
使用Sherman–Morrison–Woodbury identity求得。

sequential linear least squares:

yₙ = Aₙ xₙ

objective function:

Jₙ = sum ‖yᵢ - Aᵢ xᵢ‖²
    i=0⋯n

solution:

     ⎛              ⎞⁻¹⎛              ⎞
xₙ = ⎜  sum {AᵢᵀAᵢ} ⎟  ⎜  sum {Aᵢᵀyᵢ} ⎟
     ⎝ i=0⋯n        ⎠  ⎝ i=0⋯n        ⎠

Rₙ ≜  sum {AᵢᵀAᵢ}
     i=0⋯n
bₙ ≜  sum {Aᵢᵀyᵢ}
     i=0⋯n

xₙ = Rₙ⁻¹bₙ

Rₙ = Rₙ₋₁ + AₙᵀAₙ
bₙ = bₙ₋₁ + Aₙᵀyₙ

inverse of Gram matrix:

Pₙ ≜ Rₙ⁻¹

xₙ = Pₙbₙ

Pₙ⁻¹ = Pₙ₋₁⁻¹ + AₙᵀAₙ
bₙ = bₙ₋₁ + Aₙᵀyₙ

iteration of xₙ:

xₙ = Pₙbₙ
   = Pₙ (bₙ₋₁ + Aₙᵀyₙ)
   = Pₙ (Pₙ₋₁⁻¹xₙ₋₁ + Aₙᵀyₙ)
   = Pₙ ((Pₙ⁻¹ - AₙᵀAₙ) xₙ₋₁ + Aₙᵀyₙ)
   = xₙ₋₁ + PₙAₙᵀyₙ - PₙAₙᵀAₙxₙ₋₁
   = xₙ₋₁ + Kₙ (yₙ - Aₙxₙ₋₁)

where Kₙ = PₙAₙᵀ is pseudoinverse of Aₙ

Sherman–Morrison–Woodbury identity:

(M + UᵀV)⁻¹ = M⁻¹ - M⁻¹Uᵀ (I + VM⁻¹Uᵀ)⁻¹ VM⁻¹

Pₙ = (Pₙ₋₁⁻¹ + AₙᵀAₙ)⁻¹
Pₙ = Pₙ₋₁ - Pₙ₋₁Aₙᵀ (I + AₙPₙ₋₁Aₙᵀ)⁻¹ AₙPₙ₋₁

where M = Pₙ₋₁ and U = V = Aₙ

anonymous identity:

I - (I+M)⁻¹M = (I+M)⁻¹

Kₙ = PₙAₙᵀ
   = Pₙ₋₁Aₙᵀ (I + AₙPₙ₋₁Aₙᵀ)⁻¹

where M = AₙPₙ₋₁Aₙᵀ

iteration of Pₙ:

Pₙ = Pₙ₋₁ - KₙAₙPₙ₋₁
   = (I - KₙAₙ) Pₙ₋₁

sequential linear least squares改造版本。

專著《Optimal State Estimation: Kalman, H_∞ and Nonlinear Approaches》。

(1) minimize the squared error of y
    https://mdav.ece.gatech.edu/ece-6250-fall2019/notes/21-notes-6250-f19.pdf
(2) minimize the squared error of x
    https://web.mit.edu/kirtley/kirtley/binlustuff/literature/control/Kalman%20filter.pdf

Wiener filter使用正常版本(1)。
Kalman filter使用改造版本(2)。

兩種版本不等價。
誤差和目標函數完全不同。
P和K是完全不同的東西，只是外觀相像。

推導過程出現無法處理的項Pₙ⁻。
大家敷衍了事。Pₙ⁻當作是Pₙ₋₁。

sequential linear least squares (with noise):

⎰ yₙ = Aₙ xₙ + vₙ                  linear least squares
⎱ x̂ₙ = x̂ₙ₋₁ + Kₙ (yₙ - Aₙ x̂ₙ₋₁)    sequential LLS update rule

objective function:

Jₙ = sum ‖xₙ - x̂ₙ‖²
    i=0⋯n

initial value:
 
x̂₀ = E[x₀] = (A₀ᵀA₀)⁻¹A₀ᵀy₀        estimation mean
P₀ = E[(x₀-x̂₀)(x₀-x̂₀)ᵀ] = O        error covariance matrix
                                   (assume two means are equal)
iteration:

Pₙ⁻ = Pₙ₋₁                         敷衍了事
Rₙ = E[vₙvₙᵀ]                      noise covariance matrix
                                   (assume the mean is zero)
Kₙ = Pₙ⁻ Aₙᵀ (Aₙ Pₙ⁻ Aₙᵀ + Rₙ)⁻¹
x̂ₙ = x̂ₙ⁻ + Kₙ (yₙ - Aₙ x̂ₙ⁻)
Pₙ = (I - Kₙ Aₙ) Pₙ⁻ (I - Kₙ Aₙ)ᵀ - Kₙ Rₙ Kₙᵀ

v: zero-mean Gaussian noise
R: covariance matrix of v
x: true x (from y = Ax)
x̂: estimate of x (from sequential LLS update rule)
e: difference of x and x̂
e⁻: difference of current x and previous x̂
P: covariance matrix of e
P⁻: covariance matrix of e⁻
K: Kalman gain

complementary filter:

x̂ₙ = x̂ₙ⁻ + Kₙ (yₙ - Aₙ x̂ₙ⁻)
   = x̂ₙ⁻ + Kₙ (Aₙ xₙ - Aₙ x̂ₙ⁻)
   = (I - Kₙ Aₙ) x̂ₙ⁻ + Kₙ Aₙ xₙ

input: previous estimation x̂ₙ⁻ and true xₙ
output: current estimation x̂ₙ

sequential linear least squares:

⎰ yₙ = Aₙ xₙ
⎱ x̂ₙ = x̂ₙ₋₁ + Kₙ (yₙ - Aₙ x̂ₙ₋₁)

error:

eₙ  ≜ xₙ - x̂ₙ
eₙ⁻ ≜ xₙ - x̂ₙ₋₁

eₙ = xₙ - x̂ₙ
   = xₙ - x̂ₙ₋₁ - Kₙ (yₙ    - Aₙ x̂ₙ₋₁)
   = xₙ - x̂ₙ₋₁ - Kₙ (Aₙ xₙ - Aₙ x̂ₙ₋₁)
   = eₙ⁻       - Kₙ  Aₙ eₙ⁻
   = (I - Kₙ Aₙ) eₙ⁻

objective function:

Jₙ = e₁² + e₂² + ... + eₙ²
   = eᵀe
   = tr(eeᵀ)
   = tr(Pₙ)

error covariance matrix:

Pₙ  ≜ eₙeₙᵀ   = (xₙ - x̂ₙ)(xₙ - x̂ₙ)ᵀ
Pₙ⁻ ≜ eₙ⁻eₙ⁻ᵀ = (xₙ - x̂ₙ₋₁)(xₙ - x̂ₙ₋₁)ᵀ

Pₙ = eₙeₙᵀ
   = ((I - Kₙ Aₙ) eₙ⁻) ((I - Kₙ Aₙ) eₙ⁻)ᵀ
   = (I - Kₙ Aₙ) Pₙ⁻ (I - Kₙ Aₙ)ᵀ

solution:

dJₙ/dKₙ = 2 (I - Kₙ Aₙ) Pₙ⁻ (-Aₙ)ᵀ = 0
Kₙ = Pₙ⁻ Aₙᵀ (Aₙ Pₙ⁻ Aₙᵀ)⁻¹

sequential linear least squares (with noise):

⎰ yₙ = Aₙ xₙ + vₙ
⎱ x̂ₙ = x̂ₙ₋₁ + Kₙ (yₙ - Aₙ x̂ₙ₋₁)

zero-mean Gaussian noise:

vₙ ~ N(0,σₙ²)
E[vₙ] = 0

error:

eₙ  ≜ xₙ - x̂ₙ
eₙ⁻ ≜ xₙ - x̂ₙ₋₁

E[eₙ] = E[xₙ - x̂ₙ]
      = E[xₙ - x̂ₙ₋₁ - Kₙ (yₙ - Aₙ x̂ₙ₋₁)]
      = E[xₙ - x̂ₙ₋₁ - Kₙ (Aₙ xₙ + vₙ - Aₙ x̂ₙ₋₁)]
      = E[eₙ⁻       - Kₙ  Aₙ eₙ⁻     - Kₙ vₙ   ]
      = E[(I - Kₙ Aₙ) eₙ⁻ - Kₙ vₙ]
      = (I - Kₙ Aₙ) E[eₙ⁻] - Kₙ E[vₙ]

objective function:

Jₙ = E[e₁² + e₂² + ... + eₙ²]
   = E[eᵀe]
   = E[tr(eeᵀ)]
   = tr(E[eeᵀ])
   = tr(Pₙ)

error covariance matrix:

Pₙ  ≜ E[eₙeₙᵀ]   = E[(xₙ - x̂ₙ)(xₙ - x̂ₙ)ᵀ]
Pₙ⁻ ≜ E[eₙ⁻eₙ⁻ᵀ] = E[(xₙ - x̂ₙ₋₁)(xₙ - x̂ₙ₋₁)ᵀ]

Pₙ = E[eₙeₙᵀ]
   = E[eₙ]E[eₙ]ᵀ
   = ((I - Kₙ Aₙ) E[eₙ⁻] - Kₙ E[vₙ]) (......)ᵀ
   = (I - Kₙ Aₙ) Pₙ⁻ (I - Kₙ Aₙ)ᵀ - Kₙ Rₙ Kₙᵀ
                                    ^^^^^^^^^
*noise covariance matrix:

Rₙ ≜ E[vₙvₙᵀ]

*assume vₙ and eₙ⁻ are independent:

E[vₙeₙ⁻] = E[vₙ]E[eₙ⁻] = 0

solution:

dJₙ/dKₙ = 2 (I - Kₙ Aₙ) Pₙ⁻ (-Aₙ)ᵀ - 2 Kₙ Rₙ = 0
Kₙ = Pₙ⁻ Aₙᵀ (Aₙ Pₙ⁻ Aₙᵀ + Rₙ)⁻¹

linear quadratic observer (Kalman filter)

原始名稱Kalman filter。
雖然是observer，但是原始作者取名filter。
畢竟當時observer這個詞彙尚未流行。

stochastic LTV state-space model:
⎰ x[n+1] = Aₙ x[n] + Bₙ u[n] + w[n]
⎱ y[n]   = Cₙ x[n] + Dₙ u[n] + v[n]

一、y[n] = Cₙ x[n]
　　套用sequential linear least squares。
二、v[n]
　　追加雜訊。公式解多出一項。
三、Dₙ u[n]
　　使用中學數學的配方法。
　　簡單起見，教科書習慣省略這項。
四、x[n+1] = Aₙ x[n]
　　事先套用線性變換Aₙ，調整一下估計平均數x̂、誤差共變異矩陣P。
五、Bₙ u[n] + w[n]
　　同上。公式解多出一項。

stochastic LTV state-space model (Dₙ = 0):
⎰ x[n+1] = Aₙ x[n] + Bₙ u[n] + w[n]
⎱ y[n]   = Cₙ x[n] + v[n]

initial value:
0. estimate:    x̂₀ = E[x₀]
                P₀ = E[(x₀-x̂₀)(x₀-x̂₀)ᵀ]

iteration:
1. transform:   x̂ₙ⁻ = Aₙ x̂ₙ₋₁ + Bₙ uₙ
                Pₙ⁻ = Aₙ Pₙ₋₁ Aₙᵀ + Qₙ
2. Kalman gain: Kₙ = Pₙ⁻ Cₙᵀ (Cₙ Pₙ⁻ Cₙᵀ + Rₙ)⁻¹
3. estimate:    x̂ₙ = x̂ₙ⁻ + Kₙ (yₙ - Cₙ x̂ₙ⁻)
                Pₙ = (I - Kₙ Cₙ) Pₙ⁻ (I - Kₙ Cₙ)ᵀ

x: state
y: measurement
A: state transition matrix
C: measurement matrix
w: state transition noise
v: measurement noise
Q: covariance matrix of w
R: covariance matrix of v
x̂: mean of estimate of x
P: covariance matrix of difference of x and x̂
x̂ₙ⁻: transform of x̂ₙ₋₁ under A
Pₙ⁻: transform of Pₙ₋₁ under A

affine transformation:

y = A x + v

assumptions:

E[x] = x
E[y] = y
E[v] = 0     (since v is zero-mean white noise)

cov[x] = 0
cov[y] = 0
cov[v] = r   (r is user-defined value)

mean:

E[y] = E[Ax + v]
     = A E[x] + E[v]

E[x] = (AᵀA)⁻¹Aᵀ(E[y] - E[v])     pseudoinverse 
     = (AᵀA)⁻¹Aᵀy                 by assumptions
     = A⁺y

covariance matrix:

cov[y] = cov[Ax + v]
       = A cov[x] Aᵀ + cov[v]

cov[x] = A⁺(cov[y] - cov[v])Aᵀ⁺   pseudoinverse【尚待確認】
       = -(A⁺ cov[v] Aᵀ⁺)         by assumptions

內部狀態初始值有兩種設定方式：
一、自行設定：
　　觀察現實狀況，判斷初始值。
　　例如慣性運動測量inertial motion measurement，
　　物體從靜止狀態開始運動，
　　位置是零、速度是零，
　　內部狀態初始值就是零向量。
二、解方程式：
　　取得足夠的輸出訊號、輸入訊號，
　　利用observability matrix那道等式，求得內部狀態初始值。
　　另外如果C恰是方陣而且可逆，那麼可以直接用C來算。

⎡ y[0]   ⎤   ⎡ C     ⎤        ⎡ D    0   ⋯ ⎤ ⎡ u[0]   ⎤
⎢   :    ⎥ = ⎢ CA    ⎥ x[0] + ⎢ CB   D   ⋯ ⎥ ⎢   :    ⎥
⎢   :    ⎥   ⎢  :    ⎥        ⎢ CAB  CB  ⋯ ⎥ ⎢   :    ⎥
⎣ y[k-1] ⎦   ⎣ CAᵏ⁻¹ ⎦        ⎣ :    :     ⎦ ⎣ u[k-1] ⎦
   y⃗             𝑂                   𝐺      u⃗

x̂₀ = (𝑂ᵀ𝑂)⁻¹𝑂ᵀ(y⃗ - 𝐺u⃗)   where rank(𝑂) = k

x̂₀ = C⁻¹(y[0] - Du[0])   if rank(C) = k and C is invertible

誤差共變異矩陣初始值有兩種設定方式：
一、如果方才採用自行設定：
　　就是零矩陣。畢竟誤差是0。
二、如果方才採用解方程式：
　　雜訊的共變異矩陣，套用變異數變換公式。
　　P₀ = (𝑂ᵀR⁻¹𝑂)⁻¹

nonlinear system

專著《Kalman Filter for Beginners: with MATLAB Examples》。

Kalman filter              線性系統
extended Kalman filter     非線性系統
unscented Kalman filter    平均數和變異數變換，透過取樣點。
particle Kalman filter     取樣點變換，透過遷移機率。
error-state Kalman filter  線性、非線性，兩者分開處理。

延伸閱讀：Riccati equation / Lyapunov equation

最後順便介紹兩種方程式。
一、Pₙ的遞迴公式，將Kₙ換掉，稱作Riccati equation。
二、Pₙ的遞迴公式，趨近穩態時，稱作Lyapunov equation。

Riccati equation

Riccati equation:

Pₙ = Aₙ Pₙ₋₁ Aₙᵀ + Qₙ
   - Aₙ Pₙ₋₁ Cₙᵀ (Cₙ Pₙ₋₁ Cₙᵀ + Rₙ)⁻¹ Cₙ Pₙ₋₁ Aₙᵀ

Lyapunov equation

                   | discrete-time    | continuous-time
                   | system           | system
--------------------------------------------------------
Lyapunov equation  | AXAᵀ - X + C = 0 | AX + XAᵀ + C = 0
Sylvester equation | AXB  - X + C = 0 | AX + XB  + C = 0

stochastic LTV state-space model (B = D = 0):

⎰ x[n+1] = Aₙ x[n] + w[n]
⎱ y[n]   = Cₙ x[n] + v[n]

mean:

x̄ₙ₊₁ = E[xₙ₊₁]
     = E[Aₙ xₙ + wₙ]
     = Aₙ E[xₙ] + E[wₙ]
     = Aₙ x̄ₙ

since E[wₙ] = 0

error covariance matrix:

Pₙ₊₁ = E[(xₙ₊₁ - x̄ₙ₊₁)(xₙ₊₁ - x̄ₙ₊₁)ᵀ]
     = E[(Aₙ xₙ + wₙ - Aₙ x̄ₙ)(Aₙ xₙ + wₙ - Aₙ x̄ₙ)ᵀ]
     = E[(Aₙ (xₙ - x̄ₙ) + wₙ)(Aₙ (xₙ - x̄ₙ) + wₙ)ᵀ]
     = ......
     = Aₙ Pₙ Aₙᵀ + E[wₙ wₙᵀ]

since E[xₙ wₙ] = 0

error covariance matrix at steady state:

P = Aₙ P Aₙᵀ + E[wₙ wₙᵀ]

since Pₙ₊₁ = Pₙ

Lyapunov equation:

P = Aₙ P Aₙᵀ + Qₙ        where Qₙ = E[wₙ wₙᵀ]
Aₙ P Aₙᵀ - P + Qₙ = 0

linear quadratic controller

專著《Linear Quadratic Control: An Introduction》。

optimal control

受控廠：LTI state-space model、LTV state-space model。
參考訊號：regulator、tracker。
觀測器：Kalman filter。
控制器：P controller、PI controller。
目標函數：quadratic function、convex function。
約束條件：自訂。
最佳化演算法：公式解、梯度下降法、牛頓法、線性規劃、動態規劃。

優點：為所欲為。想到什麼條件，都可以放進來。
缺點：難以控制收斂速度，甚至不會收斂。
　　　不定時得到最佳解，甚至不確定得到最佳解。

上個世紀，最佳控制曾經搞到噴射機失事。
大家研判這種作法行不通，導致最佳控制發展停滯。

現在教科書只會介紹其特例linear quadratic control：
LTV state-space model、regulator、P controller、quadratic form。
然後推導其公式解。
畢竟只有公式解沒有上述缺點。

linear quadratic control

linear是指線性系統。包括線性非時變系統、線性時變系統。
quadratic是指二次函數。

採用二次最佳化。
當矩陣是半正定矩陣，導致二次最佳化有最小值，可以推導公式解。
請見本站文件「quadratic optimization」。

極值位於一次微分等於零的地方。推導過程需要用到矩陣微分。
請見本站文件「differential calculus」。

一些變種。
linear quadratic control:          P controller
linear quadratic integral control: PI controller
linear quadratic Gaussian control: Gaussian noise

separation principle

separation principle:

min { f(x) + g(y) } = min f(x) + min g(y)
x,y                    x          y

如果兩個目標函數各自使用不同變數、那麼可以分別最佳化。
就這麼簡單。

separation principle:

  optimal controller of state-space model
= optimal observer + optimal controller using internal state

狀態空間模型的控制器，可以拆成兩部分：
一、最佳觀測器：估計內部狀態。
二、使用內部狀態的最佳控制器：利用內部狀態實施控制。

目標函數是輸出訊號y的函數。
利用y[n] = Cₙ x[n] + Dₙ u[n]強行展開數學式子，
形成內部狀態x、輸入訊號u的函數。
於是可以分別最佳化。
就這麼簡單。

control:

     ┌────────────┐  u  ┌─────────┐
r ──→│ controller │────→│  plant  │──┬─→ y
  ╭─→└────────────┘     └─────────┘  │
  ╰──────────────────────────────────╯

separation principle:

     ┌────────────┐  u  ┌─────────┐
r ──→│ controller │──┬─→│  plant  │──┬───────────────→ y
  ╭─→└────────────┘  │  └─────────┘  ╰→┌──────────┐
  │                  ╰────────────────→│ observer │──╮
  │                                    └──────────┘  │
  ╰──────────────────────────────────────────────────╯
                          x̂

new-style name:
1. linear quadratic observer
2. linear quadratic controller

old-fashioned name:
1. linear quadratic estimator (Kalman filter)
2. linear quadratic regulator/tracker (LQR/LQT)

linear quadratic controller

(1) r[n] = 0 , Dₙ = 0 , single objective
(2) r[n] = 0 , Dₙ = 0
(3) r[n] = 0
(4) r[n] is constant    (linear quadratic regulator)
(5) r[n] is arbitrary   (linear quadratic tracker)
(6) system with noise

依序介紹六種版本，從特例到通例。
每種版本都需要引進新的公式解求解技巧。
不過我不打算仔細講解求解技巧。
甚至我不確定公式解是什麼。
主要是因為我沒找到參考文獻。
這領域的書籍都寫得有點那啥。嘛，微妙。

r[n] = 0 , Dₙ = 0 , single objective

LTV state-space model:

⎰ x[n+1] = Aₙ x[n] + Bₙ u[n]
⎱ y[n]   = Cₙ x[n] + Dₙ u[n]

objective function (r[n] = 0 and Dₙ = 0):

Jₙ = sum ‖y[n] - r[n]‖²
   = sum ‖y[n]‖²             assume r[n] = 0 for all n
   = sum ‖Cₙ x[n]‖²          assume Dₙ = 0 for all n
   = sum x[n]ᵀ Cₙᵀ Cₙ x[n]
   = sum x[n]ᵀ Qₙ x[n]       rename variable
   = Jₙ₋₁ + x[n]ᵀ Qₙ x[n]

每個時刻互不相干。我們只需要分別最佳化每一項x[n]ᵀ Qₙ x[n]。
Qₙ = Cₙᵀ Cₙ形成對稱半正定矩陣，導致二次最佳化有最小值。
最小值位於一次微分等於零的地方。

linear quadratic controller:

dJₙ/dx[n] = 2 Qₙ x[n] = 0
x[n] = 0

答案就是零。缺乏討論意義。當作是暖身。

r[n] = 0 , Dₙ = 0

objective function (r[n] = 0 and Dₙ = 0):

Jₙ = sum { x[i]ᵀ Qᵢ x[i] + u[i]ᵀ Rᵢ u[i] }
    i=0⋯n

ansatz:

Jₙ = sum { x[i]ᵀ Pᵢ x[i] }
    i=0⋯n

u[n] = -Kₙ x[n]

solution (Riccati equation):

Pₙ = AₙᵀPₙ₊₁Aₙ - (AₙᵀPₙ₊₁Bₙ)(Rₙ+BₙᵀPₙ₊₁Bₙ)⁻¹(BₙᵀPₙ₊₁Aₙ) + Qₙ
Kₙ = (BₙᵀPₙ₊₁Bₙ+Rₙ)⁻¹(BₙᵀPₙ₊₁Aₙ)

infinite-horizon solution (algebraic Riccati equation):

當無窮時間n→∞且趨近穩態Pₛₛ = Pₙ = Pₙ₊₁ = ...
Pₛₛ = AₙᵀPₛₛAₙ - (AₙᵀPₛₛBₙ)(Rₙ+BₙᵀPₛₛBₙ)⁻¹(BₙᵀPₛₛAₙ) + Qₙ
Kₙ = (BₙᵀPₛₛBₙ+Rₙ)⁻¹(BₙᵀPₛₛAₙ)

linear quadratic controller:

Pₙ = AₙᵀPₙ₊₁Aₙ - (AₙᵀPₙ₊₁Bₙ)(Rₙ+BₙᵀPₙ₊₁Bₙ)⁻¹(BₙᵀPₙ₊₁Aₙ) + Qₙ
Kₙ = (BₙᵀPₙ₊₁Bₙ+Rₙ)⁻¹(BₙᵀPₙ₊₁Aₙ)
u[n] = -Kₙ x[n]
x[n+1]​ = Aₙ x[n]​ + Bₙ u[n]​ = (Aₙ - BₙKₙ) x[n]​

linear quadratic observer有自己的P和K。
linear quadratic controller也有自己的P和K。
而且兩邊的P和K意義不同。

兩邊都是二次型最佳化求公式解，
兩邊的P和K地位相仿。
方便起見，大家將數學符號設為相同。

多目標最佳化：
一、內部狀態限制　x[n]ᵀ Qₙ x[n] for all n
二、控制訊號限制　u[n]ᵀ Rₙ u[n] for all n

正常來說Qₙ = CₙᵀCₙ、Rₙ = 0。

separation principle用在Pₙ。

(1) Lagrange multiplier
    https://math.stackexchange.com/questions/3119575/
(2) Hamilton–Jacobi–Bellman equation
    https://math.stackexchange.com/questions/3776909/

兩種推導方式等價。
中規中矩的做法是(1)。
教科書總是介紹(2)。

Hamilton–Jacobi–Bellman equation / Riccati equation

訊號學家自創一個詞彙Hamilton–Jacobi–Bellman equation。
Bellman發明動態規劃。
Jacobi發展微分方程式、動態系統。
Hamilton發明辛動態系統。

Bellman's equation源自應用數學。
它是動態規劃當中的遞迴公式。
目標函數，可以改寫成遞迴公式。

Hamilton-Jacobi equation源自物理學。
它是辛動態系統的一個關係式。
目標函數，其格式宛如作用量。
遞迴公式，其格式宛如Hamilton-Jacobi equation。

訊號學家自創一個詞彙Riccati equation。
Riccati equation源自數學。它是一種一階二次微分方程式的解。
HJE equation，右式移除遞迴項，就是一階二次微分方程式。
導致HJE equation的解就是Riccati equation。

訊號學家自創一個詞彙infinite horizon。
其實就是converge to steady state at infinity。

r[n] = 0

Dₙ≠0，目標函數從二次型變成二次函數。
常數項不影響最佳解位置，可以消除。
一次項用配方法消除。

objective function:

Jₙ = sum { x[i]ᵀ Qᵢ x[i] + 2 x[i]ᵀ Nᵢ u[i] + u[i]ᵀ Rᵢ u[i] }
    i=0⋯n

where Q = CᵀC , N = CᵀD , R_e​ = R + DᵀD

solution:

P = AᵀPA + Q - (AᵀPB+N) (R_e + BᵀPB)⁻¹ (BᵀPA+Nᵀ)
  = AᵀPA + CᵀC - (AᵀPB + CᵀD) (R+DᵀD+BᵀPB)⁻¹ (BᵀPA+DᵀC)
K = (R+DᵀD+BᵀPB)⁻¹(BᵀPA+DᵀC)​

r[n] is constant

當參考訊號是常數函數而不是零。
最佳控制器追加一個Kₛₛ，對付參考訊號。

objective function:

Jₙ = sum { (x[i] - r)ᵀ Qᵢ (x[i] - r) + u[i]ᵀ Rᵢ u[i] }
    i=0⋯n

linear quadratic controller:
yₛₛ → r

optimal controller:

u[n] = - K x[n] - Kₛₛ r

steady-state feedforward gain Kₛₛ:

x[n+1] = (A - BK) x[n] + B Kₛₛ r
xₛₛ    = (A - BK) xₛₛ  + B Kₛₛ r
xₛₛ     = (I - A + BK)⁻¹ B Kₛₛ r
yₛₛ     = [(C - DK) (I - A + BK)⁻¹ B + D] Kₛₛ r
Kₛₛ     = [(C - DK) (I - A + BK)⁻¹ B + D]⁻¹

PI controller

P controller改成PI controller。
目標函數是累進誤差z。
累進誤差z作為擴充變數。

LTI state-space model:

⎰ x[n+1] = A x[n] + B u[n]
⎱ y[n]   = C x[n] + D u[n]

accumulative error (integral action):

z[n+1] = z[n] + (r - y[n])
       = z[n] + r - C x[n] - D u[n]

augmented internal state:

⎡ x[n+1] ⎤ = ⎡  A O ⎤ ⎡ x[n] ⎤ + ⎡  B ⎤ u[n] + ⎡ O ⎤ r
⎣ z[n+1] ⎦   ⎣ -C I ⎦ ⎣ z[n] ⎦   ⎣ -D ⎦        ⎣ I ⎦
  x'[n+1]        A'    x'[n]        B'

O: zero matrix
I: identity matrix

objective function:

J[n] = sum { x'[i]ᵀ Q x'[i] + u[i]ᵀ R u[i] }
      i=0⋯n

optimal controller:

u[n] = - Kₓ x[n] - K₞ z[n]

r[n] is arbitrary

我沒有學會。
大家自己看著辦吧。
總之公式解很長一串。

system with noise

請自行翻閱專著。
總之變成二次最佳化。利用配方法。

stochastic LTV state-space model:

⎰ x[n+1] = Aₙ x[n] + Bₙ u[n] + w[n]
⎱ y[n]   = Cₙ x[n] + Dₙ u[n] + v[n]

延伸閱讀：dynamic programming

dynamic programming

linear quadratic controller的解，
只需要貪心法，不需要動態規劃。
只需要公式解，不需要遞迴公式。
即便是非線性系統，也可以用線性化硬撐過去，不需要動態規劃。
不過還是順便介紹一下動態規劃。

cost[n] = x[n]ᵀ Qₙ x[n] + u[n]ᵀ Rₙ u[n]

J[n] =   min   {  sum  cost1[i] }
      u[0]⋯u[n]  i=0⋯n

     =   min   { J[n-1] + cost1[n] }
         u[n]

每個時刻窮舉輸入訊號數值u[n]，將最小成本儲存於表格。
因此只能處理數位訊號（離散）、無法處理類比訊號（連續）。
時間複雜度O(RN)。空間複雜度O(RN)。
R是訊號數值種類，N是訊號長度。

regulator:
反向查表，以便找到每個時刻的最佳解。
每回合追加一個新訊號，過去的選擇往往不再是最佳解。
解法是改變遞推方向（這個目標函數雙向都可以遞推）。

自行設定計算範圍L。
假設接下來L個新訊號足以達到穩態。
初始值設定為穩態的目標函數數值（正常來說是零）。
時間複雜度O(RL)。

tracker:
每回合都整個重算一遍。
每回合都重新對準了穩態，朝著穩態前進。
計算過程宛如滑動視窗，視窗寬度L。
每回合時間複雜度O(RL)。總時間複雜度O(RLN)。

視窗寬度越長，每回合計算時間越久。
所幸計算時間是固定的。
儘管無法即時得到最佳解，但是可以定時得到最佳解。
Δt可以設定成計算時間O(RL)。

Luenberger observer

專著《Control System Design Guide: Using Your Computer to Understand and Diagnose Feedback Controllers》。

Luenberger observer

x̂[n+1] = f(x̂[n], u[n], l(y[n] - ŷ[n]))
找到一個幾乎一樣的state-space model。
觀測器輸入訊號：輸入訊號u[n]、輸出訊號誤差y[n] - ŷ[n]。
觀測器輸出訊號：狀態估計x̂[n+1]。

如果知道四個矩陣ABCD，即可遞推計算內部狀態x。不必繞圈子。
如果不知道四個矩陣ABCD，Luenberger observer不可靠。

output error:

         ┌──────────┐  y
u ──┬───→│  system  │──┬──┐
    │    └──────────┘  │  │
    │ ╭────────────────╯  ↓+
    │ │  ┌───────────┐    ⊕──→ y - ŷ
    │ ╰─→│   model   │    ↑-
    └───→│ estimator │────┘
         └───────────┘  ŷ

observer:

     ┌──────────┐
u ──→│ observer │──→ x̂
y ──→│          │
     └──────────┘

Luenberger observer:

                               ┌───┐
    ┌─────────────────────────→│ B │──┐
    │  ┌───────────┐           ╞═══╡  ↓+
u ──┴─→│   model   │──→ x̂ ────→│ A │─→⊕──┬─→ x̂
y ──┬─→│ estimator │──→ ŷ ──┐- ╞═══╡  ↑+ │
    │  └─────┬─────┘        ⊕─→│ L │──┘  │
    └────────↑──────────────┘+ └───┘     │
             ╰───────────────────────────╯
所以那個C和D咧？

https://commons.wikimedia.org/wiki/File:Luenberger_Observer.svg
https://www.researchgate.net/figure/Schematic-diagram-of-the-Luenberger-observer_fig1_352751281

complementary filter

Luenberger observer可以視作互補濾波器。
輸入訊號與輸出訊號灌入互補濾波器，藉以消除感測器的影響。
詳情請見專著：
https://books.google.com.tw/books?id=hTXsXTt5kewC&pg=PA199 

P controller

   ┌──────────────┐  u  ┌─────────┐
╭─→│ P controller │──┬─→│  plant  │──┬───────────────→ y
│  │ u[n] = Fx̂[n] │  │  └─────────┘  ╰→┌──────────┐
│  └──────────────┘  ╰────────────────→│ observer │──╮
│                                      └──────────┘  │
╰────────────────────────────────────────────────────╯
                          x̂

追加u[n] = Fx̂[n]形成閉迴路系統。
整體也得滿足穩定性。
特徵值變成A+BF和A+LC。

https://www.columbia.edu/~ja3451/courses/E6602/9_observers.pdf
https://control.asu.edu/Classes/MAE507/507Lecture09.pdf

control system design

專著《Control System Design》。

曾經介紹：

(1) 受控廠：LTI system
    控制器：PID controller
    數學工具：Laplace transform
    主題名稱：feedback control

此處介紹：

(1) 受控廠：LTV state-space model
    觀測器：linear quadratic observer (Kalman filter)
    數學工具：linear least squares
    主題名稱：optimal estimation

(2) 受控廠：LTV state-space model
    控制器：linear quadratic controller
    數學工具：quadratic optimization
    主題名稱：optimal control

continuous-time system

continuous-time LTV state-space model:

⎰ x′(t) = A(t) x(t) + B(t) u(t)
⎱ y(t)  = C(t) x(t) + D(t) u(t)

以上僅介紹離散時間系統。
沒有介紹連續時間系統。
若有興趣請自行尋找資料。

連續時間系統，答案完全不同。
上述所有內容都要重新推導一遍。阿娘喂。

theorem / principle / criterion

separation principle
certainty equivalence principle
internal model principle
Pontryagin's maximum principle
Kalman decomposition

intractable problem / open problem

NP-hardness of some linear control design problems
https://web.mit.edu/jnt/www/Papers/J067-97-vb-contr.pdf

Open Problems in Mathematical Systems and Control Theory
https://perso.uclouvain.be/vincent.blondel/op/