pulpとamplify の相互変換

*amplify は AWSとは関係ない方です.

TLDR

本稿ではpythonから使える最適化モデラーであるpulpとamplify間でのモデリングの相互変換コードを紹介します.

1. pulpで定義した問題をamplifyで実行できるように変換し結果をpulpの変数に格納する
2. amplifyで定義した問題をpulpで解いて解を取得する

pulp

coin-orが提供するpython用の線形計画ソルバー

amplify

株式会社フィックスターズが提供する量子アニーリング・イジングマシン
https://amplify.fixstars.com/ja/
イジング形式, QUBO, また整数線形計画のモデリングが可能です.

pulp → amplify

pulpで定義した問題をamplifyで実行できるように変換し結果をpulpの変数に格納する

# -------------------------------------------
#       pulp による問題作成
# -------------------------------------------
import pulp

x = pulp.LpVariable('x', cat='Binary')  # バイナリ変数の作成
y = pulp.LpVariable('y', cat='Binary')  # バイナリ変数の作成

prob = pulp.LpProblem()
prob += -x + 4*y + 100  # 目的関数の設定

print('-'*10, 'pulp', '-'*10)
print(prob)
print('-'*30)


# -------------------------------------------
#       amplify用に変換
# -------------------------------------------
from amplify import BinaryPoly, gen_symbols, Solver, decode_solution

vals =  prob.variables()  # 変数リスト
index2name = dict()  # 結果取り出し用
name2var = {var.name: var for var in vals} # 結果取り出し用

# イジング変数を変数数分定義
s = gen_symbols(BinaryPoly, len(vals))

# 目的関数 を定義
f = prob.objective.constant
for i in range(len(vals)):
    var_data = prob.objective.to_dict()[i]
    index2name[i] = var_data['name']
    f += var_data['value'] * s[i]

print('-'*10, 'amplify', '-'*10)
print(f"f = {f}")


# -------------------------------------------
#       求解( amplify )
# -------------------------------------------
from amplify.client import FixstarsClient
client = FixstarsClient()
# client.parameters.time_limit=1000
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"  # ローカル環境で使用する場合は、Amplify AEのアクセストークンを入力してください
solver = Solver(client)
result = solver.solve(f)  # 問題を入力してマシンを実行


# -------------------------------------------
#       結果の取り出し -> pulp変数に代入
# -------------------------------------------
sol = result.solutions[0]  # 解の取得
solution = decode_solution(s, sol.values)  #  変数配列sをsol.valuesでデコード
for i, value in enumerate(solution):
    var = name2var[ index2name[i] ]
    var.varValue = value # pulpの変数に結果を代入

print(f'result: {s} = {solution} (f = {sol.energy})')
print('objective value', pulp.value(prob.objective))
print('-'*30)

実行結果

---------- pulp ----------
NoName:
MINIMIZE
-1*x + 4*y + 100
VARIABLES
0 <= x <= 1 Integer
0 <= y <= 1 Integer

------------------------------
---------- amplify ----------
f = - q_0 + 4.000000 q_1 + 100.000000
result: [q_0, q_1] = [1, 0] (f = 99.0)
objective value 99
------------------------------

amplify → pulp

amplifyで定義した問題をpulpで解いて解を取得する

# -------------------------------------------
#       amplifyによる問題作成
# -------------------------------------------
from amplify import IsingPoly, BinaryPoly, gen_symbols

# イジング変数s_0, s_1を定義
s = gen_symbols(IsingPoly, 2)

# 目的関数 f = 1 - s_0 * s_1 を定義
f = 1 - s[0] * s[1]

print('-'*10, 'amplify', '-'*10)
print(f"f = {f}")


# -------------------------------------------
#       pulp に変換
# -------------------------------------------
import pulp

x = [ pulp.LpVariable(f'q{i}', cat='Binary')
      for i in range(len(s)) ]
y = dict()  #   y[i, j] = x[i] * x[j]

prob = pulp.LpProblem()
g = f.to_Binary()  # IsingからBinaryに変換

obj = g.constant()
for key, value in g.asdict().items():
    if len(key) == 1:
        obj += value * x[key[0]]
    else:
        i, j = key  # value * x_i * x_j 項を追加
        if (i, j) not in y:
            y[i, j] = pulp.LpVariable(f'y{i}_{j}', cat='Binary') # y_ij = x_i * x_j 変数を作成
        obj += value * y[i, j]
        prob += y[i, j] >= x[i] + x[j] - 1  # y_ij = x_i * x_j を成立させる制約
        prob += x[i] >= y[i, j]             # y_ij = x_i * x_j を成立させる制約
        prob += x[j] >= y[i, j]             # y_ij = x_i * x_j を成立させる制約

prob += obj

print('-'*10, 'pulp', '-'*10)
print(prob)


# -------------------------------------------
#       求解( pulp )
# -------------------------------------------
print('-'*10, 'pulp', '- solve -', '-'*10)
prob.solve() # 求解
print('objective value', pulp.value(prob.objective)) # 目的関数値の出力
print('-'*30)

実行結果

---------- amplify ----------
f = - s_0 s_1 + 1.000000
---------- pulp ----------
NoName:
MINIMIZE
2.0*q0 + 2.0*q1 + -4.0*y0_1 + 0.0
SUBJECT TO
_C1: - q0 - q1 + y0_1 >= -1

_C2: q0 - y0_1 >= 0

_C3: q1 - y0_1 >= 0

VARIABLES
0 <= q0 <= 1 Integer
0 <= q1 <= 1 Integer
0 <= y0_1 <= 1 Integer

---------- pulp - solve - ----------
Welcome to the CBC MILP Solver 
Version: 2.9.0 
Build Date: Feb 12 2015 

(中略)

Result - Optimal solution found

Objective value:                0.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.01

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.01

objective value 0.0
------------------------------