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  medium 

Problem Vellino

From “Constraint Programming in OPL” by L. Michel, L. Perron, and J.-C. R'egin, CP’99: this problem involves putting components of different materials (glass, plastic, steel, wood, copper) into bins of various types (identified by red, blue, green colors), subject to capacity (each bin type has a maximum capacity) and compatibility constraints. Every component must be placed into a bin and the total number of used bins must be minimized. The compatibility constraints are:

  • red bins cannot contain plastic or steel
  • blue bins cannot contain wood or plastic
  • green bins cannot contain steel or glass
  • red bins contain at most one wooden component
  • green bins contain at most two wooden components
  • wood requires plastic
  • glass excludes copper
  • copper excludes plastic

See also CSPLib-Problem 116

Vellino’s Problem. image from freesvg.org Vellino

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.

capacities = [3, 1, 4]     # for red, blue and green bins
demands = [1, 2, 1, 3, 2]  # for glass, plastic, steel, wood, and copper materials

As indicated above, capacities are orderly given for red, blue and green bins, and demands (numbers of components) are orderly given for glass, plastic, steel, wood, and copper materials.

We define a few constants:

# 0 is a special color, ’Unusable’, to be used for any empty bin
Unusable, Red, Blue, Green = BIN_COLORS = 0, 1, 2, 3
Glass, Plastic, Steel, Wood, Copper = MATERIALS = 0, 1, 2, 3, 4
nColors, nMaterials = len(BIN_COLORS), len(MATERIALS)

We also insert a special capacity (0) for the special unusable bins

capacities.insert(0,0) # unusable bins have capacity 0
maxCapacity, nBins = max(capacities), sum(demands)

We control some of our data:

print("Capacities: ", capacities)
print("Max capacity: ", maxCapacity)
print("Number of bins: ", nBins)

Capacities:  [0, 3, 1, 4]
Max capacity:  4
Number of bins:  9

We start our COP model with an array $c$ of 9 variables (one per bin), each variable having ${0,1,2,3}$ as domain.

# c[i] is the color of the ith bin
c = VarArray(size=nBins, dom=range(nColors))

We can display (the structure of) the array as well as the domain of the first variable.

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

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

Then, we introduce a two-dimensional array $p$ of variables in order to record (after computing) the quantity of each material in each bin.

# p[i][j] is the number of components of the jth material put in the ith bin
p = VarArray(size=[nBins,nMaterials], dom=lambda i, j: range(min(maxCapacity, demands[j]) + 1))

We can display (the structure of) the array as well as the domain of a few variables.

print("Array of variable p: ", p)
for j in range(nMaterials):
    print(f"Domain of any variable related to material {j}: {p[0][j].dom}")

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

For dealing with unusable bins, we post the following constraints that are constraints Intension involving constraints Sum:

satisfy (
    # every bin with a real colour must contain something, and vice versa
    (c[i] == Unusable) == (Sum(p[i]) == 0) for i in range(nBins)
);

We can display the internal representation of the posted constraints; this way, although a little bit technical, we can check that the constraints are correctly posted (note that ‘eq’ stands for ‘equal to’). Because of the “complexity” of the constraints, mixing Intension and Sum, some auxiliary variables called $aux_gb$ are automatically introduced to deal with the involved constraining parts related to Sum.

print(posted())

intension(function:iff(eq(c[0],0),eq(aux_gb[0],0)))
intension(function:iff(eq(c[1],0),eq(aux_gb[1],0)))
intension(function:iff(eq(c[2],0),eq(aux_gb[2],0)))
intension(function:iff(eq(c[3],0),eq(aux_gb[3],0)))
intension(function:iff(eq(c[4],0),eq(aux_gb[4],0)))
intension(function:iff(eq(c[5],0),eq(aux_gb[5],0)))
intension(function:iff(eq(c[6],0),eq(aux_gb[6],0)))
intension(function:iff(eq(c[7],0),eq(aux_gb[7],0)))
intension(function:iff(eq(c[8],0),eq(aux_gb[8],0)))
sum(list:[p[0][0], p[0][1], p[0][2], p[0][3], p[0][4]], condition:(eq,aux_gb[0]))
sum(list:[p[1][0], p[1][1], p[1][2], p[1][3], p[1][4]], condition:(eq,aux_gb[1]))
sum(list:[p[2][0], p[2][1], p[2][2], p[2][3], p[2][4]], condition:(eq,aux_gb[2]))
sum(list:[p[3][0], p[3][1], p[3][2], p[3][3], p[3][4]], condition:(eq,aux_gb[3]))
sum(list:[p[4][0], p[4][1], p[4][2], p[4][3], p[4][4]], condition:(eq,aux_gb[4]))
sum(list:[p[5][0], p[5][1], p[5][2], p[5][3], p[5][4]], condition:(eq,aux_gb[5]))
sum(list:[p[6][0], p[6][1], p[6][2], p[6][3], p[6][4]], condition:(eq,aux_gb[6]))
sum(list:[p[7][0], p[7][1], p[7][2], p[7][3], p[7][4]], condition:(eq,aux_gb[7]))
sum(list:[p[8][0], p[8][1], p[8][2], p[8][3], p[8][4]], condition:(eq,aux_gb[8]))

