1. 概要
-
前回、
python-contorlを用いて、ステップ応答やfeedbackループを構築した。
- 上記の考え方を少し応用して、PIDパラメーターを
Scipy.optimizeを用いて、最適化してみた。
- 最終的には、古典的な従来法のZiegler-Nicholsのステップ応答法と比較して、有効性を考察した(結果的には、
Scipy.optimizeを使った方がよかったです。 - 前提として、化学プラントで大多数を占めるPI制御ではなく、PID制御にて実施してみました。
2. 検討手法
2-1. 評価関数の作り方
まず、最小化する評価関数を作成する。基本的な考えとして、最小2乗法をベースとする。これだけでは、追従性の良し悪しが分からないので、二乗誤差に時間を掛けることで、その欠点をカバーする。つまり、
$$ (1-y(t))^2 * t$$
ここで、t:経過時間, y(t):時間tでの制御変数
とする。
$y(t)$に関しては、python-controlを用いて以下の関数T(x)の様に表現する。
また、評価する期間に関しては100秒とした。
from control import matlab
import numpy as np
import sympy
from scipy.optimize import minimize
def T(x):
"""
Loss Function(評価関数)
parameters
---------------------------------
x[0](=K);proprotional gain
x[1](=Ti);integral time
x[2](=Td);Derivarive time
return
---------------------------------
valuation of Loss Function
"""
"""process model"""
Kt = 10 #Steady-State gain
J = 25 #Time Constant
Cm = 4 #Dead time
num = [Kt] #分子の係数
dem = [J,1] #分母の係数
P = matlab.tf(num,dem)
num2, dem2 = matlab.pade(Cm, 6)
Pade = matlab.tf(num2,dem2)
Process = P * Pade
"""PID controller"""
K = x[0] #Proprotional Gain
Ti = x[1]#Integral Time
Td = x[2] #Derivarive Time
#伝達関数用の係数に変換
kp = K
ki = kp/Ti
kd = Td * kp
num = [kd, kp, ki]
den = [1, 0]
C = matlab.tf(num, den)
"""Feedback Loop"""
sys = matlab.feedback(Process*C,1,-1)
"""Step Response"""
time= 100 #step応答を算出する時間[秒]
t = np.linspace(0,time,1000) #時間を指定する。
y1,T = matlab.step(sys,t)
return sum(T*(y1*1000 - 1000)**2)
1000を最後にかけているのは、step応答なので1以下の値になることが考えられる。その為、二乗誤差がさらに小さくなることを懸念し、1000倍した。(やりすぎだった。。。)
2-2. 最小化の方法
最小化についてはScipy.optimize無いのminimize関数を用いた。
関数内の説明を見ると、methodとしては色々種類が有るらしく、以下の様な手法があった。
method : str or callable, optional
Type of solver. Should be one of
- Nelder-Mead
- Powell
- CG
- BFGS
- Newton-CG
- L-BFGS-B
- TNC
- COBYLA
- SLSQP
- dogleg
- trust-ncg
If not given, chosen to be one ofBFGS,L-BFGS-B,SLSQP,
depending if the problem has constraints or bounds.
よく分からないので、全部試してみた。
最後に、制約条件について述べる。minimizeするに当り、P,I,Dの値がそれぞれ負にならないように制約条件をつける必要がある。以下のコードにて、それを行った。
def cons(x):
return x
cons = ({'type':'ineq', 'fun':cons})
で、それぞれ最適計算を行った。
#初期値
init = [1,1,1]
#最適化
nelder_mead = minimize(T,init, method='Nelder-Mead')
powell = minimize(T,init, method='Powell')
cg = minimize(T,init, method='CG')
bfgs = minimize(T, init, method='BFGS')
newton_cg = minimize(T,init, method='Newton-CG')
l_bfgs_b = minimize(T, init, method='L-BFGS-B')
tnc = minimize(T, init, method='TNC')
cobyla = minimize(T, init, method='COBYLA')
slsqp = minimize(T, init, method='SLSQP')
dogleg = minimize(T, init, method='dogleg')
trust_ncg = minimize(T, init, method='trust-ncg')
2-3. 古典的なPID手法との比較
比較対象として、古典的でかつ今でも有用である"Ziegler-Nichols"のステップ応答法
を算出し、比較した。1
具体的な数式は以下の表の通り、
| 御則 | 比例ゲインKp | 積分時間Ti | 微分時間Td |
|---|---|---|---|
| P | T/KL | - | - |
| PI | 0.9T/KL | 3.33L | - |
| PID | 1.2T/KL | 2L | 0.5L |
ここで、上記と同じ様に$K=10, T=25,L=4$でPIDパラメータを推定してみると、
$$ K_p = \frac{(1,2)(25)}{(10)(4)} = 0.75$$ $$ T_I = (2)(4) = 8 $$ $$ T_D = (0.5)(4) = 2 $$
となる。
3. 検討結果
上記を全て実施し、以下の表とグラフに纏めた。
| method | P | I | D | result |
|---|---|---|---|---|
| Nelder-mead | 0.5391335 | 25.07839392 | 2.00272656 | 99,900,739.6 |
| Powell | -0.59740189 | 5.5945868 | 123.97611877 | int |
| CG | 0.34445442 | 1.7559693 | 1.13736957 | 95,857,663,934,531.3 |
| BFGS | 0.34445442 | 1.75596931 | 1.13736958 | 95,857,651,292,393.7 |
| L-BFGS-B | 0.53913351 | 25.07851699 | 2.00272495 | 99,900,739.6 |
| COBYLA | 0.02800935 | 2.90217303 | 2.39647272 | 3,593,250,685.4 |
| SLSQP | -1.28195490e+26 | 1.47833896e+26 | 2.68633657e+25 | 0.0 |
| Ziegler-Nichols(従来法) | 0.75 | 8 | 2 | 222,281,988.5 |
- "Newton-CG", "TNC", "dogleg", "trust-ncg"については、収束が不可能だった。
ここらへんは追加で勉強中。。。
- "Powell"法と"SLSQP"では、Pの値が負となったので、除外。
- ここから、
CGやBFGSの挙動がステップ応答に追従せず、発散している。
-
COBYLAに関しては、時定数が大きいのもありかなりゆっくりした挙動だった。100秒までで収束していない。
4. 考察
- 結果として、今回のプロセスでは従来法より良い結果が得られた。
- 今後は、プロセス定数を変化させた場合の挙動などを改めて追っていく。
- また、ガウス最適化を用いた場合についても追っていく。