Pythonを用いた最適化

Pythonを用いた最適化
In [1]:
from IPython.display import display,Math
from pulp import *

Pythonを用いた最適化
小林和博

  • モデラー PuLP, PICOSとOpenOpt
  • PuLPによる線形最適化
  • PICOSによる二次錐最適化
  • OpenOptによる非線形最適化

1 モデラーPuLPとPICOS

  • PuLPは,線形最適化のモデラ,整数変数も可
  • PICOSは,錐最適化のモデラ,整数変数も可
    • 錐最適化は,線形最適化の拡張
    • PuLPでできることは,手間を厭わなければ,PICOSでもできる
  • 線形最適化は,PuLPの方が楽に書ける.多くのORモデルを簡単に書ける.
  • PICOSは行列をそのまま扱える.行列の演算も可.制御など工学向き.

1. モデラーOpenOpt

  • OpenOptは,線形最適化,非線形最適化のモデラ
  • 錐最適化も,整数変数を可

1-1 線形最適化

  • 目的関数と制約式がすべて線形な最適化
  • 大規模な入力データでも,たいてい高速に解ける.充分に色々なORモデルを表現できる.
  • 整数変数を含んでいても,注意してモデル化すれば,かなり大きなデータでも解ける.

線形最適化の例

最小化 $2x + 3 y$

制約条件 $x + y \leq 2,$

$0 \leq x \leq 3,$

$0 \leq y \leq 1 $

PuLPで書くと,

In [2]:
prob=pulp.LpProblem("myLP",pulp.LpMinimize)
x=pulp.LpVariable("x",0,3)
y=pulp.LpVariable("y",0,1)
prob+=2*x+3*y
prob+=x+y<=2 

線形最適化

一般形

最小化 $\displaystyle \sum_{j=1}^n c_j x_j$

制約条件 $\displaystyle \sum_{j=1}^n a_{ij}x_j = b_i \ (i=1,2,\ldots,m)$

$x_j \geq 0$

1-2. 非線形最適化

  • 目的関数または制約式に非線形な式が含まれている.
  • 局所最適解を求める方法が発達している.
  • 効率のよい計算には,微分を与えるなどの処理が要ることがある.

1-3. 錐最適化

  • 変数が「錐」に入っているという条件(非線形)
  • 目的関数と(錐意外の制約式は)すべて線形
  • 内点法など,効率のよい計算手法が開発されている(ソフトウェアも利用可)

非負象限と錐

Out[3]:

2 PuLPによる線形最適化

2-1 集合分割問題

  • 集合を,いくつかの部分集合に分割する問題
  • 各部分集合には,利益(コスト)がついている
  • 利益を最大(コストを最小)にするような分け方を求める
Out[4]:

例 : 結婚式の席決め

Out[5]:

線形最適化としての定式化

Out[6]:

線形最適化としての定式化

Out[7]:

2-2 変数の定義

まず,ゲストのすべての分け方(部分集合)を求める.

In [8]:
max_tables = 5
max_table_size = 4
guests=['A','B','C','D']
possible_tables=[tuple(c) for c in pulp.allcombinations(guests,4)]
print possible_tables
[('A',), ('B',), ('C',), ('D',), ('A', 'B'), ('A', 'C'), ('A', 'D'), ('B', 'C'), ('B', 'D'), ('C', 'D'), ('A', 'B', 'C'), ('A', 'B', 'D'), ('A', 'C', 'D'), ('B', 'C', 'D'), ('A', 'B', 'C', 'D')]

possible_tablesの各要素に対して0-1変数を定義する.(先ほどのx=pulp.LpVariable("x",0,1)のかわりに)

In [9]:
x=pulp.LpVariable.dicts('table',possible_tables,0,1,cat=pulp.LpInteger) 
print x
print x[('B','C')]
{('B', 'C'): table_('B',_'C'), ('A',): table_('A',), ('C', 'D'): table_('C',_'D'), ('C',): table_('C',), ('A', 'B', 'D'): table_('A',_'B',_'D'), ('B',): table_('B',), ('D',): table_('D',), ('B', 'C', 'D'): table_('B',_'C',_'D'), ('A', 'D'): table_('A',_'D'), ('A', 'C', 'D'): table_('A',_'C',_'D'), ('A', 'B', 'C'): table_('A',_'B',_'C'), ('A', 'B'): table_('A',_'B'), ('A', 'B', 'C', 'D'): table_('A',_'B',_'C',_'D'), ('A', 'C'): table_('A',_'C'), ('B', 'D'): table_('B',_'D')}
table_('B',_'C')