Then, we have to guarantee that all components are put in the bins, by introducing a group of constraints Sum.

satisfy(
    # all components of each material are spread across all bins
    Sum(p[:,j]) == demands[j] for j in range(nMaterials)
);

Note that we use the notation $p[:,j]$ to extract the jth column of the array $p$, as in NumPy.

By calling posted(-1) we can get the constraints that have been posted at the last call to satisfy().

print(posted(-1))

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

The next constraints concern the capacity of each bin (that must not be exceeded).

Because we are going to use constraints Element, we need to slightly transform one list as follows:

capacities = cp_array(capacities)

The good news is that when the data are loaded from a file (the usual case), data have systematically the right types, and so there is no need to call the function cp_array() explicitly.

We can now post a group of constraints Sum involving constraints *Element on the left side of the expressions.

satisfy(  
   # the capacity of each bin is not exceeded
   Sum(p[i]) <= capacities[c[i]] for i in range(nBins)
);

Because of the “complexity” of the constraints, mixing Sum and Element, some auxiliary variables called $aux_gb$ are automatically introduced to deal with the involved constraining parts related to Element.

print(posted(-1))

sum(list:[p[0][0], p[0][1], p[0][2], p[0][3], p[0][4]], condition:(le,aux_gb[9]))
sum(list:[p[1][0], p[1][1], p[1][2], p[1][3], p[1][4]], condition:(le,aux_gb[10]))
sum(list:[p[2][0], p[2][1], p[2][2], p[2][3], p[2][4]], condition:(le,aux_gb[11]))
sum(list:[p[3][0], p[3][1], p[3][2], p[3][3], p[3][4]], condition:(le,aux_gb[12]))
sum(list:[p[4][0], p[4][1], p[4][2], p[4][3], p[4][4]], condition:(le,aux_gb[13]))
sum(list:[p[5][0], p[5][1], p[5][2], p[5][3], p[5][4]], condition:(le,aux_gb[14]))
sum(list:[p[6][0], p[6][1], p[6][2], p[6][3], p[6][4]], condition:(le,aux_gb[15]))
sum(list:[p[7][0], p[7][1], p[7][2], p[7][3], p[7][4]], condition:(le,aux_gb[16]))
sum(list:[p[8][0], p[8][1], p[8][2], p[8][3], p[8][4]], condition:(le,aux_gb[17]))
element(list:[0, 3, 1, 4], index:c[0], condition:(eq,aux_gb[9]))
element(list:[0, 3, 1, 4], index:c[1], condition:(eq,aux_gb[10]))
element(list:[0, 3, 1, 4], index:c[2], condition:(eq,aux_gb[11]))
element(list:[0, 3, 1, 4], index:c[3], condition:(eq,aux_gb[12]))
element(list:[0, 3, 1, 4], index:c[4], condition:(eq,aux_gb[13]))
element(list:[0, 3, 1, 4], index:c[5], condition:(eq,aux_gb[14]))
element(list:[0, 3, 1, 4], index:c[6], condition:(eq,aux_gb[15]))
element(list:[0, 3, 1, 4], index:c[7], condition:(eq,aux_gb[16]))
element(list:[0, 3, 1, 4], index:c[8], condition:(eq,aux_gb[17]))

By calling the function solve(), we can check that the problem is satisfiable (SAT). We can also display the values assigned to the variables in the found solution by calling the function value().

if solve() is SAT:
    print(values(c))
    print(values(p))

[0, 0, 0, 0, 0, 1, 1, 3, 0]
[
  [0, 0, 0, 0, 0]
  [0, 0, 0, 0, 0]
  [0, 0, 0, 0, 0]
  [0, 0, 0, 0, 0]
  [0, 0, 0, 0, 0]
  [0, 2, 0, 0, 0]
  [0, 0, 0, 1, 2]
  [1, 0, 1, 2, 0]
  [0, 0, 0, 0, 0]
]

