資料請求番号：SH43 TS53

CONTENTS

自動制御の基礎であるPID制御モデルを解くプログラムをVBAで作成

PID制御という言葉は、工学に携わったことのある人なら一度は聞いたことがあると思います。工学の目的は「適切な機械や装置を使って、その装置を意のままに操って必要なモノを作ること」です。この中で自動制御という考え方が重要になってきます。

今回は自動制御の基本やPID制御の基本を押さえながら、
実際のプロセスモデル（2槽タンクモデル）にPID制御系を組ませたときのプロセス挙動を予測するプログラムを作成していきたいと思います。

※本来、自動制御ではラプラス変換を行って微分方程式を解きますが、VBAで書かれたルンゲクッタ法のプログラムがあるので、それを使用して直接微分方程式を解きます。

制御とは？

PID制御・・・とは言いますが、そもそも、制御とはどういう意味でしょうか？身近な例で制御のお話をしたいと思います。

オマエは「制御とは何ですか？」って聞かれたときにしっかり答えられるか？

制御？せ、制御は制御でしょ？？

例えばさ、オマエが東名に乗るためにETCレーンを20km/hで通過する。それでランプを進み40km/hまで加速する。
このまま本線に入るか？

入るわけないでしょ！！？危ないでしょ！

じゃあどうする？

どうするって・・・そんなのアクセル踏んで100km/hまで加速するに決まってるじゃない。

それで、アクセルはずっと同じ踏み方でいいか？ずっとベタ踏みか？

いや、アクセルは柔らかく踏んで加速してきたと同時に強く踏み込んでいくのよ。いきなりベタ踏みはダメ。
それで100km/h超えそうだなと思ったら少しアクセル緩めて、最終的に一定の強さで踏みつづけて100㎞/hで走り続けるのよ。
当たり前でしょう？

その「あたりまえでしょう？」って言ったことが制御っていうたいね。

・・・どういうこと？

時速40km/hの状態から時速100㎞/hで運転するために、アクセルを強弱の調整を加えながら目標のスピードに到達していく。
こういう風に制御というのは、装置を自分の意のままに操る。
装置や機械を望み通りの運転状態にすることを言うたい。
今の話で言う装置や機械というのは車だな。

装置や機械を望み通りの運転状態にする・・・？

せや。オマエがランプから本線に入るまでに100km/h巡航という目標に向かって車を意のままに操っているやろ？オマエが車を制御しているたい。

確かに、ワタシが車という装置を、アクセルを使って制御して、100km/hというを望み通りの運転状態を維持した。という表現ができるわね。

せや。アクセルを踏むというのが因、100km/hのスピードが出るというのが果。因果はインプット、アウトプットという。機械や装置にインプットをしたら、何らかのアウトプットが返ってくるわけたい。

そして、これから、俺らがやろうとしている制御っていうのは、さっきオマエが言ったようなインプットの手順をアルゴリズム化してコンピュータでも動けるようにする。
つまり、装置を意のままに操る手順をアルゴリズム化してコンピュータにプログラムして自動制御を作るということになるたい。

～2槽タンクモデルの場合～

なるほど・・。さっき、ワタシが「当たり前でしょ？」って言ったことはあくまでワタシだけの「当たり前」であって、その部分をちゃんとアルゴリズムにして、自動制御にすることで、誰がやっても同じ結果になるようにするのね。

そういうこと！

制御というのは「装置や機械を意のままに操る」ことを言います。
車を時速100km/hで走らせる、槽内温度を75℃に保つ、回転数６００rpmで攪拌する、タンク内水位を8mに保つ・・・
それぞれの目的に合わせて、どのようにしたら、上記の目標を達成できるのか、その方法をアルゴリズム化してコンピュータにプログラムすることができれば、コンピュータによる自動制御ができるようになります。

PID制御はそのアルゴリズムの1種と位置付けられます。

問題設定（2槽タンクモデル）

今回、制御してみようと思うプロセスは２槽のタンクで、その液面高さを制御するたい。

バルブの開度はそれぞれで固定で、２槽目の液面高さy₂を俺らの好きなように調整したい。このとき、どうする？

