魚の体重を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)
このときの経時変化は以下のようになります。 比例制御のみだと全く目標値に追随しません。
比例 積分制御
比例制御に加えて積分制御も入れてみます。
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)
経時変化のグラフは以下のようになり、目標値に追随できていることがわかります
今回はここまでにします。