Source code for the experiment of optimizing over a circle


from gurobipy import *
from math import *
import random
import time
import sys

# Work on a circle defined on a million constraints
t0      = time.time()
n       = 1024 * 1024
m       = 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()