えっと・・・１槽目の流量、q_in1ってやつを増やしたり減らしたりすれば、できそうね。

では、今タンクの水面の高さは4mだが、水位を上げて8mにしたい。１槽目の流量はどうすればいい？

えっと・・・とりあえず増やして・・・でもどれだけ増やせばいいか分からないし・・・増やし過ぎて１槽目のタンクが溢れたら大変・・・。

せやろ？普通はそうなるたい。何回もこの操作をやったことのあるベテランだったら、流入バルブはこのくらい開ければいいかな？って感覚がついとるかもしれんがな。
でも、実際には現場行くとさ、水位8mって設定してOKボタンクリックしたら、自動で8mになるような機械があるやろ？

うん。そうだね。

あれ、自動で流入流量を調節してんねん。その自動でされる部分がどういう仕組みになっているのかを知るって言うのが今回のテーマ。

なるほど・・・。

プロセスのな、インプットをアルゴリズム化するには、プロセスのインプットとアウトプットを定量的に評価できなあかん。
まずは物質収支とエネルギー収支を立てる。インプットは１槽目の流量。アウトプットは２槽目の液面高さやな。
今回は温度変化は考えない系だから物質収支だけでよか。1槽目、2槽目の物質収支式は以下の通り。

えっと・・・。蓄積量の時間変化＝流入量―流出量よね。

せや。それで、それぞれの槽の底面積をA₁,A₂としたとき、体積は底面積×液面高さであることと、
それぞれの槽の流出流量をバルブ抵抗R[hr/m²]を使って以下の様に定義する。

これ、1槽目の流出流量と2槽目の流入流量は一緒よね。

せやせや。1槽目の出口をそのまま2槽目につないどるからな。
こうして、さっきの物質収支式はこんな風に式変形される。

最後に両辺をρAで割れば、液面高さに関する2つの連立微分方程式ができる。

おおっ、ルンゲクッタ法で解けそうな微分方程式になった！

とりあえず、系のモデル化はこれで終わりや。次からはこの系に制御系を組んで行く。

PID制御とは？

PID制御は、PはProportional、比例制御、IはIntegral、積分制御、Dはdifferental、微分制御を略してPIDと呼んでいます。
では、何に比例し、何を積分し、何を微分しているのでしょうか？そしてその操作は系の制御にどのように貢献するのでしょうか？

P制御

PIDを学ぶときはPとIとDそれぞれ分けて考えたほうがわかりやすい。
まず、Pだが、これは比例制御。

比例って・・・何を何に比例させるのよ・・？

おっ、ええ質問やなぁ～。
例えばオマエ、20km/hから100km/hに加速する時と60km/hから100km/hに加速する時とではどっちの方がアクセルを強く踏む？

20km/hの方に決まってるじゃない。

何で？

何でって・・・・。何でっていわれても・・・。

２０～１００と、６０～１００との間では前者の方が目標の速さとのギャップが大きいから、よりエンジンふかさないかんな。
ポイントは今の値と設定値との差、つまり偏差やねん。これをeで表す。

y_2setって言うのは2槽目の液面高さの設定値、y₂は現在の液面高さ。例えば液面高さ4mで定常状態になっているプロセスに対して、8mの設定値を与えると偏差は4mになるな。これに比例する量だけ流入流量を増やす。というのが比例制御たい。

流量qss、液面高さ4mで定常になってたところ、液面高さ8mに設定したら、上の式に従って流量が増えるような仕組みを表しているのね！

せや。偏差に比例する分だけ流入流量を追加しとるやろ？y=8mって設定したら上の式に従って流量を自動的に変えてくれるたい。

それで、この式を微分すると・・・

流入流量q_in1に関する微分方程式ができる。流入流量は偏差に従って動く関数。そしてその偏差と言うのは設定値と現在地の差のことやから、流入流量が2槽目のタンクの液面高さの関数になった。これが何を意味しているかわかるか？

ある目的の液面高さを達成するために必要な流入流量が一義に決まったってことでしょ？

