Step By Step Modeling Of PuzzlOr Electrifying ProblemPuzzlOr problems are nice because they are simplified versions of real world problems of interest. Last December problem is a simplified version of an interesting logistics problem. A variety of method can be used to attack them, see for instance this interesting post by Isaac Slavitt where he tried both a brute force search and simulated annealing. Not surprisingly, I will try CPLEX on it. Here is the statement of the problem
A new city is being built which will include 20 distinct neighborhoods as shown by the house icons in the map. As part of the planning process, electricity needs to be connected to each of the neighborhoods.
Let's model this problem using CPLEX and OPL. First, as is customary with OPL, we need to define a data file that contains the problem data. Separating the data file from the model file is good practice as it enables to run the same model with different data sets more easily. Here, this may be an overkill as there is only one data set, but let's keep the good habit still. The data file content is pretty straightforward, we provide
/*** In a real application this data would be fetched from an application data source. Let us now build the OPL model. We start with reading the data, uing the = ...; OPL statement
int nrows = ...; Note that we create ranges constructed from input data for the sake of clarity of the rest of the model. Let us now provide some utilities to move from a 1D array to a 2D matrix representation. These utilities use the construction of arrays via comprehensions in OPL:
int pos2row [p in positions] = p div ncols; We can now start the bulk of the problem definition. We have to define the decision variables, the objective function, and the constraints of the problem. We need to make two sets of decisions to solve the problem:
We know that the second set of decisions is constrained: each neighborhood is served by the closest generator. We will deal with this, and other constraints later on, let us first define a decision variable for each of the decisions we need to make. We model the first set of decisions by one binary decision variable per possible location of the generators. The decision variable is set to 1 if there is a generator at that location dvar int genAt[positions] in 0..1; We model the second set of decisions by one binary decision variable for each pair of possible generator positions, and neighborhood. The decision variable is set to 1 if the neighborhood is served by a generator put at that position.
dvar int serv Wait a minute, we have not defined the set of neighborhoods in the model. This set should be constructed from the given array in the input data. setof(int) neighborhoods = { pos  pos in positions : given[pos] > 0 }; The objective function is the sum of the distances of neighborhoods to the closest generator. Let us first define a utility for computing distances, again using comprehensions to create an OPL array.
float distance[n in neighborhoods, pos in positions] = The objective function is then straighforward. minimize sum (pos in positions, n in neighborhoods) distance[n, pos] * servedBy[n, pos]; Let us now look at the constraints of the problem. We have a fixed number ngens of generators to position. It means that exactly ngens of the genAt decision variables must be set to 1. We can capture this by saying that their sum is equal to ngens. sum (pos in positions) genAt[pos] == ngens; Each neighborhood is served by exactly one generator. It means that for a given neighborhood, exactly one of the decision variables servedBy can be set to 1. The usual way to state this in mathematical optimization is to state that the sum of all these variables is equal to 1. forall (n in neighborhoods) sum (pos in positions) servedBy[n, pos] == 1; A neighborhood is served by a generator at a given position only if there is actually a generator at that position. This may look silly when we read it, but we need to explicitly state the constraint. This is one of the difficulties when modeling a problem for mathematical optimization. It is easy to forget constraints that are so obvious that we don't even think to state them! Here it is not that hard, we know about the constraint. but when you model a problem that you do not know that well, you must make sure that all relevant constraints are documented and modeled. if not, then you end up producing mathematical solutions that aren't applicable to the real world problem you wanted to solve. Back to our problem: a neighborhood is served by a generator at a given position only if there is actually a generator at that position That condition is stated by a simple inequality between decision variables. forall (pos in positions, n in neighborhoods) servedBy[n, pos] <= genAt[pos]; Last, we should state that each neighborhood is served by the closest generator. How can we do that? We will show at the end of this post how to state this constraint, but let's not do it here. Why? We won't state this constraint because we don't have to! Indeed, assume we have a feasible solution where a neighborhood is not served by its closest generator. This is possible, as we do not explicitly forbid it. Let us say, for the sake of clarity, that the neighborhood n is served by a generator at position pos1 when the closest generator is at position pos2. In this solution we have servedBy[n, pos1] = 1 servedBy[n, pos2] = 0 distance[n, pos1] > distance[n, pos2] Let us now modify that solution by swapping the value of servedBy[n, pos1] and servedBy[n, pos2]. Does this yield another feasible solution of the problem? Let's check the constraints. There are two constraints involving servedBy variables. The first set is sum (pos in positions) servedBy[n, pos] == 1; It stays true after the swap, given we still ave exactly one servedBy variable set to 1 per neighborhood after the swap. The second set of constraints is a little trickier to analyze forall (pos in positions, n in neighborhoods) servedBy[n, pos] <= genAt[pos]; Our swap impacts two constraint in this set of constraint servedBy[n, pos1] <= genAt[pos1]; servedBy[n, pos2] <= genAt[pos2]; By definition of our solution, we have a generator at position pos1 and one at position pos2. genAt[pos1] = 1 genAt[pos2] = 1 . Therefore, the following constraints are met whatever the values of ServedBy[n, pos1] and servedBy[n, pos2]. servedBy[n, pos1] <= genAt[pos1]; servedBy[n, pos2] <= genAt[pos2]; It means that our swap will not change the fact that the constraints are met. In short, our swap yields a second feasible solution to the problem! In the second feasible solution, the neighborhood n is served by a generator at position pos2 while in the original feasible solution it is served by a generator at position pos1 It means that both feasible solutions have the same objective value, except for the following; The term distance[n, pos1] * servedBy[n, pos1] in the first feasible solution is replaced by distance[n, pos2] * servedBy[n, pos2] in the second feasible solution. Given the latter is smaller than the former, the second feasible solution is better. We say that the first feasible solution is dominated by the second feasible solution. The magic comes from the fact that mathematical optimization outputs a feasible solution only if it is optimal, i.e. only if there are not better other feasible solutions . In our case, it means that the first feasible solution will not be output given it is dominated by another feasible solution. Therefore, we achieve our goal of only computing solutions where neighborhoods are served by the closest generator without explicitly stating it. Let us finish our model with a little script to display the result using the coordinate names from the input data
execute { The full model is appended at the end of tis post. Running it on my laptop solve sthe problem instantaneously, and the script displays
// solution (optimal) with objective 35.6274637069481 Graphically, this yields (picture is from Isaac Slavitt post): For the sake of curiosity, let us model explicitly the fact that a neighborhood must be served by the closest generator. We must forbid solutions where there is a neighborhood n, and two positions pos1 and pos2 such that servedBy[n, pos1] = 1 servedBy[n, pos2] = 0 genAt[pos1] = 1 genAtPos[2] = 1 distance[n, pos1] > distance[n, pos2] One way to state it is as follows.
forall (n in neighborhoods, pos1 in positions, pos2 in positions: distance[n,pos1] > distance[n, pos2]) It is easy to check that it would reject the case we want to reject while accepting the case where we swap the servedBy variables values as follows servedBy[n, pos1] = 0 servedBy[n, pos2] = 1 genAt[pos1] = 1 genAtPos[2] = 1 distance[n, pos1] > distance[n, pos2] By our previous discussion, such constraint is a dominance constraint. Stating it does not change the optimal solutions of the optimization problem. Running this model is much slower, about one second, instead of an instantaneous solve. The reason is that the resulting model has many more constraints than the original one. In that case, it is better to not explicitly state these additional constraints. One should not turn this into a general rule though. It may be that for some other problems, explicitly stating dominance constraints actually help solve the problem. Knowing when we can safely omit constraints because of solutions dominance is one of the key skills for modeling problems with mathematical optimization. Best rule of thumb is to try with and without them. This concludes our discussion of this nice PuzzlOr problem. The complete model is the following. It has to be used with the data file provided above.
/***
/* This problem is a toy problem, but it is linked to real problems. For instance, a major problem retailers have is to define their distribution network. They have a set of stores, and they need to decide where to put warehouses to serve these stores. The problem is similar to the above one, when you replace neighborhoods by stores, and generators by warehouses. There are other constraints of course, for instance a given warehouse has a limited capacity. It means that the number of stores that can be served by a warehouse is limited. We could easily add that constraint to the above model, by stating that the number of the servedBy[n, pos] variables for a given pos is limited. The usual way to state this is to use the sum of these variables, as follows. forall (pos in positions) sum (n in neighborhoods) servedBy[n, pos] <= capacity[pos]; Of course, the capacity of each warehouse must be provided as input.
