Google Optimization Toolsで組み合わせ最適化
巡回セールスマン問題やナップサック問題は、組み合わせ最適化と呼ばれる最適化問題の分野のひとつで扱われているとのこと。手法の全容は色々本が出ています。
今日から使える!組合せ最適化 離散問題ガイドブック (KS理工学専門書)
- 作者: 穴井宏和,斉藤努
- 出版社/メーカー: 講談社
- 発売日: 2015/06/23
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (1件) を見る
実際の使い道として、生産工程の最適化などがあり、便利なソルバーとしてGoogleからGoogle Optimization Toolsが出ています。ちょっと使ってみました。
インストール
Python 3.X系の場合はpy3-ortools 6.3.4431 : Python Package Indexとして登録されているようなので…
pip install py3-ortools
でインストール。どうやら64bit系でしか動かない様なので、そこはご注意ください。
Python 2.7の場合は、ortools 6.3.4431 : Python Package Indexにあって、
pip install ortools
最新版はGithubにおいてあったりしますが、楽ちんな方法でやりました。
生産工程最適化問題
業界用語(?)で、ジョブショップ問題と呼ぶそうです。Google Optimization Toolsのページにもサンプルがあって、今回はこれを参考にしました。→The Job Shop Problem | Optimization | Google Developers. サンプルの中には若干間違いもあって、Feedback入れておきました。いつ反映されるかなー。
※ 生産工程の最適化についても制約条件や問題設定に応じて、オープンショップ問題、フローショップ問題、リスケジュール問題など色々種類があるそうです。
問題設定
こちらの問題設定をお借りして、3種類の機械(旋盤、ボール盤、フライス盤とかかな)の三種類を使う工程を3種類の部品(JOB)について、最適工程割り付けを行います。
機械と部品を便宜上、数字で呼ぶと扱いやすいので、
- 旋盤
- ボール盤
- フライス盤
としましょうか。工程は下記の通り(言うまでもなく、それぞれの部品の加工順序は守らなくてはなりません)
Job (部品) | 工程1 | 工程2 | 工程3 |
---|---|---|---|
部品0 | 旋盤 3分 | ボール盤 2分 | フライス盤 3分 |
部品1 | 旋盤 2分 | フライス盤 3分 | ボール盤 4分 |
部品2 | ボール盤 3分 | 旋盤 2分 | フライス盤 1分 |
これを、ソルバー用の書き方に書き換えると
# [m, t] : m 機械のID, t 加工時間(分) jobs = [ [[0,3],[1,2],[2,3]], # Job 0 / 部品0 [[0,2],[2,3],[1,4]], # Job 1 / 部品1 [[1,3],[0,2],[2,1]] # Job 2 / 部品2 ]
この問題を解かせます
結果のイメージ
ソルバーの宣言と問題の設定
横着して3次元の配列変数操作のためだけにnumpy使います。
# OR toolsのインポート from ortools.constraint_solver import pywrapcp # 変数操作を楽にするためNumpy import numpy as np # 視覚化用 import matplotlib.pyplot as plt def main(): # ソルバー定義 solver = pywrapcp.Solver('jobshop') # 機械の種類 # 0. 旋盤 # 1. ボール盤 # 2. フライス盤 machines_count = 3 all_machines = range(0, machines_count) # Jobの定義 # [m, t] : m 機械のID, t 加工時間(分) jobs = np.array([ [[0,3],[1,2],[2,3]], # Job 0 / 部品0 [[0,2],[2,3],[1,4]], # Job 1 / 部品1 [[1,3],[0,2],[2,1]] # Job 2 / 部品2 ]) jobs_count = 3 all_jobs = range(0, jobs_count) machines = jobs[:][:][0] processing_times = jobs[:][:][1] # 全てのJobを直列に実施した場合の時間 # すなわち最悪値 horizon = 0 for i in all_jobs: horizon += sum(processing_times[i])
変数と制約条件の定義
ジョブショップ問題を、制約つき離散値組み合わせ最適化問題に変換するための作業をします
# Jobをタスクの定義にばらします all_tasks = {} for i in all_jobs: for j in range(0, len(machines[i])): all_tasks[(i, j)] = solver.FixedDurationIntervalVar(0, horizon, processing_times[i][j], False, 'Job_%i_%i' % (i, j)) # 変数の設定と離接制約(機械は1度に1つのタスク) all_sequences = [] all_machines_jobs = [] for i in all_machines: machines_jobs = [] for j in all_jobs: for k in range(0, len(machines[j])): if machines[j][k] == i: machines_jobs.append(all_tasks[(j, k)]) disj = solver.DisjunctiveConstraint(machines_jobs, 'machine %i' % i) all_sequences.append(disj.SequenceVar()) solver.Add(disj)
それぞれの加工順序を守るための順序制約を追加
# 作業順序の順序制約 for i in all_jobs: for j in range(0, len(machines[i]) - 1): solver.Add(all_tasks[(i, j + 1)].StartsAfterEnd(all_tasks[(i, j)]))
最適化の目的関数設定
全行程の一番最後のタスクが終わる時間を評価関数の目的変数に。(全てのタスクの完了時刻の最大値 = 全行程の完了時刻)
# 目的変数定義 obj_var = solver.Max([all_tasks[(i, len(machines[i])-1)].EndExpr() for i in all_jobs])
ソルバーのオプション設定
全行程の完了時刻が最小値になるものを探します
# 探索的に目的変数が小さいものになるように… sequence_phase = solver.Phase([all_sequences[i] for i in all_machines], solver.SEQUENCE_DEFAULT) vars_phase = solver.Phase([obj_var], solver.CHOOSE_FIRST_UNBOUND, solver.ASSIGN_MIN_VALUE) main_phase = solver.Compose([sequence_phase, vars_phase])
ソルバー
最適化実行
# 解を保存する変数 collector = solver.LastSolutionCollector() collector.Add(all_sequences) collector.AddObjective(obj_var) for i in all_machines: sequence = all_sequences[i]; sequence_count = sequence.Size(); for j in range(0, sequence_count): t = sequence.Interval(j) collector.Add(t.StartExpr().Var()) collector.Add(t.EndExpr().Var()) solver.Solve(main_phase, [objective_monitor, collector])
結果表示
ぐちゃぐちゃしてきてしまったので、最後のソースコード全てのところで…。
ソース
from __future__ import print_function # OR toolsのインポート from ortools.constraint_solver import pywrapcp # 変数操作を楽にするためNumpy import numpy as np # 視覚化用 import matplotlib.pyplot as plt def main(): # ソルバー定義 solver = pywrapcp.Solver('jobshop') # 機械の種類 # 0. 旋盤 # 1. ボール盤 # 2. フライス盤 machines_count = 3 all_machines = range(0, machines_count) # Jobの定義 # [m, t] : m 機械のID, t 加工時間(分) jobs = np.array([ [[0,3],[1,2],[2,3]], # Job 0 / 部品0 [[0,2],[2,3],[1,4]], # Job 1 / 部品1 [[1,3],[0,2],[2,1]] # Job 2 / 部品2 ]) jobs_count = 3 all_jobs = range(0, jobs_count) machines = jobs[:,:,0].tolist() processing_times = jobs[:,:,1].tolist() # 全てのJobを直列に実施した場合の時間 # すなわち最悪値 horizon = 0 for i in all_jobs: horizon += sum(processing_times[i]) # Jobをタスクの定義にばらします all_tasks = {} for i in all_jobs: for j in range(0, len(machines[i])): all_tasks[(i, j)] = solver.FixedDurationIntervalVar(0, horizon, processing_times[i][j], False, 'Job_%i_%i' % (i, j)) # 変数の設定と離接制約(機械は1度に1つのタスク) all_sequences = [] all_machines_jobs = [] for i in all_machines: machines_jobs = [] for j in all_jobs: for k in range(0, len(machines[j])): if machines[j][k] == i: machines_jobs.append(all_tasks[(j, k)]) disj = solver.DisjunctiveConstraint(machines_jobs, 'machine %i' % i) all_sequences.append(disj.SequenceVar()) solver.Add(disj) # 作業順序の順序制約 for i in all_jobs: for j in range(0, len(machines[i]) - 1): solver.Add(all_tasks[(i, j + 1)].StartsAfterEnd(all_tasks[(i, j)])) # 目的変数定義 obj_var = solver.Max([all_tasks[(i, len(machines[i])-1)].EndExpr() for i in all_jobs]) objective_monitor = solver.Minimize(obj_var, 1) # 探索的に目的変数が小さいものになるように… sequence_phase = solver.Phase([all_sequences[i] for i in all_machines], solver.SEQUENCE_DEFAULT) vars_phase = solver.Phase([obj_var], solver.CHOOSE_FIRST_UNBOUND, solver.ASSIGN_MIN_VALUE) main_phase = solver.Compose([sequence_phase, vars_phase]) # 解を保存する変数 collector = solver.LastSolutionCollector() collector.Add(all_sequences) collector.AddObjective(obj_var) for i in all_machines: sequence = all_sequences[i]; sequence_count = sequence.Size(); for j in range(0, sequence_count): t = sequence.Interval(j) collector.Add(t.StartExpr().Var()) collector.Add(t.EndExpr().Var()) disp_col_width = 10 # 最適化計算 if solver.Solve(main_phase, [objective_monitor, collector]): print("\nOptimal Schedule Length:", collector.ObjectiveValue(0), "\n") sol_line = "" sol_line_tasks = "" print("Optimal Schedule", "\n") Result = np.zeros((3,3,3)) for i in all_machines: seq = all_sequences[i] sol_line += "Machine " + str(i) + ": " sol_line_tasks += "Machine " + str(i) + ": " sequence = collector.ForwardSequence(0, seq) seq_size = len(sequence) for j in range(0, seq_size): t = seq.Interval(sequence[j]); sol_line_tasks += t.Name() + " " * (disp_col_width - len(t.Name())) for j in range(0, seq_size): t = seq.Interval(sequence[j]); sol_tmp = "[" + str(collector.Value(0, t.StartExpr().Var())) + "," sol_tmp += str(collector.Value(0, t.EndExpr().Var())) + "] " sol_line += sol_tmp + " " * (disp_col_width - len(sol_tmp)) sol_line += "\n" sol_line_tasks += "\n" # 結果 for j in range(0, seq_size): t = seq.Interval(sequence[j]); r, c = int(t.Name()[4]),int(t.Name()[6]) Result[r,c,0]=i Result[r,c,1]=collector.Value(0, t.StartExpr().Var()) Result[r,c,2]=collector.Value(0, t.EndExpr().Var()) print(sol_line_tasks) print("Time Intervals for Tasks\n") print(sol_line) print(Result) # 結果可視化 lglabels = ["Machine %i"%i for i in range(0,machines_count)] ylabels = ["Job %i"%i for i in range(0,jobs_count)] colormap = {0:'#27d67e', 1:'#2780d6', 2:'#d62728'} color_mapper = lambda cps:[colormap[c] for c in cps] temp = np.array(processing_times).cumsum(axis=1) machines = np.array(machines) fig, ax = plt.subplots() for i in range(0,Result.shape[1]): ax.barh(all_jobs, Result[:,i,2]-Result[:,i,1], left=Result[:,i,1], color=color_mapper(Result[:,i,0]),align='center') ax.set_yticks(all_jobs) ax.set_yticklabels(ylabels) ax.invert_yaxis() ax.set_xlabel('Time') ax.set_title('Optimized Process') plt.legend(lglabels) fig, ax = plt.subplots() for i in range(0,temp.shape[1]): left_ind = temp[:,i-1] if i>0 else [0,0,0] ax.barh(all_jobs, temp[:,i], left=left_ind, color=color_mapper(machines[:,i]),align='center') ax.set_yticks(all_jobs) ax.set_yticklabels(ylabels) ax.invert_yaxis() ax.set_xlabel('Time') ax.set_title('Job Definition') plt.legend(lglabels) plt.show() if __name__ == '__main__': main()
Pandas Groupbyの練習 (チートシート?)
Pandasを使っているとGroupbyな処理をしたくなることが増えてきます。ドキュメントを読んだりしながらよく使ったりする機能の骨格をまとめました。手っ取り早く勉強するなら、本が簡単そうです。
Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
- 作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
- 出版社/メーカー: オライリージャパン
- 発売日: 2013/12/26
- メディア: 大型本
- この商品を含むブログ (19件) を見る
では、さっそくスタート。
pandasライブラリなどなどの読み込み
普段使っている感じのものをコピペ。matplotlibはinlineの方が好き。
%matplotlib inline # %matplotlib notebook import sqlite3 import pandas as pd from matplotlib import pylab as plt import numpy as np from datetime import datetime import dateutil
csvからpandasへのデータの読み込み
Climate Data Onlineから引っ張ってきた2016年1月の東京と京都に関する気象データをちょっと加工したものを使用。アメリカが出典なだけあって、単位がアメリカ風味。
df = pd.read_csv('./201601.csv',encoding="SHIFT-JIS") df.head()
Station | YEARMODA | TEMP | MAXT | MINT | WDSP | MXSPD | |
---|---|---|---|---|---|---|---|
0 | Tokyo | 20160101 | 45.1 | 53.8 | 37.2* | 4.9 | 9.9 |
1 | Tokyo | 20160102 | 45.3 | 56.1 | 36.7 | 2.5 | 6.0 |
2 | Tokyo | 20160103 | 49.6 | 61.2 | 37.4 | 3.0 | 5.1 |
3 | Tokyo | 20160104 | 49.9 | 59.4 | 39.9 | 3.4 | 6.0 |
4 | Tokyo | 20160105 | 51.5 | 59.7 | 41.2 | 3.3 | 5.1 |
- Station : 観測地点。今回は東京か京都。
- YEARMODA : YYYYMMDDフォーマットの年月日データ
- TEMP : 1日の平均気温。単位 華氏[℉]
- MAXT : 1日の最高気温。単位 華氏[℉] *印付きには事情があるのだが、今はゴミ。
- MINT : 1日の最低気温。単位 華氏[℉] *印付きは事情があるのだが、今はゴミ。
- WDSP : 1日の平均風速。単位 1/10ノット
- MXSP : 1日の最大風速。単位 1/10ノット
df.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 60 entries, 0 to 59
Data columns (total 7 columns):
Station 60 non-null object
YEARMODA 60 non-null int64
TEMP 60 non-null float64
MAXT 60 non-null object
MINT 60 non-null object
WDSP 60 non-null float64
MXSPD 60 non-null float64
dtypes: float64(3), int64(1), object(3)
memory usage: 3.4+ KB
pandasでのdatetime型への変換; 文字列→浮動小数点変換
ちょっとした前処理をします。MAXT, MINTがたまに入っている*のせいでobject扱いに。floatに直そう。ついでにYEARMODAも変換。
# YEARMODAをdatetime型にパース df['YEARMODA'] = df['YEARMODA'].apply(lambda x: dateutil.parser.parse(str(x))) # *を除去してfloat型に。 df['MAXT'] = df['MAXT'].apply(lambda x: float(str.replace(x, '*', ''))) df['MINT'] = df['MINT'].apply(lambda x: float(str.replace(x, '*', ''))) df.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 60 entries, 0 to 59
Data columns (total 7 columns):
Station 60 non-null object
YEARMODA 60 non-null datetime64[ns]
TEMP 60 non-null float64
MAXT 60 non-null float64
MINT 60 non-null float64
WDSP 60 non-null float64
MXSPD 60 non-null float64
dtypes: datetime64[ns](1), float64(5), object(1)
memory usage: 3.4+ KB
df.head()
Station | YEARMODA | TEMP | MAXT | MINT | WDSP | MXSPD | |
---|---|---|---|---|---|---|---|
0 | Tokyo | 2016-01-01 | 45.1 | 53.8 | 37.2 | 4.9 | 9.9 |
1 | Tokyo | 2016-01-02 | 45.3 | 56.1 | 36.7 | 2.5 | 6.0 |
2 | Tokyo | 2016-01-03 | 49.6 | 61.2 | 37.4 | 3.0 | 5.1 |
3 | Tokyo | 2016-01-04 | 49.9 | 59.4 | 39.9 | 3.4 | 6.0 |
4 | Tokyo | 2016-01-05 | 51.5 | 59.7 | 41.2 | 3.3 | 5.1 |
Pandasの行方向計算
華氏から摂氏への変換
[C] = ([F] − 32) * 5/ 9で変換
df["TEMP"] = (df["TEMP"]-32)*5/9 df["MAXT"] = (df["MAXT"]-32)*5/9 df["MINT"] = (df["MINT"]-32)*5/9
風速:0.1ノットからm/sへの変換
0.1 ノット = 0.0514444[m/s]だそうです
df["WDSP"] *= 0.051444 df["MXSPD"] *= 0.051444
df.head()
Station | YEARMODA | TEMP | MAXT | MINT | WDSP | MXSPD | |
---|---|---|---|---|---|---|---|
0 | Tokyo | 2016-01-01 | 7.277778 | 12.111111 | 2.888889 | 0.252076 | 0.509296 |
1 | Tokyo | 2016-01-02 | 7.388889 | 13.388889 | 2.611111 | 0.128610 | 0.308664 |
2 | Tokyo | 2016-01-03 | 9.777778 | 16.222222 | 3.000000 | 0.154332 | 0.262364 |
3 | Tokyo | 2016-01-04 | 9.944444 | 15.222222 | 4.388889 | 0.174910 | 0.308664 |
4 | Tokyo | 2016-01-05 | 10.833333 | 15.388889 | 5.111111 | 0.169765 | 0.262364 |
groupbyとグループ毎での統計量
観測地点ごとにグルーピングして統計情報を出力をする練習。東京、京都それぞれの2016年1月の平均気温の平均、2016年1月の最高気温・最低気温を出してみましょう。
個別に出力する
統計量を個別に出力する場合です。あとでまとめて出す方法やります。
平均気温の平均
df.groupby("Station")["TEMP"].mean()
Station
Kyoto 5.881226
Tokyo 6.035842
Name: TEMP, dtype: float64
最高気温
df.groupby("Station")["MAXT"].max()
Station
Kyoto 16.000000
Tokyo 16.222222
Name: MAXT, dtype: float64
最低気温
df.groupby("Station")["MINT"].min()
Station
Kyoto -4.111111
Tokyo -2.611111
Name: MINT, dtype: float64
参考:Pandas組み込みのGroup用統計関数
上で使った基本的なもの以外にも、Essential Basic Functionality — pandas 0.19.2 documentationにある通り、いろいろな関数が組み込まれています。抜粋して翻訳。
関数 | 機能 |
---|---|
count | Non-Nullなデータの数 |
sum | 合計値 |
mean | 平均値 |
mad | 平均絶対偏差 |
median | メジアン(中央値)を出力 |
min | 最小値 |
max | 最大値 |
mode | モード(最頻値) |
abs | 絶対値 |
prod | 総乗 |
std | 標準偏差. 自由度を引数にとれるがデフォルトは1. |
var | 不偏分散 |
sem | 平均値の標準誤差(cf : 推計統計学) |
skew | 歪度 |
kurt | 尖度 |
quantile | 分位数またはクォンタイル。デフォルトでは中央値を出力する。 .quantile(.25)で25パーセンタイルを出力。 [参考]pandas.quantile — pandas |
cumsum | 累積和を出力(デフォルトで行方向) |
cumprod | 累積乗を出力 (デフォルトで行方向) |
cummax | 累積最大値。axis指定が必要。 |
cummin | 累積最小値。累積最大値と似てる。 |
まとめて出力する
個別に統計量を出力させてもよいのですが、まとめてバッっと出したい場合は。.agg関数を使うと素敵になる。
aggregations = { 'TEMP': { '1月の平均気温の平均': 'mean' }, 'MAXT': { '1月の最高気温' : 'max' }, 'MINT': { '1月の最低気温' : 'min' } } df.groupby("Station").agg(aggregations)
TEMP | MAXT | MINT | |
---|---|---|---|
1月の平均気温の平均 | 1月の最高気温 | 1月の最低気温 | |
Station | |||
Kyoto | 5.881226 | 16.000000 | -4.111111 |
Tokyo | 6.035842 | 16.222222 | -2.611111 |
Groupbyとmatplotlib
Pandasでグループ毎に分けてプロットする方法。連続してプロットすると重ね描きになるのはmatplotlibと同様。
cmap = plt.get_cmap('rainbow') fig, ax = plt.subplots() #df['YEARMODA'] = df.YEARMODA.astype(np.int64) for i, (key, group) in enumerate(df.groupby(['Station']), start=1): ax = group.plot(ax=ax, x='YEARMODA', y='TEMP',color=cmap(i / 4.0),label = key+"_avg") for i, (key, group) in enumerate(df.groupby(['Station']), start=1): ax = group.plot(ax=ax, x='YEARMODA', y='MAXT',color=cmap(2*(i+1) / 4.0),label = key+"_max") plt.title(u"2016/01 Temp") fig.patch.set_facecolor('white') plt.show()
Groupレベルの処理をDataFrameに戻す
Group別の処理結果を使って、元のデータを処理したい時のお話。
DataFrameにGroup別最大値を追加
DataFrameに、観測地別2016年1月最高気温をMAXT_Monthという列名で追加
df['MAXT_Month'] = df.loc[df.groupby('Station')['MAXT'].transform('idxmax'), 'MAXT'].values ## 別の方法 # df = df.join(df.groupby('Station')['MAXT'].max(), on='Station', rsuffix='_Month') df.head()
Station | YEARMODA | TEMP | MAXT | MINT | WDSP | MXSPD | MAXT_Month | |
---|---|---|---|---|---|---|---|---|
0 | Tokyo | 2016-01-01 | 7.277778 | 12.111111 | 2.888889 | 0.252076 | 0.509296 | 16.222222 |
1 | Tokyo | 2016-01-02 | 7.388889 | 13.388889 | 2.611111 | 0.128610 | 0.308664 | 16.222222 |
2 | Tokyo | 2016-01-03 | 9.777778 | 16.222222 | 3.000000 | 0.154332 | 0.262364 | 16.222222 |
3 | Tokyo | 2016-01-04 | 9.944444 | 15.222222 | 4.388889 | 0.174910 | 0.308664 | 16.222222 |
4 | Tokyo | 2016-01-05 | 10.833333 | 15.388889 | 5.111111 | 0.169765 | 0.262364 | 16.222222 |
Group別に他の列の条件に応じて代入値を決める
例えば、観測地別2016年1月 最高気温日の最大風速をMXSPD_of_DAY_MAXTという列名で追加
df['MXSPD_of_DAY_MAXT'] = df.loc[df.groupby('Station')['MAXT'].transform('idxmax'), 'MXSPD'].values df.head()
Station | YEARMODA | TEMP | MAXT | MINT | WDSP | MXSPD | MAXT_Month | MXSPD_of_DAY_MAXT | |
---|---|---|---|---|---|---|---|---|---|
0 | Tokyo | 2016-01-01 | 7.277778 | 12.111111 | 2.888889 | 0.252076 | 0.509296 | 16.222222 | 0.262364 |
1 | Tokyo | 2016-01-02 | 7.388889 | 13.388889 | 2.611111 | 0.128610 | 0.308664 | 16.222222 | 0.262364 |
2 | Tokyo | 2016-01-03 | 9.777778 | 16.222222 | 3.000000 | 0.154332 | 0.262364 | 16.222222 | 0.262364 |
3 | Tokyo | 2016-01-04 | 9.944444 | 15.222222 | 4.388889 | 0.174910 | 0.308664 | 16.222222 | 0.262364 |
4 | Tokyo | 2016-01-05 | 10.833333 | 15.388889 | 5.111111 | 0.169765 | 0.262364 | 16.222222 | 0.262364 |