We obtain a solution, but unfortunately, not all constraints have been taken into consideration.

To manage the fact that some bins cannot contain some types of materials, we post the three following groups of constraints Intension.

satisfy(
  # red bins cannot contain plastic or steel
  [If(c[i] == Red, Then=[p[i][Plastic] == 0, p[i][Steel] == 0]) for i in range(nBins)],
  
  # blue bins cannot contain wood or plastic
  [If(c[i] == Blue, Then=[p[i][Wood] == 0, p[i][Plastic] == 0]) for i in range(nBins)],
  
  # green bins cannot contain steel or glass
  [If(c[i] == Green, Then=[p[i][Steel] == 0, p[i][Glass] == 0]) for i in range(nBins)]
);

We can display them.

print(posted(-1))

intension(function:or(ne(c[0],1),eq(p[0][1],0)))
intension(function:or(ne(c[0],1),eq(p[0][2],0)))
intension(function:or(ne(c[1],1),eq(p[1][1],0)))
intension(function:or(ne(c[1],1),eq(p[1][2],0)))
intension(function:or(ne(c[2],1),eq(p[2][1],0)))
intension(function:or(ne(c[2],1),eq(p[2][2],0)))
intension(function:or(ne(c[3],1),eq(p[3][1],0)))
intension(function:or(ne(c[3],1),eq(p[3][2],0)))
intension(function:or(ne(c[4],1),eq(p[4][1],0)))
intension(function:or(ne(c[4],1),eq(p[4][2],0)))
intension(function:or(ne(c[5],1),eq(p[5][1],0)))
intension(function:or(ne(c[5],1),eq(p[5][2],0)))
intension(function:or(ne(c[6],1),eq(p[6][1],0)))
intension(function:or(ne(c[6],1),eq(p[6][2],0)))
intension(function:or(ne(c[7],1),eq(p[7][1],0)))
intension(function:or(ne(c[7],1),eq(p[7][2],0)))
intension(function:or(ne(c[8],1),eq(p[8][1],0)))
intension(function:or(ne(c[8],1),eq(p[8][2],0)))
intension(function:or(ne(c[0],2),eq(p[0][3],0)))
intension(function:or(ne(c[0],2),eq(p[0][1],0)))
intension(function:or(ne(c[1],2),eq(p[1][3],0)))
intension(function:or(ne(c[1],2),eq(p[1][1],0)))
intension(function:or(ne(c[2],2),eq(p[2][3],0)))
intension(function:or(ne(c[2],2),eq(p[2][1],0)))
intension(function:or(ne(c[3],2),eq(p[3][3],0)))
intension(function:or(ne(c[3],2),eq(p[3][1],0)))
intension(function:or(ne(c[4],2),eq(p[4][3],0)))
intension(function:or(ne(c[4],2),eq(p[4][1],0)))
intension(function:or(ne(c[5],2),eq(p[5][3],0)))
intension(function:or(ne(c[5],2),eq(p[5][1],0)))
intension(function:or(ne(c[6],2),eq(p[6][3],0)))
intension(function:or(ne(c[6],2),eq(p[6][1],0)))
intension(function:or(ne(c[7],2),eq(p[7][3],0)))
intension(function:or(ne(c[7],2),eq(p[7][1],0)))
intension(function:or(ne(c[8],2),eq(p[8][3],0)))
intension(function:or(ne(c[8],2),eq(p[8][1],0)))
intension(function:or(ne(c[0],3),eq(p[0][2],0)))
intension(function:or(ne(c[0],3),eq(p[0][0],0)))
intension(function:or(ne(c[1],3),eq(p[1][2],0)))
intension(function:or(ne(c[1],3),eq(p[1][0],0)))
intension(function:or(ne(c[2],3),eq(p[2][2],0)))
intension(function:or(ne(c[2],3),eq(p[2][0],0)))
intension(function:or(ne(c[3],3),eq(p[3][2],0)))
intension(function:or(ne(c[3],3),eq(p[3][0],0)))
intension(function:or(ne(c[4],3),eq(p[4][2],0)))
intension(function:or(ne(c[4],3),eq(p[4][0],0)))
intension(function:or(ne(c[5],3),eq(p[5][2],0)))
intension(function:or(ne(c[5],3),eq(p[5][0],0)))
intension(function:or(ne(c[6],3),eq(p[6][2],0)))
intension(function:or(ne(c[6],3),eq(p[6][0],0)))
intension(function:or(ne(c[7],3),eq(p[7][2],0)))
intension(function:or(ne(c[7],3),eq(p[7][0],0)))
intension(function:or(ne(c[8],3),eq(p[8][2],0)))
intension(function:or(ne(c[8],3),eq(p[8][0],0)))

