二次型最優控制(一)

本文來自潘神的知乎筆記什麼是二次型最優控制

什麼是二次型

說到“最”字，我們自然需要一個标準。比如，我們說到“最高”，可以用長度來丈量；說到“最重”，一般用品質來衡量，這些都是比較容易獲得标準的。現在如果問，什麼是“最美”？我們可能就得思考一會了。同樣的道理，現在如果問什麼是“最優”控制？——在沒有“标準”的情況下讨論這個問題，就沒有太大的意義。是以，我們要先來搞一個标準。既然我們搞控制，簡單來說就是要“指東打東，指西打西”，翻譯一下就是“輸出能跟随輸入”，那什麼是“最優控制”呢？——一個很直覺的想法就是：“輸入能完全跟随輸入，在每一個時間點上都是100%一樣”，當然這是理想情況，實際上很難實作。那我們能不能拉開來看呢？——我們要有更大的視野，不要求每一點都一樣，那在整個工作時間内是不是可以更接近呢？比如我們可以這樣：把每一時刻的狀态的誤差都累加起來，隻要累加的值更小，那是不是就更接近呢？比如我們可以令：

J = ∫ 0 ∞ ∣ e ∣ d t = ∫ 0 ∞ ∣ x o ( t ) − x i ( t ) ∣ d t J=\int_{0}^{\infty}|e| d t=\int_{0}^{\infty}\left|x_{o}(t)-x_{i}(t)\right| d t J=∫0∞​∣e∣dt=∫0∞​∣xo​(t)−xi​(t)∣dt

其中 x o ( t ) x_{o}(t) xo​(t)為實際狀态， x i ( t ) x_{i}(t) xi​(t)為期望狀态。這樣顯然是可以的，如果 x o ( t ) x_{o}(t) xo​(t)和 x i ( t ) x_{i}(t) xi​(t)完全一樣， J J J則為零，達到最小值，顯然 J J J 越小越好——我們一般将 J J J稱之為代價函數（cost function）。

但是我們知道，絕對值這個操作，數學上不太好處理——因為它的導數不連續，會帶來很多麻煩，那怎麼辦？——聰明的你可能已經想到了，可以用平方啊，它也能表征距離，而且性能非常好。 這樣代價函數就變成了：

J = ∫ 0 ∞ e 2 d t = ∫ 0 ∞ ( x o ( t ) − x i ( t ) ) 2 d t J=\int_{0}^{\infty} e^{2} d t=\int_{0}^{\infty}\left(x_{o}(t)-x_{i}(t)\right)^{2} d t J=∫0∞​e2dt=∫0∞​(xo​(t)−xi​(t))2dt

為簡單起見，我們先假設期望狀态 x i ( t ) x_{i}(t) xi​(t)為零，即

J = ∫ 0 ∞ e 2 d t = ∫ 0 ∞ ( x o ( t ) ) 2 d t J=\int_{0}^{\infty} e^{2} d t=\int_{0}^{\infty}\left(x_{o}(t)\right)^{2} d t J=∫0∞​e2dt=∫0∞​(xo​(t))2dt

我們看一下被積分的部分，這隻有一個狀态，我們可以擴充一下，假如現在有 n n n個狀态，即

J = ∫ 0 ∞ x 1 2 ( t ) + x 2 2 ( t ) + … + x n 2 ( t ) d t J=\int_{0}^{\infty} x_{1}^{2}(t)+x_{2}^{2}(t)+\ldots+x_{n}^{2}(t) d t J=∫0∞​x12​(t)+x22​(t)+…+xn2​(t)dt

這樣不夠緊湊，我們換個寫法：

J = ∫ 0 ∞ x T x d t J=\int_{0}^{\infty} \boldsymbol{x}^{T} \boldsymbol{x} d t J=∫0∞​xTxdt

其中

x T x = [ x 1 , x 2 , … x n ] [ 1 0 ⋯ 0 0 1 … 0 ⋮ ⋮ ⋮ ⋮ 0 0 … 1 ] [ x 1 x 2 ⋮ x n ] = x 1 2 + x 2 2 + … + x n 2 \begin{array}{c}{\boldsymbol{x}^{T} \boldsymbol{x}=\left[x_{1}, x_{2}, \ldots x_{n}\right] \left[ \begin{array}{cccc}{1} & {0} & {\cdots} & {0} \\ {0} & {1} & {\ldots} & {0} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ {0} & {0} & {\ldots} & {1}\end{array}\right] \left[ \begin{array}{c}{x_{1}} \\ {x_{2}} \\ {\vdots} \\ {x_{n}}\end{array}\right]} \\ {=x_{1}^{2}+x_{2}^{2}+\ldots+x_{n}^{2}}\end{array} xTx=[x1​,x2​,…xn​]⎣⎢⎢⎢⎡​10⋮0​01⋮0​⋯…⋮…​00⋮1​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​x1​x2​⋮xn​​⎦⎥⎥⎥⎤​=x12​+x22​+…+xn2​​

