-
Notifications
You must be signed in to change notification settings - Fork 0
/
fold_equalities.py
75 lines (54 loc) · 2.19 KB
/
fold_equalities.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
# Given matrices A, C, and vectors b, d representing the polytope
# Ax <= b
# Cx = d, create a derived representation A'x' <= b' where any feasible
# point x' can be translated into an x so that Ax <= b, Cx = d.
# Example problem:
# A = |-1 0 0 0 | b = |0| C = |1 1 1 0| d = | 50|
# | 0 -1 0 0 | |0| |2 3 4 0| |158|
# | 0 0 -1 0 | |0| |4 3 2 0| |142|
# | 0 0 0 -1 | |0|
# | 0 0 0 1 | |1|
import numpy as np
from scipy.linalg import null_space, lstsq
def incorporate_equality(A, b, C, d):
# Now use the following strategy: Any solution to Cx=d is given by
# Fz + x_0, where x_0 is some solution to Cx=d, and F is the null
# space of C.
# Much easier than mucking about with backsubstitution.
F = null_space(C)
# This makes it easier to read the transformed matrices.
F = F / F.max(axis = 0)
# Any solution to Cx=d will do, but this is the only that scipy
# natively supports.
x_0 = lstsq(C, d)[0]
# The reduced problem is AFz <= b - Ax0.
# To turn z to x, x = Fz + x_0.
Ax0 = np.dot(A, x_0)
return np.dot(A, F), b - Ax0, F, x_0
if __name__ == "__main__":
A = np.array([[-1, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0],
[0,0,0,-1], [0, 0, 0, 1]])
b = np.array([0, 0, 0, 0, 1])
C = np.array([[1, 1, 1, 0], [2, 3, 4, 0], [4, 3, 2, 0]])
d = np.array([50, 158, 142])
AF, bn, F, x_0 = incorporate_equality(A, b, C, d)
# For verification purposes, consider the linear program
# max cTx s.t. Ax <= b, Cx=d
# with c being [0, 1, 0, 0]
# The optimal solution is x = [0, 42, 8, 0].
# The modified version is
# max cTFz s.t. AFz <= b - Ax0
# The test consists of running an LP against both the unmodified problem
# and the problem with equality constraints incorporated.
from scipy.optimize import linprog
c = np.array([0, -1, 0, 0])
x_opt = linprog(c, A, b, C, d)["x"] # (0, 42, 8, 0)
# We're no longer constrained to nonnegative solutions, so explicitly
# state that we're not because otherwise the test fails.
cF = np.dot(c, F)
z_opt = linprog(cF, AF, bn, bounds=((None, None), (None, None)))["x"]
x_sec_opt = np.dot(F, z_opt) + x_0
if np.all(np.abs(x_sec_opt - x_opt) < 1e-12):
print "Passed test."
else:
print "Failed test."