Time windows cplex

timewindow.py

# --------------------------------------------------------------------
# The MIT License (MIT)
#
# Copyright (c) 2009-2015  Paola Lizarralde-Bejarano
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
# --------------------------------------------------------------------

import sys

import cplex as cp
from cplex.callbacks import UserCutCallback, LazyConstraintCallback
from cplex.exceptions import CplexError

# --------------------------------------------------------------------

promotional       = False
reduceTimeWindows = True
twojobcuts        = True
maxRows           = 10
archivo           = 2
if archivo == 1:
  databasename  = "atw-repository"
  filecosts     = "costMatrix.txt"
  filetime      = "timeprocessing.txt"
  filewindows   = "timewindows.txt"

if archivo == 2:
  # profe No. Nodes 50
  databasename  = "baldoquin50"
  filecosts     = "datosCostos.txt"
  filetime      = "datosPrepro.txt"
  filewindows   = "datosVentanas.txt"

if archivo == 3:
  # ejemplo No. Nodes 4
  databasename  = "polis4nodes"
  filecosts     = "datosCostos1"
  filetime      = "datosPrepro1"
  filewindows   = "datosVentanas1"

# --------------------------------------------------------------------

costMatrix,timewindows, timeprocessing = [],[],[]

def x(i,j):
  return "X"+str(i)+"m"+str(j)

def y(i,j):
  return "Y"+str(i)+"m"+str(j)

def c(i,j):
  return costMatrix[i][j]

def isEdge(i,j):
  return c(i,j) > 0

def v(i,j):
  return p(i) + c(i,j)

rmemo = {}
def r(i):
  if i in rmemo: return rmemo[i]
  rmemo[i] = timewindows[i][0]
  return rmemo[i]

ddmemo = {}
def dd(i):
  if i in ddmemo: return ddmemo[i]
  ddmemo[i] = timewindows[i][1]
  return ddmemo[i]

def p(i):
  return timeprocessing[i]

# --------------------------------------------------------------------

def tightening():
  # release data adjustment
  n = numNodes
  for k in range(1, n):
    dminus = [i for i in range(n) if isEdge(i,k)]
    if len(dminus) > 0:
      rmemo[k] = max(r(k), min(r(i) + v(i,k) for i in dminus))

  for k in range(1, n):
    dplus = [j for j in range(n) if isEdge(k,j)]
    if len(dplus) > 0:
      rmemo[k] = max(r(k), min(dd(k) , min(r(j)-v(k,j) for j in dplus)))

  # due date adjustment
  for k in range(1, n):
    dminus = [i for i in range(n) if isEdge(i,k)]
    if len(dminus) > 0:
      ddmemo[k] = min(dd(k), max(r(k), max(dd(i) + v(i,k) for i in dminus)))

  for k in range(1, n):
    dplus = [j for j in range(n) if isEdge(k, j)]
    if len(dplus) > 0:
      ddmemo[k] = min(dd(k), max(dd(j) - v(k,j) for j in dplus))

  print "new time windows:"
  for i in range(n):
    print "node {2} : [{0}, {1}]".format(r(i), dd(i), i)

# --------------------------------------------------------------------

costMatrix = []
with open(filecosts) as sys.stdin:
  while True:
    try:
      row =  map(float, raw_input().replace(",", ".").split())
      if promotional:
        row = row[:maxRows]
      costMatrix.append(row)
      if promotional and len(costMatrix) == maxRows:
        break
    except Exception as e:
      break
  assert len(costMatrix) > 0
  assert len(costMatrix) == len(costMatrix[0])


numNodes  = len(costMatrix)
numVars   = numNodes * numNodes

assert numNodes > 0

my_obj = []
for i in xrange(numNodes):
  for d in costMatrix[i]:
    my_obj += [ d * 1.0 ]

assert isinstance(my_obj[0], float)
assert len(my_obj) == numVars

my_colnames = []
for i in xrange(numNodes):
  for j in xrange(numNodes):
    my_colnames += [x(i,j)] # x_{ij} = Ximj

assert len(my_colnames) == numVars

