fromgurobipyimportModel,quicksum,GRB# from mypulp import Model, quicksum, GRBfromcollectionsimportOrderedDict,defaultdictimportnetworkxasnx
様々な最適化のベンチマーク問題例を収集しているサイトとしては、昔からあるOR-Library (http://people.brunel.ac.uk/~mastjjb/jeb/info.html) が有名だ。1990年に開始したこのライブラリには、論文で実験された問題のデータが配布されている。昔は、専用のアルゴリズムを作成してやっと解けていた問題が、今では数理最適化ソルバーを使って比較的簡単に解けるようになっている。簡単な練習問題ではなく、もっと難しい問題に挑戦したい人には、このライブラリはうってつけだ。
ここでは、以下のシフト最適化問題を解いてみる.
複数のジョブを作業員に割り当てたい。ジョブには開始時刻と終了時刻が与えられており、ジョブの処理には1人の作業員が必要で、同時に2つのジョブの処理はできないものとする。作業員には固定費用があり、使用する作業員の総固定費用を最小化せよ。
もとになった論文やデータは以下ののサイトを参照されたい。
http://people.brunel.ac.uk/~mastjjb/jeb/orlib/ptaskinfo.html
最初の制約は、すべてのジョブが処理されなければならないことを表す。2番目の制約は、クリークに含まれるジョブは、高々1人の作業員にしか割り当てられないことと割り当てられた作業員は使われなければならないことを同時に表す。
まず,データを読み込んで準備をする.
fname="../data/ptask/data_10_51_111_66.dat"f=open(fname,"r")lines=f.readlines()f.close()
n_jobs=int(lines[4].split()[2])print("number of jobs=",n_jobs)start,finish=OrderedDict(),OrderedDict()forjinrange(n_jobs):start[j],finish[j]=map(int,lines[5+j].split())
number of jobs= 111
n_workers=int(lines[5+n_jobs].split()[2])print("number of workers=",n_workers)Jobs=defaultdict(set)forwinrange(n_workers):line=lines[6+n_jobs+w].split()forjinline[1:]:Jobs[w].add(int(j))
number of workers= 51
2つのジョブの開始時刻と終了時刻に重なりがあるかどうか判定する関数を作っておく.
defintersect(j,k):ifstart[j]>start[k]:j,k=k,jiffinish[j]>start[k]:returnTrueelse:returnFalse
ジョブを点とした交差グラフを生成する.このグラフのクリーク(完全部分グラフ)が,同時に処理できないジョブの集合になる.
G=nx.Graph()forjinrange(n_jobs):forkinrange(j+1,n_jobs):ifintersect(j,k):G.add_edge(j,k)
networkXを用いて極大クリークを列挙しておく.
Cliques=[]forcinnx.find_cliques(G):Cliques.append(set(c))print("Number of cliques=",len(Cliques))
Number of cliques= 32
model=Model()x,y={},{}forwinrange(n_workers):y[w]=model.addVar(vtype="B",name=f"y[{w}]")forjinJobs[w]:x[j,w]=model.addVar(vtype="B",name=f"x[{j},{w}]")model.update()forjinrange(n_jobs):model.addConstr(quicksum(x[j,w]forwinrange(n_workers)if(j,w)inx)==1,name=f"job_assign[{j}]",)forwinrange(n_workers):forj,cinenumerate(Cliques):model.addConstr(quicksum(x[j,w]forjinJobs[w].intersection(c))<=y[w],name=f"connection[{j},{w}]",)model.setObjective(quicksum(y[w]forwinrange(n_workers)),GRB.MINIMIZE)model.optimize()
Academic license - for non-commercial use only - expires 2023-06-24Using license file /Users/mikiokubo/gurobi.licGurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)Thread count: 8 physical cores, 16 logical processors, using up to 16 threadsOptimize a model with 1743 rows, 3943 columns and 48251 nonzerosModel fingerprint: 0x6315480aVariable types: 0 continuous, 3943 integer (3943 binary)Coefficient statistics: Matrix range [1e+00, 1e+00] Objective range [1e+00, 1e+00] Bounds range [1e+00, 1e+00] RHS range [1e+00, 1e+00]Found heuristic solution: objective 51.0000000Presolve removed 449 rows and 0 columnsPresolve time: 0.08sPresolved: 1294 rows, 3943 columns, 36453 nonzerosVariable types: 0 continuous, 3943 integer (3943 binary)Root relaxation: objective 4.000000e+01, 2652 iterations, 0.17 seconds Nodes | Current Node | Objective Bounds | Work Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time* 0 0 0 40.0000000 40.00000 0.00% - 0sExplored 0 nodes (3964 simplex iterations) in 0.36 secondsThread count was 16 (of 16 available processors)Solution count 2: 40 51 Optimal solution found (tolerance 1.00e-04)Best objective 4.000000000000e+01, best bound 4.000000000000e+01, gap 0.0000%
print("Optimal value=",model.ObjVal)
Optimal value= 40.0
典型的な看護婦スケジューリング問題(nurse scheduling problem)では,以下の制約が付加される.
このような付加条件が課された看護婦スケジューリング問題は,制約最適化ソルバーSCOPで簡単に定式化できる.
看護婦 $i$ の 日 $d$におけるシフトを表す変数 $X_{id}$ (領域は昼,夕,夜,休暇の集合)を用いる.
対応する値変数 $x_{ids}$ は $0$-$1$ 変数であり, 看護婦 $i$ の $d$ 日のシフトが $s$ のとき $1$ になる.
制約を式で表現する.
日 $d$ のシフト $s$ の必要人数を $L_{ds}$ とする.$$ \sum_{i} x_{ids} \geq L_{ds} \ \ \ \forall d,s$$
休暇以外のシフト集合を $W$,下限と上限を $LB,UB$ とする.$$ LB \leq \sum_{d, s \in W} x_{ids} \leq UB \ \ \ \forall i$$
看護婦 $i$ が休日を希望する日の集合を $R_i$ とする.$$\sum_{d \in R_i, s \in W} x_{ids} \leq 0 \ \ \ \forall i$$
$d$ 日から始まる連続7日間の集合を $C7_d$ とする.休日(値は $r$) の制約だけ記述する.$$ \sum_{t \in C7_d} x_{itr} \geq 1 \ \ \ \forall i, d $$
例として3連続夜勤の場合の制約を考える.夜勤の値を $n$, $d$ 日から始まる連続3日間の集合を $C3_d$ とする.
$$ \sum_{t \in C3_d} x_{itn} \leq 2 \ \ \ \forall i, d $$夜勤以外のシフトの集合を $N$ とする.$$ \sum_{s \in N} x_{ids} + x_{i,d+1,n} + \sum_{s \in N} x_{i,d+2,s} \leq 2 \ \ \ \forall i, d $$
1つのチームに含まれる看護婦の集合を $T$ とし,以下の制約を考慮制約とする.$$ \sum_{i \in T} x_{ids} = \sum_{i \not\in T} x_{ids} \ \ \ \forall d, s $$
以下の変数を追加する.
$z_{id}$: 看護婦 $i$ が日 $d$ に開始するシフト; 対応する値変数 $z_{ids}$ は,看護婦 $i$ がシフト $s$ を日 $d$ に開始するとき $1$
$w_{id}$: 看護婦 $i$ が日 $d$ に終了するシフト;対応する値変数 $w_{ids}$は,看護婦 $i$ がシフト $s$ を日 $d$ に終了するとき $1$
看護婦スケジューリングの国際コンペティション(First International Nurse Rostering Competition 2010)では,SCOPは3つの異なるトラックですべてメダル(2位,3位,4位)を獲得している.
実際のシフト最適化においては,複数の日(たとえば1ヶ月)にスタッフを割り当てる必要がある.このとき,各日の各時間帯においてどのような業務をどのスタッフに割り当てるかを決めたい.もちろん,日や時間帯によって,各業務の遂行に必要な人数は異なる.このような諸条件を考慮した最適化は,問題依存であり,一般には難しい.
実務から発生した諸条件を組み込んだシフトシフト最適化システムとして OptShift (https://www.logopt.com/optshift/) が開発されている. これは,制約最適化ソルバー SCOP を用いたものであり,大規模な実際問題に対する実用的な解を,現実的な計算時間で算出する.