Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

1
##########################################
2
# Linear Programming Scripts
3
##########################################
4
#
5
# Johann Thiel
6
# ver 11.17.24
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 var:
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