時間の入ったスケジュール最適化を試してみた

久保幹雄先生の Python言語による実務で使える100+の最適化問題の時間を元にしたスケジュール最適化を試したのでまとめます。

久保幹雄先生の解説YouTube

https://www.youtube.com/watch?v=UdgWRyStR-U

コードは扱いやすいように以下に手を加えています。

• 計算量軽減のため、time_df 50個からサンプルで30個に絞った
• 使っているデータセットを pandasのDataFrame型に変換
• ノード表示する際に分かりやすいように何番目のデータかどうかノードに重ねてテキストで表示
  
• 最適化結果も pandas の DataFrame で出力。なお、group の数字が同じものが組み合わせたスケジュール
  
• 組み合わせ可能なスケジュール候補 cost_df をロジックで作成するようにした。ただし cost は 適当に time と同じ数値にした

コードは GitHub がこちら、Google Colab がこちら

使用しているデータ

• time_df
  それぞれのタスクがいつから(start)いつまで(finish)なのかをまとめた表です。
  time は start から finish までの経過時間です。
  

• cost_df
  time_df にまとめられたタスクのうち、組み合わせ可能な候補のようです(間違ってたらすみません)。

• イメージ図
  

最適化結果

group が同じ数値のものが組み合わせられたスケジュールです。

index は 元データの time_df の index で、ネットワーク図に表示されている番号と同じものです。

上記の例ですと、

576から708まで(その間132) => 待機 => 711から765まで(その間54) => 766から791まで(その間25)となっています。

参考イメージ図

数値を秒、分や時間に置き換えれば使い勝手がよさそうですね。

パス型定式化すると組み合わさったもの以外が消えてしまったのですが、ネットワークでも表示することができます。

以上になります、最後までお読みいただきありがとうございました。

コード全文

## ライブラリのインポート
GUROBI = True
if GUROBI:
    from gurobipy import Model, quicksum, GRB
else:
    from mypulp import Model, quicksum, GRB
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from collections import OrderedDict, defaultdict
import networkx as nx
pd.set_option('display.max_rows', 200)

## シード値の固定
SEED = 42
np.random.seed(SEED)
random.seed(SEED)

## データの読込
fname = "csp50.txt"
f = open(fname,"r")
lines = f.readlines()
n_tasks, time_limit = map(int, lines[0].split())
K = 27

# 組み合わせた後の合計時間の上限を修整
time_limit = 300

## DataFrame処理
### データ① スケジュール
time_df = pd.DataFrame()

for t in range(1,n_tasks+1):
    time_df = pd.concat(
      [time_df, pd.DataFrame(lines[t].split(" ")).T.replace('\n', '', regex=True)]
      )

time_df = time_df[[1, 2]].reset_index(drop=True)
time_df.columns = ["start", "finish"]
time_df = time_df.astype(int)
time_df["time"] = time_df["finish"] - time_df["start"]

# 計算量を減らすために50行から30行に削減
time_df = time_df.sample(30, random_state=SEED)

# indexを1からスタートに変更
time_df.index = np.arange(1, len(time_df)+1)

# n_tasks を変更
n_tasks = len(time_df)

### データ② 組み合わせ候補
cost_df = pd.DataFrame()

for i in range(len(time_df)):
    for j in range(len(time_df)):
        if i==j:
            continue
        else:
             # 時間が被ってないかどうか判定
            if time_df.iloc[i:i+1]["finish"].values[0] < \
                time_df.iloc[j:j+1]["start"].values[0]:

                # 組み合わせたペアのトータル時間がtime_limit を超えない
                if time_df.iloc[j:j+1]["finish"].values[0] \
                  - time_df.iloc[i:i+1]["start"].values[0] <= time_limit:

                    cost_df = pd.concat([
                        cost_df,
                        pd.DataFrame([
                                i+1,# tail
                                j+1, # head
                                time_df.iloc[j:j+1]["start"].values[0] - \
                                   time_df.iloc[i:i+1]["finish"].values[0], # cost
                                time_df.iloc[j:j+1]["start"].values[0] - \ 
                                  time_df.iloc[i:i+1]["finish"].values[0]  # time
                        ]).T                                     
                    ])

cost_df.columns = ["tail", "head", "cost", "time"]

# indexを1からスタートに変更
cost_df.index = np.arange(1, len(cost_df)+1)


## ネットワークの作成
D = nx.DiGraph()
for i in range(len(cost_df)):
    D.add_edge(
          cost_df["tail"].iloc[i], cost_df["head"].iloc[i], 
          weight=cost_df["cost"].iloc[i],
          time=cost_df["time"].iloc[i]
          )
for i in range(len(time_df)):
    D.add_edge(0,i+1, weight=0, time=0)
for i in range(len(time_df)):
    D.add_edge(i+1, len(time_df)+1,weight=0, 
               time=time_df["time"].iloc[i])