そういうこと！インプットとアウトプットの関係がこれでハッキリしたよな！
これと、さっきの2槽タンクの物質収支式

を、連立して解けば、P制御系がどのように動くかシミュレーションできる。

問題①
底面積30m²、高さ10m、バルブ抵抗0.08の2つのタンクを用いて以下のようなプロセスを組んだ。

現在、流量q_ss = 50m³/h、両タンクの液面高さ4mで定常状態になっている。
ここで2槽目の液面高さの設定値を8m（y_2set=8m）としたら、2槽目の液面高さはどのような応答をするか。P制御系で予測せよ。ただし、比例ゲインKc=20とする。

つまり、3つの連立微分方程式をこのシートで解くVBAプログラムを組めばいいわけね。

q_in1の初期条件は？

q_in1はさっきの偏差の式に従う。

あ、そっか。t=0の時も例外なくこの式に従うから、以下の計算式で求めればいいのね。

・・・

できた！

※VBAソースコード　P制御

/* –Pcontrol　rk-method ——– */
Public A1, R1, A2, R2, qss, y2set, Kc As Double

Sub rkdoublep()

Dim y10, y20, q0, init, ed, h, h2, ky1(4), ky2(4), kq(4) As Double
Dim i, j As Integer

‘Calculate propertv
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)

‘Constants
A1 = Cells(5, 3)
R1 = Cells(6, 3)
A2 = Cells(7, 3)
R2 = Cells(8, 3)
qss = Cells(9, 3)
y2set = Cells(10, 3)
Kc = Cells(11, 3)

‘Initial Condition
y10 = Cells(4, 7)
y20 = Cells(4, 8)
Cells(4, 9) = Kc * (y2set – y20) + qss
q0 = Cells(4, 9)
y1 = y10
y2 = y20
q = q0
t = init

‘Runge-kutta roop

For i = 0 To ((ed – init) / h)

j = 6 + i / h2

Cells(1, 6) = i
Cells(2, 6) = j – 6

If i Mod h2 = 0 Then

Cells(j, 6) = t
Cells(j, 7) = y1
Cells(j, 8) = y2
Cells(j, 9) = q

End If

ky1(1) = h * F1(t, y1, y2, q)
ky2(1) = h * F2(t, y1, y2, q)
kq(1) = h * F3(t, y1, y2, q)
ky1(2) = h * F1(t + h / 2, y1 + ky1(1) / 2, y2 + ky2(1) / 2, q + kq(1) / 2)
ky2(2) = h * F2(t + h / 2, y1 + ky1(1) / 2, y2 + ky2(1) / 2, q + kq(1) / 2)
kq(2) = h * F3(t + h / 2, y1 + ky1(1) / 2, y2 + ky2(1) / 2, q + kq(1) / 2)
ky1(3) = h * F1(t + h / 2, y1 + ky1(2) / 2, y2 + ky2(2) / 2, q + kq(2) / 2)
ky2(3) = h * F2(t + h / 2, y1 + ky1(2) / 2, y2 + ky2(2) / 2, q + kq(2) / 2)
kq(3) = h * F3(t + h / 2, y1 + ky1(2) / 2, y2 + ky2(2) / 2, q + kq(2) / 2)
ky1(4) = h * F1(t + h, y1 + ky1(3), y2 + ky2(3), q + kq(3))
ky2(4) = h * F2(t + h, y1 + ky1(3), y2 + ky2(3), q + kq(3))
kq(4) = h * F3(t + h, y1 + ky1(3), y2 + ky2(3), q + kq(3))
ny1 = y1 + (ky1(1) + 2 * ky1(2) + 2 * ky1(3) + ky1(4)) / 6
ny2 = y2 + (ky2(1) + 2 * ky2(2) + 2 * ky2(3) + ky2(4)) / 6
nq = q + (kq(1) + 2 * kq(2) + 2 * kq(3) + kq(4)) / 6
nt = t + h

t = nt
y1 = ny1
y2 = ny2
q = nq

Next i

End Sub

