PIDチューニングにScipyを使ってみたら、思ったよりも良いチューニングが出来た。

#1. 概要  

- [前回](https://qiita.com/hikaruyaku/private/021769a76e72ed34a881)、`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秒とした。

```python: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の値がそれぞれ負にならないように制約条件をつける必要がある。以下のコードにて、それを行った。


```python:constraint.py
def cons(x):
    return x
cons = ({'type':'ineq', 'fun':cons})
```

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

```python: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の値が負となったので、除外。  


![図1.png](https://qiita-image-store.s3.amazonaws.com/0/171392/8f537df4-e400-3991-ca68-3ca72a01ea2f.png)

- ここから、`CG`や`BFGS`の挙動がステップ応答に追従せず、発散している。  
- `COBYLA`に関しては、時定数が大きいのもありかなりゆっくりした挙動だった。100秒までで収束していない。

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