この記事中のプログラムコードには誤りがあります。
詳細は下記記事を御覧ください。

otepipi.hatenablog.com

今回はJuliaでPIDパラメーターのオートチューニングを行ってみます。

今回、PIDのオートチューニングにはジーグラ・ニコルスのステップ応答を使用します。

ja.wikipedia.org

以下、wikiから転載します。

ジーグラ・ニコルスのステップ応答法では以下のような手順でゲイン値を決定する[17]。まず、調整器を介さずに制御対象単体にステップ入力を加える。これにより得られたステップ応答曲線から、曲線の勾配が最も急なところに接線を引き、その勾配をR=K/T（ここでKは定常ゲイン、Tは時定数）とする。この接線が横軸（時間軸）と交わる時刻と入力を加えた時刻との差を取り、この時間をLとする。これらの値から調節器の各パラメータを下の表のように決める。

f:id:Otepipi:20190406225408p:plain

今回はこの方法を利用してPIDのそれぞれのパラメーターを求めてみます。

例題

0℃の固体を100℃の雰囲気下に置いたときの固体の温度上昇を例題とします。

出力を固体の温度とします。
入力は雰囲気の温度です。
今回はt=0秒のときに雰囲気温度Ta=100℃にステップ入力したときの応答曲線を確認します。
遅れ時間として、100秒設定しました。
固体の温度上昇はニュートンの冷却法則に従うとします。

ja.wikipedia.org

この法則によると媒質中の固体から媒質に熱が伝わる速度は、固体の表面積及び固体と媒質の温度差に比例する。すなわち固体の持つ熱量Q 、時刻t 、固体の表面積S 、固体の温度T 、媒質の温度Tm の間には次の関係が成り立つ。

f:id:Otepipi:20190406225908p:plain

これを物体の温度T について解けば次の解が得られ、物体の温度変化の様子が求められる。

f:id:Otepipi:20190406230030p:plain

今回は簡単のため $$ \frac{\alpha S}{C} = 0.005 $$ としました。

実装

コードは以下のようになります。

using Plots
using LinearAlgebra

tm = 1000 #シミュレーション時間
T0 = 0 #温度初期値（℃）
td = 100 #Delay

# 空の配列の作成
u_ = zeros(tm + 1)
y_ = zeros(tm + 1)

for  t in 0:tm
      u = 100
 if t < td
     y=0
 else
     y =  (T0 - u) * exp(-0.005*(t-td)) + u
 end
  u_[t+1] = u
  y_[t+1] = y

end

# 傾斜の最も大きな値を探すため、差分をとる。

#差分行列の作成
dv = ones(tm + 1)
ev = -1 *ones(tm)
A =  Bidiagonal(dv,ev, :L)
A[1,:] = zeros(1,tm+1)

⊿y = A * y_
#最大値とそのインデックス
maxb,indexb = findmax(⊿y)

#最大接線の式
maxy = y_[indexb]
cons = maxy - maxb * indexb

#式は y = maxb * x + con

yda = maxb .* (0:250) + cons .* ones(251)

p = plot([0:tm],[y_],title="STEP",label=["y"])
 plot!([0:250],[yda])
display(p)
R = maxb
L = -1 * cons / maxb

Kp = 1.2/(R*L)
Ti = 2*L
Td = 0.5*L

println("Kp=",Kp,"Ti=",Ti,"Td=",Td)