システムとモデリング

modelica, Julia, Design Structure Matrix, SysML, 他モデリング全般について。

捕食者の個体数を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)

制御しないとどうなるか

制御しない場合の挙動を示します。繰り返しのループ構造になります。

f:id:Otepipi:20190421001829p:plain

P制御のみの場合

P制御のみの場合、振動しつつゆっくりと目標値に向かっていきます。 kp=2としています。

f:id:Otepipi:20190421002103p:plain

ただ時間を長くしてシミュレーションしてみたところ、目標値には一致しないところで振動してしまうことがわかりました。

f:id:Otepipi:20190421002333p:plain

PI制御の場合

PI制御の場合、出力は発散してしまいました。 下図はKp = Ki = 2 です。

f:id:Otepipi:20190421002448p:plain

I制御はやめたほうが良いと考えました。

PD制御の場合

PD制御で、ようやく出力は目標値と一致するようになりました。 下図はKp = 2, Kd = 50 です。

f:id:Otepipi:20190421002756p:plain

もう少し制御工学に習熟したら根軌跡など理論的な取扱も検討してみたいのですが 今回はここまでにします。