数理最適化で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])))
あなたの回答
tips
プレビュー