Mixed Integer Programming to solve a retail category price optimization problem

In this post we are going to implement the procedure given in "A fractional programming approach for retail category price optimization" by S. Subramanian and H. D. Sherali. We are going to use the SCIP optimization suite with Python bindings. I have constructed a small example with n = 2 products and m = 3 prices per product. Here is the python code:

  """
Solves the retail category price optimization in "A fractional programming approach for retail category
price optimization" from S. Subramanian and H. D. Sherali
Copyright (c) by Orges Leka (2016)
"""
from pyscipopt import Model, quicksum
import math

"""

"""
def retail(n,m,P, d, p0, l, u, phi0, phiMax, R0, elast, mu, lmbda, alpha, beta, K, a,b,c,yMin,yMax):
  model = Model("retail price")
  v, r, g, x, z = {}, {}, {}, {},{}
  y = model.addVar(vtype="CONTINUOUS",name="y")
  phi = model.addVar(vtype="CONTINUOUS",name="phi")
  w = model.addVar(vtype="CONTINUOUS",name="w")
  for i in range(1,n+1):
    for j in range(1,m+1):
      v[i,j] = math.exp(mu[i]+lmbda[i]*P[i,j])
      r[i,j] = v[i,j] * P[i,j]
      g[i,j] = v[i,j] * (P[i,j] - d[i])
      x[i,j] = model.addVar(vtype="CONTINUOUS", name="x(%s,%s)"%(i,j))
      z[i,j] = model.addVar(vtype="BINARY",name="z(%s,%s)"%(i,j))

  model.addCons(quicksum(v[i,j]*x[i,j] for (i,j) in x) == phi) # 3b
  model.addCons( w == math.log(phi0) + elast/(1.0*n) * quicksum( math.log(P[i,j]/p0[i])*z[i,j] for (i,j) in z ) ) # 3d
  model.addCons( alpha*phi0 <= phi ) #3e
  model.addCons( phi <= phiMax ) # 3e
  model.addCons( math.log(alpha*phi0) <= w ) # 3e
  model.addCons( w <= math.log(phiMax) )  # 3e
  model.addCons( quicksum( r[i,j]*x[i,j] for (i,j) in x ) >= beta * R0 ) #3f
  for k in K.keys():
    model.addCons( a[k] * quicksum( P[k,l]*z[k,l] for l in range(1,m+1) ) <= b[k] * quicksum( P[k,l] * z[k,l] for l in range(1,m+1) ) + c[k] ) # 3g
  for k in K.keys():
    model.addCons( a[k] * quicksum( P[k,l]*x[k,l] for l in range(1,m+1) ) <= b[k] * quicksum( P[k,l] * x[k,l] for l in range(1,m+1) ) + c[k]*y ) # 3h
  for i in range(1,n+1):
    model.addCons( quicksum( z[i,j] for j in range(1,m+1)) == 1) # 3i
  for i in range(1,n+1):
    model.addCons( quicksum( x[i,j] for j in range(1,m+1)) == y) # 3j
  for i in range(1,n+1):
    for j in range(1,m+1):
      model.addCons( x[i,j] <= yMax * z[i,j] ) # 3k
      model.addCons( x[i,j] >= yMin * z[i,j] ) # 3l
  model.addCons( yMin <= y ) # 3m
  model.addCons( y <= yMax ) # 3m
  model.addCons( phi >= alpha*phi0*(1+w-math.log(alpha*phi0)) ) # 4b. T=2
  model.addCons( phi >= phiMax*(1+w-math.log(phiMax)) ) # 4b .  T = 2
  model.addCons( phi <= (phiMax-alpha*phi0)/(math.log(phiMax)-math.log(alpha*phi0))*(w-math.log(alpha*phi0))+alpha*phi0 )# 4a, S= 1
  model.setObjective(quicksum(g[i,j]*x[i,j] for (i,j) in x), "maximize") # 3a
  print g
  return model

def generateExample():
  n = 2
  m = 3
  P = { (1,1):2.49, (1,2):2.59, (1,3) : 2.69,
        (2,1):8.99, (2,2):9.09, (2,3) : 9.19 }
  d = { 1 : 1.49, 2 : 7.49 }
  p0 ={ 1 : 2.50, 2 : 9.05 }
  l = { 1 : 2.00, 2 : 8.49 }
  u = { 1 : 3.00, 2 : 9.99 }
  phi0 = 500
  phiMax = 700
  R0 = 50
  elast = -2.0
  mu = { 1 : 0.1, 2: 0.2}
  lmbda = { 1: 0.02, 2: 0.01 }
  alpha = 0.8
  beta = 0.9
  K = { 1 : 1}
  a = { 1 : 0.5}
  b = { 1 : 1.5}
  c = { 1 : 0.5}
  yMin = 0.0
  yMax = 1000.0
  return n,m,P, d, p0, l, u, phi0, phiMax, R0, elast, mu, lmbda, alpha, beta, K, a,b,c, yMin, yMax


if __name__ == "__main__":

    n,m,P, d, p0, l, u, phi0, phiMax, R0, elast, mu, lmbda, alpha, beta, K, a,b,c,yMin,yMax = generateExample()
    model = retail( n,m,P, d, p0, l, u, phi0, phiMax, R0, elast, mu, lmbda, alpha, beta, K, a,b,c,yMin,yMax  )

    model.hideOutput() # silent mode
    model.optimize()
    cost = model.getObjVal()
    print()
    print("Optimized:")
    print("Optimal value:", cost)
    #model.printAttr("X")
    for v in model.getVars():
      print(v.name, "=", model.getVal(v))

Please refer to the article for the meaning of the variables.