Some bins must necessarily contain some types of materials:

satisfy( 
  # red bins contain at most one wooden component
  [If(c[i] == Red, Then=p[i][Wood] <= 1) for i in range(nBins)],
  
  # green bins contain at most two wooden components
  [If(c[i] == Green, Then=p[i][Wood] <= 2) for i in range(nBins)]
);

We can display them.

print(posted(-1))

intension(function:imp(eq(c[0],1),le(p[0][3],1)))
intension(function:imp(eq(c[1],1),le(p[1][3],1)))
intension(function:imp(eq(c[2],1),le(p[2][3],1)))
intension(function:imp(eq(c[3],1),le(p[3][3],1)))
intension(function:imp(eq(c[4],1),le(p[4][3],1)))
intension(function:imp(eq(c[5],1),le(p[5][3],1)))
intension(function:imp(eq(c[6],1),le(p[6][3],1)))
intension(function:imp(eq(c[7],1),le(p[7][3],1)))
intension(function:imp(eq(c[8],1),le(p[8][3],1)))
intension(function:imp(eq(c[0],3),le(p[0][3],2)))
intension(function:imp(eq(c[1],3),le(p[1][3],2)))
intension(function:imp(eq(c[2],3),le(p[2][3],2)))
intension(function:imp(eq(c[3],3),le(p[3][3],2)))
intension(function:imp(eq(c[4],3),le(p[4][3],2)))
intension(function:imp(eq(c[5],3),le(p[5][3],2)))
intension(function:imp(eq(c[6],3),le(p[6][3],2)))
intension(function:imp(eq(c[7],3),le(p[7][3],2)))
intension(function:imp(eq(c[8],3),le(p[8][3],2)))

The last constraints are:

satisfy(
  # wood requires plastic
  [If(p[i][Wood] > 0, Then=p[i][Plastic] > 0) for i in range(nBins)],
  
  # glass excludes copper
  [If(p[i][Glass] > 0, Then=p[i][Copper] == 0) for i in range(nBins)],
  
  # copper excludes plastic
  [If(p[i][Copper] > 0, Then=p[i][Plastic] == 0) for i in range(nBins)]
);

print(posted(-1))

intension(function:imp(gt(p[0][3],0),gt(p[0][1],0)))
intension(function:imp(gt(p[1][3],0),gt(p[1][1],0)))
intension(function:imp(gt(p[2][3],0),gt(p[2][1],0)))
intension(function:imp(gt(p[3][3],0),gt(p[3][1],0)))
intension(function:imp(gt(p[4][3],0),gt(p[4][1],0)))
intension(function:imp(gt(p[5][3],0),gt(p[5][1],0)))
intension(function:imp(gt(p[6][3],0),gt(p[6][1],0)))
intension(function:imp(gt(p[7][3],0),gt(p[7][1],0)))
intension(function:imp(gt(p[8][3],0),gt(p[8][1],0)))
intension(function:imp(gt(p[0][0],0),eq(p[0][4],0)))
intension(function:imp(gt(p[1][0],0),eq(p[1][4],0)))
intension(function:imp(gt(p[2][0],0),eq(p[2][4],0)))
intension(function:imp(gt(p[3][0],0),eq(p[3][4],0)))
intension(function:imp(gt(p[4][0],0),eq(p[4][4],0)))
intension(function:imp(gt(p[5][0],0),eq(p[5][4],0)))
intension(function:imp(gt(p[6][0],0),eq(p[6][4],0)))
intension(function:imp(gt(p[7][0],0),eq(p[7][4],0)))
intension(function:imp(gt(p[8][0],0),eq(p[8][4],0)))
intension(function:imp(gt(p[0][4],0),eq(p[0][1],0)))
intension(function:imp(gt(p[1][4],0),eq(p[1][1],0)))
intension(function:imp(gt(p[2][4],0),eq(p[2][1],0)))
intension(function:imp(gt(p[3][4],0),eq(p[3][1],0)))
intension(function:imp(gt(p[4][4],0),eq(p[4][1],0)))
intension(function:imp(gt(p[5][4],0),eq(p[5][1],0)))
intension(function:imp(gt(p[6][4],0),eq(p[6][1],0)))
intension(function:imp(gt(p[7][4],0),eq(p[7][1],0)))
intension(function:imp(gt(p[8][4],0),eq(p[8][1],0)))

