PROC OPTMODEL Statements and Output

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.

Figure 11.1: benefit Parameter

The OPTMODEL Procedure

benefit
  Brighton Bristol London
A 10 10 0
B 20 15 0
C 15 10 0
D 15 20 0
E 15 5 0


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.

Figure 11.2: comm Parameter

comm
  A B C D E
A 0.0 0.0 1.0 1.5 0.0
B 0.0 0.0 1.4 1.2 0.0
C 1.0 1.4 0.0 0.0 2.0
D 1.5 1.2 0.0 0.0 0.7
E 0.0 0.0 2.0 0.7 0.0


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.

Figure 11.3: cost Parameter

cost
  Brighton Bristol London
Brighton 5 14 9
Bristol 14 5 13
London 9 13 10


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 $\Variable{Assign[i,j]} = 1$ and $\Variable{Assign[k,l]} = 1$ together imply $\Variable{Product[i,j,k,l]} = 1$:

   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 $\Variable{Product[i,j,k,l]} = 1$ implies both $\Variable{Assign[i,j]} = 1$ and $\Variable{Assign[k,l]} = 1$:

   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.

Figure 11.4: Output from Mixed Integer Linear Programming Solver

Problem Summary
Objective Sense Maximization
Objective Function NetBenefit
Objective Type Linear
   
Number of Variables 105
Bounded Above 0
Bounded Below 0
Bounded Below and Above 105
Free 0
Fixed 0
Binary 105
Integer 0
   
Number of Constraints 278
Linear LE (<=) 273
Linear EQ (=) 5
Linear GE (>=) 0
Linear Range 0
   
Constraint Coefficients 660

Performance Information
Execution Mode Single-Machine
Number of Threads 1

Solution Summary
Solver MILP
Algorithm Branch and Cut
Objective Function NetBenefit
Solution Status Optimal
Objective Value 14.90002665
   
Relative Gap 0
Absolute Gap 0
Primal Infeasibility 5E-7
Bound Infeasibility 5E-7
Integer Infeasibility 2E-6
   
Best Bound 14.90002665
Nodes 1
Iterations 212
Presolve Time 0.00
Solution Time 0.55

TotalBenefit TotalCost
80 65.1

Assign
  Brighton Bristol London
A 5.0E-07 1.0E+00 -5.0E-07
B 1.0E+00 5.0E-07 5.0E-07
C 1.0E+00 5.0E-07 0.0E+00
D 1.0E-06 1.0E+00 0.0E+00
E 1.0E+00 0.0E+00 0.0E+00


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.

Figure 11.5: Output from Mixed Integer Linear Programming Solver (Compact Linearization)

Problem Summary
Objective Sense Maximization
Objective Function NetBenefit
Objective Type Linear
   
Number of Variables 105
Bounded Above 0
Bounded Below 0
Bounded Below and Above 105
Free 0
Fixed 0
Binary 105
Integer 0
   
Number of Constraints 68
Linear LE (<=) 3
Linear EQ (=) 65
Linear GE (>=) 0
Linear Range 0
   
Constraint Coefficients 270

Performance Information
Execution Mode Single-Machine
Number of Threads 1

Solution Summary
Solver MILP
Algorithm Branch and Cut
Objective Function NetBenefit
Solution Status Optimal
Objective Value 14.9000088
   
Relative Gap 0
Absolute Gap 0
Primal Infeasibility 1.110223E-16
Bound Infeasibility 0
Integer Infeasibility 2E-6
   
Best Bound 14.9000088
Nodes 1
Iterations 95
Presolve Time 0.00
Solution Time 0.11

TotalBenefit TotalCost
80 65.1

Assign
  Brighton Bristol London
A 0.000001 0.999999 0.000000
B 1.000000 0.000000 0.000000
C 0.999999 0.000001 0.000000
D 0.000001 0.999999 0.000000
E 0.999999 0.000001 0.000000


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.