In [10]:
print x[('A',)]
table_('A',)

2-3 制約式と目的関数の定義

(その前に問題の定義)

In [11]:
seating_model=pulp.LpProblem("Seating",pulp.LpMinimize) 

各席割の"幸福度"の定義

In [12]:
def happiness(table):
    return abs(ord(table[0]) - ord(table[-1]))

目的関数

In [13]:
seating_model += sum([happiness(table)*x[table] for table in possible_tables])

制約式

テーブル数の上限

In [14]:
seating_model += sum([x[table] for table in possible_tables]) <= max_tables

全てのゲストはどこかのテーブルに席がある

In [15]:
for guest in guests:
    seating_model += sum([x[table] for table in possible_tables  if guest in table]) == 1, "Must_seat_%s"%guest
In [16]:
seating_model.solve()
Out[16]:
1

3 PICOSによる二次錐最適化

PICOS: A Python Interface for Conic Optimization Solvers

- http://picos.zib.de/
- 錐最適化のためのインターフェイス

3-1 二次錐最適化

  • 錐最適化のひとつ
  • 線形最適化の拡張
  • 内点法で効率的に解ける
  • 二次関数が表現できる

線形最適化と二次錐最適化

線形最適化

最小化 $\displaystyle \sum_{j=1}^n c_j x_j $

制約条件 $ \displaystyle \sum_{j=1}^n a_{ij} x_j = b_i \ (i=1,2,\ldots,m), (x_1,x_2, \ldots,x_n) \geq (0,0, \ldots,0) $

をベクトルを使って

最小化 $ \mathbf{c}^{\top} \mathbf{x}$

制約条件 $ \mathbf{a}_i^{\top} \mathbf{x} = b_i \ (i=1,2,\ldots,m), \mathbf{x} \geq \mathbf{0} $

と書くことにすると,

線形最適化と二次錐最適化

線形最適化

最小化 $ \mathbf{c}^{\top} \mathbf{x}$

制約条件 $ \mathbf{a}_i^{\top} \mathbf{x} = b_i \ (i=1,2,\ldots,m), \mathbf{x} \geq \mathbf{0} $

二次錐最適化

最小化 $ \mathbf{c}^{\top} \mathbf{x}$

制約条件 $ \mathbf{a}_i^{\top} \mathbf{x} = b_i \ (i=1,2,\dots,m), \mathbf{x} \succeq_{\text{Q}} \mathbf{0}$

ただし,$\mathbf{x} \succeq_{\text{Q}} \mathbf{0} \Leftrightarrow \displaystyle \sqrt{\sum_{i=2}^n x_i^2}, x_i \geq 0 $

3-2 Weber問題

  • 久保他「あたらしい数理最適化 -Python言語とGurobiで解く」より引用

砂漠に7組の家がある.表10.1に各家の位置と人数が書かれている.この地域には井戸がなく,水は何キロも離れたところへ,毎日取りにいかねばならないのが,重労働であった.そこでみんなで出資して新しく井戸を掘ることになった.なるべく水平に井戸の位置を決めたいので,「各家が水を運ぶ距離×各家の水の消費量」の総和が最小になる場所に井戸を掘ることにした.どこを掘っても水が出るものとしたとき,どのようにして掘る場所を決めれば良いだろうか.ただし,各家の水の消費量は人数に比例するものとする.

3-2 Weber問題

x座標 y座標 人数
1 24 54 2
2 60 63 1
3 1 84 2
4 23 100 3
5 84 48 4
6 15 64 5
7 52 74 4
Out[17]:

3-2 Weber問題