At this stage, we can run again the solver.

if solve() is SAT:
    print(values(c))
    print(values(p))
    print("Number of used bins: ", sum(x.value > 0 for x in c))

[0, 0, 2, 1, 1, 3, 3, 0, 0]
[
  [0, 0, 0, 0, 0]
  [0, 0, 0, 0, 0]
  [0, 0, 1, 0, 0]
  [1, 0, 0, 0, 0]
  [0, 0, 0, 0, 2]
  [0, 1, 0, 1, 0]
  [0, 1, 0, 2, 0]
  [0, 0, 0, 0, 0]
  [0, 0, 0, 0, 0]
]
Number of used bins:  5

One can see that the number of used bins is 5 (actually, this value may depend on the version of the solver). Anyway, we now need to find an optimal solution. We can write:

minimize(
  # minimizing the number of used bins
  Sum(c[i] != Unusable for i in range(nBins))
);

We can also display the internal representation of the posted objective.

print(objective())

minimize(sum(list:[ne(c[0],0), ne(c[1],0), ne(c[2],0), ne(c[3],0), ne(c[4],0), ne(c[5],0), ne(c[6],0), ne(c[7],0), ne(c[8],0)]))

We can run again 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(values(c))
    print(values(p))
    print("Number of used bins: ", bound())

[2, 3, 1, 1, 0, 0, 0, 0, 3]
[
  [0, 0, 1, 0, 0]
  [0, 1, 0, 2, 0]
  [0, 0, 0, 0, 2]
  [1, 0, 0, 0, 0]
  [0, 0, 0, 0, 0]
  [0, 0, 0, 0, 0]
  [0, 0, 0, 0, 0]
  [0, 0, 0, 0, 0]
  [0, 1, 0, 1, 0]
]
Number of used bins:  5

We have now proved that it was not possible to use less than 5 bins.

Finally, we give below the model in one piece. Note how we have used the control function Match to gather some constraints. Here the data is expected to be given by the user (in a command line).

from pycsp3 import *

# 0 is a special color, 'Unusable', to be used for any empty bin
Unusable, Red, Blue, Green = BIN_COLORS = 0, 1, 2, 3  
Glass, Plastic, Steel, Wood, Copper = MATERIALS = 0, 1, 2, 3, 4
nColors, nMaterials = len(BIN_COLORS), len(MATERIALS)

capacities, demands = data
capacities.insert(0, 0)  # unusable bins have capacity 0
maxCapacity, nBins = max(capacities), sum(demands)

# c[i] is the color of the ith bin
c = VarArray(size=nBins, dom=range(nColors))

# p[i][j] is the number of components of the jth material put in the ith bin
p = VarArray(size=[nBins, nMaterials],
             dom=lambda i, j: range(min(maxCapacity, demands[j]) + 1))

satisfy(
    # every bin with a real colour must contain something, and vice versa
    [(c[i] == Unusable) == (Sum(p[i]) == 0) for i in range(nBins)],

    # all components of each material are spread across all bins
    [Sum(p[:, j]) == demands[j] for j in range(nMaterials)],

    # the capacity of each bin is not exceeded
    [Sum(p[i]) <= capacities[c[i]] for i in range(nBins)],

    # handling compatibility of materials
    [
        Match(
            c[i],
            Cases={
                Red: [p[i][Plastic] == 0, p[i][Steel] == 0, p[i][Wood] <= 1],
                Blue: [p[i][Wood] == 0, p[i][Plastic] == 0],
                Green: [p[i][Steel] == 0, p[i][Glass] == 0, p[i][Wood] <= 2]
            }
        ) for i in range(nBins)
    ],

    # wood requires plastic
    [If(p[i][Wood] > 0, Then=p[i][Plastic] > 0) for i in range(nBins)],

    # glass excludes copper
    [If(p[i][Glass] > 0, Then=p[i][Copper] == 0) for i in range(nBins)],

    # copper excludes plastic
    [If(p[i][Copper] > 0, Then=p[i][Plastic] == 0) for i in range(nBins)],

    # tag(symmetry-breaking)
    [LexIncreasing(p[i], p[i + 1]) for i in range(nBins - 1)]
)
  
minimize(
  # minimizing the number of used bins
  Sum(c[i] != Unusable for i in range(nBins))
)