這就是一種簡單的二次型函數，因為變量的最高次數是2。

當然，我們的工具可以更強大，比如機關矩陣 I I I可以擴充為對角矩陣 λ \lambda λ：

λ = [ λ 1 0 … 0 0 λ 2 … 0 ⋮ ⋮ ⋮ ⋮ 0 0 … λ n ] \boldsymbol{\lambda}=\left[ \begin{array}{cccc}{\lambda_{1}} & {0} & {\ldots} & {0} \\ {0} & {\lambda_{2}} & {\dots} & {0} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ {0} & {0} & {\ldots} & {\lambda_{n}}\end{array}\right] λ=⎣⎢⎢⎢⎡​λ1​0⋮0​0λ2​⋮0​……⋮…​00⋮λn​​⎦⎥⎥⎥⎤​

這相當于給各個變量誤差取個權重，然後對角矩陣 λ \lambda λ還可以擴充為更一般的矩陣 Q Q Q：

Q = [ a 11 a 12 … a 1 n a 21 a 22 … a 2 n ⋮ ⋮ ⋮ ⋮ a n 1 a n 2 … a n n ] \boldsymbol{Q}=\left[ \begin{array}{cccc}{a_{11}} & {a_{12}} & {\dots} & {a_{1 n}} \\ {a_{21}} & {a_{22}} & {\dots} & {a_{2 n}} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ {a_{n 1}} & {a_{n 2}} & {\dots} & {a_{n n}}\end{array}\right] Q=⎣⎢⎢⎢⎡​a11​a21​⋮an1​​a12​a22​⋮an2​​……⋮…​a1n​a2n​⋮ann​​⎦⎥⎥⎥⎤​

但無論怎麼擴充，展開的函數最高次數都是2次，是以我們把這種類型的函數統稱為二次型函數。

二階系統的阻尼比為何選擇0.707

一個簡單例子，二階系統，框圖如下，現在加入我們要實作最優控制，看看最佳阻尼比應該是什麼樣。

其開環傳遞函數為

C ( s ) E ( s ) = 1 s ( s + 2 ζ ) \frac{C(s)}{E(s)}=\frac{1}{s(s+2 \zeta)} E(s)C(s)​=s(s+2ζ)1​

誤差為：

E ( s ) = R ( s ) − C ( s ) E(s)=R(s)-C(s) E(s)=R(s)−C(s)

将上面兩個傳遞函數變成時域方程為：

c ¨ + 2 ζ c ˙ = e e = r − c \begin{array}{l}{\ddot{c}+2 \zeta \dot{c}=e} \\ {e=r-c}\end{array} c¨+2ζc˙=ee=r−c​

聯立以上兩式就可以得到

e ¨ + 2 ζ e ˙ + e = r ¨ + 2 ζ r ˙ \ddot{e}+2 \zeta \dot{e}+e=\ddot{r}+2 \zeta \dot{r} e¨+2ζe˙+e=r¨+2ζr˙

傳遞函數隻能描述單輸入單輸出系統，現在我們要把它改成狀态方程的形式：

x 1 = e x 2 = e ˙ \begin{aligned} x_{1} &=e \\ x_{2} &=\dot{e} \end{aligned} x1​x2​​=e=e˙​

也就是我們要研究誤差和誤差變化率，為簡單起見，我們假設輸入為階躍信号，即 r ¨ = 0 \ddot{r}=0 r¨=0， r ˙ = 0 \dot{r}=0 r˙=0， r = 1 \quad r=1 r=1，對應的 e = 1 e=1 e=1， e ˙ = 0 \quad \dot{e}=0 e˙=0。這樣狀态方程就變成

x ˙ = A x \dot{\boldsymbol{x}}=\boldsymbol{A x} x˙=Ax

其中

A = [ 0 1 − 1 − 2 ζ ] A=\left[ \begin{array}{cc}{0} & {1} \\ {-1} & {-2 \zeta}\end{array}\right] A=[0−1​1−2ζ​]

完整狀态方程為

