Source Code Examples:#
Source Code for the Experiment of Optimizing over a Circle#
from math import cos, sin, pi
import random
import time
import sys
import gurobipy as gp
# Work on a circle defined on a million constraints
t0 = time.time()
n = 1024 * 1024
m = gp.Model("Circle Optimization")
X = m.addVars(2, lb=-2, ub=2)
Wb = 0
Wc = 0
Wd = 0
maxdiff = 0
niter = 0
margin = 1.01
m.addConstrs(
X[0] * cos((2 * pi * i) / n) + X[1] * sin((2 * pi * i) / n) <= 1 for i in range(n)
)
print("Added 2 Vars and %d constraints in %.2f seconds" % (n, time.time() - t0))
m.Params.OutputFlag = 0
m.Params.Presolve = 0
# Now select random objectives and optimize. Resulting optimal
# solution must be in the circle
for i in range(4096):
theta = 2 * pi * random.random()
a = cos(theta)
b = sin(theta)
m.setObjective(X[0] * a + X[1] * b)
m.optimize()
niter += m.IterCount
# See how far is the solution from the boundary of a circle of
# radius one, if we minimize (a,b) the optimal solution should be (-a,-b)
error = (X[0].X + a) * (X[0].X + a) + (X[1].X + b) * (X[1].X + b)
# Display most inacurate solution
if (
error > margin * maxdiff
or m.BoundVio > margin * Wb
or m.ConstrVio > margin * Wc
or m.DualVio > margin * Wd
):
maxdiff = max(maxdiff, error)
Wb = max(Wb, m.BoundVio)
Wc = max(Wb, m.ConstrVio)
Wd = max(Wd, m.DualVio)
print(
"Errors: %g %g %g %g Iter %d %d Kappa %g"
% (maxdiff, Wb, Wc, Wd, i, niter, m.KappaExact)
)
sys.stdout.flush()
Source Code for the Experiment on a Thin Feasible Region#
import random
import sys
import gurobipy as gp
from gurobipy import GRB
# Test the effect of small perturbations on the optimal solutions
# for a problem with a thin feasible region
rhs = 1e3
m = gp.Model("Thin line Optimization")
x = m.addVar(obj=1)
y = m.addVar(obj=0, lb=-GRB.INFINITY, ub=GRB.INFINITY)
c1 = m.addConstr(1e-5 * y + 1e-0 * x <= rhs)
c2 = m.addConstr(-1e-5 * y + 1e-0 * x <= rhs)
m.Params.OutputFlag = 0
m.Params.Presolve = 0
m.optimize()
xval = x.X
yval = y.X
maxdiff = 0
for i in range(1024 * 1024):
c1.Rhs = rhs + 2e-6 * random.random()
c2.Rhs = rhs + 2e-6 * random.random()
x.Obj = 1 + 2e-6 * random.random()
y.Obj = 0 + 2e-6 * random.random()
m.optimize()
x2val = x.X
y2val = y.X
error = (xval - x2val) * (xval - x2val) + (yval - y2val) * (yval - y2val)
if error > 1e-5 + maxdiff:
print(
"New maxdiff %g Iter %d Kappa %g Violations: %g %g %g"
% (error, i, m.KappaExact, m.BoundVio, m.ConstrVio, m.DualVio)
)
sys.stdout.flush()
maxdiff = error
Source Code for the Experiment with Column Scalings#
import random
import argparse
import gurobipy as gp
from gurobipy import GRB
# Use parameters for greater flexibility
parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument("-f", "--infile", help="Problem File", default=None, required=True)
parser.add_argument("-s", "--scale", help="Scaling Factor", type=float, default=10000.0)
args = parser.parse_args()
# Load input problem
m = gp.read(args.infile)
# Scale domain of all columns randomly in the given domain
for var in m.getVars():
if var.vtype == GRB.CONTINUOUS:
scale = random.uniform(args.scale / 2.0, args.scale * 2.0)
flip = random.randint(0, 3)
if flip == 0:
scale = 1.0
elif flip == 1:
scale = 1.0 / scale
col = m.getCol(var)
for i in range(col.size()):
coeff = col.getCoeff(i)
row = col.getConstr(i)
m.chgCoeff(row, var, coeff * scale)
var.obj = var.obj * scale
if var.lb > -GRB.INFINITY:
var.lb = var.lb / scale
if var.ub < GRB.INFINITY:
var.ub = var.ub / scale
# Optimize
m.optimize()
if m.Status == GRB.OPTIMAL:
print("Kappa: %e\n" % m.KappaExact)