Link Search Menu Expand Document
PyCSP3
Models & Data GitHub XCSP3 About
Documentation
XCSP web site is provided by CRIL
This site uses Just the Docs, a documentation theme for Jekyll.
download notebook
COP  difficult 

Problem Open Stacks

From Steven Prestwich: ``A manufacturer has a number of orders from customers to satisfy. Each order is for a number of different products, and only one product can be made at a time. Once a customer’s order is started a stack is created for that customer. When all the products that a customer requires have been made the order is sent to the customer, so that the stack is closed. Because of limited space in the production area, the number of stacks that are simultaneously open should be minimized.’’

Note than when a product is made, the produced quantity is sufficient for all orders.

To build a COP (Constraint Optimization Problem) model, we need first to import the library PyCSP$^3$:

from pycsp3 import *

Then, we need some data. Here, we have a certain number of orders. Each row of the two-dimensional array below corresponds to a customer order indicating with 0 or 1 which products are needed: orders[i][j] is 1 if the jth product is needed for the ith order.

orders = [
    [1, 0, 1, 0, 0, 1, 0, 1, 0, 0], 
    [1, 0, 0, 0, 1, 1, 0, 0, 1, 0],
    [0, 0, 1, 0, 0, 0, 0, 0, 1, 0], 
    [0, 1, 0, 0, 0, 1, 0, 0, 0, 1], 
    [0, 0, 0, 1, 0, 1, 0, 1, 0, 0], 
    [1, 0, 0, 0, 0, 0, 1, 0, 0, 1], 
    [0, 1, 0, 0, 0, 1, 0, 1, 1, 0], 
    [1, 0, 1, 1, 1, 0, 0, 0, 0, 0], 
    [1, 1, 1, 1, 0, 0, 1, 1, 0, 1], 
    [1, 0, 0, 0, 0, 1, 0, 1, 0, 1]
]

We can define a few useful constants.

n = len(orders)     # n orders (customers)
m = len(orders[0])  # m possible products

We gently start our COP model with an array $p$ of $m$ variables, one per product.

# p[j] is the period (time) of the jth product
p = VarArray(size=m, dom=range(m))

We can display the structure of the array, as well as the domain of the first variable (note that all variables have the same domain).

print("Array of variable p: ", p)
print("Domain of any variable in p: ", p[0].dom)

Array of variable p:  [p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], p[9]]
Domain of any variable in p:  0..9

Concerning the constraints, we first post a constraint AllDifferent.

satisfy(
   # all products are scheduled at different times
   AllDifferent(p)
);

Then, we need to know at which time a stack must be open for a costumer. We then introduce an array $n$ of 10 variables.

# s[i] is the starting time of the ith stack
s = VarArray(size=n, dom=range(m))

To “compute” the starting time of each stack, we post a group of constraints Minimum.

satisfy(
   # computing starting times of stacks
   Minimum(p[j] for j in range(m) if orders[i][j]) == s[i] for i in range(n)
);

We can display the internal representation of the posted constraints; this way, although a little bit technical, we can see that the arguments to the constraints are correct (note that ‘eq’ stands for ‘equal to’).

print(posted())

allDifferent(list:[p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], p[9]])
minimum(list:[p[0], p[2], p[5], p[7]], condition:(eq,s[0]))
minimum(list:[p[0], p[4], p[5], p[8]], condition:(eq,s[1]))
minimum(list:[p[2], p[8]], condition:(eq,s[2]))
minimum(list:[p[1], p[5], p[9]], condition:(eq,s[3]))
minimum(list:[p[3], p[5], p[7]], condition:(eq,s[4]))
minimum(list:[p[0], p[6], p[9]], condition:(eq,s[5]))
minimum(list:[p[1], p[5], p[7], p[8]], condition:(eq,s[6]))
minimum(list:[p[0], p[2], p[3], p[4]], condition:(eq,s[7]))
minimum(list:[p[0], p[1], p[2], p[3], p[6], p[7], p[9]], condition:(eq,s[8]))
minimum(list:[p[0], p[5], p[7], p[9]], condition:(eq,s[9]))

For example, we have:

$min(p[0],p[2],p[5],p[7]) = s[0]$

which is coherent since the order of customer 0 requires products 0, 2, 5 and 7.

Interestingly, by calling the function solve(), we can check that the problem is satisfiable (SAT). We can also display the found solution. Here, we call the function values() that collects the values assigned to a specified list of variables.

if solve() is SAT:
    print("Production time of products: ", values(p))
    print("Starting time of stacks: ", values(s))

Production time of products:  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
Starting time of stacks:  [0, 0, 2, 1, 3, 0, 1, 0, 0, 0]

This looks coherent. For example, the starting time of a customer stack is 0 when it involves the product planned to be produced at time 0.

We also need to know at which time a stack can be closed for a costumer. We then introduce an array $e$ of $n$ variables.

# e[i] is the ending time of the ith stack
e = VarArray(size=n, dom=range(m))

To “compute” the ending time of each stack, we post a group of constraints Maximum.

