状態フィードバック制御による安定化

フィードバック制御といわれるとすぐに思いつくのはPIDに代表される出力フィードバック制御である.
しかしプラントが高次となると, 閉ループ系全体の次数が上がり, 解析するのが困難になる. また, 制御器そのものが複雑になる. そこで, プラントの内部状態をフィードバックさせることで, 見かけ上の出力フィードバックを行う. それが状態フィードバック制御である.
ここでは不安定な系を示し, 実際にPythonを用いてシミュレーションを行う.
そしてそこから, 状態フィードバックによる安定化を行う.
そして, その状態フィードバックを用いて, 安定化できていることをシミュレーションで示す.

不安定な系
状態フィードバックによる安定化
参考文献

不安定な系

ここでは参考文献[1]に書かれた不安定な系を示す. 次の通りである.

\begin{gather}
\begin{bmatrix}
\dot{x}_1\\
\dot{x}_2
\end{bmatrix}=
\begin{bmatrix}
0&1\\
0.5&0.5
\end{bmatrix}
\begin{bmatrix}
x_1\\
x_2
\end{bmatrix}
+\begin{bmatrix}
0\\
1
\end{bmatrix}
u\\
y = \begin{bmatrix}
1 &0
\end{bmatrix}
\begin{bmatrix}
x_1\\
x_2
\end{bmatrix}
\end{gather}
この系の伝達関数は$P(s) = \frac{1}{(s+0.5)(s-1)}$で不安定となる.
この系に, 入力１のステップ関数を入力した際のシミュレーション結果を, 図1に示す.

図1を見ると大きく発散していくことが確認できる.
コードを次に示す通りである.

import numpy as np
import matplotlib.pyplot as plt

A = np.array([[0,1],[0.5,0.5]])
B = np.array([[0],[1]])
C  = np.array([[1,0]])

x0 = np.array([[0], [0]])
Tsim = 30
sampleT = 0.001
Tn = int(round(Tsim / sampleT))

u = np.ones(Tn + 1)
x_s = np.zeros((2, Tn + 1))
x_s[:, 0] = x0[:, 0]
t_s = np.zeros(Tn )
t_o = np.zeros(Tn  )
x_o = np.zeros(Tn  )

for i in range(Tn):
    t_s[i] = i* sampleT
    # ルンゲクッタ法
    k1 = A @ x_s[:, i] + B[:, 0] * u[i]
    k2 = A @ (x_s[:, i] + sampleT / 2 * k1) + B[:, 0] * u[i]
    k3 = A @ (x_s[:, i] + sampleT / 2 * k2) + B[:, 0] * u[i]
    k4 = A @ (x_s[:, i] + sampleT * k3) + B[:, 0] * u[i]
    x_s[:, i + 1] = x_s[:, i] + sampleT / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
    t_o[i] = t_s[i]
    x_o[i] = C @ x_s[:, i]

# Figure
fss = 20

FONTSIZE=12
fig,axes = plt.subplots(1,1,figsize=(10,10))
axes.plot(t_o,x_o)
axes.grid()
axes.set_xlabel('time/s',fontsize=FONTSIZE)
axes.set_ylabel('x/m',fontsize=FONTSIZE)
plt.show()

状態フィードバックによる安定化

図2のようなフィードバックゲインを設計して安定化を行う.
図2の$x_2$直後の1は、C行列を表し, 出力フィードバックと区別するために書いている.

フィードバックゲインと入力および, 状態変数の関係は次のようになる.
\begin{gather}
u = v-\begin{bmatrix} K_1&K_2\end{bmatrix}
\begin{bmatrix}
x_1\\
x_2
\end{bmatrix}
\end{gather}
このとき, 閉ループ系の状態方程式は
\begin{gather}
\begin{bmatrix}
\dot{x}_1\\
\dot{x}_2
\end{bmatrix}=
\begin{bmatrix}
0&1\\
0.5&0.5
\end{bmatrix}
\begin{bmatrix}
x_1\\
x_2
\end{bmatrix}+
\begin{bmatrix}
0\\
1
\end{bmatrix}
(v-\begin{bmatrix}
K_1&K_2
\end{bmatrix}
\begin{bmatrix}
x_1\\
x_2
\end{bmatrix})\\
=\begin{bmatrix}
0&1\\
0.5-K_1&0.5-K_2
\end{bmatrix}
\begin{bmatrix}
x_1\\
x_2
\end{bmatrix}+
\begin{bmatrix}
0\\1
\end{bmatrix}v
\end{gather}

となる. この特性多項式は
\begin{gather}
\det(s\boldsymbol{I}-\begin{bmatrix}
0&1\\
0.5-K_1&0.5-K_2
\end{bmatrix})\\=
s^2+(K_2-0.5)+(K_1-0.5)\tag{1}
\end{gather}

ここでフィードバックゲイン$K_1, K_2$を適切に選ぶことにより, 任意の極を設定することができる.
例えば, 極を$-1\pm j$とすると,
\begin{gather}
(s+1-j)(s+1+j)=s^2+2s+2 \tag{2}
\end{gather}
(2)式と(1)式を計数比較すると$K_1=2.5$, $K_2=2.5$と求まる

状態フィードバックにより安定化したシミュレーション結果を図3に示す.

図3より, 安定していることが確認できる.
コードを次に示す. (もし何か間違いがあったらご指摘いただきたい.)

import numpy as np
import matplotlib.pyplot as plt

A = np.array([[0,1],[0.5,0.5]])
B = np.array([[0],[1]])
C  = np.array([[1,0]])

k1 = 2.5
k2 = 2.5

A = np.array([[0,1],[0.5-k1,0.5-k2]])

x0 = np.array([[0], [0]])
Tsim = 30
sampleT = 0.001
Tn = int(round(Tsim / sampleT))

u = np.ones(Tn + 1)
x_s = np.zeros((2, Tn + 1))
x_s[:, 0] = x0[:, 0]
t_s = np.zeros(Tn )
t_o = np.zeros(Tn  )
x_o = np.zeros(Tn  )

for i in range(Tn):
    t_s[i] = i* sampleT
    # ルンゲクッタ法
    k1 = A @ x_s[:, i] + B[:, 0] * u[i]
    k2 = A @ (x_s[:, i] + sampleT / 2 * k1) + B[:, 0] * u[i]
    k3 = A @ (x_s[:, i] + sampleT / 2 * k2) + B[:, 0] * u[i]
    k4 = A @ (x_s[:, i] + sampleT * k3) + B[:, 0] * u[i]
    x_s[:, i + 1] = x_s[:, i] + sampleT / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
    t_o[i] = t_s[i]
    x_o[i] = C @ x_s[:, i]

# Figure
fss = 20

FONTSIZE=12
fig,axes = plt.subplots(1,1,figsize=(10,10))
axes.plot(t_o,x_o)
axes.grid()
axes.set_xlabel('time/s',fontsize=FONTSIZE)
axes.set_ylabel('x/m',fontsize=FONTSIZE)
plt.show()