timeprocessing = []
with open(filetime) as sys.stdin:
  while True:
    try:
      timeprocessing.append( int(raw_input()) )
      if len(timeprocessing) == maxRows and promotional:
        break
    except:
      break
  assert len(timeprocessing) == numNodes

timewindows = []
with open(filewindows) as sys.stdin:
  while True:
    try:
      timewindows.append( map(int, raw_input().split()) )
      if len(timewindows) == maxRows and promotional:
        break
    except:
      break
  assert len(timewindows) == numNodes


names_yij = []
for i in xrange(numNodes):
  for j in xrange(numNodes):
    names_yij += [y(i,j)] # x_{ij} = Ximj

assert len(timewindows) == len(costMatrix)

# --------------------------------------------------------------------

class LazyTwoJob(LazyConstraintCallback):

    def __call__(self):

      path = [0]
      t    = {0:0}

      u = 0
      while len(path) < numNodes:
        for i in range(numNodes):
          if i != u:
            varName = x(u,i)
            varYName = y(u,i)
            if abs(self.get_values(varName)) > 1.e-10:
              path.append(i)
              if u == 0:
                t[i] = r(i)
              else:
                t[i] = max( t[u] + v(u,i), r(i) )
              u = i

      print "path:", ' > '.join(map(str,path))
      # print "parent:", parent
      for i in t:
        print "t[%d] := "%i, t[i]

      # now, we add the constraint if its violations happens
      lhs = []
      rhs = []

      for i in range(1, numNodes):
        for j in range(1, numNodes):
          if i != j:
            if r(i) < r(j) + p(j) and r(j) < r(i) + p(i):
              lhs.append([
                 [y(i,j)],
                  [(v(i,j) + r(i) - r(j))]
                #,[0.0]
                ])
              rhs.append([(
                  v(j, i) * (v(i, j) - 1) +
                  r(j) * (v(i, j) - 1) +
                  r(i) * (v(j, i) + 1)
                  #-(v(i,j) + r(i) - r(j)) * t[i]
              )])


        numCuts = len(rhs)

        for i in range(numCuts):
          # calculate activity of cut (dot product)
          nameVar, coefVar = lhs[i][0][0], lhs[i][1][0]
          val     = coefVar * self.get_values(nameVar)
          # print val, rhs[i][0]
          #if val < rhs[i][0]:
            #self.add(constraint=lhs[i], sense="G", rhs=rhs[i][0])

# --------------------------------------------------------------------

