COP easy

Problem Warehouse Location

In the Warehouse Location problem (WLP), a company considers opening warehouses at some candidate locations in order to supply its existing stores. Each possible warehouse has the same maintenance cost, and a capacity designating the maximum number of stores that it can supply. Each store must be supplied by exactly one open warehouse. The supply cost to a store depends on the warehouse. The objective is to determine which warehouses to open, and which of these warehouses should supply the various stores, such that the sum of the maintenance and supply costs is minimized. See CSPLib (Problem 034) for more information.

Palumbo Fruit Company Warehouse. Image from commons.wikimedia.org

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.

fixed_cost = 30                 # cost of opening a warehouse
capacities = [1, 4, 2, 1, 3]    # capacities of the warehouses
costs = [                       # costs of supplying stores
   [100, 24, 11, 25, 30], 
    [28, 27, 82, 83, 74],
    [74, 97, 71, 96, 70], 
    [ 2, 55, 73, 69, 61],
    [46, 96, 59, 83,  4], 
    [42, 22, 29, 67, 59],
    [ 1,  5, 73, 59, 56], 
    [10, 73, 13, 43, 96],
    [93, 35, 63, 85, 46], 
    [47, 65, 55, 71, 95]
]
nWarehouses, nStores = len(capacities), len(costs)

Here, we have 5 warehouses and 10 stores. Everytime a wahehouse is open, it costs 30, the capacity of the second warehouse (at index 1) is 4 (it can supply upt to 4 stores), and the cost of supplying the second store with the third warehouse is 82.

We start our COP model with an array $w$ of 10 variables (one per store), each variable having ${0,1,\dots,4}$ as domain.

# w[i] is the warehouse supplying the ith store
w = VarArray(size=nStores, dom=range(nWarehouses))

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

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

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

For dealing with warehouse capacities, we use the constraint Count by calling the function Count(): the number of variables in a given list (here, $w$) that take the value specified by the named parameter value must be less than a constant.

satisfy(
    # capacities of warehouses must not be exceeded
    Count(w, value=j) <= capacities[j] for j in range(nWarehouses)
);

We can display the internal representation of the posted constraints; this way, although a little bit technical, we can see that the right capacities are set (note that ‘le’ stands for ‘less than or equal to’).

print(posted())

count(list:[w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8], w[9]], values:[0], condition:(le,1))
count(list:[w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8], w[9]], values:[1], condition:(le,4))
count(list:[w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8], w[9]], values:[2], condition:(le,2))
count(list:[w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8], w[9]], values:[3], condition:(le,1))
count(list:[w[0], w[1], w[2], w[3], w[4], w[5], w[6], w[7], w[8], w[9]], values:[4], condition:(le,3))

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(w))
    for i in range(nStores):
        print(f"Warehouse supplying the store {i} is {value(w[i])} with cost {costs[i][value(w[i])]}")

[0, 1, 1, 1, 1, 2, 2, 3, 4, 4]
Warehouse supplying the store 0 is 0 with cost 100
Warehouse supplying the store 1 is 1 with cost 27
Warehouse supplying the store 2 is 1 with cost 97
Warehouse supplying the store 3 is 1 with cost 55
Warehouse supplying the store 4 is 1 with cost 96
Warehouse supplying the store 5 is 2 with cost 29
Warehouse supplying the store 6 is 2 with cost 73
Warehouse supplying the store 7 is 3 with cost 43
Warehouse supplying the store 8 is 4 with cost 46
Warehouse supplying the store 9 is 4 with cost 95

We can control that the capacities of the warehouse are respected. We also note that all warehouses are open.

Concerning optimization, in a first step, let us suppose that we just want to minimize the number of warehouses to be open. This is exactly the task of minimizing the number of distinct values taken by variables of $w$. This can be written:

minimize(
    # minimizing the number of open warehouses
    NValues(w)
);

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

print(objective())

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

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(w))
    print("Number of open warehouses: ", bound())
    print("Cost of opening warehouses: " , bound()*fixed_cost)

[0, 1, 1, 1, 1, 2, 2, 4, 4, 4]
Number of open warehouses:  4
Cost of opening warehouses:  120

In a second step, let us suppose that we just want to minimize the cost of supplying the stores. For the ith store, this cost is given by $costs[i][w[i]]$ when the variable $w[i]$ is assigned (technically, here, a constraint ElementMatrix is involved). This can be written:

minimize(
    # minimizing the cost of supplying stores
    Sum(costs[i][w[i]] for i in range(nStores))
);

Unfortunately, this is not valid (an error occurs, if we run this piece of code). This is because it is not possible to use a special object (here, an object Var) as an index of a list of integers in Python. We need to transform the list of integers into a more specific list with the function cp_array(). So we write:

costs = cp_array(costs)

And we post the new objective:

minimize(
    # minimizing the cost of supplying stores
    Sum(costs[i][w[i]] for i in range(nStores))
);

Currently, whenever we call the function minimize() or maximize(), we overwrite the previous objective if any was present.

We can run again the solver, with this new optimization task.

if solve() is OPTIMUM:
    print(values(w))
    for i in range(nStores):
        print(f"Cost of supplying the store {i} is {costs[i][value(w[i])]}")
    print("Total supplying cost: ", bound())

[3, 1, 4, 0, 4, 1, 1, 2, 1, 2]
Cost of supplying the store 0 is 25
Cost of supplying the store 1 is 27
Cost of supplying the store 2 is 70
Cost of supplying the store 3 is 2
Cost of supplying the store 4 is 4
Cost of supplying the store 5 is 22
Cost of supplying the store 6 is 5
Cost of supplying the store 7 is 13
Cost of supplying the store 8 is 35
Cost of supplying the store 9 is 55
Total supplying cost:  258

Finally, we can combine both parts of the objective as follows:

minimize(
    # minimizing the overall cost
    Sum(costs[i][w[i]] for i in range(nStores)) + NValues(w) * fixed_cost
);

We run the solver, with this combined optimization task.

if solve() is OPTIMUM:
    print(values(w))
    for i in range(nStores):
        print(f"Cost of supplying the store {i} is {costs[i][value(w[i])]}")
    print("Overall cost: ", bound())

[4, 1, 4, 0, 4, 1, 1, 2, 1, 2]
Cost of supplying the store 0 is 30
Cost of supplying the store 1 is 27
Cost of supplying the store 2 is 70
Cost of supplying the store 3 is 2
Cost of supplying the store 4 is 4
Cost of supplying the store 5 is 22
Cost of supplying the store 6 is 5
Cost of supplying the store 7 is 13
Cost of supplying the store 8 is 35
Cost of supplying the store 9 is 55
Overall cost:  383

Finally, we give below the model in one piece. Here the data is expected to be given by the user (in a command line). Note that when the data is loaded from a file (which is the usual case), all lists of integers have the specific type of list returned by cp_array(), and so, it is very rare to need to call this function explicitly.

from pycsp3 import *

fixed_cost, capacities, costs = data  # fixed_cost is the fixed cost when opening a warehouse
nWarehouses, nStores = len(capacities), len(costs)

# w[i] is the warehouse supplying the ith store
w = VarArray(size=nStores, dom=range(nWarehouses))

satisfy(
    # capacities of warehouses must not be exceeded
    Count(w, value=j) <= capacities[j] for j in range(nWarehouses)
)

minimize(
    # minimizing the overall cost
    Sum(costs[i][w[i]] for i in range(nStores)) + NValues(w) * fixed_cost
)