プラントの理論モデルに使用する式には注意しよう（戒め）

以下の記事に使用した計算モデルに間違いがあったため修正と、教訓のために記事を残しておきます。

otepipi.hatenablog.com

何を間違ったか

プラントモデルに使用する数式を間違えました 上記記事では以下の例題を扱っておりました。

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

ニュートンの冷却法則の基礎式は以下になります。

f:id:Otepipi:20190414195757p:plain

上記記事ではこれを積分した以下の数式を用いていましたが、これが間違いでした。

f:id:Otepipi:20190414195919p:plain

今回の例題では制御入力が$T{m}$なので、$T{m}$は時間変化することになるためこの積分は使用できません。微分形のまま用いいるべきでした。 つまり、物体の熱容量をCとすると以下の関係が成り立ちますので

f:id:Otepipi:20190414200327p:plain

以下の式で各計算ステップごとの⊿Tを計算して、前ステップの値に足していくという積分処理をするべきでした。　 f:id:Otepipi:20190414200445p:plain

入力と出力の関係が比例か、微分か、積分かという制御工学の基礎的な部分を疎かにしておりました。

　修正したコード

修正したコードは以下になります。

using Plots

tm = 5000 #シミュレーション時間
ie = 0 #誤差積分値初期値

function PID(err,prerr)
    kp = 1 #比例ゲイン
    Ti=20 #積分時間
    Td=5　#微分時間
    ki = kp/Ti #積分ゲイン
    kd = kp*Td #微分ゲイン
    global ie = err + ie 
    ⊿err = err - prerr
    kp * err + ki * ie + kd * ⊿err
end

#目標値の変更
function setpoint(simtime)
 
    return 50

end

function Plant(ypre,input,simtime)
   td = 100 #Delay
   T0 = 0 #温度初期値（℃）
   temp_max = 100 #温度最大値
   temp_min = 0 #温度最小値
   u = max(input,temp_min)
   u = min(u,temp_max)
   
    if simtime < td
        ⊿T = 0
    else
        ⊿T =  (u - ypre)*0.005*1
    end

    return ypre + ⊿T

end

# 空の配列の作成
u = zeros(tm + 1)
y = zeros(tm + 1)
r = zeros(tm + 1)
er = zeros(tm + 1)

for  t in 1:tm
      prer = r[t] - y[t]
      y[t+1] = Plant(y[t],u[t],t)
      r[t+1] = setpoint(t)
      er[t+1] = r[t+1] - y[t+1]
      u[t+1] = PID(er[t+1],prer)            
end

p = plot([0:tm],[y,r],title="kp=1 Ti=20 Td=5",label=["y","setpoint"])
display(p)