p5.jsとDjangoでインタラクティブシミュレーション(2) 簡単な離散事象シミュレーションの実装からアニメーションまで

はじめに

筆者のラボでは，人を含むシステムの挙動を分析するために，またそうしたシステムの中に置かれた人の挙動を分析するために，ユーザが途中で介入できる形式のコンピュータシミュレーション（インタラクティブシミュレーション）をウェブ上のゲームとして開発し，よく利用している．

このドキュメントは，そうしたインタラクティブシミュレーションのアプリケーションをp5.jsとDjangoを用いて開発するための基礎を身につけてもらうことを狙いとしたもので，あまり一般的なニーズはないかもしれないが，もし多少でもどなたかの参考になれば幸いだ．

今回は全4回中の2回目で，全体のコードはまとめてGitHubに置いた．

簡単な離散事象シミュレーションの実装からアニメーションまで

今回はごく簡単な在庫管理モデルを例にして，離散事象シミュレーションの実装とそのアニメーション化に取り組んでみよう．最初に到達点のイメージを貼り付けておく（いつか動画に差替えたい）．

離散事象シミュレーションとは

ここでいうシミュレーションとは，何らかのシステムを対象として取りあげ，その挙動をコンピュータ上で再現するものである．対象とするシステムとしては，例えば，工場などの生産システム，輸送や物流のシステム，サプライチェーン，サービス業の店舗などが考えられ，挙動を再現するとは，対象システムの状態が時間とともにどのように変化していくかを追いかけることを指す．

ここでは特に，離散事象シミュレーションというフレームワークを考える．これは，対象システムの（興味の対象となる本質的な）状態は，何かのイベントが生じたときに（のみ）変化し，何もイベントが生じない間は不変であると考えるシミュレーションのフレームワークである．

このフレームワークに沿ってシミュレーションモデルを開発するためには，

1. 対象システムの状態をどのように表現するか
2. 状態に変化をもたらすイベントにはどのようなものがあるか
3. 個々のイベントはどのようなタイミングで生起するか
4. 個々のイベントが生起した際に状態はどのように変化するか

の4つを明確にしておく必要がある．

簡単な例として，ある在庫管理モデルを考えよう．ランダムにやってくる顧客に在庫の中から商品を1個ずつ販売していく．このとき，もし在庫が空だと欠品で機会損失となる．在庫量が減り，ある値（発注点という）になれば，所定の量の商品を発注する．そして，そこからある所定の時間（リードタイムという）が経過すると発注した商品が届き，その分だけ在庫量が増える．

まず1については，在庫量のほか，発注（してまだ納品されていない）量，欠品の回数などが挙げられる．2は，顧客の来店，発注，納品の3種類である．3は，顧客の来店はランダムなので，ある顧客が来店してから次の顧客が来店するまでの時間間隔を指数分布に従う確率変数としてモデル化すればよい（一般に，ランダム事象の発生間隔は指数分布に，単位時間内の発生回数はポワソン分布に，それぞれ従う）．

発注は，在庫量が減っていき，それが発注点に等しくなった時点で生起する．在庫量が減るのは顧客が来店したときなので，発注は，来店と同時に（条件付きで）生起することになる．また，納品は，発注からリードタイム分の時間が経過した時点で生起する．

最後は4についてである．まず来店では，在庫が空でなければ在庫量が1個だけ減少し，在庫が空であれば欠品回数が1回増えることになる．発注では，発注量が所定の量だけ増える．納品では，逆に発注量が所定の量だけ減少し，同時に在庫量が同じだけ増える．

シミュレーションモデルのコーディング

続いて，このモデルに基づいて，p5.jsでシミュレーションのコードを書いていこう．離散事象シミュレーションでは，これから生起するイベントを時刻順に並べたイベントカレンダを用意しておくと便利である．これは，例えば，次のように作成できる．

// event calendar
function Calendar(initial) {
  this.events = initial;
}

// add a new event to the calendar
Calendar.prototype.extend = function(e) {
  for(var i = 0; i < this.events.length; i ++) {
    if(this.events[i].time > e.time) {
      this.events.splice(i, 0, e);
      return;
    }
  }
  this.events.push(e);
}

