前回 : はじめての現代制御理論講義13 サーボ系の設計 - ma-hi-ro Hobby

講義14 最適制御

現代制御のひと段落が「最適制御」

講義14 最適制御

評価関数と最適制御則

ここまでの状態フィードバックでは、極配置を行うことで、制御特性（収束具合）を調整できた
しかし、現実問題では、フィードバックゲインを大きくしすぎると、
アクチュエータがその出力を出せなかったり、負荷がかかったりで、実現できない

最適制御では以下の評価関数 $J$ を設定して、その値が最小となる制御を実現する

$\displaystyle \begin{aligned} J &= \int_0^\infty \lbrace x^T(t) Q x(t) + r u^2(t) \rbrace dt \\ &= \int_0^\infty \lbrace x^T(t) \frac{Q}{r} x(t) + u^2(t) \rbrace dt \end{aligned}$

ここで $Q$ は対角行列、 $r$ はスカラーとして、設定するパラメータになる
第１項が、計測値（状態量）の大きさ（2乗）を表し、
第２項が、出力（操作量）の大きさを表す
上記２つのパラメータの重みを割り振ることで、どちらを優先するかを設定できる

最適制御問題

上記の $J$ を最小とする制御パラメータ $u(t)$ を求める

まずは、次のリカッチ代数方程式から行列 $P$ を求める

$\displaystyle A^T P + P A - P b r^{-1} b^T P + Q = 0 \\$

また [tex: P = P^T > 0] を満たす正定行列である

求まった $P$ から、 $u(t)$ を下記で求める

$\displaystyle u(t) = - r^{-1} b^T P x(t)$

これで、最適制御が実現できる

ちなみに、この時に評価関数は

$\displaystyle J_{min} = x^T(0) P x(0)$

になる

Pythonによる数値シミュレーション

2x2の以下のシステムでの最適制御を考える

$\displaystyle \begin{aligned} \frac{d}{dt} \begin{bmatrix} x_1 (t) \\ x_2 (t) \end{bmatrix} &= \begin{bmatrix} 0 & 1 \\ -6 & -5 \end{bmatrix} \begin{bmatrix} x_1 (t) \\ x_2 (t) \end{bmatrix} - \begin{bmatrix} 0 \\ 1 \end{bmatrix} u(t) \\ Q &= \begin{bmatrix} 13 & 0 \\ 0 & 1 \end{bmatrix} \\ r &= 1 \end{aligned}$

この時、

$\displaystyle P = \begin{bmatrix} p_1 & p_2 \\ p_2 & p_3 \end{bmatrix}$

とすると、リカッチ代数方程式は、

$\displaystyle -p_2^2 - 12 p_2 + 13 = 0 \\ p_1 - 5 p_2 - 6 p_3 - p_2 p_3 = 0 \\ -p_3^2 + 2 p_2 - 10 p_3 + 1 = 0$

となり、結局、

$\displaystyle P = \begin{bmatrix} -30 + 14 \sqrt{7} & 1 \\ 1 & -5 + 2 \sqrt{7} \end{bmatrix}$

となる

フィードバックは、

$\displaystyle u(t) = - \lbrack 1 \quad -5 + 2 \sqrt{7} \rbrack x(t) \quad = - \lbrack 1 \quad 0.2915 \rbrack x(t)$

となる

t = [0.0]
x1_arr = [-1.0]
x2_arr = [ 0.0]
u_arr  = []

dt = 0.001
for i in range(10000):
    u_arr  += [-(1 * x1_arr[-1] + 0.2915 * x2_arr[-1])]
    x1_arr += [x1_arr[-1] + (0  * x1_arr[-1] + 1  * x2_arr[-1] + 0.0 * u_arr[-1]) * dt]
    x2_arr += [x2_arr[-1] + (-6 * x1_arr[-1] - 5  * x2_arr[-1] + 1.0 * u_arr[-1]) * dt]
    t      += [t[-1] + dt]

極配置法との比較

極を-3の重根にした場合、 $- \lbrack 3 \quad 1 \rbrack x(t)$ となる

t = [0.0]
x1_arr2 = [-1.0]
x2_arr2 = [ 0.0]
u_arr2  = []

dt = 0.001
for i in range(10000):
    u_arr2  += [-(3 * x1_arr2[-1] + 1 * x2_arr2[-1])]
    x1_arr2 += [x1_arr2[-1] + (0  * x1_arr2[-1] + 1  * x2_arr2[-1] + 0.0 * u_arr[-1]) * dt]
    x2_arr2 += [x2_arr2[-1] + (-6 * x1_arr2[-1] - 5  * x2_arr2[-1] + 1.0 * u_arr[-1]) * dt]
    t      += [t[-1] + dt]

fig = plt.figure(figsize=(12, 4))
ax1 = fig.add_subplot(1, 3, 1)
ax2 = fig.add_subplot(1, 3, 2, sharex=ax1)
ax3 = fig.add_subplot(1, 3, 3, sharex=ax1)

ax1.plot(t, x1_arr,  '-',   label='LQ')
ax1.plot(t, x1_arr2, '--',  label='pole')
ax2.plot(t, x2_arr,  '-',   label='LQ')
ax2.plot(t, x2_arr2, '--',  label='pole')
ax3.plot(t[1:], u_arr,  '-',   label='LQ')
ax3.plot(t[1:], u_arr2, '--',  label='pole')
ax1.grid()
ax1.legend()
ax1.set_ylabel('$x_1(t)$')
ax1.set_ylim(-1.2, 2.2)
ax2.grid()
ax2.legend()
ax2.set_ylabel('$x_2(t)$')
ax2.set_ylim(-1.2, 2.2)
ax3.grid()
ax3.legend()
ax3.set_ylabel('$u(t)$')
ax3.set_ylim(-0.1, 3.2)
ax1.set_xlim(0, 5)
fig.tight_layout()

最適制御では、出力 $u(t)$ が小さくできている

さらに出力を小さくする

パラメータを下記にする

$\displaystyle \begin{aligned} Q &= \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \\ r &= 3 \end{aligned}$

フィードバックゲインは $u(t) = - \lbrack 0.0277 \quad 0.0387 \rbrack x(t)$ になる

t = [0.0]
x1_arr3 = [-1.0]
x2_arr3 = [ 0.0]
u_arr3  = []

dt = 0.001
for i in range(10000):
    u_arr3  += [-(0.0277 * x1_arr3[-1] + 0.0387 * x2_arr3[-1])]
    x1_arr3 += [x1_arr3[-1] + (0  * x1_arr3[-1] + 1  * x2_arr3[-1] + 0.0 * u_arr3[-1]) * dt]
    x2_arr3 += [x2_arr3[-1] + (-6 * x1_arr3[-1] - 5  * x2_arr3[-1] + 1.0 * u_arr3[-1]) * dt]
    t      += [t[-1] + dt]

fig = plt.figure(figsize=(5, 4))
ax1 = fig.add_subplot(1, 1, 1)

ax1.plot(t, x1_arr, '-',  label='$x_1$')
ax1.plot(t, x2_arr, '--', label='$x_2$')
ax1.grid()
ax1.legend()
ax1.set_xlim(0, 5)
ax1.set_ylim(-1.2, 3.2)

fig.tight_layout()