捕食者の個体数をPID制御する - システムとモデリング

今回は前回に引き続き、PID制御とはあまり関係のなさそうな分野で無理やりPID制御を導入してみます。

前回の記事は以下

otepipi.hatenablog.com

今回は生物の捕食ー非捕食関係下の個体数の変動を取り扱います。捕食者、非捕食者の個体数の関係はロトカ・ヴォルテラ方程式により以下のように示されます。

$$ \frac{dx}{dt} = (a - by)x \\ \frac{dy}{dt}=(cx - d)y $$ ここで非捕食者の個体数を$x$、捕食者の個体数を$y$としています。

さて、今回は非捕食者の数をPIDで調整することによって、捕食者の個体数を一定に留めるようにしてみます。このときの調整数を$u$として、また$a,b,c,d$のパラメーターを適当に入れた以下の式を今回は使用します。

$$ \frac{dx}{dt} = (1 - 0.01y)x + u \\ \frac{dy}{dt}=(0.02x - 1)y \\ x(0) = y(0) = 20 \\ -10≦u≦10　の範囲でPIDにより調整 $$

コード

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

using Plots

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

function PID(err,prerr)
         kp =1 #比例ゲイン
         ki=0 #積分ゲイン
         kd=100　#微分ゲイン
       
    global ie = err + ie 
    ⊿err = err - prerr
    kp * err + ki * ie + kd * ⊿err
end

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

end

function Plant(y,x,input,simtime)
   u_max = 10 #非捕食者の最大投入数
   u_min = -10 #非捕食者の最大屠殺数
   u = max(input,u_min)
   u = min(u,u_max)
   println(u)
    ⊿x = ((1-0.01*y)*x + u)*0.01 
    ⊿y = ((-1 + 0.02*x)*y)*0.01

    hosyoku = y + ⊿y
    hihosyoku = x + ⊿x

    if y + ⊿y <0
        hosyoku = 0
    end

    if x + ⊿x < 0
        hihosyoku = 0
    end

    return hosyoku, hihosyoku

end

# 空の配列の作成
u = zeros(tm + 1);
x = zeros(tm + 1);  #非捕食者の個体数
y = zeros(tm + 1);　#捕食者の個体数
r = zeros(tm + 1);
er = zeros(tm + 1);

y[1] = 20 #捕食者の初期値
x[1] = 20 #非捕食者の初期値

for  t in 1:tm
      prer = r[t] - y[t];
      y[t+1],x[t+1] = Plant(y[t],x[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,x,r],label=["y","x","setpoint"])
display(p)