def csp_callback(model, where):
    if where !=GRB.Callback.MIPSOL:
        return
    SolGraph = nx.DiGraph()
    for (i,j) in x:
        if model.cbGetSolution(x[i,j])  >  0.1:
            SolGraph.add_edge(i,j)

    for path in nx.all_simple_paths(SolGraph, 
                                    source=0, target=n_tasks+1):
        i = path[1]
        t_sum = time_df["finish"].iloc[i-1] - time_df["start"].iloc[i-1]
        for j in  path[2:-1]:
            t_sum += time_df["finish"].iloc[j-1] - time_df["start"].iloc[j-1]
            t_sum += max(time_df["start"].iloc[j-1] - time_df["finish"].iloc[j-1], 0)
            i = j
        if t_sum > time_limit:
            edges =[]
            i = path[0]
            for j in path[1:]: 
                edges.append( (i,j) )
                i = j
            model.cbLazy(  quicksum( x[i,j] 
           for (i,j) in edges) <=len(edges) -1)   
    return

## 最適化の実行
model = Model()
x = {}
for (i,j) in D.edges():
    x[i,j] = model.addVar(vtype="B",name=f"x[{i},{j}]")
model.update()
for i in range(1,n_tasks+1):
    model.addConstr( quicksum( x[i,j] for j in D.successors(i) ) 
                    == 1, name=f"task_execution[{i}]" )
for i in range(1,n_tasks+1):
    model.addConstr( quicksum( x[j,i] for j in D.predecessors(i) )
                    == quicksum( x[i,j] for j in D.successors(i) ),
                    name=f"flow_conservation[{i}]" )   
model.addConstr(quicksum( x[0,j]  for j in D.successors(0) )==K)
model.setObjective( quicksum(D[i][j]["weight"]*x[i,j] 
                 for (i,j) in x), GRB.MINIMIZE) 
if GUROBI:
    model.Params.DualReductions = 0
    model.params.LazyConstraints = 1
    model.optimize(csp_callback)


## ネットワークの表示
plt.figure(figsize=(15, 8))
alpha = 0.1

SolGraph = nx.DiGraph()
for (i,j) in x:
    if x[i,j].X > 0.1:
        SolGraph.add_edge(i,j)
SolGraph.remove_node(0)
SolGraph.remove_node(n_tasks+1)

pos = nx.spring_layout(SolGraph)

nx.draw(SolGraph, pos=pos)

# 末尾に ; を入れて jupyter にログが出ないようにしている
nx.draw_networkx_labels(SolGraph, pos=pos, font_color='w');

## パス型定式化
def k_th_sp(G, source, sink, k, weight="weight", time_limit=time_limit):
    """
    Find k-th shortest paths and returns path and its cost
    """
    time_list, path_list = [], [] 
    for i, p in enumerate(nx.shortest_simple_paths(G,
           source, sink, weight=weight)):
        if i >= k:
            break
        v = p[0]
        time = 0.
        for w in p[1:]:
            time += G[v][w][weight]
            v = w
        if time>time_limit:
            break
        time_list.append( time )
        path_list.append( tuple(p)) 
    return time_list, path_list

## ネットワークの作成
time_list, path_list = k_th_sp(D,0,n_tasks+1,100000,weight="time")

C = {}
Paths = defaultdict(set)
for p, path in enumerate(path_list):
    cost_ = 0
    for j in range(1,len(path)-1):
        cost_ += D[path[j]][path[j+1]]["weight"]
        Paths[path[j]].add(p)
    C[p] = cost_


## 最適化の実行
model = Model()
X = {}
for p in C:
    X[p] = model.addVar(vtype="B",name=f"x({p})")
model.update()  
for t in range(1,n_tasks+1):
    model.addConstr( quicksum(X[p] for p in Paths[t]) >=1 )
model.addConstr( quicksum( X[p] for p in X) ==K)
model.setObjective( quicksum(C[p]*X[p] for p in X), 
                   GRB.MINIMIZE) 

model.optimize()
print(model.ObjVal)


## ネットワークの表示
plt.figure(figsize=(15, 8))
pair_list = []

SolGraph = nx.DiGraph()
for p in X:
    if X[p].X >0.01:
        pair_list.append(list(path_list[p][1:-1]))
        for j in range(1,len(path_list[p])-2):
            SolGraph.add_edge( path_list[p][j], path_list[p][j+1] )

pos = nx.spring_layout(SolGraph)
nx.draw(SolGraph, pos=pos)

# 末尾に ; を入れて jupyter にログが出ないようにしている
nx.draw_networkx_labels(SolGraph, pos=pos, font_color='w');

## 最適化結果をDataFrameで出力
pair_df = pd.DataFrame()

for i in range(len(pair_list)):
    pair = pair_list[i]

    tmp_df = time_df.iloc[np.array(pair)-1]
    tmp_df["group"]=i
    pair_df = pd.concat([pair_df, tmp_df])

pair_df = pair_df[["group", "start", "finish", "time"]]
pair_df = pair_df.reset_index()