system🚧

system

「系統」。多變量函數,輸入多道訊號,輸出多道訊號。

已經發展成熟,累積大量懸案,最近三十年沒有太大變化。

電力、機械相關科系基礎知識。暫且先做簡介,有空再來展開。

生成器:系統沒有輸入訊號。輸入一些參數,輸出一道訊號。
濾波器:系統輸出之後,追加一個系統,用來調整輸出訊號。
控制器:系統輸入之前,追加一個系統,用來調整輸入訊號暨輸出訊號。
觀測器:系統輸出之後,追加一個系統,用來觀測系統狀態。

課程關係圖

 多變量微積分  微分方程  動態系統  訊號與系統   聲學
        ↘       ↘        ↘          ↘       ↘
力學 → 電磁學 → 電路學 → 化學動力學 → 自動控制 → 電聲學

generator🚧

generator

「生成器」。函數,輸入參數,輸出訊號。

filter🚧

filtering

「濾波」。函數,輸入訊號,輸出訊號。

system

finite impulse response / infinite impulse response

系統有兩種款式。
有限脈衝響應:沒迴圈=>輸入的加權總和=>輸入只取幾項=>有限
無限脈衝響應:有迴圈=>輸入與輸出的加權總和=>展開迴圈,變成取所有項=>無限

moving average model / autoregressive model

系統有三種款式。
移動平均數模型:輸入訊號的加權總和=輸出訊號 (all zero) 
自迴歸模型:輸出訊號的加權總和=輸出訊號 (all pole)
兩種都用:輸入訊號的加權總和=輸出訊號的加權總和

z-transform / convolution theorem / zero / pole

z-transform:訊號=>多項式
convolution theorem:輸入訊號加權總和=>多項式乘法
convolution theorem:輸出訊號加權總和=>多項式除法
zero:分子多項式的根(多項式變成0)。
pole:分母多項式的根(多項式變成±∞)。
大家用示波器看zero pole,然後反向推理系統是什麼。

Fourier transform / Laplace transform

離散傅立葉轉換: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。

很多人將拉普拉斯轉換解讀成z轉換的連續版本,但是這種解讀不太對。
嚴格來說,z轉換是普通生成函數,拉普拉斯轉換是指數生成函數。性質截然不同。

傅立葉轉換是特例,拉普拉斯轉換是通例,導致教科書很喜歡用拉普拉斯轉換。
但是現實世界無法取得數值、無法拿來計算啦。
實務上都是使用傅立葉轉換。

time domain / frequency domain

因為真實世界的訊號幾乎都是一堆波,
所以大家用傅立葉轉換,把訊號分解成波。
傅立葉轉換:振幅為1、頻率為自然數/各種數值。
拉普拉斯轉換:振幅為各種數值、頻率為各種數值。
原本訊號稱作時域(座標軸是時間),傅立葉轉換之後稱作頻域(座標軸是頻率)。

spectrum / Bode plot / magnitude spectrum / phase spectrum

訊號實施傅立葉轉換(時域轉頻域)。
一串數列的傅立葉轉換是一串數列,每個數值都是複數。
複數長度是強度。頻譜每個數值的強度,形成強度頻譜。
複數角度是相位。頻譜每個數值的相位,形成相位頻譜。
兩者合稱頻譜。
強度頻譜橫軸與縱軸都取log,相位頻譜橫軸取log,兩者合稱波特圖。

frequency response / transfer function / gain / phase shift

訊號實施傅立葉轉換(時域轉頻域)。
一串數列的傅立葉轉換是一串數列,
每個數值代表每種頻率(frequency)的波(wave)的強度(magnitude)和相位(phase)。

值得一提的是,除了訊號可以分解成波,其實系統也可以分解成波。
系統弄成加權總和的形式,然後把權重當作訊號,拿去分解。
(因為z-transform和convolution theorem)

兩種方式可以得到系統的頻譜:
一、系統的頻譜(權重做傅立葉轉換)。
二、輸出訊號頻譜除以輸入訊號頻譜(訊號做傅立葉轉換,再相除)。