[ x ˙ 1 x ˙ 2 ] = [ 0 1 − 1 − 2 ζ ] [ x 1 x 2 ] \left[ \begin{array}{c}{\dot{x}_{1}} \\ {\dot{x}_{2}}\end{array}\right]=\left[ \begin{array}{cc}{0} & {1} \\ {-1} & {-2 \zeta}\end{array}\right] \left[ \begin{array}{l}{x_{1}} \\ {x_{2}}\end{array}\right] [x˙1​x˙2​​]=[0−1​1−2ζ​][x1​x2​​]

我們定義代價函數為：

J = ∫ 0 ∞ ( e 2 + e ˙ 2 ) d t J=\int_{0}^{\infty}\left(e^{2}+\dot{e}^{2}\right) d t J=∫0∞​(e2+e˙2)dt

即兩個變量之間同等權重，用變量表示

J = ∫ 0 ∞ ( x 1 2 + x 2 2 ) d t J=\int_{0}^{\infty}\left(x_{1}^{2}+x_{2}^{2}\right) d t J=∫0∞​(x12​+x22​)dt

現在我們想要 J J J最小，具體應該怎麼做呢？先将代價函數寫成緊湊一點的形式：

J = ∫ 0 ∞ x T x d t J=\int_{0}^{\infty} \boldsymbol{x}^{T} \boldsymbol{x} d t J=∫0∞​xTxdt

這是個積分的問題，是以我們要先找到被積函數 x T x \boldsymbol{x}^{T} \boldsymbol{x} xTx的原函數，這個數學問題不是我們的重點，假設我們現在已經知道原函數的長相為 x T P x \boldsymbol{x}^{T} \boldsymbol{P} \boldsymbol{x} xTPx，即

d d t ( x T P x ) = − x T x \frac{d}{d t}\left(\boldsymbol{x}^{T} \boldsymbol{P} \boldsymbol{x}\right)=-\boldsymbol{x}^{T} \boldsymbol{x} dtd​(xTPx)=−xTx

其中 P \boldsymbol{P} P是一個待定的實對稱矩陣，我們把它看成待定系數，現在要研究一下這系數需要滿足什麼條件，我們把上式的左邊微分展開：

d d t ( x T P x ) = x ˙ T P x + x T P x ˙ \frac{d}{d t}\left(\boldsymbol{x}^{T} \boldsymbol{P} \boldsymbol{x}\right)=\dot{\boldsymbol{x}}^{T} \boldsymbol{P} \boldsymbol{x}+\boldsymbol{x}^{T} \boldsymbol{P} \dot{\boldsymbol{x}} dtd​(xTPx)=x˙TPx+xTPx˙

把狀态方程 x ˙ = A x \dot{\boldsymbol{x}}=\boldsymbol{A x} x˙=Ax代入上式，于是我們就得到：

d d t ( x T P x ) = x T ( A T P + P A ) x \frac{d}{d t}\left(\boldsymbol{x}^{T} \boldsymbol{P} \boldsymbol{x}\right)=\boldsymbol{x}^{T}\left(\boldsymbol{A}^{T} \boldsymbol{P}+\boldsymbol{P} \boldsymbol{A}\right) \boldsymbol{x} dtd​(xTPx)=xT(ATP+PA)x

如果我們保證 A T P + P A = − I \boldsymbol{A}^{T} \boldsymbol{P}+\boldsymbol{P} \boldsymbol{A}=-\boldsymbol{I} ATP+PA=−I，由于 A A A 是已知的，這樣我們就能确定 P P P，也就能确定原函數 x T P x \boldsymbol{x}^{T} \boldsymbol{P} \boldsymbol{x} xTPx ，此時

J = ∫ 0 ∞ x T x d t = − x T P x ∣ 0 ∞ = x T ( 0 ) P x ( 0 ) − x T ( ∞ ) P x ( ∞ ) J=\int_{0}^{\infty} \boldsymbol{x}^{T} \boldsymbol{x} d t=-\boldsymbol{x}^{T} \boldsymbol{P}\left.\boldsymbol{x}\right|_{0} ^{\infty}=\boldsymbol{x}^{T}(0) \boldsymbol{P} \boldsymbol{x}(0)-\boldsymbol{x}^{T}(\infty) \boldsymbol{P} \boldsymbol{x}(\infty) J=∫0∞​xTxdt=−xTPx∣0∞​=xT(0)Px(0)−xT(∞)Px(∞)

我們的系統是穩定的，而且輸入為零，這樣我們就可以得到 x ∞ = 0 \boldsymbol{x}_{\infty}=0 x∞​=0 ，是以代價函數可以進一步簡化為：

