The first two READ DATA statements populate the DEPTS and CITIES index sets:

proc optmodel;
   set <str> DEPTS;
   read data dept_data into DEPTS=[dept];

   set <str> CITIES;
   read data city_data into CITIES=[city];

The following READ DATA statement reads dense two-dimensional data, as in previous examples, with an initial value of 0 for benefit to account for the possibility of staying in London:

   num benefit {DEPTS, CITIES} init 0;
   read data benefit_data into [city] {dept in DEPTS} 
      <benefit[dept,city]=col(dept)>;
   print benefit;

Note that this READ DATA statement does not repopulate CITIES using the following:

   CITIES=[city]

Doing so would have removed London from the CITIES index set, because London does not appear in the benefit_data data set. Figure 11.1 shows the resulting values of benefit.

The following NUM and assignment statements read an upper triangular matrix and use these values to populate the lower triangular part. The INIT option initializes the parameter comm to a missing value, and the body of the FOR loop replaces each missing value with the value obtained by reflection across the main diagonal.

   num comm {DEPTS, DEPTS} init .;
   read data comm_data into [i j] comm;
   for {i in DEPTS, j in DEPTS} do;
      if i = j then comm[i,j] = 0;
      else if comm[i,j] = . then comm[i,j] = comm[j,i];
   end;
   print comm;

Figure 11.2 shows the resulting values of comm.

Similar statements are used to populate the cost parameter, but instead the main diagonal is read from the cost_data data set:

   num cost {CITIES, CITIES} init .;
   read data cost_data into [i j] cost;
   for {i in CITIES, j in CITIES: cost[i,j] = .}
      cost[i,j] = cost[j,i];
   print cost;

Figure 11.3 shows the resulting values of cost.

The following declarations are straightforward:

   var Assign {DEPTS, CITIES} binary;

   set IJKL = {i in DEPTS, j in CITIES, k in DEPTS, l in CITIES: i < k};
   var Product {IJKL} binary;

   impvar TotalBenefit 
      = sum {i in DEPTS, j in CITIES} benefit[i,j] * Assign[i,j];
   impvar TotalCost 
      = sum {<i,j,k,l> in IJKL} comm[i,k] * cost[j,l] * Product[i,j,k,l];
   max NetBenefit = TotalBenefit - TotalCost;

   con Assign_dept {dept in DEPTS}:
      sum {city in CITIES} Assign[dept,city] = 1;

   con Cardinality {city in CITIES}:
      sum {dept in DEPTS} Assign[dept,city] <= &max_num_depts;

The following CON statement enforces the rule that and together imply :

   con Product_def {<i,j,k,l> in IJKL}:
      Assign[i,j] + Assign[k,l] - 1 <= Product[i,j,k,l];

The following two CON statements enforce the converse rule that implies both and :

   con Product_def2 {<i,j,k,l> in IJKL}:
      Product[i,j,k,l] <= Assign[i,j];
   con Product_def3 {<i,j,k,l> in IJKL}:
      Product[i,j,k,l] <= Assign[k,l];

Because of the maximization objective in this problem, these two constraints could be omitted and would still be satisfied by an optimal solution. But including them can sometimes reduce the number of simplex iterations and branch-and-bound nodes.

   solve;

   print TotalBenefit TotalCost;
   print Assign;

Figure 11.4 shows the output from the mixed integer linear programming solver.

An alternative “compact linearization” formulation involves fewer constraints (Liberti 2007). The following DROP statement removes three families of constraints from the original formulation:

   drop Product_def Product_def2 Product_def3;

The following two CON statements declare two new constraints that enforce the desired relationship between the Product and Assign variables:

   con Product_def4 {i in DEPTS, k in DEPTS, l in CITIES: i < k}:
      sum {<(i),j,(k),(l)> in IJKL} Product[i,j,k,l] = Assign[k,l];
   con Product_def5 {k in DEPTS, i in DEPTS, j in CITIES: i < k}:
      sum {<(i),(j),(k),l> in IJKL} Product[i,j,k,l] = Assign[i,j];

   solve;

   print TotalBenefit TotalCost;
   print Assign;
quit;

Figure 11.5 shows the output from the mixed integer linear programming solver for the compact linearization formulation.

As expected, this solution agrees with the solution reported in Figure 11.4 for the original formulation.

In both formulations, the Product variable can be relaxed to be nonnegative instead of binary. The integrality of Assign, together with the various Product_def* constraints, automatically implies integrality of Product. For real-world problems, you should try both ways to determine which alternative performs better in specific cases.

Similarly, the compact formulation is weaker but contains fewer constraints than the original formulation. The net effect on total solve time is difficult to predict. For real-world problems, you should try both formulations.