「制御のしやすさ」を定量化した話

　はじめまして、[りょーつ](https://twitter.com/Gomachan_2021)といいます。高専出身の大学院1年生です。研究の専門は力学や機構学で、主に制御工学や数学に関する記事を書いています。
　今回はシステムの「制御のしやすさ」について考えてみようと思います。

# 目次
[1.はじめに](#1-はじめに)
[2.時間の無い方へ](#2-時間の無い方へ)
[3.レギュレータ制御の時間応答](#3-レギュレータ制御の時間応答)
[4.レギュレータ制御の入力最適化](#4-レギュレータ制御の入力最適化)
[5.制御のしやすさと可制御性グラミアン](#5-制御のしやすさと可制御性グラミアン)
[6.おわりに](#6-おわりに)
[7.参考文献](#7-参考文献)


# 1. はじめに
　本稿では線形なシステムの「制御のしやすさ」を定量的に表す方法を紹介し、制御しやすいシステムを作る方法を考えてみます。
　制御工学の講義では、あるシステムが与えられて、そのシステムをうまく動かす方法を学ぶことがメインとなりがちですが、**「良い制御は良いハードから」** と言われるように、制御しやすいシステムを作ることは非常に大事です。内容は制御工学っぽいですが、良いハードを作るためにも、**機構を作っている方に是非読んでみてほしい内容** となってます。
　大学の一般教養程度の数学と、工学部で学ぶ現代制御の知識を使います。

# 2. 時間の無い方へ
　時間が無い方のために冒頭で結論を述べておきます。「制御のしやすさ」は「入力$u$の大きさを制限した場合に、より多様な偏差$\boldsymbol{x_e}$を無効化できるシステム」と定義しました。数式にすると以下のように表され、(1)式の右辺が大きければ大きいほど制御しやすくなるはずです。
```math
制御のしやすさ
= 
\mathrm{tr} (W_c)
\tag{1}
```
```math
W_c = \int_{0}^{t_f} e^{-At}\boldsymbol{b} \boldsymbol{b}^\mathrm{T} e^{-A^\mathrm{T}t} \ dt
\tag{2}
```
ただし、状態空間を
```math
\dfrac{d}{dt} \boldsymbol{x}  = A\boldsymbol{x} +\boldsymbol{b} u 
\tag{3}
```
とし、$\boldsymbol{x}$を$n$次元状態ベクトル、$A$を$n \times n$行列、$\boldsymbol{b}$を$n$次元ベクトル、$u$をスカラー入力としました。

# 3. レギュレータ制御の時間応答
　まずは制御対象を定義します。(3)式と同様に制御対象の状態空間を
```math
\dfrac{d}{dt} \boldsymbol{x}(t)  = A\boldsymbol{x}(t) +\boldsymbol{b} u(t) 
\tag{4}
```
とします。$\boldsymbol{x}(t)$を$n$次元状態ベクトル、$A$を$n \times n$行列、$\boldsymbol{b}$を$n$次元ベクトル、$u(t)$をスカラー入力としました。※簡単のために1入力系を考えており、$u$はベクトルではなくスカラーとしています。
　(4)式を変形すると、時刻$t_f$におけるシステムの状態$\boldsymbol{x}(t_f)$は以下のように表現されます$^{1) \ 2)}$。
```math
\boldsymbol{x}(t_f)
=
e^{At_f} \biggl(\boldsymbol{x}(0) + \int_{0}^{t_f}e^{-At}\boldsymbol{b} u(t) \ dt\biggr)
\tag{5}
```
ここで時刻$t = 0$のときにシステムが何らかの偏差$\boldsymbol{x}(0) = \boldsymbol{x_e}$を持っていたとします。そして$t = t_f$において$\boldsymbol{x}(t_f) = \boldsymbol{0}$となるような制御(レギュレータ制御)を考えます。これらを(5)式に代入すると
```math
\boldsymbol{0}
=
e^{At_f} \biggl(\boldsymbol{x_e} + \int_{0}^{t_f}e^{-At}\boldsymbol{b} u(t) \ dt \biggr)
```
```math
\therefore
\boldsymbol{x_e} + \int_{0}^{t_f}e^{-At}\boldsymbol{b} u(t) \ dt = \boldsymbol{0}
\tag{6}
```
の関係が得られます。(6)式は、式を満たすように入力$u(t)$を調整することで、$t=t_f$において偏差$\boldsymbol{x_e}$が無効化されるということを表しています。

# 4. レギュレータ制御の入力最適化
　(6)式を満たす入力$u(t)$には様々なものが考えられ、$u(t)$は一意には定まりません。ですが、$u(t)$の解がたくさんあるというのも変な感じがします。モヤモヤを解消するために、ここでは数ある$u(t)$のうち、
```math
J = \int_{0}^{t_f} {u(t)}^2 dt
\tag{7}
``` 
で定義される評価関数$J$を最小化するような入力$u_{opt}(t)$を求めてみましょう。(6)式は$\boldsymbol{x_e}$を無効化するために必要な入力$u(t)$に関する陰関数形式の「制約条件」と捉えることができます。(7)式の評価関数と(6)式の制約条件が見つかったので、**ラグランジュの未定乗数法**が使えそうです。ここで注意しないといけないのは、(6)式が$n$次元ベクトルに関する式だということです。そのため、未定乗数が$n$個あることに気をつけてください。$n$個の未定乗数をまとめた$n$次元ベクトル$\boldsymbol{\lambda} $を以下のように定義します。
```math
\boldsymbol{\lambda} 
=
\begin{pmatrix}
 \lambda_1 \\
 \lambda_2 \\
 \vdots\\
 \lambda_n
\end{pmatrix}
\tag{8}
``` 
少しややこしいですが、未定乗数$\boldsymbol{\lambda}$を使った関数$\mathcal{L}$を
```math
\mathcal{L} 
=
\int_{0}^{t_f} {u(t)}^2 dt
- 
\boldsymbol{\lambda}^\mathrm{T} \biggl( \boldsymbol{x_e}　+ \int_{0}^{t_f}e^{-At}\boldsymbol{b} u(t)  \ dt  \biggr)
\tag{9}
``` 
とします。$\mathcal{L}$の極値を探すために$u$で偏微分してみましょう。特に深いことは考えず、積分の中身を微分します。
```math
\dfrac{\partial \mathcal{L}}{\partial u}
=
\int_{0}^{t_f} 2u(t) \ dt
- 
\boldsymbol{\lambda}^\mathrm{T} \int_{0}^{t_f}e^{-At}\boldsymbol{b}\ dt
``` 
```math
=
\int_{0}^{t_f} \bigl( 2u(t) 
- 
\boldsymbol{\lambda}^\mathrm{T} e^{-At}\boldsymbol{b}\bigr) \ dt
\tag{10}
``` 
$\mathcal{L}$が極値を取るとき(10)式の右辺は$0$になるので最適な入力$u_{opt}$は以下のように書き表せます。
```math
2u_{opt}(t) - \boldsymbol{\lambda}^\mathrm{T} e^{-At}\boldsymbol{b} = 0
``` 
```math
\therefore
u_{opt}(t) = \dfrac{1}{2}\boldsymbol{\lambda}^\mathrm{T} e^{-At}\boldsymbol{b} 
=
\dfrac{1}{2}\boldsymbol{b}^\mathrm{T}e^{-A^\mathrm{T}t}\boldsymbol{\lambda}
\tag{11}
``` 
次に$\boldsymbol{\lambda}$をもとめます。$u_{opt}(t)$は$u(t)$の1つなので、(6)式を満たします。そこで(11)式を(6)式に代入すると
```math
\boldsymbol{x_e} +\dfrac{1}{2} \int_{0}^{t_f}e^{-At}\boldsymbol{b} \boldsymbol{b}^\mathrm{T}e^{-A^\mathrm{T}t}\boldsymbol{\lambda} \ dt = \boldsymbol{0}
\tag{12}
```
となります。ここで(12)式中の積分について以下の行列$W_c$を定義すると
```math
W_c = \int_{0}^{t_f} e^{-At}\boldsymbol{b} \boldsymbol{b}^\mathrm{T} e^{-A^\mathrm{T}t} \ dt
\tag{13}
```
(12)式は次のように簡略化されます。なお、この$n\times n$行列$W_c$を可制御性グラミアンと呼びます$^{3)}$。
```math
\boldsymbol{x_e} +\dfrac{1}{2} W_c \boldsymbol{\lambda} = \boldsymbol{0}
\tag{14}
```
(14)式の両辺を等しくするためには
```math
\boldsymbol{\lambda} = -2 {W_C}^{-1}\boldsymbol{x_e} 
\tag{15}
```
である必要があり、(15)式を(11)式に代入することで最適入力$u_{opt}(t)$が次のように決定されます。
```math
u_{opt}(t) = - \boldsymbol{b}^\mathrm{T}e^{-A^\mathrm{T}t}{W_C}^{-1}\boldsymbol{x_e} 
\tag{16}
```

# 5. 制御のしやすさと可制御性グラミアン

では最適入力$u_{opt}(t)$を使ったときに、(7)式の評価関数はどうなるでしょうか？(16)式を(7)式に代入して計算してみましょう。
```math
J_{opt} = \int_{0}^{t_f} {u_{opt}(t)}^2 \ dt
``` 
```math
 = \int_{0}^{t_f} \boldsymbol{x_e}^\mathrm{T}{W_C}^{-1} e^{-At}\boldsymbol{b} \boldsymbol{b}^\mathrm{T}e^{-A^\mathrm{T}t}{W_C}^{-1}\boldsymbol{x_e} \ dt
``` 
```math
 = \boldsymbol{x_e}^\mathrm{T}{W_C}^{-1} 
 \biggl( \int_{0}^{t_f}  e^{-At}\boldsymbol{b} \boldsymbol{b}^\mathrm{T}e^{-A^\mathrm{T}t} \ dt \biggr){W_C}^{-1}\boldsymbol{x_e}
``` 
```math
 = \boldsymbol{x_e}^\mathrm{T}{W_C}^{-1} 
 W_c {W_C}^{-1}\boldsymbol{x_e}
``` 
```math
\therefore
J_{opt} = \boldsymbol{x_e}^\mathrm{T}{W_C}^{-1} \boldsymbol{x_e}
 \tag{17}
``` 
驚くことに、最適入力$u_{opt}(t)$を蓄積した量である$J_{opt}$は(17)式の右辺に示すような二次形式で表されることが分かりました。そこで
```math
J_{opt} = \int_{0}^{t_f} {u_{opt}(t)}^2 \ dt = 1
 \tag{18}
``` 
という制約を与えてみると、可制御性グラミアンが正定値行列であることから無効化可能な偏差$\boldsymbol{x_e}$のはんいは以下の楕円体の式で表されます$^{4) \ 5)}$。
```math
\boldsymbol{x_e}^\mathrm{T}{W_C}^{-1} \boldsymbol{x_e} = 1
 \tag{19}
``` 
楕円体が大きいということは偏差$\boldsymbol{x_e}$が$\boldsymbol{0}$からより大きな値を取ったとしても無効化することができ、「制御しやすい」ことを意味しています。楕円体の大きさを特徴づけるのは可制御性グラミアン$W_c$の固有値であるため、楕円体の体積は$W_c$の行列式に対応して大きくなります$^{6) \ 7)}$。具体的には楕円体の各軸の長さは可制御グラミアン$W_c$の固有値の平方根になるため、楕円体の体積は$W_c$の行列式の平方根に比例します。
　まとめると、「制御のしやすさ」は以下のように定量化できそうです。
```math
制御のしやすさ = \det{W_c}
 \tag{20}
``` 
(20)式のほかにも楕円体の各軸の長さの二乗を足し合わせたものとして
```math
制御のしやすさ = \mathrm{tr} (W_c)
 \tag{21}
``` 
も「制御のしやすさ」の指標として使えそうです。

# 6. おわりに
　本稿では可制御性グラミアン$W_c$の行列式とトレースを使って「制御のしやすさ」を定義し、紹介しました。今後は倒立振子のモデルで実際に可制御性グラミアン$W_c$を計算してみて、**「振子は長いほうが倒立させやすい」という直観的な予想を数値的に示してみたい** と考えています。
　かなりゴリゴリした内容でしたが、さいごまで読んでいただきありがとうございました！

# 7. 参考文献
1)中溝高好, 小林伸明. [システム制御の講義と演習](https://www.amazon.co.jp/%E3%82%B7%E3%82%B9%E3%83%86%E3%83%A0%E5%88%B6%E5%BE%A1%E3%81%AE%E8%AC%9B%E7%BE%A9%E3%81%A8%E6%BC%94%E7%BF%92-%E5%AE%9F%E7%94%A8%E7%90%86%E5%B7%A5%E5%AD%A6%E5%85%A5%E9%96%80%E8%AC%9B%E5%BA%A7-%E4%B8%AD%E6%BA%9D-%E9%AB%98%E5%A5%BD/dp/4817301619). 日新出版, 2015, pp.12-13

2)中野道雄, 美多勉. [制御基礎理論 -古典から現代まで-](https://www.coronasha.co.jp/np/isbn/9784339032130/). コロナ社, 2016, pp.121-123

3)川田昌克. [MATLAB/Simulinkによる現代制御入門](https://www.morikita.co.jp/books/mid/092041). 森北出版, 2015, pp.70-71, pp.200-202

4)佐藤一宏. 可制御性グラミアンと可観測性グラミアン. 初級Mathマニアの寝言, [https://ogyahogya.hatenablog.com/entry/2017/07/04/](https://ogyahogya.hatenablog.com/entry/2017/07/04/), (参照 2025-02-15)

5)日本機械学会. [機械工学便覧](https://shop.jsme.or.jp/shopdetail/000000000018/mecheng/page1/order/)$\beta \ 6$ 編, 2024, pp. 6-9

6)岩見正見, 川又政征, 樋口龍雄. [可制御性と可観測性に着目したロボットアームの性能評価](https://www.jstage.jst.go.jp/article/sicetr1965/23/2/23_2_149/_pdf). 計測自動制御学会論文集. 1986, Vol. 23, No. 2, pp.49-54

7)坂野幾海, 東俊一, 有泉亮, 浅井徹, 井村順一. [可制御性グラミアンのデータ駆動型感度制御](https://www.jstage.jst.go.jp/article/jacc/64/0/64_1201/_pdf/-char/ja). 第64回自動制御連合講演会講演予稿集. 2021, pp.1201-1205