#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秒とした。

T(x).py

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 of BFGS, L-BFGS-B, SLSQP,
depending if the problem has constraints or bounds.

よく分からないので、全部試してみた。

最後に、制約条件について述べる。minimizeするに当り、P,I,Dの値がそれぞれ負にならないように制約条件をつける必要がある。以下のコードにて、それを行った。

constraint.py

def cons(x):
    return x
cons = ({'type':'ineq', 'fun':cons})

で、それぞれ最適計算を行った。

minimize(PID).py

#初期値
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. 考察

結果として、今回のプロセスでは従来法より良い結果が得られた。
今後は、プロセス定数を変化させた場合の挙動などを改めて追っていく。
また、ガウス最適化を用いた場合についても追っていく。