| Hosted by CoCalc | Download
1
##########################################
2
# Linear Programming Scripts
3
##########################################
4
#
5
# Johann Thiel
6
# ver 12.18.18
7
# Functions to solve linear
8
# programming problems.
9
#
10
##########################################
11
12
##########################################
13
# Necessary modules
14
##########################################
15
from sage.numerical.backends.generic_backend import get_solver
16
import sage.numerical.backends.glpk_backend as backend
17
from tabulate import tabulate
18
##########################################
19
20
##########################################
21
# Generic linear programming solver that
22
# produces a sensitivity report
23
##########################################
24
# var = list of variable names
25
# con = list of constraint names
26
# ob = list of coefficients of objective
27
# functions
28
# M = matrix of constraint coefficients
29
# inq = list of inequality direction
30
# 1 for '<=', -1 for '>=' and
31
# 0 for '='.
32
# bnd = list of constraint bounds
33
# mx = Boolean to determine if the
34
# problem is a maximization (True)
35
# or minimization (False) problem.
36
# Default is set to maximization.
37
##########################################
38
def lp(var, con, ob, M, inq, bnd, mx=True):
39
# loads the solver
40
q = get_solver(solver = 'GLPK')
41
# sets solver to min, max otherwise
42
if not mx:
43
q.set_sense(-1)
44
# adds all variables
45
for v in va:
46
q.add_variable(name=v)
47
# adds all constraints
48
for i in range(len(M)):
49
if inq[i] == 1:
50
q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), None, bnd[i], con[i])
51
elif inq[i] == -1:
52
q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], None, con[i])
53
else:
54
q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], bnd[i], con[i])
55
# sets objective
56
q.set_objective(ob)
57
q.solver_parameter(backend.glp_simplex_or_intopt, backend.glp_simplex_only)
58
# solves the problem
59
q.solve()
60
# produces sensitivity report
61
q.print_ranges()
62
##########################################
63
64
##########################################
65
# Generic linear programming solver for
66
# integer programming problems (produces
67
# no sensitivity report)
68
##########################################
69
# var = list of variable names
70
# con = list of constraint names
71
# ob = list of coefficients of objective
72
# functions
73
# M = matrix of constraint coefficients
74
# inq = list of inequality direction
75
# 1 for '<=', -1 for '>=' and
76
# 0 for '='.
77
# bnd = list of constraint bounds
78
# mx = Boolean to determine if the
79
# problem is a maximization (True)
80
# or minimization (False) problem.
81
# Default is set to maximization.
82
##########################################
83
def lpInt(var, con, ob, M, inq, bnd, mx=True):
84
# loads the solver
85
q = get_solver(solver = 'GLPK')
86
# sets solver to min, max otherwise
87
if not mx:
88
q.set_sense(-1)
89
# adds all variables
90
for v in var:
91
q.add_variable(integer=True, name=v)
92
# adds all constraints
93
for i in range(len(M)):
94
if inq[i] == 1:
95
q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), None, bnd[i], con[i])
96
elif inq[i] == -1:
97
q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], None, con[i])
98
else:
99
q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], bnd[i], con[i])
100
# sets objective
101
q.set_objective(ob)
102
q.solver_parameter(backend.glp_simplex_then_intopt)
103
# solves the problem
104
q.solve()
105
# The following lines all produce an output report
106
if mx:
107
st = ' (max) '
108
else:
109
st = ' (min) '
110
sol = [q.get_variable_value(i) for i in range(q.ncols())]
111
print 'Optimal'+st+'value: ', q.get_objective_value()
112
print ''
113
results = [[a,b] for a,b in zip(var,sol)]
114
print tabulate(results, headers=['Variable', 'Value'])
115
print ''
116
slack = [[q.row_name(j), round(max(0,abs(bnd[j]-sum([a*b for a,b in zip(sol,M[j])]))),3)] for j in range(q.nrows())]
117
print tabulate(slack, headers=['Constraint', 'Slack'])
118
##########################################
119
120
##########################################
121
# Generic linear programming solver for
122
# binary integer programming problems
123
# (produces no sensitivity report)
124
##########################################
125
#
126
# var = list of variable names
127
# con = list of constraint names
128
# ob = list of coefficients of objective
129
# functions
130
# M = matrix of constraint coefficients
131
# inq = list of inequality direction
132
# 1 for '<=', -1 for '>=' and
133
# 0 for '='.
134
# bnd = list of constraint bounds
135
# mx = Boolean to determine if the
136
# problem is a maximization (True)
137
# or minimization (False) problem.
138
# Default is set to maximization.
139
##########################################
140
def lpBin(var, con, ob, M, inq, bnd, mx=True):
141
# loads the solver
142
q = get_solver(solver = 'GLPK')
143
# sets solver to min, max otherwise
144
if not mx:
145
q.set_sense(-1)
146
# adds all variables
147
for v in var:
148
q.add_variable(binary=True, name=v)
149
# adds all constraints
150
for i in range(len(M)):
151
if inq[i] == 1:
152
q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), None, bnd[i], con[i])
153
elif inq[i] == -1:
154
q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], None, con[i])
155
else:
156
q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], bnd[i], con[i])
157
# sets objective
158
q.set_objective(ob)
159
q.solver_parameter(backend.glp_simplex_then_intopt)
160
# solves the problem
161
q.solve()
162
# The following lines all produce an output report
163
if mx:
164
st = ' (max) '
165
else:
166
st = ' (min) '
167
sol = [q.get_variable_value(i) for i in range(q.ncols())]
168
print 'Optimal'+st+'value: ', q.get_objective_value()
169
print ''
170
results = [[a,b] for a,b in zip(va,sol)]
171
print tabulate(results, headers=['Variable', 'Value'])
172
print ''
173
slack = [[q.row_name(j), round(max(0,abs(bnd[j]-sum([a*b for a,b in zip(sol,M[j])]))),3)] for j in range(q.nrows())]
174
print tabulate(slack, headers=['Constraint', 'Slack'])
175
##########################################
176