// get the next event from the calendar
Calendar.prototype.fire = function() {
  var e = this.events[0];
  this.events.shift();
  return e;
}

Calendar()は，カレンダオブジェクトのコンストラクタである．オブジェクト指向でプログラミングする場合，javaScriptでは，クラスは書かずに，コンストラクタのみを定義する．そして，この例だと，

var my_cal = new Calendar();

のようにしてオブジェクトを生成することができる．

なお，eventsはイベントを生起時刻の順に並べたリストで，initialはそのリストの初期値である．リスト内の個々のイベントには少なくとも生起時刻の情報（time）をもたせておく．

オブジェクトのメソッドは，通常は上記のようにprototypeの中に定義していけばよい．extend()は新しいイベントをカレンダに追加するメソッド，fire()は次に生起するイベントをカレンダから取り出すメソッドである．

次に，このカレンダも利用しつつ，対象とする在庫管理モデルのオブジェクトを作る．

function Model() {
  this.par = {
    MTB: 2,  //mean time between shipments
    LT: 6,  // lead time to replenishment
    OQ: 20,  // order quantity
    OP: 4  // replenishment point
  };
  this.state = {
    time: 0,  // what time is it now?
    vol: this.par.OQ,  // stock volume at hand
    ordered: [],  // quantities ordered
    outs: 0  // number of stockouts
  };
  this.calendar = new Calendar([
    {time:exp_rand(1 /this.par.MTB), type:"ship_out"},
    {time:900, type:"over"}
  ]);
  this.stateLog = [this.state];
}

Model()は，モデルオブジェクトのコンストラクタである．parはモデルのパラメータで，来店の平均間隔（MTB），リードタイム（LT），1回あたりの発注量（OQ），発注点（OP）の値を定めている．stateは，対象システムの状態を表す変数群である．時刻（time），在庫量（vol），発注量（ordered），欠品回数（outs）を含んでいることがわかる．

これらのほか，上で作成したカレンダを1つ生成し，状態推移の履歴をリストとして保持するためのstateLogも定義してある．カレンダには，最初の来店イベントと，シミュレーション終了のイベントが初期値として登録してあることもわかる．なお，指数分布に従う乱数を生成する関数が見当たらなかったので，下記のように自作した．

function exp_rand(lambda) {
  var x = -log(1 - random()) /lambda;
  return x;
}

メソッドもいくつか定義しておこう．

Model.prototype.reduce = function() {
  if(this.state.vol > 0) {
    this.state.vol --;
    var total_ordered = 0;
    for(o of this.state.ordered) {
      total_ordered += o;
    }
    if(this.state.vol +total_ordered <= this.par.OP) {
      this.calendar.extend({  // a new order is issued
        time: this.state.time +this.par.LT,
        type: "refill"
      });
      this.state.ordered.push(this.par.OQ);
    }
  } else {
    this.state.outs ++;
  }
  this.calendar.extend({
    time: this.state.time +exp_rand(1 /this.par.MTB),
    type: "ship_out"
  });
}

Model.prototype.raise = function() {
  this.state.vol += this.state.ordered[0];
  this.state.ordered.shift();
}

Model.prototype.update = function() {
  var e = this.calendar.fire();
  this.state = Object.assign({}, this.state);
  this.state.ordered = this.state.ordered.concat();
  this.state.time = e.time;
  if(e.type == "over") {
    noLoop();
  } else if(e.type == "ship_out") {
      this.reduce();
  } else if(e.type == "refill") {
    this.raise();
  }
  this.stateLog.push(this.state);
}

update()で次のイベントを生起させている．イベントのタイプがoverであればnoLoop()でdraw()のループを止め，ship_outであればreduce()を，refillであればraise()をそれぞれ呼んでいることがみてとれる．

なお，update()メソッドの3行目と4行目はそれぞれオブジェクトと配列をdeep copyするためのトリックである．

reduce()では，もし在庫量（vol）が0より大きければ，volの値を1つ減らし，発注するかどうかの条件判定を行っている．そして，発注することになった場合は，それに対応する補充イベントをカレンダに追加している．

