Volo di notte

お勉強の成果メモや日常のこと

夜間飛行へ移行中です

 

Google Optimization Toolsで組み合わせ最適化

巡回セールスマン問題やナップサック問題は、組み合わせ最適化と呼ばれる最適化問題の分野のひとつで扱われているとのこと。手法の全容は色々本が出ています。

今日から使える!組合せ最適化 離散問題ガイドブック (KS理工学専門書)

今日から使える!組合せ最適化 離散問題ガイドブック (KS理工学専門書)

実際の使い道として、生産工程の最適化などがあり、便利なソルバーとして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)について、最適工程割り付けを行います。

機械と部品を便宜上、数字で呼ぶと扱いやすいので、

  1. 旋盤
  2. ボール盤
  3. フライス盤

としましょうか。工程は下記の通り(言うまでもなく、それぞれの部品の加工順序は守らなくてはなりません)

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
    ]

この問題を解かせます

結果のイメージ f:id:Chachay:20170805170527p:plain

ソルバーの宣言と問題の設定

横着して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読み込み~整形まで

本記事は日経平均オプション理論価格等情報のPandas読み込み~整形まで - 夜間飛行へ更新とともに移転しました

Pandas Groupbyの練習 (チートシート?)

Pandasを使っているとGroupbyな処理をしたくなることが増えてきます。ドキュメントを読んだりしながらよく使ったりする機能の骨格をまとめました。手っ取り早く勉強するなら、本が簡単そうです。

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

では、さっそくスタート。

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
  1. Station : 観測地点。今回は東京か京都。
  2. YEARMODA : YYYYMMDDフォーマットの年月日データ
  3. TEMP : 1日の平均気温。単位 華氏[℉]
  4. MAXT : 1日の最高気温。単位 華氏[℉] *印付きには事情があるのだが、今はゴミ。
  5. MINT : 1日の最低気温。単位 華氏[℉] *印付きは事情があるのだが、今はゴミ。
  6. WDSP : 1日の平均風速。単位 1/10ノット
  7. 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()

f:id:Chachay:20170212151243p:plain

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