質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

0回答

340閲覧

数理最適化、分岐限定法で0-1整数2次計画

NORMSENSEI

総合スコア11

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

0クリップ

投稿2020/06/29 09:01

数理最適化で0-1整数2次計画問題を解きかたを調べています。
下記のような問題の数理最適化です。
イメージ説明

どうやら調べるとGroubiOptimizer,CPLEXといった有償のソルバーで解くのが一般的なようですが、お金がないのでフリーのものがないか調査しました。
2次計画のフリーのライブラリーにcvxoptというものを見つけましたが、
実数値の2次計画は扱えても、0-1整数の2次計画を扱えないようでした。
結局pythonのライブラリーでもよいのが見つけられなかったので、
ナップザック問題で使われている分岐限定法というアルゴリズムを改造して作ってみたのですが、
なかなか問題が解き終わらず、そもそもうまく解けているのかよくわからない状態です。

1.下記のソース(0-1整数2次計画を分岐限定法で解く)の間違っている箇所がないか
2.フリーで0-1整数2次計画を扱えるpythonライブラリーはないか。
(実数の2次計画ではなく)

python

1import cvxpy as cvx 2import copy 3from heapq import * 4import numpy as np 5import itertools 6counter = itertools.count() 7 8class BBTreeNode(): 9 def __init__(self, vars = set(), constraints = [], objective=0, bool_vars=set()): 10 self.vars = vars 11 self.constraints = constraints 12 self.objective = objective 13 self.bool_vars = bool_vars 14 self.children = [] 15 def buildProblem(self): 16 prob = cvx.Problem(cvx.Minimize(self.objective), self.constraints) #i put Minimize, just so you know that I'm assuming it 17 return prob 18 def is_integral(self): 19 return all([abs(v.value - 1) <= 1e-3 or abs(v.value - 0) <= 1e-3 for v in self.bool_vars]) 20 def branch(self): 21 children = [] 22 for b in [0,1]: 23 n1 = copy.deepcopy(self) #yeesh. Not good performance wise, but is simple implementation-wise 24 v = n1.heuristic() #dangerous what if they don't do the same one? I need to do it here though because I need access to copied v. 25 n1.constraints.append( v == b ) # add in the new binary constraint 26 n1.children = [] 27 n1.bool_vars.remove(v) #remove binary constraint from bool var set 28 n1.vars.add(v) #and add it into var set for later inspection of answer 29 #self.children.append(n1) # eventually I might want to keep around the entire search tree. I messed this up though 30 children.append(n1) 31 return children 32 def heuristic(self): 33 # a basic heuristic of taking the ones it seems pretty sure about 34 return min([(min(1 - v.value, v.value) , i, v) for i, v in enumerate(self.bool_vars)])[2] 35 def bbsolve(self): 36 root = self 37 res = root.buildProblem().solve() 38 heap = [(res, next(counter), root)] 39 bestres = 1e20 # 1e20 a big arbitrary initial best objective value 40 bestnode = root # initialize bestnode to the root 41 print(heap) 42 nodecount = 0 43 while len(heap) > 0: 44 nodecount += 1 # for statistics 45 print("Heap Size: ", len(heap)) 46 _, _, node = heappop(heap) 47 prob = node.buildProblem() 48 res = prob.solve() 49 print("Result: ", res) 50 if prob.status not in ["infeasible", "unbounded"]: 51 if res > bestres - 1e-3: #even the relaxed problem sucks. forget about this branch then 52 print("Relaxed Problem Stinks. Killing this branch.") 53 pass 54 elif node.is_integral(): #if a valid solution then this is the new best 55 print("New Best Integral solution.") 56 bestres = res 57 bestnode = node 58 else: #otherwise, we're unsure if this branch holds promise. Maybe it can't actually achieve this lower bound. So branch into it 59 new_nodes = node.branch() 60 for new_node in new_nodes: 61 heappush(heap, (res, next(counter), new_node ) ) # using counter to avoid possible comparisons between nodes. It tie breaks 62 print("Nodes searched: ", nodecount) 63 return bestres, bestnode 64 65 66# Generate a random non-trivial quadratic program. 67m = 2 68n = 100 69p = 10 70 71np.random.seed(1) 72P = np.random.randn(n, n) 73P = P.T @ P 74q = np.random.randn(n) 75G = np.random.randn(m, n) 76h = G @ np.random.randn(n) 77A = np.random.randn(p, n) 78b = np.random.randn(p) 79 80# Define and solve the CVXPY problem. 81x = cvx.Variable(n) 82objective = (1/2)*cvx.quad_form(x, P) + q.T @ x 83bool_vars = {x[i] for i in range(n)} 84root = BBTreeNode(constraints = [A@x==b,G@x<= h,x<=1,0<= x], objective= objective, bool_vars = bool_vars) 85res, sol = root.bbsolve() 86print(sorted(list([(v.name(), v.value) for v in sol.bool_vars] + [(v.name(), v.value) for v in sol.vars])))

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

まだ回答がついていません

会員登録して回答してみよう

アカウントをお持ちの方は

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問