CPLEXサンプル 航空交通管理 shed_aftm を読み解く

# 概要
CPLEXのサンプルの[航空交通管理](https://www.ibm.com/docs/ja/icos/20.1.0?topic=library-transportation#usroplsamples.uss_opl_modlib.1016594__section1205855103301)（CPLEX OPL CPLEXサンプル CP）を読み解きます。

オリジナルのサンプルは
<サンプル導入ディレクトリ> examples/opl/models/sched_atfm
にあります。

## コード
https://github.com/kuru-to/air_traffic_flow_management

`docplex` を用いて python で実行できるようにしています。

### OPL
https://www.ibm.com/docs/ja/icos/20.1.0?topic=SSSA5P_20.1.0/ilog.odms.ide.help/examples/html/opl/models/sched_atfm/sched_atfm.mod.htm

### サンプルデータ
https://github.com/kuru-to/air_traffic_flow_management/blob/main/data/example/sched_atfm.dat

#### 実行環境
- CPLEX 22.1.1
- Windows 11 64bit
- Python 3.10.11
- docplex 2.28.240

# 問題
米国では、重要な航空区域での混雑を避けるため、地域を区画に分け、各区画に入ることができる航空機の数に制限を設けている。
（[参考論文](https://www.mit.edu/~dbertsim/papers/AirTransportation/The%20air%20traffic%20flow%20management%20problem%20with%20enroute%20capacities.pdf)より抜粋）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3873809/e1840972-4c9d-26f0-7df8-33be6eaf85e9.png)

また、航空機（今回の問題では旅客機/貨物機両方を想定）には離発着空港とルートが決まっており、出発予定時刻をもとに、何時にどの区画へ進入するか（以下、進入イベントと表記する）が決まっている。
（[参考論文](https://www.mit.edu/~dbertsim/papers/AirTransportation/The%20air%20traffic%20flow%20management%20problem%20with%20enroute%20capacities.pdf)より抜粋）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3873809/771df24e-9012-2dbd-bcbc-928c1e8738d7.png)

そのため、航空機が進入イベントの予定時刻において区画の航空機数上限に達している場合、その航空機は対象区画に入るのを待たなければならない。ここで、空中で待っていると燃料コストがかかったり、天候変動のようなリスクがあったりするので、地上で出発を遅らせることで対応する。

最終的な目的は、各区画での航空機数制限を満たしたうえで、到着予定からの遅延を少なくするような航空機の進入計画を作成すること。

## 問題の種類
制約プログラミング（スケジューリング）

## データ
OPLコードとともに詳しく解説していきます。

### 時刻
```
// times are specified in hours, minutes, seconds 
// (in general year, month, day are needed as well)
tuple Time {
    int hours;
    int minutes;
    int seconds;
};
```
本問題において時刻を表す。本来なら年月日も入るが、今回のモデルにおいては使用しない。

### period
```:sched_aftm.mod
// limitations of air capacity are expressed in terms of regulated
// periods, i.e. an interval with an hourly capacity rate
tuple Period {
   Time start;
   Time end;
   int rate;
};

{Period} periods[SectorNames] = ...;
```

区画ごとに、ある期間における1時間あたりの進入可能な航空機機数を表す。つまり、`start`~`end` において 1時間あたり `rate`機が、対象の区画において進入可能な航空機数となる。
今回のモデルは10分区切りでどのフライトがどの区画へ進入するかを決定するため、10分あたりの機数に変換する必要がある（後述）。

period が設定されていない時間帯は上限がないものとする（つまり、19:00~20:00 で `rate`=10だった場合、19:00~20:00 は10機までしか進入できないが、00:00~19:00, 20:00~24:00 は何機進入してもok）。

サンプルデータでは以下のように与えられている。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3873809/ec10b5d2-a869-f69e-4b44-71db0021293b.png)


### 進入イベント
```:sched_aftm.mod
// an enter event specifies that a given flight will enter a given sector
// at an expected time (called expected time over)
tuple Enter {
   int flight;
   string sector;
   Time eto;
};
```

どのフライトが、どの区画に、どの時刻で進入する予定(`eto`= expected time over)か。

サンプルデータでは以下のように与えられている（一部抜粋）。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3873809/e2e4634f-866c-539a-470b-59d68a53bff7.png)


## 各種定数
1. 最大遅延時間
1. 何分区切りで最適化を行うか

### 最大遅延時間
```:shed_atfm.mod
// flight delays will be limited to 2 hours
int maxDelay = 120;
```

航空機を遅らせることができる時間（分）。例題のモデルでは 120分（2時間）。

### 何分区切りで最適化を行うか
```:shed_atfm.mod
// capacity of the resource will be made available by time steps of 10 minutes
int timeStep = 10;
```

今回のモデルでは、10分区切りで航空機の進入計画を作成するものとする。

## 決定変数
1. 各フライトの遅延
1. 各進入イベントにおける計画上の進入時刻・時間

また、決定式として与えられる「各区画における各 timeStep での進入イベント数」についても本項で解説する。

### 各フライトの遅延
```:shed_atfm.mod
// flight delays are expressed by integer variables
dvar int delay[Flights] in 0 .. maxDelay;
```

各フライトの遅延（分）を表す。0 以上 `maxDelay` （最大遅延時間）以下の整数変数。

### 各進入イベントにおける計画上の進入時刻・時間
```:shed_atfm.mod
// each enter event is modelled by an activity of duration 1
dvar interval a[Enters] size 1;
```

各進入イベント（どのフライトがどの区画へいつ進入する予定か）において、計画上で進入開始する時刻とその時間を表す、間隔変数型。

「間隔変数型」[interval](https://www.ibm.com/docs/ja/icos/20.1.0?topic=keywords-interval)は、スケジュールの制約プログラミングに特化した変数です。間隔変数は、開始値、終了値、サイズなどの情報をもったクラスのような変数です。CPLEXではこのようなスケジューリングに特化した変数を用意することで、スケジュールの制約プログラミングをわかりやすく、かつ速く解けるようにしています。

本モデルでは、 `size 1` としているため、各航空機は 1 timeStep(=10分)の進入イベントとなる。

サンプルデータを解いてみると以下のような解が出力される。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3873809/7400ec4f-9b7d-7e2b-70a8-7308c3a5ff5c.png)
横軸が何 timeStep 目かを表し、サンプルデータで与えられた進入イベントが順に縦に並んでいる（表示1行目では、サンプルデータの進入イベントのうち、1番目は 33 timeStep目 = 05:30 に進入する計画が出力されている）。

### 各区画における各 timeStep での進入イベント数
```:sched_atfm.mod
// each sector is modelled by a resource
cumulFunction r[i in SectorNames] = sum(en in Enters : e[en].sector == i) pulse(a[en], 1);
```

[pulse 関数](https://www.ibm.com/docs/ja/icos/20.1.0?topic=functions-pulse)の説明は以下の通り

> Pulse は、個々の間隔変数または固定の時間間隔のコントリビューションを表す基本的な累積関数式です。Pulse は、アクティビティーの開始時にリソース使用量関数が増加し、その終了時にリソースが解放されて使用量が減少する際の、累積リソースまたは再生可能リソースの使用量をカバーします。
pulse 関数の間隔は、a によって、または開始時点 u と終点時点 v で表されます。関数の高さは h で表されるか、hmin と hmax で境界されます。

今回の例では、`pulse`関数によって、進入イベントにおける計画上の進入開始時刻~終了時刻において 1 が立つようになっている。
そして、`sum` 関数で合計値をとっているため、上記 `r` は各区画の各 timeStep において何機の進入が計画されているか、を表す決定式となっている。

サンプルデータを解いてみると以下のような解が出力される。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3873809/a920d17a-daf4-068d-77e8-1cee411fd1f6.png)
各 timeStep で何機の進入が計画されたかを、各区画で色分けしている。

## 目的関数
```:sched_atfm.mod
dexpr int totalDelay = sum(i in Flights) delay[i];

minimize totalDelay;
```
合計遅延分数の最小化

## 制約
1. 10分間隔での進入航空機上限を下回るようにする
1. 進入開始時刻は進入イベントの予定時刻+遅延

### 10分間隔での進入航空機上限を下回るようにする
```:sched_atfm.mod
// the capacity rate is adapted to intervals of 10 minutes;
// the time scale of a resource is divided by the time step
forall (i in SectorNames)
    forall (p in periods[i])
        alwaysIn(
            r[i],
            (p.start.hours * 60 + p.start.minutes) div timeStep,
            (p.end.hours * 60 + p.end.minutes) div timeStep,
            0,
            (p.rate * timeStep + 59) div 60
        );
```
(読みやすくするため、適宜改行をいれております)

[alwaysIn 関数](https://www.ibm.com/docs/ja/icos/20.1.0?topic=functions-alwaysin)の説明は以下の通り。

> この制約は、cumulFunction または stateFunction f の値を、負でない範囲 [hmin, hmax] に限定するために使用されます。この関数は、間隔 [u, v) のすべての範囲、または間隔 a の開始から終了までの間 (存在する場合) で、この範囲に制限されます。

また、`div` は割り算の商を返す演算子（Python でいう`//`）であるので、
```
(p.start.hours * 60 + p.start.minutes) div timeStep
```
は `period` の開始時刻が 何`timeStep`目なのかを表す（07:10 であれば `(7 * 60 + 10) div 10 = 43`. 07:19 も切り捨てられて 43番目となる）。`period` の終了時刻についても同様。
また、
```
(p.rate * timeStep + 59) div 60
```
は10分あたりの進入可能な航空機機数を表す。`p.rate` が1時間あたりの航空機上限機数なので、`timeStep` で積をとったあとに60（分）で割るが、この時に上限機数を過小評価しないように切り上げのための59を足している（例えば、`rate=10`（1時間あたり10機が進入可能）とした場合10分あたりの進入可能な機数は、切り上げがなければ `10 * 10 dev 60 = 1`. しかしこの場合、60分あたりの機数に換算すると6機になるため、本来の1時間あたり10機より少ない。そのため、59を足して切り上げして `(10 * 10 + 59) div 60 = 2` としている）。
ただしこの場合、**1時間あたりの航空機上限機数を過大評価することになる**（先の例では10分あたり2機となったので、各 `timeStep` で2機進入する計画が出力されうる。これは1時間当たり12機進入する計画となるため、本来の上限機数10機を超えてしまう）。

以上より本制約は、`period`が設定されている区画において、「各区画の各 timeStep において何機の進入イベントがあるか」を表す `r` が `alwaysIn` 関数によって、10分あたりの進入可能な航空機機数までに収める制約となる。

### 進入開始時刻は進入イベントの予定時刻+遅延
```:sched_atfm.mod
// a flight enters a sector at its expected time-over plus its delay;
// since the time scale of a resource is divided by the time step,
// we do the same for the start time of the activity
forall (i in Enters)
    startOf(a[i]) == (delay[e[i].flight] + e[i].eto.hours * 60 + e[i].eto.minutes) div timeStep;
```
本制約は、間隔変数 `a` の開始点（何`timeStep`目か）が、進入イベント`e`の expected time over に遅延分数 `delay` を足した場合の `timeStep`目と一致することを表す。
予定時刻より早く進入することはない（予定時刻が 07:19 であれば`timeStep`の区切りである 07:10 に進入する計画となるかもしれないが、07:00 に進入とはならない）。

# サンプルデータの求解結果
## 諸情報
区画数: 17
フライト数: 1991
進入イベント数: 3783

詳細情報は [sched_atfm.dat](https://github.com/kuru-to/air_traffic_flow_management/blob/main/data/example/sched_atfm.dat) 参照。

## 結果
### 各 period での enter_event 数
![num_flows_and_upper.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3873809/2ddc368f-b654-7c31-ce74-98a5fbeb501f.png)
`num_enter`が最適化の結果得られた、各`period`の期間における進入航空機数。`upper`は、`period`期間全体で何機まで進入可能なのかを`rate`をもとに算出した値（小数点以下含む）。
`upper` は常に `num_enter` より大きくなければならないが、[10分間隔での進入航空機上限を下回るようにする](#10分間隔での進入航空機上限を下回るようにする)で述べたように、`num_enter`が`upper`をちょっと超える `period` が存在する。

### delay に関するヒストグラム
![delay_histogram.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3873809/eaecc632-d542-f501-baac-0034262ebb30.png)
maxDelay=120 なので2時間以内に収まっていればok、ちゃんと収まっている


## 考察
### モデルの改良案
- 10分間隔での進入可能な航空機数を計算する際に切り上げをしているため、各 `period` 期間での進入イベント数が若干超えてしまうモデルになっている。進入可能な航空機数を遵守するのであれば、period の `start`~`end` の全体において `rate`機数を超えない、という制約を追加する必要があると考える
- 本モデルは各区画において進入可能な機数に上限を与えているが、航空機が進入した後もその区画に残り続ける状況を考えていないように思う。航空機がその区画に残り続けたらどんどん他の航空機が進入してきて混雑してしまうのでは？

# 参考
#### 輸送: 航空交通管理 - IBM Documentation
https://www.ibm.com/docs/ja/icos/20.1.0?topic=library-transportation#usroplsamples.uss_opl_modlib.1016594__section1205855103301

#### CPLEX サンプルを読み解く記事一覧
https://qiita.com/spssfun2017/items/bdb85f66e5f1de6039d6

#### 参考論文
https://www.mit.edu/~dbertsim/papers/AirTransportation/The%20air%20traffic%20flow%20management%20problem%20with%20enroute%20capacities.pdf