一方，在庫量（vol）が0のときは，欠品回数（outs）の値を1つ増やしている．また，いずれの場合にも，次の来店イベントをカレンダに追加していることもがわかる．

raise()では，未補充のうち最も古い発注に対応する発注量をorderedのリストから取り出し，その量だけ在庫量（vol）を増加させている．

シミュレーションの実行

これでモデルの骨格ができたので，実際に離散事象シミュレーションを走らせてみよう．setup()とdraw()の部分を次のように書いてみる．

var my_model;

function setup() {
  var my_element = select("#mysketch");
  var my_width = min(my_element.width, 1000);
  var my_canvas = createCanvas(my_width, my_width *0.6);
  my_canvas.parent("mysketch");

  my_model = new Model();
}

function draw() {
  my_model.update();
  console.log(my_model.state);
}

draw()のループで毎回，モデルをupdate()し，その状態（state）をconsole.log()で書き出している．コンソールを開くと，シミュレーションで状態が推移していく様子が（数値として）見えるはずだ．

ここで，setup()やdraw()の外でモデルオブジェクトの変数を

var my_model;

と宣言していることに気をつけてほしい．これは，setup()とdraw()の双方で参照したい変数は大域変数として位置付けておく必要があるからである．もし，この宣言を省いて，setup()の中に（のみ），

var my_model = new Model();

と記述したとすると，draw()からmy_modelにアクセスできなくなる．一方で，今の場所のまま，

var my_model = new Model();

でオブジェクトの生成まで済ませてしまえばいいと思うかもしれないが，これも得策ではない．p5.js固有の関数は，setup()やdraw()の外では使えないからである．

シミュレーション結果のアニメーション化

さて，最後に，結果がわかりやすいようにアニメーションにしてみよう．描画用のメソッドをモデルに追加する．

Model.prototype.show_history = function() {
  var my_ratio = width /1000;
  push();
  scale(my_ratio);
  translate(50, 550);
  textSize(20);
  text("0", -5, 20);
  text("200", 180, 20);
  text("400", 380, 20);
  text("600", 580, 20);
  text("800", 780, 20);
  text("time", 890, 20);
  stroke(0, 0, 255);
  for(var i = 0; i < this.stateLog.length -1; i ++) {
    var from = this.stateLog[i];
    var to = this.stateLog[i +1];
    line(from.time, -5 *from.vol, to.time, -5 *from.vol);
    line(to.time, -5 *from.vol, to.time, -5 *to.vol);
  }
  stroke(255, 0, 0);
  line(frameCount, 0, frameCount, -150);
  stroke(0);
  strokeWeight(2);
  line(0, 0, 900, 0);
  line(0, 0, 0, -200);
  pop();
}

Model.prototype.show_stock = function() {
  var my_ratio = width /1000;
  push();
  scale(my_ratio);
  translate(50, 200);
  textSize(20);
  text("stock at hand", 100, 70);
  for (var i = 0; i < this.state.vol; i ++) {
    rect((i %10) *40, -floor(i /10) *40, 40, 40);
  }
  translate(500, 0);
  text("stockouts", 100, 70);
  fill(255, 0, 0);
  for (var i = 0; i < this.state.outs; i ++) {
    rect((i %10) *40, -floor(i /10) *40, 40, 40);
  }
  pop();
}

描画の詳細の説明は省略するが，show_history()では，在庫量の推移を折れ線グラフで描画しており，show_stock()では，在庫量と欠品回数をそれぞれ白と赤の四角形の個数で可視化している．

アニメーションは時間の流れに沿って動いた方がわかりやすいため，draw()にも少し修正を加えよう．

function draw() {
  background(200);
  my_model.show_history();
  my_model.show_stock();
  while(frameCount > my_model.calendar.events[0].time) {
    my_model.update();
  }
}

frameCountで時間の流れを捉え，モデルをupdate()するタイミングをwhileループで調整している．

setup()の方で，frameRate()もゆっくり目に調整したほうが見やすいだろう（例えば，frameRate(6)ぐらいに設定するとよい）．

これでアニメーションまで完成した．

おわりに

第2回はここで終了．第3回に続く．