見出し画像

実際に使える制御工学:寄り道編(3)「状態方程式。それは現代制御理論の入口。」

<この記事で伝えたい事> 

  • 状態方程式とは、伝達関数とは異なるシステム表現方法です
    状態方程式を使うと、伝達関数(古典制御)とは異なる数学的アプローチによる制御設計(現代制御)が可能です。




1. はじめに

自己紹介はこちらです。

この記事は、状態方程式の紹介・解説のために書いています。
応用編(4)で状態方程式に少し触れたので、その補足という意味合いもあります)


2. 状態方程式とは?

システムの表現方法の1つです。

行列とベクトルを駆使して、システムを1階の常微分方程式と出力方程式で表したものが、状態方程式です。
その基本形を示したものが下式です。

$$
\left \{
\begin{array}{l}
\dot{x}=Ax+Bu\text{ :1階常微分方程式} \\
y=Cx+Du\text{ :出力方程式}
\end{array}
\right.
$$

$$
\left(
\begin{array}{l}
x\text{:状態ベクトル}\\
u\text{:制御入力}\\
y\text{:制御出力}\\
A\text{:システム行列}\\
B\text{:入力行列}\\
C\text{:出力行列}\\
D\text{:伝達行列}
\end{array}
\right)
$$


3. 状態方程式の具体的な作り方

まずは、表現したいシステムを決めることから始まります。

今回は、基礎編(2)で取り上げたDCモータシステムの状態方程式を作ることにします。

DCモータ全体の動きを表す連立微分方程式は下式です。

$$
\left \{
\begin{array}{cl}
e_{(t)}&=L\cdot \cfrac{di_{(t)}}{dt}+R\cdot i_{(t)}+V_{e(t)} \\
\tau_{(t)}&=J\cdot \cfrac{d\omega_{(t)}}{dt} \\
\tau_{(t)}&=K_t\cdot i_{(t)} \\
V_{e(t)}&=K_e\cdot \omega_{(t)}
\end{array}
\right.
$$

制御入力を電圧$${e_{(t)}}$$、制御出力を角速度$${\omega_{(t)}}$$と選んで状態方程式を作ってみましょう。
1階常微分方程式を作る都合上、微分項を左辺に持っていきます。

$$
\left \{
\begin{array}{cl}
L\cdot \cfrac{di_{(t)}}{dt}&=-R\cdot i_{(t)}-V_{e(t)}-e_{(t)} \\
J\cdot \cfrac{d\omega_{(t)}}{dt}&=\tau_{(t)}
\end{array}
\right.
$$

$${\tau=K_t\cdot i}$$、$${V_e=K_e\cdot\omega}$$から、連立方程式を以下のように変形できます。

$$
\left \{
\begin{array}{cl}
L\cdot \cfrac{di_{(t)}}{dt}&=-R\cdot i_{(t)}-K_e\cdot \omega_{(t)}-e_{(t)} \\
J\cdot \cfrac{d\omega_{(t)}}{dt}&=K_t\cdot i_{(t)}
\end{array}
\right.
$$

そして、基本形どおり左辺にかかる係数がない形にします。

$$
\left \{
\begin{array}{cl}
\cfrac{di_{(t)}}{dt}&=-\cfrac{R}{L}\cdot i_{(t)}-\cfrac{K_e}{L}\cdot \omega_{(t)}-\cfrac{1}{L}e_{(t)} \\
\cfrac{d\omega_{(t)}}{dt}&=\cfrac{K_t}{J}\cdot i_{(t)}
\end{array}
\right.
$$

この形にできたら、連立方程式左辺の時間微分されている変数をベクトル化します。
(これが状態ベクトル$${\bold{x}}$$、状態ベクトルの構成要素を状態変数と呼びます)

$$
\bold{x}=\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}
\end{array}
\right]
$$

方程式中の状態変数の係数を行列形式でまとめます。

$$
\cfrac{d}{dt}\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}
\end{array}
\right]
=
\left[
\begin{array}{cc}
-\cfrac{R}{L}&-\cfrac{K_e}{L}\\
\cfrac{K_t}{J}&0
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}
\end{array}
\right]
+
\left[
\begin{array}{c}
\cfrac{1}{L}\cdot e_{(t)}\\
0
\end{array}
\right]
$$