頻率響應:系統的頻譜的數值,針對特定頻率。
傳遞函數:系統的頻譜的數學式,變數是頻率。
增益:系統的強度頻譜。各種頻率的波的強度變化倍率。
相位偏移:系統的相位頻譜。各種頻率的波的相位變化差距。

filter

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。
曲線越陡價格越貴。

controller🚧

control

「控制」。一個系統,一道輸入訊號,一道輸出訊號。調整輸入訊號,得到特別的輸出訊號。

control應用十分廣泛,是世上最實用的演算法之一。

https://www.mathworks.com/solutions/control-systems/feedback-control-systems.html
https://ctms.engin.umich.edu/CTMS/index.php

system

causal system / linear system / time-invariant system

因果系統:遞迴公式的變數都是過去時刻。
線性系統:遞迴公式由變數的加法和倍率組成(加權總和)。
     也就是說,輸入相加導致輸出相加、輸入倍率導致輸出倍率。
非時變系統:遞迴公式不因時刻而變。
      也就是說,輸入移位導致輸出移位。
              中譯 意譯
cause (noun)  原因 原因             
cause (verb)  導致 因此而           
cause (conj)  因為 because的精簡講法,原因是
causal        因果 原因的           

causal LTI system / BIBO stability / Routh array

大家習慣省略字首causal。
LTI system可以寫成微分方程式(隱式)、寫成卷積(顯式)。
LTI system = constant-coefficient linear ODE = convolution
假設系統是線性非時變系統(例如物理運動、電路運作)。
一、最初是微分方程式。
二、改寫成transfer function(時域轉頻域)。形成多項式分式。
三、改寫成因式分解。分母的根pole、分子的根zero。
四、改寫成部分分式。形成分式連加。
五、改寫成符號解(頻域轉時域)。形成exp()連加。

輸入受限則輸出受限:可以改寫成訊號L¹-norm受限。
LTI system的情況下:所有pole都在左半複平面。
訊號學家還希望收斂至零:所有pole都在左半複平面,但不含虛軸。
滿足穩定性,才能形成穩態,才能控制。
不滿足穩定性,輸出可能正負無限大。
導致電路過載燒掉、動力機械暴衝、反應槽爆炸。

判斷穩定性,
直覺的方法是使用多項式函數求根演算法,求得所有pole。
經典演算法是companion matrix求特徵值,時間複雜度O(N³T)。
特殊的方法是使用特殊數學公式,檢查正負號。
經典演算法是Routh array,時間複雜度O(N²)。
https://tutorial.math.lamar.edu/classes/de/IVPWithLaplace.aspx
https://lpsa.swarthmore.edu/LaplaceXform/InvLaplace/InvLaplaceXformPFE.html
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

dynamic response

已知輸入、系統,求得輸出。
輸入:已知函數(教科書習慣討論下述三種)
系統:已知函數(教科書習慣討論一次微分方程、二次微分方程)
輸出:未知函數。

1. impulse response:脈衝波。隔壁棚結構分析很常用,控制系統則不使用。
2. step response:單位步進函數。主角。例如啟動馬達至定速。
3. frequency response:複弦波。得到波特圖其中一個數值。

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

controller

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
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。極X零O。
畫出poles和zeros的位置,以便判斷穩定性。
如果極X都在左半複平面(開區間、不含虛軸),則穩定。
如果右半複平面沒有極X(閉區間、包含虛軸),則穩定。

transfer function的poles,就是分母的zeros。
畫出分母的極零圖,亦可判斷穩定性。
如果右半複平面沒有零O,則穩定。
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
lag compensator   ≈ low-pass filter
notch compensator ≈ band-pass filter

lead compensator 左X右O  K(s) = kp (Ts+1) / (aTs+1)  a<1
lag compensator  左O右X  K(s) = kp (Ts+1) / (aTs+1)  a>1

non-minimum phase system