J = x T ( 0 ) P x ( 0 ) J=\boldsymbol{x}^{T}(0) \boldsymbol{P} \boldsymbol{x}(0) J=xT(0)Px(0)

其中 P P P滿足 A T P + P A = − I \boldsymbol{A}^{T} \boldsymbol{P}+\boldsymbol{P} \boldsymbol{A}=-\boldsymbol{I} ATP+PA=−I(這其實也是李雅普諾夫第二法)

注意：所有的二次型最優控制，形式可能會很複雜，但是處理思路和方式都是這樣。

P P P使用待定系數法求解。得到：

P = [ ζ + 1 2 ζ 1 2 1 2 1 2 ζ ] \boldsymbol{P}=\left[ \begin{array}{cc}{\zeta+\frac{1}{2 \zeta}} & {\frac{1}{2}} \\ {\frac{1}{2}} & {\frac{1}{2 \zeta}}\end{array}\right] P=[ζ+2ζ1​21​​21​2ζ1​​]

将 x ( 0 ) = [ 1 0 ] T \boldsymbol{x}(0)=\left[ \begin{array}{ll}{1} & {0}\end{array}\right]^{T} x(0)=[1​0​]T帶入 J = x T ( 0 ) P x ( 0 ) J=\boldsymbol{x}^{T}(0) \boldsymbol{P} \boldsymbol{x}(0) J=xT(0)Px(0)

得到代價函數：

J = ζ + 1 2 ζ J=\zeta+\frac{1}{2 \zeta} J=ζ+2ζ1​

現在我們要計算 J J J的極小值，先對變量微分

d J d ζ = 1 − 1 2 ζ 2 \frac{d J}{d \zeta}=1-\frac{1}{2 \zeta^{2}} dζdJ​=1−2ζ21​

當 d J d ζ = 0 \frac{d J}{d \zeta}=0 dζdJ​=0時，可以計算得到：

ζ = 2 2 ≈ 0.707 \zeta=\frac{\sqrt{2}}{2} \approx 0.707 ζ=22 ​​≈0.707

也就是說，當阻尼比取0.707時， J = ∫ 0 ∞ ( e 2 + e ˙ 2 ) d t J=\int_{0}^{\infty}\left(e^{2}+\dot{e}^{2}\right) d t J=∫0∞​(e2+e˙2)dt這個損失函數取極小值——也就是說從整體來看，誤差和誤差變化率累計最小，這就是為什麼我們一般把二階系統的阻尼比設定在0.707左右的原因。

代價函數矩陣需要由 I I I變成 Q Q Q

前面選擇的代價函數為

J = ∫ 0 ∞ ( x 1 2 + x 2 2 ) d t J=\int_{0}^{\infty}\left(x_{1}^{2}+x_{2}^{2}\right) d t J=∫0∞​(x12​+x22​)dt

我們認為狀态變量的權重都一樣，沒有誰更重要。實際上，并不都是如此，比如還是前面的例子，有時候我們可能更關心誤差，對誤差的變化率容忍度高一些，這時候我們就可以給它們配置設定一個權重：

J = ∫ 0 ∞ ( x 1 2 + μ x 2 2 ) d t μ > 0 J=\int_{0}^{\infty}\left(x_{1}^{2}+\mu x_{2}^{2}\right) d t \quad \mu>0 J=∫0∞​(x12​+μx22​)dtμ>0

此時代價函數可以演變為

J = ∫ 0 ∞ ( e 2 + μ e ˙ 2 ) d t = ∫ 0 ∞ ( x 1 2 + μ x 2 2 ) d t = ∫ 0 ∞ [ x 1 x 2 ] [ 1 0 0 μ ] [ x 1 x 2 ] d t = ∫ 0 ∞ x T Q x d t \begin{aligned} J &=\int_{0}^{\infty}\left(e^{2}+\mu \dot{e}^{2}\right) d t=\int_{0}^{\infty}\left(x_{1}^{2}+\mu x_{2}^{2}\right) d t \\ &=\int_{0}^{\infty}\left[x_{1} \quad x_{2}\right] \left[ \begin{array}{cc}{1} & {0} \\ {0} & {\mu}\end{array}\right] \left[ \begin{array}{l}{x_{1}} \\ {x_{2}}\end{array}\right] d t \\ &=\int_{0}^{\infty} \boldsymbol{x}^{T} \boldsymbol{Q} \boldsymbol{x} d t \end{aligned} J​=∫0∞​(e2+μe˙2)dt=∫0∞​(x12​+μx22​)dt=∫0∞​[x1​x2​][10​0μ​][x1​x2​​]dt=∫0∞​xTQxdt​

其中

