system🚧
system
「系統」。多變量函數，輸入多道訊號，輸出多道訊號。
已經發展成熟，累積大量懸案，最近三十年沒有太大變化。
電力、機械相關科系基礎知識。暫且先做簡介，有空再來展開。
課程關係圖
多變量微積分 微分方程 動態系統 訊號與系統 聲學
↘ ↘ ↘ ↘ ↘
力學 → 電磁學 → 電路學 → 化學動力學 → 自動控制 → 電聲學
signal
digital signal / analog signal
訊號的數值，離散的和連續的。
SISO system / MIMO system
輸入輸出只有一串訊號 輸入輸出有許多串訊號
spectrum
Fourier transform / time domain / frequency domain
因為真實世界的訊號幾乎都是一堆波， 所以大家用傅立葉轉換，把訊號分解成波。 原本訊號稱作時域（座標軸是時間），傅立葉轉換之後稱作頻域（座標軸是頻率）。
spectrum / Bode plot / amplitude spectrum / phase spectrum
訊號實施傅立葉轉換（時域轉頻域）。 一串數列的傅立葉轉換是一串數列，每個數值都是複數。 複數長度是振幅。頻譜每個數值的振幅，形成振幅頻譜。 複數角度是相位。頻譜每個數值的相位，形成相位頻譜。 兩者合稱頻譜。 振幅頻譜橫軸與縱軸都取log，相位頻譜橫軸取log，兩者合稱波特圖。
frequency response / transfer function / gain / phase shift
訊號實施傅立葉轉換（時域轉頻域）。 一串數列的傅立葉轉換是一串數列， 每個數值代表每種頻率(frequency)的波(wave)的振幅(amplitude)和相位(phase)。 值得一提的是，除了訊號可以分解成波，其實系統也可以分解成波。 系統弄成加權總和的形式，然後把權重當作訊號，拿去分解。 （因為z-transform和convolution theorem） 兩種方式可以得到系統的頻譜： 一、系統的頻譜（權重做傅立葉轉換）。 二、輸出訊號頻譜除以輸入訊號頻譜（訊號做傅立葉轉換，再相除）。 頻率響應：系統的頻譜的數值，針對特定頻率。 轉移函數：系統的頻譜的數學式，變數是頻率。 增益：系統的振幅頻譜。各種頻率的波的振幅變化倍率。 相位偏移：系統的相位頻譜。各種頻率的波的相位變化差距。
generating function
z-transform / Laplace transform / z-domain / s-domain
實務上， 差分方程（離散訊號）／微分方程（連續訊號） 實施離散傅立葉轉換／連續傅立葉轉換， 改寫成複弦波，形成frequency domain。 理論上， 差分方程（離散訊號）／微分方程（連續訊號） 實施z-transform / Laplace transform， 改寫成生成函數，形成z-domain / s-domain。 上述「理論上」觀點其實是錯的。很不幸地，教科書習慣介紹這種觀點。 很多人將拉普拉斯轉換解讀成z轉換的連續版本，但是這種解讀不太對。 拉普拉斯轉換其實也有離散版本、應當先講離散版本。 嚴格來說，z轉換是普通生成函數，拉普拉斯轉換是指數生成函數。性質截然不同。 錯誤的解讀方式： 一、z轉換和拉普拉斯轉換兩者對比，宛如離散和連續兩者對比。 正確的解讀方式： 一、z轉換（普通生成函數）。 二、z轉換的重要特例是傅立葉轉換。可以推廣為連續版本。 三、z轉換的重要變種是拉普拉斯轉換（指數生成函數）。可以推廣成連續版本。 四、z轉換改寫成拉普拉斯轉換，透過bilinear transform。 五、可以推廣為連續版本，其前提是積分結果受限（不會是正負無限大）。
https://math.stackexchange.com/questions/105800/
Fourier transform / Laplace transform
傅立葉轉換：振幅為1、頻率為自然數／各種數值。 拉普拉斯轉換：振幅為各種數值、頻率為各種數值。 離散傅立葉轉換：z轉換當中，多項式代入特定數值。 分別代入N種不同數值，z = exp(𝑖(2π/N)f)，f從0到N-1。 離散時間傅立葉轉換：z轉換當中，多項式代入特定數值。 分別代入∞種不同數值，z = exp(𝑖2πω)，ω從-∞到+∞。 exp的指數從自然數f推廣成實數ω。 完整名稱應該是離散時間連續頻率傅立葉轉換。重點在於連續頻率。 連續傅立葉轉換：訊號從離散數列推廣成連續函數。 分別代入∞種不同數值，z = exp(𝑖2πω)，ω從-∞到+∞。 拉普拉斯轉換：訊號從離散數列推廣成連續函數。 分別代入∞×∞種不同數值，z = exp(s)，s是各種複數a+b𝑖。 exp的指數從實數ω推廣成複數s。 傅立葉轉換是特例，拉普拉斯轉換是通例，導致教科書很喜歡用拉普拉斯轉換。 但是現實世界無法取得數值、無法拿來計算啦。 實務上都是使用傅立葉轉換。
linear time-invariant system
causal system / linear system / time-invariant system
因果系統：遞迴函數的變數都是過去時刻。 線性系統：遞迴函數由變數的加法和倍率組成（加權總和）。 也就是說，輸入相加導致輸出相加、輸入倍率導致輸出倍率。 非時變系統：遞迴函數不因時刻而變。 也就是說，輸入移位導致輸出移位。
中譯 意譯
cause (noun) 原因 原因
cause (verb) 導致 因此而
cause (conj) 因為 because的精簡講法，原因是
causal 因果 原因的
causal LTI system
大家習慣省略字首causal。
moving average model / autoregressive model
系統有三種款式。 移動平均數模型：輸入訊號的加權總和＝輸出訊號 (all zero) 自迴歸模型：輸出訊號的加權總和＝輸出訊號 (all pole) 兩種都用：輸入訊號的加權總和＝輸出訊號的加權總和
zero / pole
z-transform：訊號=>多項式 convolution theorem：訊號加權總和=>多項式乘法 transfer function：輸出多項式作為分子，輸入多項式作為分母。 zero：分子多項式的根（多項式變成0）。 pole：分母多項式的根（多項式變成±∞）。 大家用示波器看zero pole，然後反向推理系統是什麼。
implicit form / explicit form
LTI system = constant-coefficient linear ODE = convolution LTI system可以寫成微分方程式（隱式）、寫成卷積（顯式）。 (1) implicit form (differentiation) (2) explicit form (convolution) 假設系統是線性非時變系統（例如物理運動、電路運作）。 一、最初是微分方程式。 二、改寫成transfer function（時域轉頻域）。形成多項式分式。 三、改寫成因式分解。分母的根pole、分子的根zero。 四、改寫成部分分式。形成分式連加。 五、改寫成符號解（頻域轉時域）。形成exp()連加。
system (explicit form):
y(t) = F(x(t))
system (implicit form):
L(x(t), y(t)) = 0
LTI system (implicit form) (expressed by differentiation):
L(x(t), x′(t), x″(t), ..., y(t), y′(t), y″(t), ..., ) = 0
linear differential equation (with zero initial condition):
a₀ x(t) + a₁ x′(t) + a₂ x″(t) + ...
= b₀ y(t) + b₁ y′(t) + b₂ y″(t) + ...
where x(0) = x′(0) = x″(0) = ... = 0
Laplace transform (time domain -> s-domain):
a₀ s⁰ X(s) + a₁ s¹ X(s) + a₂ s² X(s) + ...
= b₀ s⁰ Y(s) + b₁ s¹ Y(s) + b₂ s² Y(s) + ...
(a₀ s⁰ + a₁ s¹ + a₂ s² + ...) X(s)
= (b₀ s⁰ + b₁ s¹ + b₂ s² + ...) Y(s)
transfer function:
Y(s) a₀ s⁰ + a₁ s¹ + a₂ s² + ...
F(s) = ———— = ———————————————————————————
X(s) b₀ s⁰ + b₁ s¹ + b₂ s² + ...
polynomial factorization:
Y(s) (s-z₀)(s-z₁)(s-z₂)...
F(s) = ———— = —————————————————————
X(s) (s-p₀)(s-p₁)(s-p₂)...
partial fraction expansion (when all poles are distinct):
C₀ C₁ C₂
F(s) = ———— + ———— + ———— + ...
s-p₀ s-p₁ s-p₂
where C₀ = ... , C₁ = ... , C₂ = ...
inverse Laplace transform (s-domain -> time domain):
C₀ C₁ C₂
f(t) = ————————— + ————————— + ————————— + ...
exp(-p₀t) exp(-p₁t) exp(-p₂t)
= C₀ exp(p₀t) + C₁ exp(p₁t) + C₂ exp(p₂t) + ...
LTI system (explicit form) (expressed by convolution):
y(t) = x(t) ∗ f(t)
BIBO stability (and steady state is zero):
if |x(t)| < ∞ then |y(t)| < ∞ for all t≥0
=> if |x(t)| < ∞ then |Cᵢ exp(pᵢt)| < ∞ for all t≥0 and i≥0
=> |exp(pᵢt)| < ∞ for all t≥0 and i≥0
=> Re[pᵢt] < 0 for all t≥0 and i≥0
=> Re[pᵢ] < 0 for all i≥0
https://tutorial.math.lamar.edu/classes/de/IVPWithLaplace.aspx https://lpsa.swarthmore.edu/LaplaceXform/InvLaplace/InvLaplaceXformPFE.html
BIBO stability
輸入受限則輸出受限：可以改寫成訊號L¹-norm受限。
LTI system的情況下：所有pole都在左半複平面。
訊號學家還希望收斂至零：所有pole都在左半複平面，但不含虛軸。
滿足穩定性，才能形成穩態，才能控制。
不滿足穩定性，輸出可能正負無限大。
導致電路過載燒掉、動力機械暴衝、反應槽爆炸。
判斷穩定性，
(1) 直覺的方法是使用多項式函數求根演算法，求得所有pole。
經典演算法是companion matrix求特徵值，時間複雜度O(N³T)。
(2) 特殊的方法是使用特殊數學公式，檢查正負號。
經典演算法是Routh array，時間複雜度O(N²)。
dynamic response
已知輸入、系統，求得輸出。 輸入：已知函數（教科書習慣討論下述三種） 系統：已知函數（教科書習慣討論一次微分方程、二次微分方程） 輸出：未知函數。 1. impulse response：脈衝波。隔壁棚結構分析很常用，控制系統則不使用。 2. step response：單位步進函數。主角。例如啟動馬達至定速。 3. frequency response：複弦波。得到波特圖其中一個數值。
frequency response理論上是輸入複數exp波，實務上是輸入實數cos波。 必須自行推導cos波經過傅立葉轉換之後的振幅和相位。
steady state
已知輸入、系統，求得輸出最終數值。 輸入：教科書習慣討論步進函數（step response） 系統：教科書習慣討論一次微分方程、二次微分方程 輸出：求得穩態 步進函數：傅立葉轉換是1/s。 最終值定理：時域穩態（時間趨近無限大）＝頻域乘上s後頻率趨近無限大。 系統改寫成transfer function形成分式。觀察分母： 一、實根（一次多項式）：輸出指數衰減。 步進函數恰好跟最終值定理互相抵銷，剩下系統。 二、共軛複根（二次多項式）：輸出振盪。兩種表達方式。 甲、decay rate σ and damped frequency ωd 乙、damping ratio ζ and natural frequency ωn 其他情況： 一、連乘積：實係數多項式因式分解，總是得到一次暨二次多項式連乘積。 因此只需討論實根、共軛複根兩種情況。 最後讓transfer function相乘。 二、重根：一次暨二次多項式的次方。
block diagram
for transfer function of causal LTI system, 1. series ──→ G₁ ──→ G₂ ──→ = ──→ G₁G₂ ──→ 2. parallel ──┬─→ G₁ ──┬──→ = ──→ G₁+G₂ ──→ └─→ G₂ ──┘+ 3. feedback u e y ──┬─→ G₁ ──┬──→ = ──→ G₁/(1+G₁G₂) ──→ -└── G₂ ←─┘ y = G₁e = G₁(u-G₂y) = G₁u - G₁G₂y u = (y + G₁G₂y)/G₁ = y(1 + G₁G₂)/G₁ y/u = G₁/(1+G₁G₂)
4.
──→ G₁ ─┬─→ G ──→ = ──→ G₁ ──→ G ─┬─→
──→ G₂ ─┘ ──→ G₂ ──→ G ─┘
5.
──→ G ──┬─→ G₁ ──→ = ─┬─→ G ──→ G₁ ──→
└─→ G₂ ──→ └─→ G ──→ G₂ ──→
6.
──┬─→ G₁ ──┬──→ = ──→ 1/G₂ ──┬─→ G₂ ──→ G₁ ──┬──→
-└── G₂ ←─┘ -└───────←───────┘
signal-flow graph
Mason's rule https://en.wikipedia.org/wiki/Mason's_gain_formula
FIR system / IIR system
專著《Modeling, Identification and Simulation of Dynamical Systems》
FIR system / IIR system
脈衝響應分成兩種。 1. finite impulse response 有限脈衝響應。輸入脈衝函數，輸出很快歸零。時長有限。好人不長壽。 2. infinite impulse response 無限脈衝響應。輸入脈衝函數，輸出永不歸零。時長無限。禍害遺千年。
系統有兩種款式。 1. FIR system = MA model 開迴路。系統的數學式不含輸出訊號。 開迴路=>輸入的加權總和=>輸入只取幾項=>有限 2. IIR system = ARMA model 閉迴路。系統的數學式包含輸出訊號。 閉迴路=>輸入與輸出的加權總和=>強行展開=>輸入取所有項=>無限 FIR system和IIR system詞不達意。 唯一優點是後綴有system，可以跟LTI system並列。 那麼為啥不叫做MA system和ARMA system呢？我哪知道。
FIR system property:
(1) convergence: converge to zero at infinite
(2) stability: all poles at origin
(3) polarization: 2D plot
[u[k]]ᵀ[ |G| * ] [u[k]] = |G|² sin²(∠G)
[y[k]] [ -|G|cos(∠G) 1 ] [y[k]]
IIR system property:
(1) stability: initial condition matters
(2) two amplification factor
(3) s-plane vs. z-plane
跑去偷吃飯／偷放風，結果沒有聽全。
ARMA model
http://feldman.faculty.pstat.ucsb.edu/174-03/lectures/l7.pdf
AR(1) model:
u ───┬─────────────┬────→ y
+↑ ┌───────┐
│ │ delay │ amplifier should be a triangle
│ ┌────┐ └───────┘
└──│ ×a │←────┘
└────┘
y[k] = a y[k-1] + u[k]
^^^^^^
autoregression
impulse response:
┌──────────┐
│ 1 │
u ───→│ ———————— │───→ y
│ 1 - az⁻¹ │
└──────────┘
1
y[k] = ———————— u[k]
1 - az⁻¹
ARMA(1,1) model:
u ─────┬─────────────┬──┬─────────────┬────→ y
┌───────┐ +↑ +↑ ┌───────┐
│ delay │ │ │ │ delay │
└───────┘ ┌────┐ │ │ ┌────┐ └───────┘
└────→│ ×b │──┘ └──│ ×a │←────┘
└────┘ └────┘
y[k] = a y[k-1] + b u[k-1]
impulse response:
┌──────────┐
│ bz⁻¹ │
u ───→│ ———————— │───→ y
│ 1 - az⁻¹ │
└──────────┘
y[k] = aᵏ⁻¹ b where u[0] = 1
b bz⁻¹
G(z) = sum { aᵏ⁻¹ b z⁻ᵏ } = ——— = ————————
k=1⋯∞ z-a 1 - az⁻¹
ARX(1) model:
e
╷
+↓
u ───┬──────┴──────┬─────→ y
+↑ ┌───────┐
│ │ delay │
│ ┌────┐ └───────┘
└──│ ×a │←────┘
└────┘
y[k] = a y[k-1] + u[k] + e[k]
impulse response (with zero-mean white noise):
┌──────────┐ │ e ~ WN(0,σ²)
│ 1 │ ↓
u ───→│ ———————— │─────→ y
│ 1 - az⁻¹ │
└──────────┘
1
y[k] = ———————— u[k]
1 - az⁻¹
y[k] = a y[k-1] + u[k]
autocorrelation (with zero-mean white noise):
Ry[t] = σ² pow(a,|t|) (0 < a < 1)
σ² (1 - a²) σ² (1 - a²)
Φy[ω] = ———————————————— = ——————————————————
|1 - a exp(-𝑖ω)|² 1 + a² - 2a cos(ω)
(Φy is white iff a = 0)
Ry[t] Φy[ω]
| | __
1 | /\ | / \
| / \ | _- -_
| ____-- --_____ | __-- --___
|--------|----------> t |--------|----------> ω
0 0
(the peak is slightly higher than 1)
example:
respiration (low frequency) (slow)
ECG (high frequency) (fast)
Rx[t]
| __
1 | __--/\--__
| __-- / \ --__slow
| _____-- --_____fast
0 |---------|----------> t
0
system with internal state
state-space model / Kalman filter
未知因：autoregression model (LTI filter) 已知業：LTI filter 已知果：一道訊號 請見本站文件「filter」。 未知輸入（因）已知函數（業）已知輸出（果），求因。 精髓：果的誤差，通過反濾波器，用以修正因。
hidden Markov model / Bayes filter
未知因：Markov chain 未知業：hidden Markov model 已知果：許多道訊號 請見本站文件「hidden Markov model」。 未知輸入（因）未知函數（業）已知輸出（果），求因。 精髓：已知果，梯度下降法（反向傳播法）找最大值，得到機率最大的業。 未知果，動態規劃找最佳路線，得到機率最大的因。
representation / realization
表示法 表示法轉換
polynomial representation / matrix representation
工程師依照經驗命名，名稱不太精準。 polynomial representation = differential equation representation matrix representation = state-space representation polynomial representation還有一些細類，但是名稱不太精準。 time domain representation = impulse response representation frequency domain representation = transfer function representation
shift operator / Hankel realization
q
專著《Filtering and System Identification: A Least Squares Approach》
system with conditional function
ReLU
branch (if statement)從二元邏輯推廣成連續函數
數學家從未討論條件函數、從未建立數學理論。 導致深度學習只有工程沒有科學，只有經驗沒有理論。 即便已經發現許多應用，卻無法進行論證。 有人形容為煉金術。大家都在碰運氣、都在唬爛。 光靠唬爛就可以得到圖靈獎和諾貝爾獎。你們城裡人真會玩。 年輕時候缺錢，老婆出錢幫忙車庫創業。 創業成功發大財，開始為所欲為，偷吃員工、逛蘿莉島。 搞到高階主管紛紛離職，CEO從高加索人變印度人，金流轉向皮衣客。 為了讓公司再次偉大，搞起置入性行銷。找人來當模特，頒獎給他們。 三立鄉土劇演到一半出現業配，就是這個概念。
fuzzy logic
把布林數的AND和OR運算，變成函數的min和max運算。 Karnik–Mendel algorithm
system manipulation🚧
system manipulation
不知道該下什麼標題才好。
生成器：沒有輸入訊號。輸入一些參數，輸出一道訊號。 濾波器：系統輸出之後，追加一個系統，用來調整輸出訊號。 控制器：系統輸入之前，追加一個系統，用來調整輸入訊號暨輸出訊號。 觀測器：系統輸出之後，追加一個系統，用來觀測系統狀態。
verb | action noun | agent noun ---------|-------------|----------- generate | generation | generator filter | filter | filter control | control | controller observe | observation | observer
system manipulation — generator🚧
generator
「生成器」。函數，輸入參數，輸出訊號。
pulse generator / oscillator
脈衝生成器：產生三角形函數或者鐘形函數，寬度極窄，自訂寬度。 振盪器：產生週期函數。例如方波、三角波、鋸齒波、弦波。 實作方式採用程式語言，事情非常簡單。將函數直接寫成程式碼。 實作方式採用電子電路、生化反應，事情非常複雜。此處省略。
system manipulation — filter🚧
filter
「濾波器」。函數，輸入訊號，輸出訊號。
linear time-invariant system
low-pass filter / high-pass filter
濾波器只保留低頻波、只保留高頻波。
band-pass filter / band-stop filter
濾波器只保留中頻波、只刪除中頻波。
Shelving filter / Butterworth filter
濾波器保留低頻波或高頻波；其餘的波，頻率相差越遠、保留越少。比較平滑柔順啦。 濾波器保留中頻波；其餘的波，頻率相差越遠、保留越少。比較平滑柔順啦。 Shelving filter其實有時域公式喔！ http://www.cs.cf.ac.uk/Dave/CM0268/PDF/10_CM0268_Audio_FX.pdf
peak filter / notch filter
濾波器頻譜呈現一個尖峰、頻譜呈現一個尖谷。
difference filter / feedforward comb filter / feedback comb filter
延遲1刻 yn = xn + a xn-1 延遲d刻 yn = xn + a xn-d （頻譜的振幅呈現連綿圓峰） 改成回饋 yn = xn + a yn-d （頻譜的振幅呈現梳子）（連綿圓丘上下顛倒） 2nd-order其實就是延遲時刻有小數點，需要做線性內插。 https://thewolfsound.com/allpass-filter/
all-pass filter
濾波器頻譜，振幅不變（常數1）、相位改變。 例如feedforward comb filter與feedback comb filter串聯。 https://ccrma.stanford.edu/~jos/pasp/Allpass_Two_Combs.html
moving average filter
k點平均 yn = (xn + xn-1 + ... + xn-k+1) / k （頻譜的振幅呈現連綿縮小圓丘）
window function / spectral leakage
傅立葉轉換，只有整數倍頻率波。 如果訊號不是整數倍頻率波所組成，那就完蛋了。 非整數倍頻率波，將分散到各個整數倍頻率波，漏的到處都是。 訊號預先乘以窗函數，才做傅立葉轉換，稍微有點療效。 窗函數也可以想成是一種濾波器：連綿圓丘，消滅非整數倍頻率波。
sampling / aliasing
取樣，連續波變離散波。 根據取樣定理， 頻率超過取樣頻率兩倍的連續波（高頻連續波、取樣間隔太大）， 將得到頻率稍小的離散波（波長稍長）。 固定取樣頻率時，若上述連續波頻率增大，則上述離散波頻率減小。 從頻譜來看，差不多是往左鏡射、往左翻書，中譯疊頻。 解法是連續波事先做lowpass filter。 去除高頻連續波，讓它沒有東西用來鏡射翻書。 但是濾波器無法做到完美矩形，只能陡降斜下。 稱作sharp cutoff lowpass filter。 曲線越陡價格越貴。
system manipulation — controller🚧
controller
「控制器」。一個系統，一道輸入訊號，一道輸出訊號。調整輸入訊號，得到特別的輸出訊號。
control應用十分廣泛，是世上最實用的演算法之一。
https://www.mathworks.com/solutions/control-systems/feedback-control-systems.html https://ctms.engin.umich.edu/CTMS/index.php
linear time-invariant system
control
in practice:
plant
╭┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄╮
reference ─┬─→ controller ─→┆actuator ─→ process┆──┬──→ output
↑ ╰┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄┄╯ │
└────────────── sensor ←────────────────┘
in theory:
reference ─┬─→ controller ──→ plant ──┬──→ output
↑ ↓
└──────────────←───────────┘
時域 頻域
reference 參考訊號 r(t) R(s)
output 輸出訊號 y(t) Y(s)
controller 控制器 k(t) K(s)
plant 受控廠 g(t) G(s)
actuator 致動器
process 程序
sensor 感測器 h(t) H(s)
controller / plant
追加一個系統（控制器），用來控制原本系統（受控廠）。 控制器的輸出訊號，作為受控廠的輸入訊號。 控制器調整了輸入訊號，受控廠得到了特別的輸出訊號。
regulation / tracking
輸出訊號趨近（固定不變的）固定數值 輸出訊號趨近（即時改變的）參考訊號
open-loop controller / closed-loop controller
open-loop controller: r u y ──→ controller ──→ plant ──→ closed-loop controller: r e=r-y u y ───┬──────→ controller ──→ plant ──┬──→ -└─────────────←─────────────────┘
開迴路控制：控制器的輸入是參考訊號。 閉迴路控制：控制器的輸入是誤差。 誤差是參考訊號減輸出訊號（跟統計學家的習慣相反）。 用膝蓋想也知道，誤差訊息量更多，於是效果更好。 其實兩者可以一併使用，尤其是非線性系統。 甚至沒有必要計算誤差，直接使用參考訊號與輸出訊號，尤其是多變數系統。 效果更好的正式說法是靈敏度較低： (dT/T)/(dG/G)，控制系統變化與受控廠變化的比值。 教科書談靈敏度，喜歡以P controller的穩態誤差當作範例。只是在誤導大眾。
PID controller
e(t) = r(t) - y(t) error
u(t) = kp e(t) + ki ∫ e(t) dt + kd d/dt e(t) control
^^^^^^^ ^^^^^^^^^^^^ ^^^^^^^^^^^^ signal
proportional integral derivative
controller controller controller
y(t) = u(t) ∗ g(t) output
proportion 比例。乘上某個百分比的結果。 integral 積分。積分運算的結果。 derivative 導數。微分運算的結果。
P controller：原值，再乘上權重kp。誤差當前數值。 I controller：積分，再乘上權重ki。誤差前綴和。過往累計。 D controller：微分，再乘上權重kd。誤差相鄰差。瞬間變化。
closed-loop PID controller (in theory): ───┬──→ kp + ki/s + kd*s ──→ G(s) ──┬──→ -└───────────────←────────────────┘ rate feedback PID controller (in practice): ───┬──→ kp + ki/s ───────┬─→ G(s) ──┬──→ -│ -└── kd*s ←─┤ └────────────────←───────────────┘
PID tuning
此處討論tracking。參考訊號是步進函數（目標數值是1，當前數值是0）。 此時大家觀察下述指標，評定優劣。 rising time tr 到達1的時間（從0.1到0.9的時間，避免計入緩速時段） settling time ts 到達穩態的時間（保持在0.99到1.01之內） peak time tp 到達第一個局部極值的時間（最高峰） peak overshoot Mp 超過1的部分（最高峰減1，單位百分比） 此時PID controller調整參數，輸出訊號會有下述效果。 kp 調整斜率、改變頻率。增益（效果是放大縮小）。 ki 調整平均、改變穩態。高通濾波器（效果是鋸齒化） kd 調整曲率、改變振幅。低通濾波器（效果是平滑化）。 Ziegler–Nichols tuning：古聖先賢發明了經驗公式。兩種調參方法。 1. quarter decay ratio 2. ultimate sensitivity method https://my.ece.utah.edu/~ece3510/Notes_PID_Tuning_long.pdf
程式範例，受控廠是二次微分方程、0 zero 2 poles、-1 -2。
詳細圖解，位於影片後段。
stability analysis
pole–zero plot
極零圖：transfer function的poles與zeros。極Ｘ零Ｏ。 畫出poles和zeros的位置，以便判斷穩定性。 如果極Ｘ都在左半複平面（開區間、不含虛軸），則穩定。 如果右半複平面沒有極Ｘ（閉區間、包含虛軸），則穩定。 transfer function的poles，就是分母的zeros。 畫出分母的極零圖，亦可判斷穩定性。 如果右半複平面沒有零Ｏ，則穩定。
PID tuning的四個指標，可以畫成圖形，併入極零圖。
參考訊號是步進函數，R(s) = 1/s。
參考訊號R(s)、閉／開迴路系統T(s)，串聯就是乘法，
得到輸出訊號Y(s) = R(s)T(s)。
極零圖基本不變，只多了一個pole位於原點。
教科書只針對開迴路，受控廠G(s)是二次微分方程，沒有控制器K(s) = 1。
（教科書討論開迴路，但是照理應該討論閉迴路。因此以下結論沒有實用價值。）
（討論閉迴路，結果相當複雜，缺乏美感。）
（作者故意將閉迴路改成開迴路，然後挪到前面章節，我猜是為了美化。）
R(s) = 1/s unit step function
K(s) = 1 no controller
ωn²
G(s) = ——————————————————— second-order causal LTI system
s² + 2 ζ ωn s + ωₙ²
ωn²
Y(s) = R(s)K(s)G(s) = —————————————————————— open-loop
s(s² + 2 ζ ωn s + ωₙ²) controller
y(t) = 1 - exp(-ζ ωn t) sin(ωd t + ϕ) / sqrt(1 - ζ²)
where ωd = ωn(1 - ζ²)
ϕ = tan⁻¹(sqrt(1 - ζ²) / ζ)
根據y(t)，推導tr/ts/tp/Mp的滿足條件（不等式），畫在複平面。
可行解位於交集（四張圖片左半複平面重疊區域）。
https://cf.ppt-online.org/files/slide/c/C4t2nPmwWsDjy96FoSxvVrU7Iq3gMHBY5ziX1l/slide-32.jpg
root locus plot
https://control.asu.edu/Classes/MAE318/318Lecture12.pdf
根軌跡圖：transfer function附帶一個參數（例如kp）。
窮舉參數值，畫出poles的變化軌跡，以便判斷穩定性。
1. open-loop controller:
r u y
───→ controller ──→ plant ────→
K(s) G(s)
transfer function:
K(s)G(s)
2. closed-loop controller:
r e=r-y u y
───┬──────→ controller ──→ plant ──┬──→
-↑ K(s) G(s) ↓
└─────────────←─────────────────┘
transfer function:
K(s)G(s)
————————————
1 + K(s)G(s)
3. 改寫成分式
K(s) = nᴋ(s) / dᴋ(s)
G(s) = nɢ(s) / dɢ(s)
transfer function:
K(s)G(s) nᴋ(s)nɢ(s)
———————————— = ———————————————————————
1 + K(s)G(s) dᴋ(s)dɢ(s) + nᴋ(s)nɢ(s)
4. P controller:
K(s) = kp
nᴋ(s) = kp
dᴋ(s) = 1
transfer function:
K(s)G(s) kp G(s) kp nɢ(s) nɢ(s)
———————————— = ——————————— = ———————————————— = ————————————————
1 + K(s)G(s) 1 + kp G(s) dɢ(s) + kp nɢ(s) dɢ(s)/kp + nɢ(s)
transfer function的poles，就是分母的根。
1 + kp G(s) = 0 或 dɢ(s) + kp nɢ(s) = 0 或 dɢ(s)/kp + nɢ(s) = 0
注意到，如果今天不是採用P controller，那麼需要重新推導。
現在要畫出transfer function的poles軌跡。kp = 0⋯∞。
當kp = [0,∞)，根軌跡是分母1 + kp G(s)的根。
當kp = 0，根軌跡起點恰是dɢ(s)的根，即是G(s)的poles。
當kp → ∞，根軌跡終點恰是nɢ(s)的根，即是G(s)的zeros。
採用P controller的情況下，K(s) = kp只是一個倍率。
此時G(s)的zeros/poles，
恰是open-loop transfer function K(s)G(s)的zeros/poles。
導致大家認為root locus的起點和終點就是開迴路的poles/zeros。
然而一般情況下根本無法牽扯到開迴路。成為歷史共業。
古聖先賢發明了手工製圖方法。好幾條規則。
https://www.mit.edu/people/klund/weblatex/node8.html
5. closed-loop controller with sensor:
r e=r-y u y
───┬──────→ controller ──→ plant ──┬──→
-↑ K(s) G(s) │
│ │
└─────────── sensor ←───────────┘
H(s)
transfer function:
K(s)G(s)
————————————————
1 + K(s)G(s)H(s)
教科書定義K(s)G(s)H(s) = k L(s)，但是內文根本沒有用到。來亂的。
Nyquist plot
https://lpsa.swarthmore.edu/Nyquist/NyquistStability.html https://ocw.mit.edu/courses/18-04-complex-variables-with-applications-spring-2018/44f1db513a6a17d655abe0b6ff7748fc_MIT18_04S18_topic11.pdf
Nyquist圖：實施下述變換。 輸入：圍線（封閉路徑）（點集合），順時針圍住右半複平面。s = (-𝑖∞,+𝑖∞) 函數：逐點對應，s -> 1+K(s)G(s)。 輸出：新圍線。稱作Nyquist圖。 總結：右半複平面圍線，每一點s計算1+K(s)G(s)，逐點描繪新圍線。 性質：s = (-𝑖∞,0]與s = [0,+𝑖∞)的新圍線呈上下鏡面對稱。 取巧：加一就是圍線往右位移。大家習慣畫K(s)G(s)，再用-1取代原點。 取巧：當K(s) = kp，大家習慣畫G(s)，再用-1/kp取代-1。 Cauchy's integral theorem： 複變函數f(z)，圍線積分路徑圍住零個洞，圍線積分是0。 Cauchy's residue theorem： 複變函數f(z)，圍線積分路徑圍住多個洞，圍線積分是2π𝑖乘上留數和。 Cauchy's argument principle： 複變函數f′(z)/f(z)，圍線積分路徑圍住P個極、Z個零，圍線積分是2π𝑖(Z-P)。 援引winding number，上述圍線積分重新視作逆時針繞圈(Z-P)次。 Cauchy's argument principle： 複變函數f(z)，有多個極零。 一條圍線，逆時針繞圈一次，圍住P個poles、Z個zeros， 該條圍線實施f(z)變換之後， 一條圍線，逆時針繞圈(Z-P)次，圍住原點。 Nyquist stability criterion： 閉迴路系統極零圖，右半複平面不含閉迴路系統的pole，則穩定。 閉迴路系統K(s)G(s)/(1+K(s)G(s))的pole，就是分母1+K(s)G(s)的zero。 分母1+K(s)G(s)，有多個極零。 分母極零圖，右半複平面不含分母1+K(s)G(s)的zero，則穩定。 分母極零圖，右半複平面圍線，沒有圍住分母1+K(s)G(s)的zero，則穩定。 分母極零圖，令右半複平面有P個pole、Z個zero，令圍線是逆時針。 Nyquist圖，逆時針繞圈(-P)次，且圍住原點，則Z=0，則穩定。 （必須事先知道P是多少。因此此定理不實用。） 分母極零圖，右半複平面圍線，習慣畫順時針。 Nyquist圖，習慣畫K(s)G(s)而非1+K(s)G(s)，用-1取代原點。 Nyquist圖，逆時針繞圈P次，且圍住-1，則穩定。 採用P controller的情況下，K(s) = kp只是一個倍率。 Nyquist圖，習慣畫G(s)而非K(s)G(s)，用-1/kp取代-1。 Nyquist圖，逆時針繞圈P次，且圍住-1/kp，則穩定。
兩種製圖方式。 一、已知系統，以紙筆計算： 先畫波特圖（頻譜），再依此畫Nyquist圖。 s = (-𝑖∞,+𝑖∞)恰好對應傅立葉轉換的每種頻率的波。 二、未知系統，以儀器測量： 大家假設K(s)G(s)的pole比zero數量多、分母比分子次方高， 當s → ∞，則分母1+K(s)G(s) → 1。Nyquist圖可以畫得出來。 即便系統不穩定，Nyquist圖在下述情況還是畫得出來： 分母極零圖，右半複平面圍線，沒有途經分母1+K(s)G(s)的zero。 也就是說，閉迴路系統的pole不在虛軸上面。
Bode plot
https://lpsa.swarthmore.edu/Bode/BodeReviewRules.html
波特圖：系統的頻譜。分為振幅頻譜和相位頻譜。 振幅頻譜：採用log-log plot。橫軸頻率取log、縱軸振幅取log。 相位頻譜：採用semi-log plot。橫軸頻率取log、縱軸角度。 Bode stability criterion： if open-loop K(s)G(s) is stable and |K(s)G(s)| < 1 for all s: ∠K(s)G(s) ≡ 180° (mod 360°) then closed-loop (K(s)G(s))/(1+K(s)G(s)) is stable. 先看相位頻譜，-180°是哪幾個頻率。 再看振幅頻譜，這幾個頻率的振幅均小於1，則穩定。 這是利用開迴路來看閉迴路是否穩定。 實務上恰恰相反。大家利用閉迴路來讓開迴路變得穩定。
兩種製圖方式。 一、已知系統，以紙筆計算： 針對LTI system、並且已知zero/pole。 一、根是零：振幅一段：過原點45°降線。（原點取log之後是1） 相位一段：-90°水平線。 二、實根：振幅兩段：0°水平線、45°降線。 相位三段：0°水平線、45°降線、-90°水平線。 分裂點：彎曲過渡，其截距3dB。 三、共軛複根：振幅兩段：0°水平線、分裂點隆起、45°降線。 相位兩段：0°水平線、分裂點漸變、-180°水平線。 四、重根：振幅：水平線高度乘上倍率。降線斜率乘上倍率。 相位：水平線高度乘上倍率。 倍率是重根次數。 五、上述都是極。極零升降相反。 然而現在大家都用電腦軟體製圖。上述手法只能用來人工驗算。 二、未知系統，以儀器測量： 針對LTI system、不知zero/pole。 系統輸入：特定頻率的弦波，振幅一、相位零。 系統輸出：以儀器測量其振幅和相位，描出波特圖一點。 （即是frequency response。） 如果系統不穩定，系統輸出無限大，波特圖有些頻率畫不出來。
compensator
compensator = filter
補償器用來追加poles或zeros。用途如同filter。 根據transfer function串聯乘法原理，補償器接在控制器前面或後面都行。 PID controller + lead compensator是常見組合。 lead compensator ≈ high-pass filter ≈ PD controller lag compensator ≈ low-pass filter ≈ PI controller notch compensator ≈ band-pass filter lead compensator K(s) = (s-z)/(s-p) and |z| < |p| 左Ｘ右Ｏ lag compensator K(s) = (s-z)/(s-p) and |z| > |p| 左Ｏ右Ｘ
non-minimum phase system
右半複平面出現zeros。 當參考訊號是步進函數，則輸出訊號是先蹲後跳、聯結車轉彎。一開始衝向負值。 沒救了。compensator沒有辦法解決這種情況。 一種直覺的方式是追加poles抵銷zeros，分母分子約分之後一起消失不見。 然而實務上無法完全對準。 誤差、設備老化，都會導致zeros偏移。 稍有差池，輸出訊號就會偶然出現正負無限大。 導致電路過載燒掉、動力機械暴衝、反應槽爆炸。 非常危險。 實務上不能追加右半複平面poles抵銷右半複平面zeros。 我不知道有沒有其他解法。也許根本不需要解，順其自然就好。
gain margin / phase margin
兩個指標，用來粗略判斷前述四個調參指標以及穩定性。 增益邊界：相位為-180°=-π的頻率（波特圖相位曲線穿越橫軸之處）的振幅，減去0dB=1。 相位邊界：振幅為0dB=1的頻率（波特圖振幅曲線穿越橫軸之處）的相位，減去-180°=-π。 lead/lag compensator直接影響這兩個指標。
type 0/1/2 system
有0/1/2個pole等於零。 兩種出現情況： 一、補償器追加pole。 二、輸入訊號是constant/unit step/ramp function。 如果是情況二，穩態定義必須隨之改變， 例如零函數／零次常數函數／一次直線函數／二次拋物線函數。
system with internal state
optimization
(1) PID control: regression by formula (2) optimal control：regression by optimization
Hamilton–Jacobi–Bellman equation / Riccati equation
動態規劃，找到最好的動作方式。
Sylvester equation / Lyapunov equation
矩陣方程式，常見的最佳化問題的解。
system with parameters
(1) robust control
system with additional known variables
(2) adaptive control
system with additional unknown time-invariant variables
(3) model predictive control
system with additional constraints
nonlinear system
sliding mode controller
system manipulation — observer🚧
observer
「觀測器」。已知系統、輸入訊號、輸出訊號，求得系統狀態。
state transition state estimation (state observer) https://control.asu.edu/Classes/MAE507/507Lecture09.pdf
一階微分 = 取前一個時刻 二階微分 = 取前兩個時刻 markov process就是引入機率的連續動態系統 markov chain就是引入機率的離散動態系統
Luenberger observer
sliding mode observer
https://medium.com/@saurav310304/f9fc11c0177a
system indentification🚧
system indentification
system identification / system estimation
給你輸入訊號、輸出訊號，請你找到系統。 換句話說，就是迴歸！確定訊號用前者，隨機訊號用後者。
mean squared error / Kullback–Leibler divergence
兩道訊號的距離。 換句話說，就是誤差！迴歸的時候拿來用吧。
LTI system
frequency invariance
frequency invariance of LTI system: 針對LTI system，輸入弦波，輸出也是弦波。 而且頻率不變，僅振幅和相位改變。 【待補證明】
estimate
系統沒有誤差項，事情非常單純。 時域：輸出與輸入的反卷積 頻域：輸出頻譜除以輸入頻譜 系統擁有誤差項，才有討論意義。 衍生許多演算法。
time domain method:
u[k] = cos(ωk)
|G(ω)| = max {y[k]} / max {u[k]} magnitude
∠G(ω) = argmax Ruy[k] phase
= argmax { xcorr(u[k], y[k]) }
【what is the sampling rate?】
【how to do fast DTFT?】
frequency domain method:
G(ω) = Y(ω)/U(ω)
input signal
輸入特殊訊號，直接量頻譜。 (1) chirp: 餘弦波頻率漸增，依序得到輸出訊號頻譜每個bin。 (2) white noise: 訊號頻譜是常數函數，方便計算系統的gain。 (3) pseudorandom binary sequence：針對數位訊號。功能類似white noise。
sampling / filtering
數位訊號有spectral leakage。 解法：輸出平滑化／輸入窗函數。 數位訊號要考慮sampling rate。
FIR system with error (MAX model)
Gaussian noise / white noise
Gaussian: all samples of signal form Gaussian distribution
(time domain)
white: all bins of power spectrum form uniform distribution
(frequency domain)
Gaussian：每一個訊號都是浮動數值，呈高斯分布。 white：能量／功率頻譜是常數函數。
system identification
time domain method:
(1) correlation: least squares estimation
input signal has time-invariant autocorrelation.
e.g. weakly stationary process
(2) correlation spectrum: H₁ estimate and H₂ estimate
input signal has time-invariant autocorrelation.
e.g. weakly stationary process
frequency domain method:
(3) frequency response: ???
input signal is sinusoid.
e.g. cosine wave
(4) spectrum: empirical transfer function estimate
input signal is purpose-built.
e.g. white noise
correlation
correlation
兩道訊號，取兩個時間點相乘(SISO)／點積(MIMO) autocorrelation function: Rxx[t1,t2] = E[x[t1] x[t2]] cross-correlation function: Rxy[t1,t2] = E[x[t1] y[t2]]
time-invariant correlation
當time-invariant，可以簡化成位移量t = t2 - t1。位置k = t2。 autocorrelation function: Rxx[t] = E[x[k] x[k-t]] cross-correlation function: Rxy[t] = E[x[k] y[k-t]] where k is arbitrary property: Rxy[t] ≠ Ryx[t] cross-correlation is not commute Rxy[t] = Ryx[-t] useless. negative input is not defined.
dynamic response
impulse response 針對FIR系統，輸入訊號採用脈衝函數，輸出訊號很快歸零。 輸出訊號歸零之後，形成time-invariant correlation，即可使用此方法。 不覺得這很奇怪嗎？零有什麼好算。 step response 針對FIR系統，輸入訊號採用步進函數，輸出訊號很快變成常數，稱作DC gain。 輸出訊號收斂之後，形成time-invariant correlation，即可使用此方法。 不覺得這很奇怪嗎？常數有什麼好算。 輸入訊號到底該用什麼函數才對呢？ 白雜訊嗎？
ergodic process
(1) ensemble average
calculate expected value E[] with multiple experiments
(2) time average
calculate expected value E[] with 1 experiments
by sliding window
ergodic process: ensemble average = time average
多次實驗的平均值，改成滑動視窗的平均值，節省計算時間。
so far, so good.
sliding window
t1 t2
|--> |-->
experiment 1: ~~~~~~~~~~~~~~~ |
experiment 2: ~~~~~~~~~~~~~~~ | ensemble average
: V
-------------->
time average
Rx[t] = 1/N sum { x[k] x[k-t] }
k=0⋯N-1
Rxy[t] = 1/N sum { x[k] y[k-t] } (not commute)
k=0⋯N-1
least squares estimation
FIR system with zero-mean white noise
│ e (zero-mean white noise)
┌─────┐ ↓
u ───→│ g │─────→ y
└─────┘
y[t] = sum g[k] u[t-k] + e[t]
k=0⋯∞
assume e and u are independent
independent => uncorrelated E[e[t] u[k-t]] = 0
用correlation消除雜訊影響
Ryu[t] = E[y[k] u[k-t]]
= E[(......) u[k-1]]
= sum g[k] u[t-k] + E[e[t] u[k-t]]
^^^^^^^^^^^^^^
= 0 since uncorrelated
解一次方程組得到系統 given Ryu and Ru, solve g
assume g[t] = 0 when t > N
Ryu[t] = sum g[k] Ru[t-k] (N+1 points)
k=0⋯N
[ Ryu[0] ] [ Ru[0] ... Ru[N] ] [ g[0] ]
[ : ] = [ : : ] [ : ]
[ Ryu[N] ] [ Ru[-N] ... Ru[0] ] [ g[N] ]
Ryu Ru g
Toeplitz matrix
g = Ru⁻¹ Ryu
I prefer Cyu. However, Ryu is better than Cyu for memorizing, I guess.
correlation spectrum
correlation spectrum
autocorrelation spectrum Φx[ω] = sum { Rx[t] exp(-𝑖ωt) }
cross-correlation spectrum Φxy[ω] = sum { Rxy[t] exp(-𝑖ωt) }
t=-∞⋯+∞
正式名稱如下。名稱太長了，很少人使用。 autocorrelation spectral density function cross-correlation spectral density function 統計學與訊號學當中， spectral density function是指定義域是頻域的函數（統計而得）。 大家習慣簡稱spectral density，省略後綴function。 其中spectral是指定義域是頻域。 其中density function是指機率密度函數（積分等於一、已經正規化）。 但是大家不一定有做正規化。 spectrum也是定義域是頻域的函數。 因此大家常常將spectral density改稱為spectrum。 畢竟大家不一定有做正規化。無傷大雅。 甚至有人省略主角correlation。喪心病狂。縮短後名稱如下： auto spectral density / cross-spectral density autospectrum / cross-spectrum
cross-correlation spectrum 是指dtft(xcorr(x,y)) spectral cross-correlation 是指xcorr(dtft(x),dtft(y)) bilinear的情況下，兩者應該相等吧？【尚待確認】 dtft(xcorr(x,y)) = xcorr(dtft(x),dtft(y))
power spectrum
energy spectrum = squared amplitude spectrum power spectrum = squared amplitude spectrum / N energy／power差別在於沒做正規化／有做正規化。
Wiener–Khinchin theorem: autocorrelation spectrum = power spectrum periodogram: 使用autocorrelation spectrum，作為power spectrum的估計值。 time-invariant autocorrelation: 此定理的前提是兩個時間點t1 t2可以簡化成一個時間差。 完全沒有必要使用weakly stationary process。 但是教科書喜歡使用weakly stationary process。成為歷史共業。 integrability and boundness: 此定理的前提是傅立葉轉換必須存在。 傅立葉轉換是積分變換，必須滿足可積性。 可積性細分許多條件，其中包含受限性。 受限性：輸出受限（不會是正負無限大）。 積分變換的情況下：輸入受限，而且在無限遠處趨近零。輸出亦然。 即是autocorrelation approaches zero at infinite lag。 autocorrelation spectrum亦然。 stable system: 穩定的FIR system的impulse response會變成零。 穩定的IIR system的impulse response會趨近零。 進而滿足autocorrelation approaches zero at infinite lag。
H₁ estimate and H₂ estimate
此方法只能得到振幅頻譜，不能得到相位頻譜。
│ e (zero-mean white noise)
┌─────┐ ↓
u ───→│ g │─────→ y
└─────┘
Re[t] = σ² δ[t] since e is white
Φe[ω] = sum { σ² δ[t] exp(-𝑖ωt) } = σ² since e is white
Ry[t] = g[t] ∗ Ruy[t] + σ² δ[t]
Φy[ω] = G(ω) Φuy[ω] + Φe[ω] where Φe[ω] = σ²
│ d │ e (zero-mean white noise)
↓ ┌─────┐ ↓
u ─────→│ g │─────→ y
└─────┘
H₁ estimate: Ĝ₁ = Φyu / Φu = (G Φu) / Φu + Φd
H₂ estimate: Ĝ₂ = Φy / Φuy = (G Φuy + Φe) / Φuy
inequality: |Ĝ₁| ≤ |G| ≤ |Ĝ₂|
https://dsp.stackexchange.com/questions/71811/
bias–variance tradeoff
(1) time smoothing 相鄰N窗做N次，結果取平均。 (2) frequency filtering 窗函數 窗越窄窗越多，bias越大variance越小
frequency response
frequency response
cosine wave:
理論上是輸入複數exp波，實務上是輸入實數cos波。
必須自行推導cos波經過傅立葉轉換之後的振幅和相位。
手法是取real part。
impulse representation:
y[k] = sum g[i] u[k-i] = G(z) u[k]
^^^^ ^^^^
polynomial scalar
that is not related to k???
cosine wave:
let u[k] = Re( exp(𝑖ωk) ) = cos(ωk)
y[k] = sum g[i] Re( exp(𝑖ω(k-i)) )
= G(ω) Re( exp(𝑖ωk) )
= |G(ω)| cos(ωk + ∠G(ω)) 詭異的不行
輸入餘弦波，頻率ω、振幅α、相位0。
調整輸入振幅α，避免輸出訊號太弱太強而測量不到。
α cos(ωk) ┌─────┐
u ──────────→│ g │──→ y
└─────┘
u[k] = α cos(ωk)
y[k] = α |G(ω)| cos(ωk + ∠G(ω)) + e[k] + transient
看到transient忽然冒出來，就知道推導過程一定有問題。
donno the name of this estimate
如果沒有誤差，那麼很容易計算系統頻譜。
為了應付誤差，輸出做amplitude modulation。分別用餘弦波和正弦波。
cos(ωk)
╷
e ×↓ ┌───────────────┐
╷ ┌──→│N-point average│──→ Ic(ω)
α cos(ωk) ┌─────┐ +↓ │ └───────────────┘
u ──────────→│ g │────→ y ──┤
└─────┘ │ ┌───────────────┐
└──→│N-point average│──→ Is(ω)
×↑ └───────────────┘
╵
sin(ωk)
u[k] = α cos(ωk)
y[k] = α |G(ω)| cos(ωk + ∠G(ω)) + e[k] + transient
Ic(ω) = sum { y[k] cos(ωk) } = +½ α |G(ω)| cos(ωk)
k=1⋯N
Is(ω) = sum { y[k] sin(ωk) } = -½ α |G(ω)| sin(ωk)
k=1⋯N
積化和差公式
cos(a) cos(b) = ½ cos(a-b) + ½ cos(a+b)
推導過程
Ic(ω) = sum { y[k] cos(ωk) }
= sum { (......) cos(ωk) }
= ½ α |G(ω)| cos(ωk) ①
+ ½ α |G(ω)| (1/N) sum { cos(2ωk + ∠G(ω)) } ②
+ (1/N) sum { e[k] cos(ωk) } ③
②→0 as n→∞. since cos() has zero mean.
③→0 as n→∞. since e and u are independent by assumption.
hence only ① remains.
estimate Â(ω) = sqrt(Ic²(ω) + Is²(ω)) / (α/2) amplitude φ̂(ω) = -tan⁻¹(Is(ω) / Ic(ω)) phase
spectrum
empirical transfer function estimate
系統頻譜＝輸出頻譜／輸入頻譜 Ĝ = Y/U 頻譜相除＝振幅相除＆相位相減 複數相除＝長度相除＆角度相減 |Ĝ| = |Y| / |U| amplitude ∠Ĝ = ∠Y - ∠U phase ETFE，重劍無鋒大巧不工 故意造個長一點的句子，聽起來更有中二感
bias–variance tradeoff???
(1) bias (X) / absolute error (O)
|E[Ĝ-G]| < C₁ / (sqrt(N)|U|) where C₁ is constant
^^^^^^^
N變大可以減少估計誤差（平均數）
(2) variance (X) / squared error (O)
E[|Ĝ-G|²] < (Φ / U) + (C₂ / N) where C₂ is constant
^^^^^^^
N再大也減少不了估計誤差（變異數）
解法是correlation
spectral coherence???
(3) G(ω₁) and G(ω₂) are asyptocially uncorrelated.
即便ω取樣間距變小，函數曲線也不會變得比較連續平滑。
解法是frquency filtering
G
| ______
G₁ | __--- ---____
| __------__
G₂ | --- -----
------+-----+------> ω
ω₁ ω₂
IIR system with error (ARMAX model)
system identification
選個系統，跑個迴歸／最佳化。 分析一下bias和variance。
1. system model (1) type of model (2) order of model 2. least squares estimation (prediction filter) (1) linear regression (ARX model) (2) optimization
system model
AutoRegressive Moving Average model with eXogenous inputs 自迴歸移動平均模型ARMA附帶外生輸入X （此處的外生輸入是指white noise）
output error model y = (B/A)u + e ARX model y = (B/A)u + (1/A)e ARMAX model y = (B/A)u + (C/A)e Box–Jenkins model y = (B/A)u + (C/D)e
output error model:
│ e (zero-mean white noise)
┌─────┐ ↓
u ───→│ B/A │─────→ y
└─────┘
ARX model:
│ e
┌─────┐ ↓ ┌─────┐
u ───→│ B │─────→│ 1/A │───→ y
└─────┘ └─────┘
ARMAX model:
┌─────┐
e ───→│ C │──┐
└─────┘ │
┌─────┐ ↓ ┌─────┐
u ───→│ B │─────→│ 1/A │───→ y
└─────┘ └─────┘
Box–Jenkins model:
┌─────┐
e ───→│ C/D │──┐
└─────┘ │
┌─────┐ ↓
u ───→│ B/A │─────→ y
└─────┘
least squares estimation (prediction filter)
這次沒用correlation。不要問我為什麼。 總之沒用correlation，好像就叫做prediction filter。
linear regression
linear regression:
針對output model和ARX model，得以使用一次迴歸。
公式解是normal equation。
add–scale function:
output model和ARX model是加倍函數（狹義線性函數），屬於一次多項式函數。
ARMAX model和Box–Jenkins model是線性函數，不屬於一次多項式函數。
一次迴歸的迴歸函數，必須是一次多項式函數。
theroem:
(1) linear regression
(2) linear regression with zero-mean Gaussian white error
兩者公式解恰巧相同。
公式解是normal equation。
因此系統模型可以省略誤差項。
linear regression:
ĝ = argmin ε(g)
ε(g) = sum ‖y[k] - ŷ(k;g)‖²
k=1⋯N
ŷ(k;g) = g[1] y[k-1] + g[2] y[k-2] + ... + g[n] y[k-n]
linear regression (in matrix formulation):
[ y[k] ] [ y[k-1] y[k-2] ... y[k-n] ]
[ y[k+1] ] [ : : : ] [ g[1] ]
[ : ] = [ ] [ : ]
[ : ] [ ] [ g[n] ]
[ y[k+N] ] [ ]
thin matrix A
let k = 1 and N ≥ n
也可以改成k = n，甚至更大，這樣一來矩陣右上角不必補零。
solution:
g = (Aᵀ A)⁻¹ Aᵀ y
validation:
得到正確答案之後，驗證系統模型是否正確。
檢查各個時刻的誤差值，照理來說，整體呈現常態分布。
如果不是常態分布，那麼表示系統模型不正確。
optimization
gradient descent / Newton's method隨便你用 請見本站文件「optimization」。
system with internal state
subspace identification
rename
LTI system with internal state ------> subspace identification
找到迴歸函數之後，就可以預測訊號啦。 常見方式有Kalman filter和hidden Markov model。 http://www.princeton.edu/~stengel/MAE546Seminars.html
stochastic process🚧
stochastic process
stochastic process
訊號的數值，從固定的改成浮動的。
deterministic process / random process
deterministic process = evolution rule is fixed random process = evolution rule is random 訊號的數值的數學公式，確定的和隨機的。 確定的：訊號數值有明確數學公式。甚至使用了前後數值。 隨機的：訊號數值亂七八糟。甚至看不出是否使用了前後數值。 數學公式使用了過往數值，此時稱作causal。 數學公式使用了過往數值，而且每個數值的數學公式都相同，此時形成recurrence。
constant function / time-invariant function
in K-12 mathematics,
(1) constant = value remains unchanged
(2) constant/variable = value is known/unknown
in mathematics,
(1) constant function = output remains unchanged
by different input
(2) time-invariant function = function remains unchanged
by different time
in signals and systems,
(1) constant input/output = signal remains unchanged
(2) time-invariant system = system remains unchanged
stationary process / mutually independent process
stationary = looks like non-changing (aperiodic) mutually independent = looks like non-regressive (acausal) 等號兩側其實不一樣，只是表面上看起來很可能是那樣。
stationary: variables and time are independent. mutually independent: not depend on others.
stationary => identically distributed mutually independent => pairwise independent
identically distributed process = time-invariant distribution
X₁ = X₂ = X₃ = ...
stationary process = time-invariant joint distribution
X₁ = X₂ = X₃ = ...
(X₁, X₂) = (X₂, X₃) = (X₃, X₄) = ...
(X₁, X₂, X₃) = (X₂, X₃, X₄) = (X₃, X₄, X₅) = ...
:
pairwise independent process = constant marginal distribution
in 2-joint distribution
X₁ = (X₁, X₂ = x) for all x
X₁ = (X₁, X₃ = x) for all x
X₁ = (X₁, X₄ = x) for all x
: :
X₂ = (X₂, X₁ = x) for all x
X₂ = (X₂, X₃ = x) for all x
X₂ = (X₂, X₄ = x) for all x
: :
mutually independent process = constant marginal distribution
in joint distribution
X₁ = (X₁, X₂ = x₂) for all x₂
X₁ = (X₁, X₃ = x₃) for all x₃
X₁ = (X₁, X₄ = x₄) for all x₄
:
X₁ = (X₁, X₂ = x₂, X₃ = x₃) for all x₂, x₃
X₁ = (X₁, X₂ = x₂, X₄ = x₄) for all x₂, x₄
X₁ = (X₁, X₃ = x₃, X₄ = x₄) for all x₃, x₄
:
https://math.stackexchange.com/questions/3269920/ https://math.stackexchange.com/questions/1920473/
independent and identically distributed process
原本定義是 i.i.d. = identically distributed & mutually independent 但是有一個定理是 i.i.d. => stationary 因此實際上是 i.i.d. = stationary & mutually independent
https://www.statlect.com/glossary/stationary-sequence
moment / coherence
moment: statistics of one random variable coherence: statistics of two random varibles
indentical => time-invariant moment stationary => time-invariant moment and coherence pairwise independent => no 1st coherence mutually independent => no higher-order coherence
kth moment
E[Xᵏ]
kth central moment
E[(X-X̄)ᵏ] where X̄ = E[X] is mean
1st coherence
E[XY]
1st standardized coherence
E[(X-X̄)(Y-Ȳ)]
———————————————————————
(E[(X-X̄)²]E[(Y-Ȳ)²])¹⸍²
higher-order coherence
E[(XY,YZ,XZ)] C(3,2)
E[(AB,AC,AD,BC,BD,CD)] C(4,2)
moment: mean central moment: zero, variance, skewness coherence: correlation standardized coherence: correlation coefficient
weakly stationary process
weakly stationary process
教科書習慣採用weakly stationary達成time-invariant autocorrelation。 但是(1)是多餘的，成為歷史共業。
identically distributed process = time-invariant distribution stationary process = time-invariant joint distribution
stationary process: (1) time-invariant moment (mean/variance/skewness/...) (2) time-invariant coherence (autocorrelation/...) stationary process example: 1. i.i.d. process 2. steady states
weakly stationary process: (1) time-invariant 1st moment (mean) (2) time-invariant 1st coherence (autocorrelation) 性質更弱。失去了高階動差和高階相干。 property: (3) time-invariant 2nd central moment (variance) (1) & (2) => (3)。證明省略。也有人將(3)直接納入定義。 weakly stationary process example: 1. white noise 2. random phase cosine wave
white noise
white noise:
power spectrum is uniform
white noise in time domain:
never discovered
a sufficient but not necessary condition:
weakly stationary process with Wiener–Khinchin theorem
white noise in theory:
(1) time-invariant 1st moment (mean)
(2) time-invariant 1st coherence (autocorrelation)
(3) autocorrelation approaches zero at infinite lag
其中(2)和(3)足以實施逆向傅立葉轉換。
white noise in practice:
(1) constant mean
E[x[k]] = μ for all k
(1&2) constant variance
Var[x[k]] = σ² for all k
(2&3) impulse autocorrelation (uncorrelated)
Rx[t] = σ² δ[t]
實務上無視高階動差，只處理mean和variance。
實務上無視高階相干，只處理autocorrelation，並且視作不相關。
white noise example:
1. white Gaussian noise N(μ,σ²)
2. white uniform noise U(a,b)
stable system
stable system
definition: 1. stable: difference of two sequences is bounded 2. convergent: sequence approaches constant (steady state) 3. stationary: sequence is constant theorem: convergent => stationary at steady states (talking about a trivial thing like a theorem)
for linear time-invariant system, theorem: 1. stable: sequence is bounded 2. stable <=> Re[pᵢ] < 0 <=> convergent 3. convergent: sequence approaches zero (steady state) 大家習慣省略主角convergent。喪心病狂。 1. stable => stationary at steady states 2. stable <=> Re[pᵢ] < 0 (people pretend not to see 'convergent') (textbooks replace all 'convergent' with 'stable')
system design🚧
system design
「系統設計」。系統是函數、甚至是函數網路。工程師設計函數網路，達成特定任務，例如生成／濾波／控制／觀測。
我目前只知道三種流派：
一、電子電路：風靡全世界。雖然台灣是地球上最大的電子零件生產基地，台灣也有專門設計電子電路的公司，但是我不太確定台灣是否有這方面的專家。有言道：十萬青年十萬肝，ＧＧ輪班救台灣。關鍵字：電路學、積體電路設計、自動控制、控制理論。
二、生化反應：發展中。雖然台灣之前打算成立藥物代工廠、學名藥設計實驗室，但是以失敗告終。有言道：一日生科，終生科科。關鍵字：化學動力學、酵素動力學、藥物動力學、系統生物學。
三、深度學習：最近十年才剛萌芽。已經做到生成／濾波，目前做不到控制／觀測。正在經歷大浪潮，準備迎接大泡沫。有言道：站在風口上，豬都會飛。關鍵字：人工智慧、機器學習、深度學習。
system design — electronic circuit🚧
electric circuit
RC circuit / resistor–capacitor circuit
low-pass filter high-pass filter
RLC circuit / resistor–inductor–capacitor circuit
sine wave
electronic circuit
oscillator
https://en.wikipedia.org/wiki/Electronic_oscillator
synchronization
phase-locked loop
Phase Detector (PD): Compares the phase of the input signal with that of the output signal. Low-Pass Filter (LPF): Filters the output of the phase detector to produce a control voltage. Voltage-Controlled Oscillator (VCO): Adjusts its frequency based on the control voltage to lock onto the input signal.
VCO, Phase comparator and PID make a PLL.
electrical–mechanical–acoustical analogy
電流、機動、聲響，三者物理量相仿，可以類比。
https://en.wikipedia.org/wiki/Mechanical–electrical_analogies https://en.wikipedia.org/wiki/Impedance_analogy https://en.wikipedia.org/wiki/Harmonic_oscillator
mechanical impedance https://www.bksv.com/media/doc/17-179.pdf
system design — chemical kinetics🚧
chemical kinetics
反應速度取決於當前濃度。
production and decay
production: k ───→ A ȧ = k a(t) = a(0) + kt
decay:
k
A ───→
ȧ = -k a
a(t) = a(0) exp(-kt)
steady state:
ȧₛₛ = -k aₛₛ = 0
aₛₛ = 0
production and decay: k₀ k₁ ───→ A ───→ ȧ = k₀ - k₁ a a(t) = (a(0) - k₀/k₁) exp(-k₁t) + k₀/k₁ steady state: ȧₛₛ = k₀ - k₁ aₛₛ = 0 aₛₛ = k₀/k₁
closed system
irreversible conversion: k A ───→ B ⎰ ȧ = -k a ⎱ ḃ = k a closed system: ȧ + ḃ = 0
reversible conversion: k₁ A ⇌⇌⇌ B k₋₁ ⎰ ȧ = -k₁ a + k₋₁ b ⎱ ḃ = k₁ a - k₋₁ b closed system: ȧ + ḃ = 0
open system
interaction with decay: k₁ k₂ A ⇌⇌⇌ B ───→ (k₁ + k₋₁ >> k₂) k₋₁ ⎰ ȧ = -k₁ a + k₋₁ b ⎱ ḃ = k₁ a - k₋₁ b - k₂ b rapid equilibrium assumption: a and b rapidly become and remain steady states. b̃/ã = k₁/k₋₁
equivalent model: k₂k₁/(k₁+k₋₁) C ──────────────→ c̃ = ã + b̃
interaction with production and decay:
k₀ k₁ k₂
───→ A ⇌⇌⇌ B ───→ (k₁ + k₋₁ >> k₂)
k₋₁
⎰ ȧ = k₀ - k₁ a + k₋₁ b
⎱ ḃ = k₁ a - k₋₁ b - k₂ b
quasi-steady state assumption:
a is always steady state respect to b
aₛₛ/bₛₛ = k₁/(k₋₁+k₂)
cooperativity
Michaelis–Menten equation:
k₁
P + X ⇌⇌⇌ PX
k₋₁
number of occupied binding sites
Y = ————————————————————————————————
total number of binding sites
[PX]
= ——————————
[P] + [PX]
[X]
= ————————— where Kₛₛ = k₋₁/k₁ = [P][X]/[PX]
Kₛₛ + [X]
Adair equation:
4k₁
P + X₁ ⇌⇌⇌ PX₁
k₋₁
3k₂
P + X₂ ⇌⇌⇌ PX₂
2k₋₂
2k₃
P + X₃ ⇄⇄⇄ PX₃
3k₋₃
k₄
P + X₄ ⇌⇌⇌ PX₄
4k₋₄
number of occupied binding sites
Y = ————————————————————————————————
total number of binding sites
1[PX₁] + 2[PX₂] + 3[PX₃] + 4[PX₄]
= ——————————————————————————————————————————
4([P] + 1[PX₁] + 2[PX₂] + 3[PX₃] + 4[PX₄])
[X]/K₁ + 3[X]²/K₁K₂ + 3[X]³/K₁K₂K₃ + [X]⁴/K₁K₂K₃K₄
= ———————————————————————————————————————————————————————
1 + 4[X]/K₁ + 6[X]²/K₁K₂ + 4[X]³/K₁K₂K₃ + [X]⁴/K₁K₂K₃K₄
Hill equation:
[X]⁴/K₁K₂K₃K₄
Y ≈ ————————————————— when K₄ << K₁,K₂,K₃
1 + [X]⁴/K₁K₂K₃K₄
([X]/K)⁴ [X]⁴
Y ≈ ———————————— = —————————
1 + ([X]/K)⁴ K⁴ + [X]⁴
interconversion
Goldbeter–Koshland kinetics:
4a₁ k₁
Zp + Ex ⇌⇌⇌ C₁ ───→ Z + Ex
a₋₁
4a₂ k₂
Zp + Ey ⇌⇌⇌ C₂ ───→ Z + Ey
a₋₂
X
k₁ ↓
Zp ⇌⇌⇌⇌⇌⇌ Z
k₂ ↑
Y
transport: diffusion
diffusion: A ←─→ B ⎰ ȧ = -D(a-b)/V₁ ⎱ ḃ = D(a-b)/V₂ V: volume D: diffusion rate
inhibition
time series plot / phase portrait
stichiometry matrix
antagonistic network
inhibition:
n₁
k₁↓├──── ↓k₃
A ───→ B
k₂↓ k₅ ↓k₄
⎰ ȧ = k₁/(1 + (b/K₂)ⁿ¹) - k₃a - k₅a
⎱ ḃ = k₂ - k₄b + k₅a
n: inhibition rate
cross-inhibition: k₁↓ n₁ ↓k₃ A |===| B k₂↓ n₂ ↓k₄ ⎰ ȧ = k₁/(1 + (b/K₂)ⁿ¹) - k₃a ⎱ ḃ = k₂/(1 + (a/K₁)ⁿ²) - k₄b n: inhibition rate
gene toggle switch
bistable
Collins toggle switch [Gardner et al. 2000]:
┌─→ protein1---┐├--inducer1
│ ┴
gene1 gene2
┬ │
inducer2--┤└---protein2 ←─┘
⎰ ȧ = α₁/(1+aⁿ¹) - δ₁a
⎱ ḃ = α₂/(1+aⁿ²) - δ₂b
⎰ ȧ = F₁(a,b) - δ₁a
⎱ ḃ = F₂(a,b) - δ₂b
oscillator
Goodwin oscillator:
┌── Y (protein) ←─┐
↓ │
Z (metabolite) X (mRNA)
╎ ↑
└----┤ gene ──────┘
⎧ ẋ = a/(kⁿ+zⁿ) - bx
⎨ ẏ = cx - dy
⎩ ż = ey - fz
Elowitz–Leibler repressilator:
┌─→ repressor1---┐
│ ┴
gene1 gene2
┬ │
╎ ↓
repressor3 repressor2
↑ ╎ ╎
└─── gene3 ├----─┘ └-----┤reporter gene
⎧ ṁ₁ = α₀ + α/(1+p₃ⁿ) - m₁
⎪ ṁ₂ = α₀ + α/(1+p₁ⁿ) - m₂
⎨ ṁ₃ = α₀ + α/(1+p₂ⁿ) - m₃
⎪ ṗ₁ = β(m₁ - p₁)
⎪ ṗ₂ = β(m₂ - p₂)
⎩ ṗ₃ = β(m₃ - p₃)
feedback
autocatalysis (self-activation):
δ
┌─--------A ──→
↓ ↑
gene ──────┘
a/K
ȧ = α ——————— - δ a
1 + a/K ^^^ degrade
autoinhibition (self-repression):
δ
┌─--------A ──→
┴ ↑
gene ──────┘
1
ȧ = α ——————— - δ a
1 + a/K
system design — deep learning🚧
deep learning
「深度學習」。可以視作迴歸，迴歸函數是各種函數串聯並聯。
https://github.com/leemengtaiwan/deep-learning-resources
deep learning是新興理論，尚在發展當中，無法整理成篇。我將發展現況分成兩種方向，僅供參考。
邏輯推理觀點
我認為深度學習可以做到邏輯推理的歸納法。
邏輯推論logical inference宛如求解equation solving。根據一些邏輯運算式子的計算結果，判斷一個邏輯運算式子是否有唯一解。
邏輯推理logical reasoning宛如迴歸regression。根據一些邏輯運算式子的計算結果，調整一個邏輯運算式子，盡量達成唯一解。
機率分布即是二元邏輯的推廣版本。convolution即是加法運算。ReLU即是分支判斷。各種函數即是各種邏輯運算。
事先網羅足夠多的函數，事後練出足夠長的邏輯運算式子。
模組編程觀點
我認為深度學習就是在寫程式。
專案即是迴歸模型。編程即是拼湊函數，除錯即是降低誤差。運算子即是初階函數，函式庫即是進階函數。物件導向精簡函數架構，演算法形成函數單元。自行拼湊適當函數，考慮越多，效果越好。
以前是撰寫程式碼，經過編譯器，自動生成機器碼。現在是拼湊函數，經過deep learning，自動生成程式碼。封建升城堡。
與其說deep learning可以取代程式設計師，不如說deep learning是一種更高層次的編程方式。大家只是換個方式寫程式。
近況
近年來許多政商人士提倡人工智慧，呼籲大眾學習如何運用深度學習。他們認為深度學習可以彌補電子電路的不足。大家可以待在家中設計神經網路，透過電腦與網路，輕鬆控制生活周遭的事物──不必事事都得委託聯發科設計電子電路、委託台積電製造電子電路。
然而深度學習目前仍是煉金術。大家憑感覺發明許多結構，但是無法利用數學來論證其功能。
metaheuristics也是類似的新流派。既有的最佳化演算法無法投入實用，大家重新發明基因演算法、蟻群演算法、粒子群演算法。
machine learning也是類似的新流派，既有的模式識別演算法無法投入實用，大家重新發明梯度下降技巧、神經網路結構。
陽春白雪變成下里巴人，引經據典變成信口開河。但是它實用。