魚の体重をPID制御する - システムとモデリング

今回はPID制御を用いて魚の体重をコントロールしてみます。

魚の体重の成長曲線を示すモデルとしてフォン・ベルタンフィーモデルがよく知られています。 $$ \frac{dw}{dt} = \alpha w^{2/3} - \beta w $$ ここで$w$は魚の体重を示します。また右辺第一項は栄養分による体重の増加を、第2項は呼吸による体重ロスを表します。

今回はエサの量≒係数$\alpha$を入力とし、魚の体重$w$をコントロールします。

比例制御のみ

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

using Plots

     tm = 50 #シミュレーション時間
     ie = 0 #誤差積分値初期値
     beta = 1 # βの値

function PID(err,prerr)
         kp =1 #比例ゲイン
         ki=0 #積分ゲイン
         kd=0　#微分ゲイン
    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)
   food_max = 100 #エサの量最大値
   food_min = 0.001 #エサの量最小値
   u = max(input,food_min)
   u = min(u,food_max)
   
    ⊿w = u * ypre^(2/3) - beta * ypre
    case = [u ⊿w]
    println(case)
    return ypre + ⊿w

end

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

y[1] = 10 #体重初期値

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",label=["y","setpoint"])
display(p)

このときの経時変化は以下のようになります。比例制御のみだと全く目標値に追随しません。

f:id:Otepipi:20190421082739p:plain

比例　積分制御

比例制御に加えて積分制御も入れてみます。

using Plots

     tm = 50 #シミュレーション時間
     ie = 0 #誤差積分値初期値
     beta = 1 # βの値

function PID(err,prerr)
         kp =0.01 #比例ゲイン
         Ti=1 #積分時間
         Td=0　#微分時間
         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)
   food_max = 100 #エサの量最大値
   food_min = 0.001 #エサの量最小値
   u = max(input,food_min)
   u = min(u,food_max)
   
    ⊿w = u * ypre^(2/3) - beta * ypre
    case = [u ⊿w]
    println(case)
    return ypre + ⊿w

end

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

y[1] = 10 #体重初期値

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=0.01, ki=1",label=["y","setpoint"])
display(p)