制御入力は電圧$${e_{(t)}}$$、制御出力は角速度$${\omega_{(t)}}$$なことを考慮すると、以下の状態方程式が得られます。

$$
\begin{array}{l}
\cfrac{d}{dt}
\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}
\end{array}
\right]
=
\left[
\begin{array}{cc}
-\cfrac{R}{L}&-\cfrac{K_e}{L}\\
\cfrac{K_t}{J}&0
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}
\end{array}
\right]
+
\left[
\begin{array}{c}
\cfrac{1}{L}\\
0
\end{array}
\right]\cdot e_{(t)}\\
y=\begin{array}{c}
\left[
\begin{array}{cc}
1&0
\end{array}
\right]\cdot
\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}
\end{array}
\right]
+0\cdot e_{(t)}
\end{array}
\end{array}
$$

説明のため細かくステップを刻みましたが、慣れれば連立微分方程式からほぼ直接に状態方程式を作れるようになります。


なお、拡張性の高さも状態方程式の特徴の1つです。
例えば、制御出力を角度$${\theta_{(t)}}$$に変えたい場合、以下の微分方程式を連立させて状態方程式を拡張できます。

$$
\frac{d\theta_{(t)}}{dt}=\omega_{(t)}
$$

連立後の微分方程式は以下になります。

$$
\left \{
\begin{array}{cl}
\cfrac{di_{(t)}}{dt}&=-\cfrac{R}{L}\cdot i_{(t)}-\cfrac{K_e}{L}\cdot \omega_{(t)}-\cfrac{1}{L}e_{(t)} \\
\cfrac{d\omega_{(t)}}{dt}&=\cfrac{K_t}{J}\cdot i_{(t)} \\
\cfrac{d\theta_{(t)}}{dt}&=\omega_{(t)}
\end{array}
\right.
$$

このとき、状態ベクトルは下式になります。

$$
\bold{x}=\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}\\
\theta_{(t)}
\end{array}
\right]
$$

制御出力は、角速度$${\omega_{(t)}}$$のときと同じ要領で、方程式中の状態変数の係数を行列形式でまとめれば良いです。

$$
\cfrac{d}{dt}\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}\\
\theta_{(t)}
\end{array}
\right]
=
\left[
\begin{array}{ccc}
-\cfrac{R}{L}&-\cfrac{K_e}{L}&0\\
\cfrac{K_t}{J}&0&0\\
0&1&0
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}\\
\theta_{(t)}
\end{array}
\right]
+
\left[
\begin{array}{c}
\cfrac{1}{L}\cdot e_{(t)}\\
0\\
0
\end{array}
\right]
$$

制御入力は電圧$${e_{(t)}}$$、制御出力は角度$${\theta_{(t)}}$$なことを考慮すると、以下の状態方程式が得られます。

$$
\begin{array}{l}
\cfrac{d}{dt}
\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}\\
\theta_{(t)}
\end{array}
\right]
=
\left[
\begin{array}{ccc}
-\cfrac{R}{L}&-\cfrac{K_e}{L}&0\\
\cfrac{K_t}{J}&0&0\\
0&1&0
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}\\
\theta_{(t)}
\end{array}
\right]
+
\left[
\begin{array}{c}
\cfrac{1}{L}\\
0\\
0
\end{array}
\right]\cdot e_{(t)}\\
y=\begin{array}{c}
\left[
\begin{array}{ccc}
1&0&0
\end{array}
\right]\cdot
\left[
\begin{array}{c}
i_{(t)}\\
\omega_{(t)}\\
\theta_{(t)}
\end{array}
\right]
+0\cdot e_{(t)}
\end{array}
\end{array}
$$

