Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
0
Help us understand the problem. What are the problem?

posted at

updated at

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 による問題作成
# -------------------------------------------
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による問題作成
# -------------------------------------------
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 

command line - /usr/local/lib/python3.9/site-packages/pulp/apis/../solverdir/cbc/osx/64/cbc /var/folders/5z/fpt4ysfj2rv0bdwxpt5sfd880000gn/T/c01e3f64ef5745aa8e8754062688112e-pulp.mps branch printingOptions all solution /var/folders/5z/fpt4ysfj2rv0bdwxpt5sfd880000gn/T/c01e3f64ef5745aa8e8754062688112e-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 8 COLUMNS
At line 25 RHS
At line 29 BOUNDS
At line 33 ENDATA
Problem MODEL has 3 rows, 3 columns and 7 elements
Coin0008I MODEL read with 0 errors
Continuous objective value is 0 - 0.00 seconds
Cgl0004I processed model has 3 rows, 3 columns (3 integer (3 of which binary)) and 7 elements
Cutoff increment increased from 1e-05 to 1.9999
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 0
Cbc0038I Before mini branch and bound, 3 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.00 seconds)
Cbc0038I After 0.00 seconds - Feasibility pump exiting with objective of 0 - took 0.00 seconds
Cbc0012I Integer solution of 0 found by feasibility pump after 0 iterations and 0 nodes (0.00 seconds)
Cbc0001I Search completed - best objective 0, took 0 iterations and 0 nodes (0.00 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 0 to 0
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

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
------------------------------
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
0
Help us understand the problem. What are the problem?