Function F1(ByVal t As Double, ByVal y1 As Double, ByVal y2 As Double, ByVal q As Double) As Double
F1 = (1 / A1) * q – (1 / (A1 * R1)) * y1
End Function

Function F2(ByVal t As Double, ByVal y1 As Double, ByVal y2 As Double, ByVal q As Double) As Double
F2 = (1 / (A2 * R1)) * y1 – (1 / (A2 * R2)) * y2
End Function

Function F3(ByVal t As Double, ByVal y1 As Double, ByVal y2 As Double, ByVal q As Double) As Double
F3 = -Kc * F2(t, y1, y2, q)
End Function

けど、このままでいいのかしら・・・。

設定値は8mのはずやが・・。

そうなのよ。それなのに6.5mくらいで安定しちゃって・・・・。

比例制御やとこういう定常偏差、つまりオフセットが起こるたい。流量のグラフは？

これ。

あ～。これ、流入流量80m³/h,液面高さ6.5mで安定しちゃったんやな。

比例制御じゃ制御できないのかしら・・・。一応、Kcを200まで大きくしたんだけど・・・

これはイカン。振動が大きくて危険や！それに、オフセットは結局残る。もっと言えば、850m³/hなんて流量が現存設備で達成できるか、怪しいやろ？

そうね。比例以外の他の関数を入れないといけないってことなのかしら？

せやな。目標値に達するまでは、流入流量は常に変化させるという工夫が必要になるたいね。

PI制御

比例制御に積分制御を追加したのがPI制御。
元々の流量q_ssに、偏差に比例した値とこれまでの偏差の積分値を足して流入流量q_in1を決定する。

偏差の積分を追加することによって、偏差がある限り、q_in1の値は一定にならないということ？

せや！仮にさっきみたいなオフセットができるような状況になったとしても、時間が経つにつれ、積分項の数値は大きくなっていくから、q_in1の値は大きくなっていく。それが、偏差がなくなるまで続くたい。

先ほどの偏差の式を微分して・・・

これを2槽タンクの物質収支式と連立させて解けばPI制御の挙動がわかる。

問題②
底面積30m²、高さ10m、バルブ抵抗0.08の2つのタンクを用いた、問題①と同じプロセスにおいて、
現在、流量q_ss = 50m³/h、両タンクの液面高さ4mで定常状態になっている。
ここで2槽目の液面高さの設定値を8m（y2set=8m）としたら、2槽目の液面高さはどのような応答をするか。PI制御系で予測せよ。ただし、比例ゲインKc=20、積分時間TI=2とする。

このエクセルシートで・・・

このプログラムを動かせばできるわよね。

※VBAソースコード　PI制御

/* –PIcontrol　rk-method ——– */Public A1, R1, A2, R2, qss, y2set, Kc, TI As Double

Sub rkdoublepi()

Dim y10, y20, q0, init, ed, h, h2, ky1(4), ky2(4), kq(4) As Double
Dim i, j As Integer

‘Calculate propertv
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)

‘Constants
A1 = Cells(5, 3)
R1 = Cells(6, 3)
A2 = Cells(7, 3)
R2 = Cells(8, 3)
qss = Cells(9, 3)
y2set = Cells(10, 3)
Kc = Cells(11, 3)
TI = Cells(12, 3)

‘Initial Condition
y10 = Cells(4, 7)
y20 = Cells(4, 8)
Cells(4, 9) = Kc * (y2set – y20) + qss
q0 = Cells(4, 9)
y1 = y10
y2 = y20
q = q0
t = init

‘Runge-kutta roop

For i = 0 To ((ed – init) / h)

j = 6 + i / h2

Cells(1, 6) = i
Cells(2, 6) = j – 6

If i Mod h2 = 0 Then

Cells(j, 6) = t
Cells(j, 7) = y1
Cells(j, 8) = y2
Cells(j, 9) = q

End If

t = nt
y1 = ny1
y2 = ny2
q = nq

Next i

End Sub

Function F1(ByVal t As Double, ByVal y1 As Double, ByVal y2 As Double, ByVal q As Double) As Double
F1 = (1 / A1) * q – (1 / (A1 * R1)) * y1
End Function