行列やベクトルを使うので、はじめのうちはとっつきにくいかもしれませんが、慣れれば拡張性の高さを便利に感じると思います。


4. 状態方程式のメリット

制御入力や制御出力が複数あるシステム(MIMO:Multi Input Multi Output)を扱えるのが一番のメリットです。
伝達関数は、制御入力・制御出力共に1つ(SISO:Single Input Single Output)であることが前提です。
 
現代制御は、システムを状態方程式で表現することを前提に理論が組まれています。
このため、伝達関数を基本とする古典制御より扱えるシステムが多いです。

基本的に、劣駆動系(関節数よりアクチュエータが少ないシステムのこと)は現代制御でないと扱えません。
 
 
デメリットを挙げるとすれば、制御設計に遅延要素(特にムダ時間)の影響を考慮しにくい点です。

ただ、それも不可能というわけではありません。
パデの近似法という方法を使えば、ムダ時間の高精度な近似モデルが得られます。
 
得られた近似モデルを状態方程式に組み込めば、問題なく現代制御理論を適用できます。


真偽は不明ですが、戦闘機の制御設計には、20次かそれ以上のパデ近似モデルを使うと聞いたことがあります。
まぁ、それくらいはやっていても不思議ではないとは思います。


5. おわりに

この記事では、状態方程式を紹介・解説しました。
 
通常、産業用ロボットは関節と同じ数のアクチュエータを備えるため、状態方程式(≒現代制御)が開発現場で使われることは稀ですが、知っておいて損はないと思います。
 
例えば、カルマンフィルタ設計には、状態方程式の知識が必要です。
カルマンフィルタは様々な分野で使われるので、そのためだけに状態方程式を知っておいて良いくらいの価値があると思います。

前回:寄り道編(2)「PWM制御」はこちら
次回:寄り道編(4)「システム・パラメータ同定」はこちら


付録1. 行列の指数関数計算

これはディジタル制御でも使う数学テクニックです。

制御をやっていると
「アレ?状態方程式を伝達関数化できるのでは?」
と思うことがあるはずです。
実際、それは可能です。


状態方程式をラプラス変換したものが下式です。