Q = [ 1 0 0 μ ] \boldsymbol{Q}=\left[ \begin{array}{ll}{1} & {0} \\ {0} & {\mu}\end{array}\right] Q=[10​0μ​]

同樣前面的處理方式一樣，代價函數的表達式仍然為

J = x ( 0 ) T P x ( 0 ) J=\boldsymbol{x}(0)^{T} \boldsymbol{P} \boldsymbol{x}(0) J=x(0)TPx(0)

不過此時 P P P需要滿足不再是 A T P + P A = − I \boldsymbol{A}^{T} \boldsymbol{P}+\boldsymbol{P} \boldsymbol{A}=-\boldsymbol{I} ATP+PA=−I，而是 A T P + P A = − Q \boldsymbol{A}^{T} \boldsymbol{P}+\boldsymbol{P} \boldsymbol{A}=-\boldsymbol{Q} ATP+PA=−Q

還是前面的狀态方程：

[ x 1 ˙ x 2 ˙ ] = [ 0 1 − 1 − 2 ζ ] ⎵ A [ x 1 x 2 ] \left[ \begin{array}{c}{\dot{x_{1}}} \\ {\dot{x_{2}}}\end{array}\right]=\underbrace{\left[ \begin{array}{cc}{0} & {1} \\ {-1} & {-2 \zeta}\end{array}\right]}_{A} \left[ \begin{array}{l}{x_{1}} \\ {x_{2}}\end{array}\right] [x1​˙​x2​˙​​]=A [0−1​1−2ζ​]​​[x1​x2​​]

此時可計算滿足 A T P + P A = − Q \boldsymbol{A}^{T} \boldsymbol{P}+\boldsymbol{P} \boldsymbol{A}=-\boldsymbol{Q} ATP+PA=−Q 的 P P P為

P = [ ζ + 1 + μ 4 ζ 1 2 1 2 1 + μ 4 ζ ] \boldsymbol{P}=\left[ \begin{array}{cc}{\zeta+\frac{1+\mu}{4 \zeta}} & {\frac{1}{2}} \\ {\frac{1}{2}} & {\frac{1+\mu}{4 \zeta}}\end{array}\right] P=[ζ+4ζ1+μ​21​​21​4ζ1+μ​​]

将初始條件 x ( 0 ) = [ 1 0 ] T \boldsymbol{x}(0)=\left[ \begin{array}{ll}{1} & {0}\end{array}\right]^{T} x(0)=[1​0​]T帶入 J = x ( 0 ) T P x ( 0 ) J=\boldsymbol{x}(0)^{T} \boldsymbol{P} \boldsymbol{x}(0) J=x(0)TPx(0)，即可得到：

J = [ 1 0 ] [ ζ + 1 + μ 4 ζ 1 2 1 2 1 + μ 4 ζ ] [ 1 0 ] = ζ + 1 + μ 4 ζ J=\left[ \begin{array}{cc}{1} & {0}\end{array}\right] \left[ \begin{array}{cc}{\zeta+\frac{1+\mu}{4 \zeta}} & {\frac{1}{2}} \\ {\frac{1}{2}} & {\frac{1+\mu}{4 \zeta}}\end{array}\right] \left[ \begin{array}{l}{1} \\ {0}\end{array}\right]=\zeta+\frac{1+\mu}{4 \zeta} J=[1​0​][ζ+4ζ1+μ​21​​21​4ζ1+μ​​][10​]=ζ+4ζ1+μ​

可見，此時代價函數取決于阻尼比 ζ \zeta ζ和輸入的權重 μ \mu μ ，為了計算 J J J的極值，我們同樣要計算 J J J 對 ζ \zeta ζ 的微分

∂ J ∂ ζ = 1 − 1 + μ 4 ζ 2 \frac{\partial J}{\partial \zeta}=1-\frac{1+\mu}{4 \zeta^{2}} ∂ζ∂J​=1−4ζ21+μ​

當 ∂ J ∂ ζ = 0 \frac{\partial J}{\partial \zeta}=0 ∂ζ∂J​=0時， ζ = 1 + μ 2 \zeta=\frac{\sqrt{1+\mu}}{2} ζ=21+μ ​​

顯然，當 μ = 1 \mu=1 μ=1時， ζ = 2 2 \zeta=\frac{\sqrt{2}}{2} ζ=22 ​​，這與我們前面的計算結果一緻。如果我們現在增加 μ \mu μ ，也就意味着誤差變化率 e ˙ \dot{e} e˙的權重更大，此時也就意味着我們需要增加阻尼比 ζ \zeta ζ。