satisfy(
   # computing ending times of stacks
   Maximum(p[j] for j in range(m) if orders[i][j]) == e[i] for i in range(n)
);

We run again the solver.

if solve() is SAT:
    print("Production time of products: ", values(p))
    print("Starting time of stacks: ", values(s))
    print("Ending time of stacks: ", values(e))

Production time of products:  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
Starting time of stacks:  [0, 0, 2, 1, 3, 0, 1, 0, 0, 0]
Ending time of stacks:  [7, 8, 8, 9, 7, 9, 8, 4, 9, 9]

At any time, we need to know if a stack is open (so as to be able to compute the number of open stacks at any time). We introduce a two-dimensional array $o$ of $n \times m$ 0/1 variables.

# o[i][t] is 1 iff the ith stack is open at time t
o = VarArray(size=[n, m], dom={0, 1})

To compute the value of these variables from some constraints, we shall use the following function:

def table(t):
    return ([(ANY, te, 0) for te in range(t)] + 
            [(ts, ANY, 0) for ts in range(t + 1, m)] +
            [(ts, te, 1) for ts in range(t + 1) for te in range(t, m)])

This function returns a set (list) of 3-tuples where, for each tuple, the first value is a possible starting time of a stack, the second value is a possible end time of a task, the third value is 0 if considering the starting and ending times (two first values), the stack is not open at time t, 1 otherwise. For example, if we print the table for t=5, we can see that :

  • when the ending time is less than or equal to 4, whatever is the starting time, the stack is not open at time t=5
  • when the starting time is greater than or equal to 6, whatever is the ending time, the stack is not open at time t=5
  • if the starting time is 0 end the ending time is 5, the stack is open at time t=5
  • etc.
print(table(5))

[(*, 0, 0), (*, 1, 0), (*, 2, 0), (*, 3, 0), (*, 4, 0), (6, *, 0), (7, *, 0), (8, *, 0), (9, *, 0), (0, 5, 1), (0, 6, 1), (0, 7, 1), (0, 8, 1), (0, 9, 1), (1, 5, 1), (1, 6, 1), (1, 7, 1), (1, 8, 1), (1, 9, 1), (2, 5, 1), (2, 6, 1), (2, 7, 1), (2, 8, 1), (2, 9, 1), (3, 5, 1), (3, 6, 1), (3, 7, 1), (3, 8, 1), (3, 9, 1), (4, 5, 1), (4, 6, 1), (4, 7, 1), (4, 8, 1), (4, 9, 1), (5, 5, 1), (5, 6, 1), (5, 7, 1), (5, 8, 1), (5, 9, 1)]

We can now post a group of constraints Extension to be able to compute the values of variables in $o$.

satisfy(
   # inferring when stacks are open
   (s[i], e[i], o[i][t]) in table(t) for i in range(n) for t in range(m)
);

We can display the 3 first constraints of this group.

print(posted(-1)[:4])

extension(list:[s[0], e[0], o[0][0]], supports:(0,0,1)(0,1,1)(0,2,1)(0,3,1)(0,4,1)(0,5,1)(0,6,1)(0,7,1)(0,8,1)(0,9,1)(1,*,0)(2,*,0)(3,*,0)(4,*,0)(5,*,0)(6,*,0)(7,*,0)(8,*,0)(9,*,0))
extension(list:[s[0], e[0], o[0][1]], supports:(0,1,1)(0,2,1)(0,3,1)(0,4,1)(0,5,1)(0,6,1)(0,7,1)(0,8,1)(0,9,1)(1,1,1)(1,2,1)(1,3,1)(1,4,1)(1,5,1)(1,6,1)(1,7,1)(1,8,1)(1,9,1)(2,*,0)(3,*,0)(4,*,0)(5,*,0)(6,*,0)(7,*,0)(8,*,0)(9,*,0)(*,0,0))
extension(list:[s[0], e[0], o[0][2]], supports:(0,2,1)(0,3,1)(0,4,1)(0,5,1)(0,6,1)(0,7,1)(0,8,1)(0,9,1)(1,2,1)(1,3,1)(1,4,1)(1,5,1)(1,6,1)(1,7,1)(1,8,1)(1,9,1)(2,2,1)(2,3,1)(2,4,1)(2,5,1)(2,6,1)(2,7,1)(2,8,1)(2,9,1)(3,*,0)(4,*,0)(5,*,0)(6,*,0)(7,*,0)(8,*,0)(9,*,0)(*,0,0)(*,1,0))
extension(list:[s[0], e[0], o[0][3]], supports:(0,3,1)(0,4,1)(0,5,1)(0,6,1)(0,7,1)(0,8,1)(0,9,1)(1,3,1)(1,4,1)(1,5,1)(1,6,1)(1,7,1)(1,8,1)(1,9,1)(2,3,1)(2,4,1)(2,5,1)(2,6,1)(2,7,1)(2,8,1)(2,9,1)(3,3,1)(3,4,1)(3,5,1)(3,6,1)(3,7,1)(3,8,1)(3,9,1)(4,*,0)(5,*,0)(6,*,0)(7,*,0)(8,*,0)(9,*,0)(*,0,0)(*,1,0)(*,2,0))