家の集合をHとする.各家の位置を$(x_i,y_i) \ (i \in H) $とする.井戸の位置を$(X,Y)$とすれば,家iから井戸までの距離は,$\sqrt{(x_i-X)^2+(y_i-Y)^2}$である。 

である.家iが1日に必要とする水の量を$w_i$としたとき,$\displaystyle \sum_{i \in H} w_i \sqrt{(x_i-X)^2+(y_i-Y)^2} $を最小にする$(x_i,y_i)$を求めよ.

3-2 Weber問題

二次錐最適化問題として書くと,

最小化 $\displaystyle \sum_{i \in H} w_i z_i $

制約条件 $\sqrt{(x_i-X)^2+(y_i-Y)^2} \leq z_i \ (\forall i \in H)$

これは,目的関数は線形で,制約条件は

ベクトル$(z_i,x_i-X,y_i-Y)$が二次錐に入っている,という条件なので,二次錐最適化問題である.

3-3 問題とデータ(PICOS)

(その前に問題の定義)

In [18]:
import picos as pic
socp=pic.Problem()

(それからデータの定義)

In [19]:
  H = [0,1,2,3,4,5,6]
  xy=[[24,54],[60,63],[1,84],[12,100],[84,48],[15,64],[52,74]]
  w=[2,1,2,3,4,5,4] 

3-4 変数の定義(PICOS)

最小化 $\displaystyle \sum_{i \in H} w_i z_i$

制約条件 $\sqrt{(x_i-X)^2+(y_i-Y)^2} \leq z_i \ (i \in H)$

2次元のベクトルとして変数を定義する(XY). 各家に変数を定義する(z).

In [20]:
XY=socp.add_variable('XY',2)
z=[socp.add_variable('z['+str(i)+']',1) for i in H] 
In [21]:
print XY
# variable XY:(2 x 1),continuous #