Function F2(ByVal t As Double, ByVal y1 As Double, ByVal y2 As Double, ByVal q As Double) As Double
F2 = (1 / (A2 * R1)) * y1 – (1 / (A2 * R2)) * y2
End Function

Function F3(ByVal t As Double, ByVal y1 As Double, ByVal y2 As Double, ByVal q As Double) As Double
F3 = -Kc * F2(t, y1, y2, q) + (Kc / TI) * (y2set – y2)
End Function

できた！すごい！ちゃんと高さ8mに制御できてる！

せやな！ええ感じになってきた。

・・・でもな～y₂の値が大きく振動しているのが気になるねん。
オーバーシュートしていて危険やし・・・。

確かに。一応これ、計算上大丈夫だけど、実際何が起こるか分からないもんね・・・。

せや。現実には何らかの外乱があって、このオーバーシュートが拡大してしまった・・・なんてことが起きたら怖いやろ？
オーバーシュートはなくさないかん。

なるほど・・・。

比例ゲインKcを大きくしたり、積分時間TIを小さくしたりすると急速に目標値に近づけることができるかもやけど、
応答に振動をもたらすことがあるけん、それはプロセス制御の観点からして好いとうない。

じゃあ、どうすればいいのかしら・・・。

I制御を取り入れることによって、「P制御だけではオフセットが残る」という問題を解決することができました。
I制御はこれまでの偏差の積分値を流入流量の値q_in1に取り入れることによって、偏差がある限り、q_in1は更新され続ける。というアルゴリズムを作ることができました。
しかし、PI制御ではゲインを大きくすると応答特性（今回の場合は液面高さ）が振動してしまうという問題点が出てきます。
応答特性の振動を抑えるために微分制御、D制御という考え方を取り入れます。次は比例・積分・微分全ての制御がそろったPID制御のお話です。

PID制御

P制御ではプロセスの現在を制御する。I制御ではプロセスの過去を制御する。この感覚がわかるか？

えっと、P制御では偏差にそのまま比例定数をかけているだけ、I制御では偏差を積分している。その積分の範囲って言うのは、偏差が発生してから現在までの時間だから、過去を見ているっていう感じなのかしら？

せやせや。PID制御って言うのは今までPI制御にDifferential、つまり微分の要素が加わった制御たい。

微分ってことは・・・未来を見ているってこと？

その通り！その部分の関数の傾きを見るということは、今後、その間数がどんな風に変化していくのかを予測しているってことたい。

そして、偏差の動きが激しい時を検出して、これを修正する働きをするのがD制御たい。

例によってq_in1を微分して・・・微分すると2階微分が出来るが、dy₂/dtは既に定式化されてて、これを微分することができることに注意して、q_in1に関する微分方程式を立てる。

これを2槽タンクの物質収支式と連立させて解けばPID制御の挙動がわかる。

問題③
底面積30m²、高さ10m、バルブ抵抗0.08の2つのタンクを用いた、問題②と同じプロセスにおいて、
現在、流量qss = 50m³/h、両タンクの液面高さ4mで定常状態になっている。
ここで2槽目の液面高さの設定値を8m（y2set=8m）としたら、2槽目の液面高さはどのような応答をするか。PID制御系で予測せよ。ただし、比例ゲインKc=20、積分時間TI=2、微分時間はTD=1.2とする。

さっきと同じようにエクセルシートで・・・

このプログラムを動かせばできるわよね。

※VBAソースコード　PID制御

/* –PIDcontrol　rk-method ——– */
Public A1, R1, A2, R2, qss, y2set, Kc, TI, TD As Double

Sub rkdoublepid()

Dim y10, y20, q0, init, ed, h, h2, ky1(4), ky2(4), kq(4) As Double
Dim i, j As Integer

‘Calculate propertv
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)

‘Constants
A1 = Cells(5, 3)
R1 = Cells(6, 3)
A2 = Cells(7, 3)
R2 = Cells(8, 3)
qss = Cells(9, 3)
y2set = Cells(10, 3)
Kc = Cells(11, 3)
TI = Cells(12, 3)
TD = Cells(13, 3)

