##########################################1# Linear Programming Scripts2##########################################3#4# Johann Thiel5# ver 11.25.176# Functions to solve linear7# programming problems.8#9##########################################1011##########################################12# Necessary modules13##########################################14from sage.numerical.backends.generic_backend import get_solver15import sage.numerical.backends.glpk_backend as backend16from tabulate import tabulate17##########################################1819##########################################20# Generic linear programming solver that21# produces a sensitivity report22##########################################23# var = list of variable names24# con = list of constraint names25# ob = list of coefficients of objective26# functions27# M = matrix of constraint coefficients28# inq = list of inequality direction29# 1 for '<=' and 0 for '>='30# bnd = list of constraint bounds31# mx = Boolean to determine if the32# problem is a maximization (True)33# or minimization (False) problem.34# Default is set to maximization.35##########################################36def lp(var, con, ob, M, inq, bnd, mx=True):37# loads the solver38q = get_solver(solver = 'GLPK')39# sets solver to min, max otherwise40if not mx:41q.set_sense(-1)42# adds all variables43for v in var:44q.add_variable(name=v)45# adds all constraints46for i in range(len(M)):47if inq[i]:48q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), None, bnd[i], con[i])49else:50q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], None, con[i])51# sets objective52q.set_objective(ob)53q.solver_parameter(backend.glp_simplex_or_intopt, backend.glp_simplex_only)54# solves the problem55q.solve()56# produces sensitivity report57q.print_ranges()58##########################################5960##########################################61# Generic linear programming solver that62# produces a sensitivity report63##########################################64# var = list of variable names65# con = list of constraint names66# ob = list of coefficients of objective67# functions68# M = matrix of constraint coefficients69# inq = list of inequality direction70# 1 for '<=', -1 for '>=' and71# 0 for '='.72# bnd = list of constraint bounds73# mx = Boolean to determine if the74# problem is a maximization (True)75# or minimization (False) problem.76# Default is set to maximization.77##########################################78def lp(var, con, ob, M, inq, bnd, mx=True):79# loads the solver80q = get_solver(solver = 'GLPK')81# sets solver to min, max otherwise82if not mx:83q.set_sense(-1)84# adds all variables85for v in var:86q.add_variable(name=v)87# adds all constraints88for i in range(len(M)):89if inq[i] == 1:90q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), None, bnd[i], con[i])91elif inq[i] == -1:92q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], None, con[i])93else:94q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], bnd[i], con[i])95# sets objective96q.set_objective(ob)97q.solver_parameter(backend.glp_simplex_or_intopt, backend.glp_simplex_only)98# solves the problem99q.solve()100# produces sensitivity report101q.print_ranges()102##########################################103104##########################################105# Generic linear programming solver for106# integer programming problems (produces107# no sensitivity report)108##########################################109# var = list of variable names110# con = list of constraint names111# ob = list of coefficients of objective112# functions113# M = matrix of constraint coefficients114# inq = list of inequality direction115# 1 for '<=', -1 for '>=' and116# 0 for '='.117# bnd = list of constraint bounds118# mx = Boolean to determine if the119# problem is a maximization (True)120# or minimization (False) problem.121# Default is set to maximization.122##########################################123def lpInt(var, con, ob, M, inq, bnd, mx=True):124# loads the solver125q = get_solver(solver = 'GLPK')126# sets solver to min, max otherwise127if not mx:128q.set_sense(-1)129# adds all variables130for v in var:131q.add_variable(integer=True, name=v)132# adds all constraints133for i in range(len(M)):134if inq[i] == 1:135q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), None, bnd[i], con[i])136elif inq[i] == -1:137q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], None, con[i])138else:139q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], bnd[i], con[i])140# sets objective141q.set_objective(ob)142q.solver_parameter(backend.glp_simplex_then_intopt)143# solves the problem144q.solve()145# The following lines all produce an output report146if mx:147st = ' (max) '148else:149st = ' (min) '150sol = [q.get_variable_value(i) for i in range(q.ncols())]151print 'Optimal'+st+'value: ', q.get_objective_value()152print ''153results = [[a,b] for a,b in zip(var,sol)]154print tabulate(results, headers=['Variable', 'Value'])155print ''156slack = [[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())]157print tabulate(slack, headers=['Constraint', 'Slack'])158##########################################159160##########################################161# Generic linear programming solver for162# binary integer programming problems163# (produces no sensitivity report)164##########################################165#166# var = list of variable names167# con = list of constraint names168# ob = list of coefficients of objective169# functions170# M = matrix of constraint coefficients171# inq = list of inequality direction172# 1 for '<=', -1 for '>=' and173# 0 for '='.174# bnd = list of constraint bounds175# mx = Boolean to determine if the176# problem is a maximization (True)177# or minimization (False) problem.178# Default is set to maximization.179##########################################180def lpBin(var, con, ob, M, inq, bnd, mx=True):181# loads the solver182q = get_solver(solver = 'GLPK')183# sets solver to min, max otherwise184if not mx:185q.set_sense(-1)186# adds all variables187for v in var:188q.add_variable(binary=True, name=v)189# adds all constraints190for i in range(len(M)):191if inq[i] == 1:192q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), None, bnd[i], con[i])193elif inq[i] == -1:194q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], None, con[i])195else:196q.add_linear_constraint(list(zip(range(len(M[0])), M[i])), bnd[i], bnd[i], con[i])197# sets objective198q.set_objective(ob)199q.solver_parameter(backend.glp_simplex_then_intopt)200# solves the problem201q.solve()202# The following lines all produce an output report203if mx:204st = ' (max) '205else:206st = ' (min) '207sol = [q.get_variable_value(i) for i in range(q.ncols())]208print 'Optimal'+st+'value: ', q.get_objective_value()209print ''210results = [[a,b] for a,b in zip(var,sol)]211print tabulate(results, headers=['Variable', 'Value'])212print ''213slack = [[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())]214print tabulate(slack, headers=['Constraint', 'Slack'])215##########################################216217218219220221222223224225226227228