In [22]:
print z
[# variable z[0]:(1 x 1),continuous #, # variable z[1]:(1 x 1),continuous #, # variable z[2]:(1 x 1),continuous #, # variable z[3]:(1 x 1),continuous #, # variable z[4]:(1 x 1),continuous #, # variable z[5]:(1 x 1),continuous #, # variable z[6]:(1 x 1),continuous #]

3-5 制約式と目的関数の定義

目的関数

In [23]:
objective = sum(w[i] * z[i] for i in H)
socp.set_objective('min', objective)

最小化$\displaystyle \sum_{i \in H} w_i z_i$

制約条件 $\sqrt{(x_i-X)^2+(y_i-Y)^2} \leq z_i \ (\forall i \in H)$

3-5 制約式と目的関数の定義

制約式

In [26]:
socp.add_list_of_constraints([abs(xy[i]-XY)<z[i] for i in H])

3-6 二次錐最適化

最小化$\displaystyle \sum_{i \in H} w_i z_i$

制約条件 $\sqrt{(x_i-X)^2+(y_i-Y)^2} \leq z_i \ (\forall i \in H)$

In [27]:
socp=pic.Problem()
XY=socp.add_variable('XY',2)
z=[socp.add_variable('z['+str(i)+']',1) for i in H]
objective=sum(w[i]*z[i] for i in H)
socp.set_objective('min',objective)
socp.add_list_of_constraints([abs(xy[i]-XY) < z[i] for i in H]) 
socp.solve(solver='cvxopt') 
--------------------------
  cvxopt CONELP solver
--------------------------
     pcost       dcost       gap    pres   dres   k/t
 0: -1.0658e-14 -4.0079e-14  1e+03  6e-01  2e-16  1e+00
 1:  5.0224e+02  5.0811e+02  2e+02  1e-01  9e-16  6e+00
 2:  6.2367e+02  6.2433e+02  2e+01  1e-02  3e-16  7e-01
 3:  6.3571e+02  6.3574e+02  9e-01  5e-04  5e-16  3e-02
 4:  6.3628e+02  6.3628e+02  3e-02  2e-05  5e-15  1e-03
 5:  6.3629e+02  6.3629e+02  3e-03  2e-06  3e-14  9e-05
 6:  6.3630e+02  6.3630e+02  4e-04  2e-07  3e-14  2e-05
 7:  6.3630e+02  6.3630e+02  1e-04  6e-08  4e-13  5e-06
 8:  6.3630e+02  6.3630e+02  4e-05  3e-08  2e-13  3e-06
 9:  6.3630e+02  6.3630e+02  3e-06  2e-09  7e-13  3e-07
Optimal solution found.
cvxopt status: optimal

Out[27]:
{'cvxopt_sol': {'dual infeasibility': 7.194144966587666e-13,
  'dual objective': 636.2960965135727,
  'dual slack': 2.992548875369039e-09,
  'gap': 2.808734890897799e-06,
  'iterations': 9,
  'primal infeasibility': 1.724097848441206e-09,
  'primal objective': 636.2960962587284,
  'primal slack': 3.3002926613789896e-08,
  'relative gap': 4.4141947534922315e-09,
  'residual as dual infeasibility certificate': None,
  'residual as primal infeasibility certificate': None,
  's': <21x1 matrix, tc='d'>,
  'status': 'optimal',
  'x': <9x1 matrix, tc='d'>,
  'y': <0x1 matrix, tc='d'>,
  'z': <21x1 matrix, tc='d'>},
 'obj': 636.2960963861506,
 'status': 'optimal',
 'time': 0.10760807991027832}

4 OpenOptによる非線形最適化

4-1 制約付きWeber問題

最小化$\displaystyle \sum_{i \in H} w_i \sqrt{(x_i-X)^2+(y_i-Y)^2}$

制約条件 $(50-X)^2+(50-Y)^2 \geq 40^2$

4-2 問題とデータ

(その前に問題の定義)

In [28]:
from FuncDesigner import *
from openopt import NLP
import numpy as np 

(それからデータの定義)

In [29]:
  H = [0,1,2,3,4,5,6]
  x=np.array([24,60,1,12,84,15,52])
  y=np.array([54,63,84,100,48,64,74])
  w=np.array([2,1,2,3,4,5,4]) 

4-3 変数の定義(OpenOpt)

最小化$\displaystyle \sum_{i \in H} w_i \sqrt{(x_i-X)^2+(y_i-Y)^2}$

制約条件 $(50-X)^2+(50-Y)^2 \geq 40^2$

In [30]:
X,Y=oovars('X','Y') 

4-4 制約式と目的関数の定義

In [31]:
objective=sum(w[i]*sqrt( (x[i]-X)**2+(y[i]-Y)**2) for i in H)

constraints=[(50-X)**2+(50-Y)**2>=40**2] 

4-5 非線形最適化

In [32]:
from FuncDesigner import *
from openopt import NLP
import numpy as np
x=np.array([24,60,1,12,84,15,52])
y=np.array([54,63,84,100,48,64,74])
w=np.array([2,1,2,3,4,5,4])
X,Y=oovars('X','Y')
objective=sum(w[i]*sqrt((x[i]-X)**2+(y[i]-Y)**2) for i in H)
constraints=[(50-X)**2+(50-Y)**2>=40**2]
startPoint = {X:0, Y:0}
p=NLP(objective,startPoint,constraints=constraints)
p.minimize('ralg')

------------------------- OpenOpt 0.5404 -------------------------
solver: ralg   problem: unnamed    type: NLP   goal: minimum
 iter   objFunVal   log10(maxResidual)   
    0  1.753e+03            -100.00 
   10  6.774e+02            -100.00 
   20  6.764e+02            -100.00 
   30  6.764e+02            -100.00 
   40  6.764e+02            -100.00 
   44  6.764e+02            -100.00 
istop: 3 (|| X[k] - X[k-1] || < xtol)
Solver:   Time Elapsed = 0.43 	CPU Time Elapsed = 0.43
objFunValue: 676.40668 (feasible, MaxResidual = 0)

Out[32]:
<runProbSolver.OpenOptResult instance at 0x108f421b8>

まとめ

  • Pythonでは,数理最適化問題を簡単に記述し,解くことができます.
  • モデラー PuLP,PICOSとOpenOpt
  • PuLPによる線形最適化
  • PICOSによる二次錐最適化
  • OpenOptによる非線形最適化