$$
\left \{
\begin{array}{l}
s\cdot x=Ax+Bu \\
y=Cx+Du
\end{array}
\right.
$$

1階常微分方程式を変形すると、下式が得られます。

$$
(s\cdot I-A)=Bu\\
x=(s\cdot I-A)^{-1}Bu
$$

これを出力方程式に代入します。

$$
\begin{array}{cl}
y&=C\cdot(s\cdot I-A)^{-1}Bu+Du\\
&=(C\cdot(s\cdot I-A)^{-1}B+D)u
\end{array}
$$

これで状態方程式を伝達関数化できます。

この結果が得られると、ほどなく
「これを逆ラプラス変換できるのでは?」
と思うことがあるはずです。
実際、それは可能です。


行列化しているといえ、$${B,C,D}$$は単なる係数の集まりです。
逆ラプラス変換前後で変わりません。
入力$${u}$$は自分で決めるので、特に問題ないと思います。

となると、残る疑問は「$${(s\cdot I-A)^{-1}}$$をどう扱うか?」です。

ラプラス変換表を見れば$${\mathcal{L}^{-1}\{(s-\alpha)^{-1}\}=e^{\alpha\cdot t}}$$だと分かります。
「なんだ、$${\bm{\mathcal{L}^{-1}\{(s-A)^{-1}\}=e^{A\cdot t}}}$$か。楽勝だな。」$${\\}$$
と思うことでしょう。
($${\mathcal{L}^{-1}}$$は逆ラプラス変換を意味します)


そして、その後気付くでしょう。
「行列が指数とは...?」と。
それをこれから説明します。


仮に、A行列が下式とします。

$$
A=\left[
\begin{array}{cc}
-2&0\\
0&-3
\end{array}
\right]
$$

すると、$${s\cdot I-A}$$は以下になります。

$$
\begin{array}{cl}
s\cdot I-A&=\left[\begin{array}{cc}
s&0\\
0&s
\end{array}
\right]-
\left[\begin{array}{cc}
-2&0\\
0&-3
\end{array}
\right]\\
&=\left[\begin{array}{cc}
s+2&0\\
0&s+3
\end{array}
\right]
\end{array}
$$

この逆行列は以下のように求められます。

$$
\begin{array}{cl}
(s\cdot I-A)^{-1}&=
\cfrac{1}{(s+2)\cdot(s+3)}
\left[\begin{array}{cc}
s+3&0\\
0&s+2
\end{array}
\right]\\&=
\left[
\begin{array}{cc}
\cfrac{1}{s+2}&0
\\0&\cfrac{1}{s+3}
\end{array}
\right]
\end{array}
$$

ここまで来ればあと一歩です。
行列の各要素を逆ラプラス変換すれば良いのです。

$$
\begin{array}{cl}
e^{A\cdot t}&=
\mathcal{L}^{-1}{(s-\alpha)^{-1}}\\
&=
\left[
\begin{array}{cc}
\mathcal{L}^{-1}\cfrac{1}{s+2}&0\\
0&\mathcal{L}^{-1}\cfrac{1}{s+3}
\end{array}
\right]\\
&=
\left[
\begin{array}{cc}
e^{-2\cdot t}&0\\
0&e^{-3\cdot t}
\end{array}
\right]
\end{array}
$$

このように、ラプラス変換を応用するとシンプルな計算で$${e^{A\cdot t}}$$を求められます。
なお、マクローリン展開を使っても$${e^{A\cdot t}}$$は求められます。


付録2. 状態方程式の拡張例

ここでは、外乱オブザーバ設計方法の一部を具体例として挙げます。

まず、モータの運動方程式を下式とします。
($${\tau_{d(t)}}$$は外乱トルク)

$$
J\cdot \frac{d\omega_{(t)}}{dt}=\tau_{(t)}+\tau_{d(t)}
$$

基本的に、外乱はあらかじめ分かるものではありません。
「そんなものを状態方程式に組み込むのはムリでは?」
と思っても仕方ないと思います。

たしかに、完璧に組み込むのはムリです。
しかし、ある仮定の元、それなりの精度で組み込むのは可能です。

その仮定とは、
外乱の1階時間微分値=0とすること
(外乱波形がステップ状だと仮定すること)
です。

このとき、以下の連立方程式が成り立ちます。

$$
\begin{cases}
J\cdot \frac{d\omega_{(t)}}{dt}=\tau_{(t)}+\tau_{d(t)}\\
\frac{d\tau_{d(t)}}{dt}=0
\end{cases}
$$

ここで、$${\tau_{d(t)}}$$を状態変数の1つとみなして、本文中と同じ要領で状態方程式化したものが以下です。

$$
\cfrac{d}{dt}\left[
\begin{array}{c}
\omega_{(t)}\\
\tau_{d(t)}
\end{array}
\right]
=
\left[
\begin{array}{cc}
0&\cfrac{1}{J}\\
0&0
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
\omega_{(t)}\\
\tau_{d(t)}
\end{array}
\right]
+
\left[
\begin{array}{c}
\cfrac{1}{J}\\
0
\end{array}
\right]
\cdot \tau_{(t)}\\
\omega_{(t)}=
\left[
\begin{array}{cc}
1&0
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
\omega_{(t)}\\
\tau_{d(t)}
\end{array}
\right]
$$

この状態方程式に基づきオブザーバを設計すれば、外乱トルクの推定が可能です。
(オブザーバの設計方法は別の機会に紹介しようと思います)


「状態方程式はかなり自由にいじれるものだ」
というのが分かってもらえたと思います。


前回:寄り道編(2)「PWM制御」はこちら
次回:寄り道編(4)「システム・パラメータ同定」はこちら


いいなと思ったら応援しよう!