右半複平面出現zeros。
當參考訊號是步進函數,則輸出訊號是先蹲後跳、聯結車轉彎。一開始衝向負值。

沒救了。compensator沒有辦法解決這種情況。
一種直覺的方式是追加poles抵銷zeros,分母分子約分之後一起消失不見。
然而實務上無法完全對準。
誤差、設備老化,都會導致zeros偏移。
稍有差池,輸出訊號就會偶然出現正負無限大。
導致電路過載燒掉、動力機械暴衝、反應槽爆炸。
非常危險。
實務上不能追加右半複平面poles抵銷右半複平面zeros。
我不知道有沒有其他解法。也許根本不需要解,順其自然就好。

nonlinear system

sliding mode controller

nonlinear

model predictive control

optimization
state space

Hamilton–Jacobi–Bellman equation

動態規劃,找到最好的動作方式。

Sylvester equation / Lyapunov equation

矩陣方程式,常見的最佳化問題的解。

fuzzy logic

把布林數的AND和OR運算,變成函數的min和max運算。
Karnik–Mendel algorithm

observer🚧

observation

「觀測」。已知系統、輸入訊號、輸出訊號,求得系統狀態。

state transition
state estimation (state observer)
一階微分 = 取前一個時刻
二階微分 = 取前兩個時刻

markov process就是引入機率的連續動態系統
markov chain就是引入機率的離散動態系統

system

stochastic process

訊號的數值,從固定的改成浮動的。甚至前後項有某種特殊性關係。

system identification / system realization / system estimation

給你輸入訊號、輸出訊號,請你找到系統。
換句話說,就是迴歸!

mean squared error / Kullback–Leibler divergence

兩道訊號的距離。
換句話說,就是誤差!迴歸的時候拿來用吧。

observer

找到迴歸函數之後,就可以預測訊號啦。
常見方式有Kalman filter和hidden Markov model。
http://www.princeton.edu/~stengel/MAE546Seminars.html

Kalman filter

未知因:autoregression model (LTI filter)
已知業:LTI filter             
已知果:一道訊號

請見本站文件「filter」。
未知輸入(因)已知函數(業)已知輸出(果),求因。
精髓:果的誤差,通過反濾波器,用以修正因。

hidden Markov model (Bayes filter)

未知因:Markov chain
未知業:hidden Markov model           
已知果:許多道訊號

請見本站文件「hidden Markov model」。
未知輸入(因)未知函數(業)已知輸出(果),求因。
精髓:已知果,梯度下降法(反向傳播法)找最大值,得到機率最大的業。
   未知果,動態規劃找最佳路線,得到機率最大的因。

state observer

有點不切實際
https://control.asu.edu/Classes/MAE507/507Lecture09.pdf

system design🚧

system design

「系統設計」。系統是函數、甚至是函數網路。工程師設計函數網路,達成特定任務,例如生成/濾波/控制/觀測。

我目前只知道三種流派:

一、電子電路:風靡全世界。雖然台灣是地球上最大的電子零件生產基地,台灣也有專門設計電子電路的公司,但是我不太確定台灣是否有這方面的專家。有言道:十萬青年十萬肝,GG輪班救台灣。關鍵字:電路學、積體電路設計、自動控制、控制理論。

二、生化反應:發展中。雖然台灣之前打算成立藥物代工廠、學名藥設計實驗室,但是以失敗告終。有言道:一日生科,終生科科。關鍵字:化學動力學、酵素動力學、藥物動力學、系統生物學。

三、深度學習:最近十年才剛萌芽。已經做到生成/濾波,目前做不到控制/觀測。正在經歷大浪潮,準備迎接大泡沫。有言道:站在風口上,豬都會飛。關鍵字:人工智慧、機器學習、深度學習。

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

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

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也是類似的新流派,既有的模式識別演算法無法投入實用,大家重新發明梯度下降技巧、神經網路結構。

陽春白雪變成下里巴人,引經據典變成信口開河。但是它實用。