def main():

  model = cp.Cplex()
  model.set_problem_name("ATSP-TW")
  model.set_problem_type(cp.Cplex.problem_type.MILP)

  # sys.stdout is the default output stream for log and results
  # so these lines may be omitted
  model.set_log_stream(sys.stdout)
  model.set_results_stream(sys.stdout)

  model.parameters.mip.strategy.search.set(
    model.parameters.mip.strategy.search.values.traditional
  )


  # vamos a minimizar
  model.objective.set_name("tourCost")
  model.objective.set_sense(model.objective.sense.minimize)

  # descripcion del modelo
  model.variables.add(
      obj   = my_obj                 # objective function cx^T
    , lb    = [0.0] * numVars        # lb <= xij
    , ub    = [1.0] * numVars        #       xij <= ub
    , names = my_colnames            #    xij
    , types = ["B"] * numVars
    )
  if reduceTimeWindows:
    tightening()

  # Horizontal constraints over Xij
  for i in range(numNodes):
    ind     = [x(i,j) for j in range(numNodes)]
    val     = [1] * numNodes
    row     = [[ind, val]]
    model.linear_constraints.add(
        lin_expr  = row
      , senses    = "E"
      , rhs       = [1]
      , names     = ["dminus" + str(i)]
    )

  # Vertical constraints over Xij
  for j in range(numNodes):
    ind   = [x(i,j)for i in range(numNodes)]
    val   = [1] * numNodes
    row   = [[ind, val]]
    model.linear_constraints.add(
        lin_expr  = row
      , senses    = "E"
      , rhs       = [1]
      , names     = ["dplus" + str(j)]
    )

  # Xii = 0 forall i in V
  for i in range(numNodes):
    model.linear_constraints.add(
        lin_expr  = [[[x(i,i)], [1]]]
      , senses    = "E"
      , rhs       = [0]
      , names     = ["trivial" + str(i)]
    )

  # adding the yij constraints
  model.variables.add(
      obj   = [0] * numVars
    , lb    = [0] * numVars
    , ub    = [cp.infinity] * numVars
    , names = names_yij
    , types = ['I'] * numVars
    )

  for i in range(1, numNodes):
    for j in range(numNodes):
      if i != j:
        model.linear_constraints.add(
            lin_expr  = [
            [
              [x(i,j), y(i,j)]
            , [r(i), -1]
            ]
            ]
          , senses    = "L"
          , rhs       = [0]
          , names     = ["r%sm%sConstr" % (str(i), str(j))]
        )

  for i in range(1, numNodes):
    for j in range(numNodes):
      if i != j:
        di = dd(i)
        model.linear_constraints.add(
            lin_expr  = [
            [
              [x(i,j), y(i,j)]
            , [-1*di, 1.0]
            ]
            ]
          , senses    = "L"
          , rhs       = [0]
          , names     = ["d%sm%sConstr" % (str(i), str(j))]
        )

  # contraints item 3
  n = numNodes
  for j in range(1, n):
    lhs = [y(i,j) for i in range(1,n) if i != j]\
        + [x(i,j) for i in range(0, n) if i != j]\
        + [y(j,k) for k in range(0, n) if k != j]
    val = [1 for i in range(1,n) if i != j]\
        + [v(i,j) for i in range(n) if i != j ]\
        + [-1 for k in range(n) if k != j]
    assert len(lhs) == len(val)

    model.linear_constraints.add(
          lin_expr  = [ [lhs, val] ]
        , senses    = "L" # change to "E" if you want equality
        , rhs       = [0]
        , names     = ["item3m" + str(j)]
      )

  # construction of precedence relation R
  R = {}
  for i in range(1,n):
    for j in range(1,n):
      if r(j) + v(j,i) > dd(i):
        R[(i,j)] = True
      else:
        R[(i,j)] = False

  for k in range(1,n):
    for i in range(1,n):
      for j in range(1,n):
        if R[(i,k)] and R[(k,j)]:
          R[(i,j)] = True

  for i in range(1,n):
    for j in range(1,n):
      if R[(i,j)]:
        model.linear_constraints.add(
          lin_expr  = [ [[x(j,i)], [1.0]] ]
        , senses    = "E" # change to "E" if you want equality
        , rhs       = [0]
        , names     = ["elimArcRrelation"+str(j)+"m"+str(i)]
        )



  model.write("atsptw-{0}.lp".format(databasename)) # save the model
  model.write("atsptw.lp") # save the model
  if twojobcuts:
    model.register_callback(LazyTwoJob) # register the cut lazy

  try:
    model.solve()
    # solution.get_status() returns an integer code
    print "Solution status = ", model.solution.get_status(), ":"
    # the following line prints the corresponding string
    print model.solution.status[model.solution.get_status()]

    # print non-zero solution values
    path  = [0]
    times = []

    u = 0
    while len(path) < numNodes:
      for i in range(numNodes):
        if i != u:
          varName = x(u,i)
          varYName = y(u,i)
          if abs(model.solution.get_values(varName)) > 1.e-10:
            path.append(i)
            times.append(model.solution.get_values(varYName))
            u = i

    ans = zip(path, times)
    ans.sort(key=lambda x: x[1])
    # assert len(path) == numNodes
    print "------------------------------------------------"
    print "Tour ATSP-TW"
    print "Tour Cost: ", model.solution.get_objective_value()
    for nodo, time in ans:
      print "node = {0}\t time = {1} \ttw={2}".format(nodo, time, [r(nodo),dd(nodo)])

    print
    print "path:", ' > '.join(map(str,path))
    print "data sources:"
    print "file costs (", filecosts, ")"
    print "file preprocessing (", filetime, ")"
    print "file time windows (", filewindows, ")"

  except CplexError as cplexerror:
    print cplexerror


if __name__ == "__main__":
  main()
{"mode":"full","isActive":false}