LoginSignup
2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

@___monta___

イジングモデル（Quantum Annealer）で不等式制約をやっつける話

Last updated at Posted at 2023-03-04

最初にお願い事項。
コードに不備があったら教えてください。

最初に結論

  • 不等式を含む数理最適化問題をD-Waveなどイジングソルバで解く論文を実装した。
  • ほどほど解けているが、目的関数の値はあまりよくなさそう、ばらつきが大きい
  • 大規模問題はあまり解けない。現時点でのQA上での実用的な利用は難しそう

論文との出会い

ある日、Twitterにこんな情報が
＊キャルちゃんさんありがとうございます。

これはよいものでは？ ということで実装してみました。

問題設定

不等式を含む対象となる問題は、"set cover problem"です。
日本語で言うと集合被覆問題その中でも重み付き集合被覆問題ですね。

結果

イジングモデルでも不等式をやり過ごす方法がありますが、この方法比較的楽そうです。
QAを用いた場合、目的関数の値としてはあまりよくなさそうです。
今回試していないスラック変数とは比較していないので、現時点ではよい方なのかも。

手法 時間 制約違反 目的関数値 備考
従来手法 1秒以下 なし 10.0 最適解
本手法 5 ~ 20秒 なし 26 ~ 108

確認方法詳細 & コード

解法 1. CBC

普通にPyomo＋CBCで解きます。
*そもそも問題に解があるかさぼってチェックしてないんので、ダメな場合はSeedを変えましょう。

Solve by Pyomo with CBC
# import numpy library
import numpy as np

# import pyomo library
import pyomo.environ as pyo

# define element size as n
n = 25
# define subset size as m
m = 50

#define subset and its cost(range 1to50)
np.random.seed(1)
subset = np.random.randint(2, size=(m, n))
subset_cost = np.random.randint(1, 50, size=m)

# define pyomo model
model = pyo.ConcreteModel()

# define variables list x which show whether subset is selected or not
model.x = pyo.Var(list(range(m)), domain=pyo.Binary)

# define objective function
model.obj = pyo.Objective(expr=sum(model.x[i]*subset_cost[i] for i in range(m)), sense=pyo.minimize)

# define constraints
model.con = pyo.ConstraintList()
for j in range(n):
    model.con.add(sum(model.x[i]*subset[i][j] for i in range(m)) >= 1)

# execute optimization
solver = pyo.SolverFactory('cbc')
solver.solve(model)

# print result
print('Optimal value: ', pyo.value(model.obj))

# print selected subset
for i in range(m):
    if pyo.value(model.x[i]) == 1:
        print('subset: ', i , 'cost: ', subset_cost[i] , 'elements: ', subset[i])
        
# print selected count of elements
for j in range(n):
    print('element: ', j, 'selected count: ', sum(pyo.value(model.x[i])*subset[i][j] for i in range(m)))

結果 1. CBC

出力
Optimal value:  10.0
subset:  3 cost:  1 elements:  [0 1 1 1 0 1 0 0 1 1 0 1 1 0 1 0 0 1 1 1 0 1 1 0 1]
subset:  21 cost:  4 elements:  [0 1 1 0 0 1 1 0 1 1 0 0 0 0 1 0 0 0 0 1 0 1 0 0 1]
subset:  22 cost:  2 elements:  [1 1 0 0 1 0 0 0 0 0 1 1 1 1 1 0 0 1 1 0 0 0 0 1 0]
subset:  38 cost:  3 elements:  [0 1 1 0 0 1 0 1 1 0 0 1 0 0 1 1 1 1 1 1 1 1 0 1 1]
element:  0 selected count:  1.0
element:  1 selected count:  4.0
element:  2 selected count:  3.0
element:  3 selected count:  1.0
element:  4 selected count:  1.0
...
element:  21 selected count:  3.0
element:  22 selected count:  1.0
element:  23 selected count:  2.0
element:  24 selected count:  3.0

難なく答えが出ます。しかも1秒以下です。 さすが（従来手法がすごすぎてこの時点で萎えますが…）

解法 2．論文の手法

前半の問題設定は解法1と同じです。
PlaceHolderなどを使うともうすこし高速化できるんですが、このテストの問題規模では
サンプリング周りに時間がかかっていたので、手抜きしています。

論文手法
# import numpy library
import numpy as np

# import pyomo library
import pyomo.environ as pyo

# define element size as n
n = 25
# define subset size as m
m = 50

#define subset and its cost(range 1to50)
np.random.seed(1)
subset = np.random.randint(2, size=(m, n))
subset_cost = np.random.randint(1, 50, size=m)

# subset number is m, element number is n
# initialize µ.λ,ρ
µ = 0.5
λ = [0.0 for i in range(n)]
ρ = 1.1

# define variables list x which show whether subset is selected or not
x = Array.create('x', shape=m, vartype='BINARY')
# define sampler
sampler = EmbeddingComposite(DWaveSampler())

# define function to get model of AL
def get_model_AL():
    AL1 = sum([subset_cost[i]*x[i] for i in range(m)])
    AL2 = 0
    for i in range(n):
        AL2 += (-λ[i]-μ/2.0)*sum([x[j] for j in range(m) if subset[j][i] == 1])
        AL2 += μ*sum([x[j]*x[k] for k in range(m) for j in range(k) if subset[j][i] == 1 and subset[k][i] == 1])
    AL = AL1 + AL2
    return  AL.compile()

# main loop
for l in range(40):
    # get model of AL if we use placeholder , we can skip this part
    model_AL = get_model_AL()
    # create bqm from model
    bqm = model_AL.to_bqm()

    # solve problem by D-Wave and decode sampleset
    response = sampler.sample(bqm, num_reads=2000)
    decoded_response = model_AL.decode_sampleset(response)

    # get best solution
    best_x = [decoded_response[0].sample[f"x[{i}]"] for i in range(m)]

    # define function of get c(x)
    cix_break_flag = True
    for i in range(n):
        cix = 1 - sum(subset[j][i] for j in range(m) if best_x[j] == 1)
        if cix > 0:
            λ[i] = λ[i] + μ*cix
            cix_break_flag = False
            
    μ = μ*ρ

    # print result
    print('iteration: ', l)
    print('Optimal value: ', sum(subset_cost[i] for i in range(m) if best_x[i] == 1))
    # all elements are selected or not
    for i in range(n):
        if sum(best_x[j]*subset[j][i] for j in range(m)) == 0:
            print('element: ', i, 'is not selected')
    
    #check break condition
    if cix_break_flag:
        print('break condition is satisfied')
        break

結果 2. D-Wave

処理時間 8.9秒でした。

出力
iteration:  0
Optimal value:  21
element:  9 is not selected
iteration:  1
Optimal value:  26
break condition is satisfied

なお、結果にはかなりばらつきがあります。

num_reads=2000での結果

試行回 Itr回数 時間 制約違反 目的関数値
1 2回 8.9秒 なし 26
2 3回 16.0秒 なし 63
3 1回 5 秒 なし 70
4 2回 9.8 秒 なし 70
5 1回 5.2 秒 なし 108

num_reads=500 での結果

試行回 Itr回数 時間 制約違反 目的関数値
1 2回 9.5秒 なし 75
2 4回 18.7秒 なし 84
3 1回 13.3秒 なし 68
4 4回 17.7秒 なし 77
5 1回 13.9秒 なし 29
2
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?