Interstingly, when we display the first solution obtained by the solver, we can directly see that it is a valid solution.

if solve() is SAT:
    print("Production time of products: ", values(p))
    print("Starting time of stacks: ", values(s))
    print("Ending time of stacks: ", values(e))
    print("Stacks: ", values(o))

Production time of products:  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
Starting time of stacks:  [0, 0, 2, 1, 3, 0, 1, 0, 0, 0]
Ending time of stacks:  [7, 8, 8, 9, 7, 9, 8, 4, 9, 9]
Stacks:  [
  [1, 1, 1, 1, 1, 1, 1, 1, 0, 0]
  [1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
  [0, 0, 1, 1, 1, 1, 1, 1, 1, 0]
  [0, 1, 1, 1, 1, 1, 1, 1, 1, 1]
  [0, 0, 0, 1, 1, 1, 1, 1, 0, 0]
  [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
  [0, 1, 1, 1, 1, 1, 1, 1, 1, 0]
  [1, 1, 1, 1, 1, 0, 0, 0, 0, 0]
  [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
  [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
]

The last thing to do is to post an objective function. Here, we need to minimize the maximum number of 1 in the columns of $o$. Note that we use the notation $o[:, t]$ to extract the column t of the array $o$, as in NumPy.

minimize(
   # minimizing the number of stacks that are simultaneously open
   Maximum(Sum(o[:, t]) for t in range(m))
);

If we display the internal representation of the objective, one can see the presence of new (auxiliary) variables, called $aux_gb$. They are introduced to deal with the complexity of the expression of the objective, as intern sums are involved.

print(objective())

minimize(maximum(list:[aux_gb[0], aux_gb[1], aux_gb[2], aux_gb[3], aux_gb[4], aux_gb[5], aux_gb[6], aux_gb[7], aux_gb[8], aux_gb[9]]))

We can display the 3 constraints that have indirectly (automatically) introduced.

print(posted()[-3:])

sum(list:[o[0][7], o[1][7], o[2][7], o[3][7], o[4][7], o[5][7], o[6][7], o[7][7], o[8][7], o[9][7]], condition:(eq,aux_gb[7]))
sum(list:[o[0][8], o[1][8], o[2][8], o[3][8], o[4][8], o[5][8], o[6][8], o[7][8], o[8][8], o[9][8]], condition:(eq,aux_gb[8]))
sum(list:[o[0][9], o[1][9], o[2][9], o[3][9], o[4][9], o[5][9], o[6][9], o[7][9], o[8][9], o[9][9]], condition:(eq,aux_gb[9]))

We can run the solver, with this optimization task. Note that we need to check that the status returned by the solver is now OPTIMUM.

if solve() is OPTIMUM:
    print("Production time of products: ", values(p))
    print("Stacks: ", values(o))
    print(f"maximum number of open stacks : {bound()}")

Production time of products:  [0, 1, 4, 9, 6, 7, 3, 8, 5, 2]
Stacks:  [
  [1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
  [1, 1, 1, 1, 1, 1, 1, 1, 0, 0]
  [0, 0, 0, 0, 1, 1, 0, 0, 0, 0]
  [0, 1, 1, 1, 1, 1, 1, 1, 0, 0]
  [0, 0, 0, 0, 0, 0, 0, 1, 1, 1]
  [1, 1, 1, 1, 0, 0, 0, 0, 0, 0]
  [0, 1, 1, 1, 1, 1, 1, 1, 1, 0]
  [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
  [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
  [1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
]
maximum number of open stacks : 8

Finally, we give below the model in one piece. Here the data is expected to be given by the user (in a command line).

from pyscp3 import *

orders = data
n, m = len(orders), len(orders[0])  # n orders (customers), m possible products

def table(t):
    return ({(ANY, te, 0) for te in range(t)} |
            {(ts, ANY, 0) for ts in range(t + 1, m)} |
            {(ts, te, 1) for ts in range(t + 1) for te in range(t, m)})
    
# p[j] is the period (time) of the jth product
p = VarArray(size=m, dom=range(m))

# s[i] is the starting time of the ith stack
s = VarArray(size=n, dom=range(m))

# e[i] is the ending time of the ith stack
e = VarArray(size=n, dom=range(m))

# o[i][t] is 1 iff the ith stack is open at time t
o = VarArray(size=[n, m], dom={0, 1})

satisfy(
   # all products are scheduled at different times
   AllDifferent(p),

   # computing starting times of stacks
   [Minimum(p[j] for j in range(m) if orders[i][j]) == s[i] for i in range(n)],

   # computing ending times of stacks
   [Maximum(p[j] for j in range(m) if orders[i][j]) == e[i] for i in range(n)],
   
   # inferring when stacks are open
   [(s[i], e[i], o[i][t]) in table(t) for i in range(n) for t in range(m)],
)

minimize(
   # minimizing the number of stacks that are simultaneously open
   Maximum(Sum(o[:, t]) for t in range(m))
)