Lineare Optimierung in Python
Contents
Lineare Optimierung in Python#
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
SciPy#
open source solver and interface
Links:
Install under Anaconda:
conda install scipy
Example:
\[\begin{split}\begin{align}
\text{max.}\;
5 x_1 + 4 x_2 & \\
\text{s.t.}\;
x_1 & \leq 6 \\
\frac{1}{4}x_1 + x_2 & \leq 6 \\
3x_1 + 2x_2 & \leq 22 \\
x_1 & \geq 0 \\
x_2 & \geq 0
\end{align}\end{split}\]
c = np.array([-5, -4])
A = np.array([[ 1, 0],
[0.25, 1],
[ 3, 2]])
b = np.array([6, 6, 22])
from scipy.optimize import linprog
res = linprog(c, A_ub=A, b_ub=b, bounds= (0, None))
print("Resultate gesamt:\n", res)
print("optimaler Wert:", -res.fun)
print("optimaler Punkt:", res.x)
Resultate gesamt:
con: array([], dtype=float64)
crossover_nit: 0
eqlin: marginals: array([], dtype=float64)
residual: array([], dtype=float64)
fun: -40.0
ineqlin: marginals: array([-0. , -0.8, -1.6])
residual: array([2., 0., 0.])
lower: marginals: array([0., 0.])
residual: array([4., 5.])
message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
nit: 2
slack: array([2., 0., 0.])
status: 0
success: True
upper: marginals: array([0., 0.])
residual: array([inf, inf])
x: array([4., 5.])
optimaler Wert: 40.0
optimaler Punkt: [4. 5.]
CVXOPT#
open source interface to own solver and other solvers
Links:
Install under Anaconda:
conda install cvxopt
Example 1:
import cvxopt as cx
A = cx.matrix([ [-1.0, -1.0, 0.0, 1.0], [1.0, -1.0, -1.0, -2.0] ])
b = cx.matrix([ 1.0, -2.0, 0.0, 4.0 ])
c = cx.matrix([ 2.0, 1.0 ])
sol = cx.solvers.lp(c, A, b)
print(sol['x'])
pcost dcost gap pres dres k/t
0: 2.6471e+00 -7.0588e-01 2e+01 8e-01 2e+00 1e+00
1: 3.0726e+00 2.8437e+00 1e+00 1e-01 2e-01 3e-01
2: 2.4891e+00 2.4808e+00 1e-01 1e-02 2e-02 5e-02
3: 2.4999e+00 2.4998e+00 1e-03 1e-04 2e-04 5e-04
4: 2.5000e+00 2.5000e+00 1e-05 1e-06 2e-06 5e-06
5: 2.5000e+00 2.5000e+00 1e-07 1e-08 2e-08 5e-08
Optimal solution found.
[ 5.00e-01]
[ 1.50e+00]
Example 2:
import cvxopt.modeling as cxm
x = cxm.variable(size=3, name='a')
c1 = (x[0] <= 1)
c2 = (x >= 0)
c3 = (sum(x) == 2)
print(c1)
scalar inequality
constraint function:
affine function of length 1
constant term:
[-1.00e+00]
linear term: linear function of length 1
coefficient of variable(3,'a'):
[ 1.00e+00 0 0 ]
x = cxm.variable()
y = cxm.variable()
c1 = ( 2*x+y <= 3 )
c2 = ( x+2*y <= 3 )
c3 = ( x >= 0 )
c4 = ( y >= 0 )
lp1 = cxm.op(-4*x-5*y, [c1,c2,c3,c4])
lp1.solve()
lp1.status
pcost dcost gap pres dres k/t
0: -8.1000e+00 -1.8300e+01 4e+00 0e+00 8e-01 1e+00
1: -8.8055e+00 -9.4357e+00 2e-01 2e-16 4e-02 3e-02
2: -8.9981e+00 -9.0049e+00 2e-03 6e-16 5e-04 4e-04
3: -9.0000e+00 -9.0000e+00 2e-05 2e-16 5e-06 4e-06
4: -9.0000e+00 -9.0000e+00 2e-07 1e-16 5e-08 4e-08
Optimal solution found.
'optimal'
print(lp1.objective.value())
print(x.value)
print(y.value)
print(c1.multiplier.value)
print(c1.value())
[-9.00e+00]
[ 1.00e+00]
[ 1.00e+00]
[ 1.00e+00]
[-1.10e-07]
Optlang#
open source interface only, install e. g. scipy or glpk as solver
Links:
Install under Anaconda:
conda install -c conda-forge optlang
or
pip install optlang
Install Solvers:
see below
Example:
import optlang
optlang.available_solvers
{'GUROBI': True,
'GLPK': True,
'MOSEK': False,
'CPLEX': False,
'COINOR_CBC': False,
'SCIPY': True,
'OSQP': True}
# from optlang.scipy_interface import Model, Variable, Constraint, Objective
from optlang.glpk_interface import Model, Variable, Constraint, Objective
# from optlang.gurobi_interface import Model, Variable, Constraint, Objective
# from optlang.cplex_interface import Model, Variable, Constraint, Objective
# All the (symbolic) variables are declared, with a name and optionally a lower
# and/or upper bound.
x1 = Variable('x1', lb=0)
x2 = Variable('x2', lb=0)
x3 = Variable('x3', lb=0)
my_vars = [x1, x2, x3]
# A constraint is constructed from an expression of variables and
# a lower and/or upper bound (lb and ub).
# c1 = Constraint(x1 + x2 + x3, ub=100)
c1 = Constraint( sum( my_var for my_var in my_vars), ub=100)
c2 = Constraint(10 * x1 + 4 * x2 + 5 * x3, ub=600)
c3 = Constraint(2 * x1 + 2 * x2 + 6 * x3, ub=300)
# An objective can be formulated
obj = Objective(10 * x1 + 6 * x2 + 4 * x3, direction='max')
# Variables, constraints and objective are combined in a Model object,
# which can subsequently be optimized.
model = Model(name='Simple model')
model.objective = obj
model.add([c1, c2, c3])
status = model.optimize()
print("status:", model.status)
print("objective value:", model.objective.value)
print("----------")
for var_name, my_var in model.variables.iteritems():
print(var_name, "=", my_var.primal)
status:Scaling...
optimal
objective value: 733.3333333333333
----------
x1 = 33.33333333333333
x2 = 66.66666666666667
x3 = 0.0
A: min|aij| = 1.000e+00 max|aij| = 1.000e+00 ratio = 1.000e+00
Problem data seem to be well scaled
print("shadow prices:", model.shadow_prices.values())
shadow prices: odict_values([3.3333333333333335, 0.6666666666666666, 0.0])
print("reduced costs:",model.reduced_costs)
reduced costs: OrderedDict([('x1', 0.0), ('x2', 0.0), ('x3', -2.6666666666666665)])
PuLP#
open source interface only, install e. g. glpk as solver
Links:
Install under Anaconda:
pip install pulp
Example: https://coin-or.github.io/pulp/CaseStudies/a_blending_problem.html
import pulp
# Create the 'prob' variable to contain the problem data
prob = pulp.LpProblem("The Whiskas Problem", pulp.LpMinimize)
# The 2 variables Beef and Chicken are created with a lower limit of zero
x1 = pulp.LpVariable("ChickenPercent", 0, None, pulp.LpInteger)
x2 = pulp.LpVariable("BeefPercent", 0)
# The objective function is added to 'prob' first
prob += 0.013*x1 + 0.008*x2, "Total Cost of Ingredients per can"
# The five constraints are entered
prob += x1 + x2 == 100, "PercentagesSum"
prob += 0.100*x1 + 0.200*x2 >= 8.0, "ProteinRequirement"
prob += 0.080*x1 + 0.100*x2 >= 6.0, "FatRequirement"
prob += 0.001*x1 + 0.005*x2 <= 2.0, "FibreRequirement"
prob += 0.002*x1 + 0.005*x2 <= 0.4, "SaltRequirement"
# The problem is solved using PuLP's or your choice of Solver
prob.solve()
#prob.solve(pulp.GLPK())
#prob.solve(pulp.GUROBI())
# The status of the solution is printed to the screen
print("Status:", pulp.LpStatus[prob.status])
# Each of the variables is printed with it's resolved optimum value
for v in prob.variables():
print(v.name, "=", v.varValue)
# The optimised objective function value is printed to the screen
print("Total Cost of Ingredients per can = ", pulp.value(prob.objective))
Welcome to the CBC MILP Solver
Version: 2.10.3
Build Date: Dec 15 2019
command line - /home/kr/.local/lib/python3.10/site-packages/pulp/apis/../solverdir/cbc/linux/64/cbc /tmp/b595477c1df54abcab51f6dc5631a773-pulp.mps timeMode elapsed branch printingOptions all solution /tmp/b595477c1df54abcab51f6dc5631a773-pulp.sol (default strategy 1)
At line 2 NAME MODEL
At line 3 ROWS
At line 10 COLUMNS
At line 25 RHS
At line 31 BOUNDS
At line 33 ENDATA
Problem MODEL has 5 rows, 2 columns and 10 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 0.966667 - 0.00 seconds
Cgl0004I processed model has 0 rows, 0 columns (0 integer (0 of which binary)) and 0 elements
Cbc3007W No integer variables - nothing to do
Cuts at root node changed objective from 0.97 to -1.79769e+308
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Result - Optimal solution found
Objective value: 0.97000000
Enumerated nodes: 0
Total iterations: 1
Time (CPU seconds): 0.00
Time (Wallclock seconds): 0.00
Option for printingOptions changed from normal to all
Total time (CPU seconds): 0.00 (Wallclock seconds): 0.01
Status: Optimal
BeefPercent = 66.0
ChickenPercent = 34.0
Total Cost of Ingredients per can = 0.97
/home/kr/.local/lib/python3.10/site-packages/pulp/pulp.py:1352: UserWarning: Spaces are not permitted in the name. Converted to '_'
warnings.warn("Spaces are not permitted in the name. Converted to '_'")
PYMPROG#
open source interface only, needs glpk as solver
Links:
Install under Anaconda:
pip install pymprog
Example:
import pymprog as mp
c = (10, 6, 4)
A = [ ( 1, 1, 1),
( 9, 4, 5),
( 2, 2, 6) ]
b = (10, 60, 30)
mp.begin('basic') # begin modelling
x = mp.var('x', 3) # create 3 variables
mp.maximize(sum(c[i]*x[i] for i in range(3)))
for i in range(3):
sum(A[i][j]*x[j] for j in range(3)) <= b[i]
mp.solve() # solve the model
mp.sensitivity() # sensitivity report
mp.end() #Good habit: do away with the model
PyMathProg 1.0 Sensitivity Report Created: 2022/08/29 Mon 14:51PM
================================================================================
Variable Activity Dual.Value Obj.Coef Range.From Range.Till
--------------------------------------------------------------------------------
*x[0] 4 0 10 6 13.5
*x[1] 6 0 6 4.44444 10
x[2] 0 -2.8 4 -inf 6.8
================================================================================
================================================================================
Constraint Activity Dual.Value Lower.Bnd Upper.Bnd RangeLower RangeUpper
--------------------------------------------------------------------------------
R1 10 2.8 -inf 10 6.66667 15
R2 60 0.8 -inf 60 40 90
*R3 20 0 -inf 30 20 20
================================================================================
model('basic') is not the default model.
x[0].primal
4.0
GLPK#
open source solver
Links:
Install under Anaconda:
conda install -c conda-forge glpk
Gurobi, gurobipy#
commercial solver, academic license available
gurobipy is an interface for python
Links:
Install:
After installation of Gurobi, one can install the Python package gurobipy under Anaconda with:
conda config --add channels http://conda.anaconda.org/gurobi
conda install gurobi
Example:
import gurobipy as gp
model = gp.Model("Test")
x = model.addVar(lb=0,name = "x")
y = model.addVar(lb=0,name = "y")
model.setObjective(3*x + 4*y, gp.GRB.MAXIMIZE)
model.addConstr(x + y <= 10)
model.addConstr(2*x + 3*y <= 20)
model.addConstr(y<=5)
model.optimize()
for v in model.getVars():
print('%s = %g' % (v.varName, v.x))
print('Obj = %g' % model.objVal)
Set parameter Username
Academic license - for non-commercial use only - expires 2023-07-13
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (linux64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 3 rows, 2 columns and 5 nonzeros
Model fingerprint: 0x9dadc558
Coefficient statistics:
Matrix range [1e+00, 3e+00]
Objective range [3e+00, 4e+00]
Bounds range [0e+00, 0e+00]
RHS range [5e+00, 2e+01]
Presolve removed 1 rows and 0 columns
Presolve time: 0.01s
Presolved: 2 rows, 2 columns, 4 nonzeros
Iteration Objective Primal Inf. Dual Inf. Time
0 4.0000000e+01 7.500000e+00 0.000000e+00 0s
2 3.0000000e+01 0.000000e+00 0.000000e+00 0s
Solved in 2 iterations and 0.01 seconds (0.00 work units)
Optimal objective 3.000000000e+01
x = 10
y = 0
Obj = 30
Pi … Dual value (also known as the shadow price)
List of all attributes.
model.printAttr(['ConstrName', 'Pi', 'Slack'])
Constraint ConstrName Pi Slack
---------------------------------------------------
R0 R0 1 0
R1 R1 1 0
R2 R2 0 5
CPLEX#
commercial solver and interface, academic license available
Links:
Cplex: https://www.ibm.com/products/ilog-cplex-optimization-studio
academic license available under https://www.ibm.com/academic/topic/data-science