/
workspace.py
433 lines (356 loc) · 13.3 KB
/
workspace.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
from __future__ import print_function
import numpy as np
# Import osqp solver
import osqp
# Miosqp files
from miosqp.node import Node
from miosqp.constants import MI_UNSOLVED, \
MI_PRIMAL_INFEASIBLE, \
MI_DUAL_INFEASIBLE, \
MI_SOLVED, \
MI_MAX_ITER_FEASIBLE, \
MI_MAX_ITER_UNSOLVED
class Workspace(object):
"""
Workspace class
Attributes
----------
data: class
miqp problem data
settings: dictionary
branch and bound settings
qp_settings: dictionary
branch and bound settings
solver: class
QP solver class
leaves: list
leaves in the tree
Other internal variables
------------------------
iter_num: int
number of iterations
osqp_iter: int
number of osqp admm iterations
osqp_iter_avg: int
average number of osqp admm iterations
osqp_solve_time: float
total time required to solve OSQP problems
run_time: double
runtime of the algorithm
first_run: int
has the problem been solved already?
upper_glob: double
global upper bound
lower_glob: double
global lower bound
status: int
MIQP solver status
x: numpy array
current best solution
"""
def __init__(self, data, settings, qp_settings=None):
self.data = data
self.settings = settings
# Setup OSQP solver instance
self.solver = osqp.OSQP()
self.qp_settings = qp_settings
if self.qp_settings is None:
self.qp_settings = {}
self.solver.setup(self.data.P, self.data.q, self.data.A,
self.data.l, self.data.u, **qp_settings)
# Define other internal variables
self.first_run = 1
self.iter_num = 1
self.osqp_solve_time = 0.
self.osqp_iter = 0
self.osqp_iter_avg = 0
self.lower_glob = -np.inf
self.status = MI_UNSOLVED
# Define root node
root = Node(self.data, self.data.l, self.data.u, self.solver)
# Define leaves at the beginning (only root)
self.leaves = [root]
# Set initial solution and objective value
self.upper_glob = np.inf
self.x = np.empty(self.data.n)
# Define timings
self.setup_time = 0.
self.solve_time = 0.
self.run_time = 0.
def set_x0(self, x0):
"""
Set initial condition
"""
# Suppose that algorithm hasn't been run yet.
# Tree has only one leaf (root)
root = self.leaves[0]
# Add initial solution and objective value
if self.satisfies_lin_constraints(x0, root.l, root.u) and \
self.is_int_feas(x0, root):
self.x = x0
self.upper_glob = self.data.compute_obj_val(x0)
else:
print('Invalid initial solution!\n')
self.upper_glob = np.inf
self.x = np.empty(self.data.n)
def can_continue(self):
"""
Check if the solver can continue
"""
# Check if there are any leaves left in the list
if not any(self.leaves):
return False
# Check if the number of iterations is within the limit
if not self.iter_num < self.settings['max_iter_bb']:
return False
return True
def choose_leaf(self, tree_explor_rule):
"""
Choose next leaf to branch from the ones that can still be expanded
depending on branch_rule
"""
if tree_explor_rule == 0:
# Depth first: Choose leaf with highest depth
leaf_idx = np.argmax([leaf.depth for leaf in self.leaves])
elif tree_explor_rule == 1:
# Two-phase method.
# - First perform depth first
# - Then choose leaves with best bound
if np.isinf(self.upper_glob):
# First phase
leaf_idx = np.argmax([leaf.depth for leaf in self.leaves])
else:
# Second phase
leaf_idx = np.argmax([leaf.lower for leaf in self.leaves])
else:
raise ValueError('Tree exploring strategy not recognized')
# Get leaf
leaf = self.leaves[leaf_idx]
# Remove leaf from the list of leaves
self.leaves.remove(leaf)
return leaf
def add_left(self, leaf):
"""
Add left node from the current leaf
"""
# Get updated QP interval bounds
l_left = np.copy(leaf.l)
u_left = np.copy(leaf.u)
# x <= floor(x_relaxed)
u_left[leaf.constr_idx] = \
np.floor(leaf.x[leaf.nextvar_idx])
# DEBUG:
if any(l_left > u_left):
import ipdb; ipdb.set_trace()
# Create new leaf
new_leaf = Node(self.data, l_left, u_left, self.solver,
depth=leaf.depth + 1, lower=leaf.lower,
x0=leaf.x, y0=leaf.y)
# Add leaf to the leaves list
self.leaves.append(new_leaf)
def add_right(self, leaf):
"""
Add right node from the current leaf
"""
# Get updated QP interval bounds
l_right = np.copy(leaf.l)
u_right = np.copy(leaf.u)
# ceil(x_relaxed) <= x
l_right[leaf.constr_idx] = \
np.ceil(leaf.x[leaf.nextvar_idx])
# DEBUG:
if any(l_right > u_right):
import ipdb; ipdb.set_trace()
# Create new leaf
new_leaf = Node(self.data, l_right, u_right, self.solver,
depth=leaf.depth + 1, lower=leaf.lower,
x0=leaf.x, y0=leaf.y)
# Add leaf to the leaves list
self.leaves.append(new_leaf)
def pick_nextvar(self, leaf):
"""
Pick next variable to branch upon
"""
# Part of solution that is supposed to be integer
x_frac = leaf.x[self.data.i_idx[leaf.frac_idx]]
if self.settings['branching_rule'] == 0:
# Get vector of fractional parts
fract_part = abs(x_frac - np.round(x_frac))
# Get index of max fractional part
next_idx = np.argmax(fract_part)
# Get index of next variable as position in the i_idx vector
nextvar = leaf.frac_idx[next_idx]
else:
raise ValueError('No variable selection rule recognized!')
# Get next variable constraint index
leaf.constr_idx = self.data.m + nextvar
# Get next variable idx
leaf.nextvar_idx = self.data.i_idx[nextvar]
def satisfies_lin_constraints(self, x, l, u):
"""
Check if solution x is within linear constraints
"""
# Check if it satisfies current l and u bounds
z = self.data.A.dot(x)
if any(z < l - self.qp_settings['eps_abs']) | \
any(z > u + self.qp_settings['eps_abs']):
return False
# If we got here, it is integer feasible
return True
def is_int_feas(self, x, leaf):
"""
Check if current solution is integer feasible
"""
# Part of solution that is supposed to be integer
x_int = x[self.data.i_idx]
# Part of the solution that is still fractional
int_feas_false = abs(x_int - np.round(x_int)) >\
self.settings['eps_int_feas']
# Index of frac parts
leaf.frac_idx = np.where(int_feas_false)[0].tolist()
# Store number of fractional elements (integer infeasible)
leaf.intinf = np.sum(int_feas_false)
if leaf.intinf > 0:
return False
# If we got here, it is integer feasible
return True
def get_integer_solution(self, x):
"""
Round obtained solution to integer feasibility
"""
x_int = np.copy(x)
x_int[self.data.i_idx] = np.round(x[self.data.i_idx])
return x_int
def prune(self):
"""
Prune all leaves whose lower bound is greater than current upper bound
"""
for leaf in self.leaves:
if leaf.lower > self.upper_glob:
self.leaves.remove(leaf)
def bound_and_branch(self, leaf):
"""
Analize result from leaf solution and bound
"""
# Update total number of OSQP ADMM iteration
self.osqp_iter += leaf.num_iter
# Update total time to solve OSQP problems
self.osqp_solve_time += leaf.osqp_solve_time
# 1) If infeasible or unbounded, then return (prune)
if leaf.status == osqp.constant('OSQP_PRIMAL_INFEASIBLE') or \
leaf.status == osqp.constant('OSQP_DUAL_INFEASIBLE'):
return
# 2) If lower bound is greater than upper bound, then return (prune)
if leaf.lower > self.upper_glob:
return
# 3) If integer feasible, then
# - update best solution
# - update best upper bound
# - prune all leaves with lower bound greater than best upper bound
if (self.is_int_feas(leaf.x, leaf)):
# Update best solution so far
self.x = leaf.x
# Update upper bound
self.upper_glob = leaf.lower
# Prune all nodes
self.prune()
return
# 4) If fractional, get integer solution using heuristic.
# If integer solution satisfies linear constraints, then:
# - compute objective value at integer x
# If objective value improves the upper bound
# - update best upper bound
# - prune all leaves with lower bound greater than current one
x_int = self.get_integer_solution(leaf.x)
if self.satisfies_lin_constraints(x_int,
self.data.l, self.data.u):
obj_int = self.data.compute_obj_val(x_int)
if obj_int < self.upper_glob:
self.upper_glob = obj_int
self.x = x_int
self.prune()
# 5) If we got here, branch current leaf producing two children
self.branch(leaf)
# 6) Update lower bound with minimum between lower bounds
self.lower_glob = min([leaf.lower for leaf in self.leaves])
def branch(self, leaf):
"""
Branch current leaf according to branching_rule
"""
# Branch obtaining two children and add to leaves list
# Pick next variable to branch upon
self.pick_nextvar(leaf)
# Add left node to the leaves list
self.add_left(leaf)
# Add right node to the leaves list
self.add_right(leaf)
def get_return_status(self):
"""
Get return status for MIQP solver
"""
if self.iter_num < self.settings['max_iter_bb']: # Finished
if self.upper_glob != np.inf:
self.status = MI_SOLVED
else:
if self.upper_glob >= 0: # +inf
self.status = MI_PRIMAL_INFEASIBLE
else: # -inf
self.status = MI_DUAL_INFEASIBLE
else: # Hit maximum number of iterations
if self.upper_glob != np.inf:
self.status = MI_MAX_ITER_FEASIBLE
else:
if self.upper_glob >= 0: # +inf
self.status = MI_MAX_ITER_UNSOLVED
else: # -inf
self.status = MI_DUAL_INFEASIBLE
def get_return_solution(self):
"""
Get exact mixed-integer solution
"""
if self.status == MI_SOLVED or \
self.status == MI_MAX_ITER_FEASIBLE:
# Part of solution that is supposed to be integer
self.x[self.data.i_idx] = \
np.round(self.x[self.data.i_idx])
def print_headline(self):
"""
Print headline
"""
print(" Nodes | Current Node | Objective Bounds | Cur Node")
print("Explr\tUnexplr\t| Obj\tDepth\tIntInf | Lower\t Upper\t Gap | Iter")
def print_progress(self, leaf):
"""
Print progress at each iteration
"""
if self.upper_glob == np.inf:
gap = " --- "
else:
gap = "%8.2f%%" % \
((self.upper_glob - self.lower_glob)/abs(self.lower_glob)*100)
if leaf.status == osqp.constant('OSQP_PRIMAL_INFEASIBLE') or \
leaf.status == osqp.constant('OSQP_DUAL_INFEASIBLE'):
obj = np.inf
else:
obj = leaf.lower
if leaf.intinf is None:
intinf = " ---"
else:
intinf = "%5d" % leaf.intinf
print("%4d\t%4d\t %10.2e\t%4d\t%s\t %10.2e\t%10.2e\t%s\t%5d" %
(self.iter_num, len(self.leaves), obj,
leaf.depth, intinf, self.lower_glob,
self.upper_glob, gap, leaf.num_iter), end='')
if leaf.status == osqp.constant('OSQP_MAX_ITER_REACHED'):
print("!")
else:
print("")
def print_footer(self):
"""
Print final statistics
"""
print("\n")
print("Status: %s" % self.status)
if self.status == MI_SOLVED:
print("Objective bound: %6.3e" % self.upper_glob)
print("Total number of OSQP iterations: %d" % self.osqp_iter)