‘Initial Condition
y10 = Cells(4, 7)
y20 = Cells(4, 8)
Cells(4, 9) = Kc * (y2set – y20) + qss
q0 = Cells(4, 9)
y1 = y10
y2 = y20
q = q0
t = init

‘Runge-kutta roop

For i = 0 To ((ed – init) / h)

j = 6 + i / h2

Cells(1, 6) = i
Cells(2, 6) = j – 6

If i Mod h2 = 0 Then

Cells(j, 6) = t
Cells(j, 7) = y1
Cells(j, 8) = y2
Cells(j, 9) = q

End If

t = nt
y1 = ny1
y2 = ny2
q = nq

Next i

End Sub

Function F1(ByVal t As Double, ByVal y1 As Double, ByVal y2 As Double, ByVal q As Double) As Double
F1 = (1 / A1) * q – (1 / (A1 * R1)) * y1
End Function

Function F2(ByVal t As Double, ByVal y1 As Double, ByVal y2 As Double, ByVal q As Double) As Double
F2 = (1 / (A2 * R1)) * y1 – (1 / (A2 * R2)) * y2
End Function

Function F3(ByVal t As Double, ByVal y1 As Double, ByVal y2 As Double, ByVal q As Double) As Double
F3 = -Kc * F2(t, y1, y2, q) + (Kc / TI) * (y2set – y2) – Kc * TD * ((1 / (A2 * R1)) * F1(t, y1, y2, q) – (1 / (A2 * R2)) * F2(t, y1, y2, q))
End Function

オーバーシュートがなくなることはないけど、さっきよりはマシになったね。

そうだな。最初の傾きが緩くなって2周期目からの振動がなくなったな。

D制御とは、偏差の微分値を求め、それをインプット（q_in1）へ反映する仕組みのことです。微分を求めるということは、未来のプロセスを予想することに相当します。
ルンゲクッタ法でも、関数の傾きから、h/2だけ進んだところのyの値を求めていました。
傾きを出すことによって、今後関数がどのように変化していくのか予測し、偏差が大きい方向に傾いていたらそれを修正することによって、応答特性を良くするのがD制御の役割です。

チューニング

この問題、オーバーシュートをなくして、しっかりとy=8mに制御することってできるのかしら？

できるぜ。例えば、Kc=20, TI=5, TD=1でシミュレーションしてみな。

あっ！できた！すご～い！！

ちょうど良い比例ゲイン、積分時間、微分時間に調整することをチューニングっていうたい。
チューニングは昔から制御においては大きな課題でな。現場の作業員が上手いことチューニングしている場合も多いんや。
水俣事業所で働く俺の知り合いにもおるけん、めっちゃチューニングプロってるベテランがな。
最近は、限界感度法のような理論に基づくオートチューニングって言うのが普及し始めているけど、プロセスが複雑になるとやっぱり難しい。

チューニングミスるとこんなカーブすることもあるでしょ？これ実際にやったら始末書よ。

※Kc = 10, TI = 1 TD＝0

だから、自動制御とはいえ、全て機械任せ、ブラックボックスにしてたらいけんくて、
プロセスの動特性の理論をしっかりさせておかないといけんよな。理論がしっかりわかっていたらいたずらにKcを大きくしたり、TIを小さくしたりしないよな。

そうね。

まとめ

今回は、PID制御モデルを解くVBAマクロの作成ということで、
まず制御とは何なのかを説明したのち、P,I,D制御それぞれのアルゴリズムを説明しました。そのアルゴリズムに基づくVBAプログラムを作成し、
2槽タンクの水位モデルにおいて、液面がどのように動くのかをシミュレーションしてみました。

PID制御など、自動制御は普通、ラプラス変換を行い、微分方程式を代数方程式に変換して解を求めるのが一般的な解析法ですが、微分方程式を計算できるVBAプログラムを持っていたので、それを活用しました。
さらにExcelのグラフシートを使うことによって、リアルタイムでシミュレーションが可能になるのもこのVBAプログラムの利点です。