equiv: Location of income tax offices

The example description in the following sections is taken from Section 15.5 `Location of income tax offices' of the book `Applications of optimization with Xpress-MP'.

The income tax administration is planning to restructure the network of income tax offices in a region. The number of inhabitants of every city and the distances between each pair of cities are known (Table Distance matrix and population of cities). The income tax administration has determined that offices should be established in three cities to provide sufficient coverage. Where should these offices be located to minimize the average distance per inhabitant to the closest income tax office?

Model formulation

Let CITIES be the set of cities. For the formulation of the problem, two groups of decision variables are necessary: a variable build_c that is one if and only if a tax office is established in city c, and a variable depend_c that takes the number of the office on which city c depends. For the formulation of the constraints, we further introduce two sets of auxiliary variables: depdist_c, the distance from city c to the office indicated by depend_c, and numdep_c, the number of cities depending on an office location.

The following relations are required to link the build_c with the depend_c variables:
(1) numdep_c counts the number of occurrences of office location c among the variables depend_c.
(2) numdep_c ≥ 1 if and only if the office in c is built (as a consequence, if the office in c is not built, then we must have numdep_c = 0).

Since the number of offices built is limited by the given bound NUMLOC

it would actually be sufficient to formulate the second relation between the build_c and depend_c variables as the implication `If numdep_c ≥ 1 then the office in c must be built, and inversely, if the office in c is not built, then we must have numdep_c = 0'.

The objective function to be minimized is the total distance weighted by the number of inhabitants of the cities. We need to divide the resulting value by the total population of the region to obtain the average distance per inhabitant to the closest income tax office. The distance depdist_c from city c to the closest tax office location is obtained by a discrete function, namely the row c of the distance matrix DIST_cd indexed by the value of depend_c:

We now obtain the following CP model:

Implementation

To solve this problem, we define a branching strategy with two parts, one for the build_c variables and a second strategy for the depdist_c variables. The latter are enumerated using the split_domain branching scheme that divides the domain of the branching variable into several disjoint subsets (instead of assigning a value to the variable). We now pass an array of type cpbranching as the argument to procedure cp_set_branching. The different strategies will be applied in their order in this array. Since our enumeration strategy does not explicitly include all decision variables of the problem, Xpress Kalis will enumerate these using the default strategy if any unassigned variables remain after the application of our search strategy.

model "J-5 Tax office location (CP)"
 uses "kalis"

 forward procedure calculate_dist

 setparam("KALIS_DEFAULT_LB", 0)

 declarations
  NC = 12
  CITIES = 1..NC                        ! Set of cities

  DIST: array(CITIES,CITIES) of integer ! Distance matrix
  POP: array(CITIES) of integer         ! Population of cities
  LEN: dynamic array(CITIES,CITIES) of integer ! Road lengths
  NUMLOC: integer                       ! Desired number of tax offices
  D: array(CITIES) of integer           ! Auxiliary array used in constr. def.

  build: array(CITIES) of cpvar         ! 1 if office in city, 0 otherwise
  depend: array(CITIES) of cpvar        ! Office on which city depends
  depdist: array(CITIES) of cpvar       ! Distance to tax office
  numdep: array(CITIES) of cpvar        ! Number of depending cities per off.
  totDist: cpvar                        ! Objective function variable
  Strategy: array(1..2) of cpbranching  ! Branching strategy
 end-declarations

 initializations from 'Data/j5tax.dat'
  LEN POP NUMLOC
 end-initializations

! Calculate the distance matrix
 calculate_dist

 forall(c in CITIES) do
  build(c) <= 1
  1 <= depend(c); depend(c) <= NC
  min(d in CITIES) DIST(c,d) <= depdist(c)
  depdist(c) <= max(d in CITIES) DIST(c,d)
  numdep(c) <= NC
 end-do

! Distance from cities to tax offices
 forall(c in CITIES) do
  forall(d in CITIES) D(d):=DIST(c,d)
  element(D, depend(c)) = depdist(c)
 end-do

! Number of cities depending on every office
 forall(c in CITIES) occurrence(c, depend) = numdep(c)

! Relations between dependencies and offices built
 forall(c in CITIES) equiv( build(c) = 1, numdep(c) >= 1 )

! Limit total number of offices
 sum(c in CITIES) build(c) <= NUMLOC

! Branching strategy
 Strategy(1):= assign_and_forbid(KALIS_MAX_DEGREE, KALIS_MAX_TO_MIN, build)
 Strategy(2):= split_domain(KALIS_SMALLEST_DOMAIN, KALIS_MIN_TO_MAX,
                            depdist, true, 5)
 cp_set_branching(Strategy)

! Objective: weighted total distance
 totDist = sum(c in CITIES) POP(c)*depdist(c)

! Solve the problem
 if not cp_minimize(totDist) then
  writeln("Problem is infeasible")
  exit(1)
 end-if

! Solution printing
 writeln("Total weighted distance: ", getsol(totDist),
        " (average per inhabitant: ",
	   getsol(totDist)/sum(c in CITIES) POP(c), ")")
 forall(c in CITIES) if(getsol(build(c))>0) then
  write("Office in ", c, ": ")
  forall(d in CITIES) write(if(getsol(depend(d))=c, " "+d, ""))
  writeln
 end-if

!-----------------------------------------------------------------

! Calculate the distance matrix using Floyd-Warshall algorithm
 procedure calculate_dist
 ! Initialize all distance labels with a sufficiently large value
  BIGM:=sum(c,d in CITIES | exists(LEN(c,d))) LEN(c,d)
  forall(c,d in CITIES) DIST(c,d):=BIGM

  forall(c in CITIES) DIST(c,c):=0    ! Set values on the diagonal to 0

! Length of existing road connections
  forall(c,d in CITIES | exists(LEN(c,d))) do
   DIST(c,d):=LEN(c,d)
   DIST(d,c):=LEN(c,d)
  end-do

! Update shortest distance for every node triple
  forall(b,c,d in CITIES | c<d )
   if DIST(c,d) > DIST(c,b)+DIST(b,d) then
    DIST(c,d):= DIST(c,b)+DIST(b,d)
    DIST(d,c):= DIST(c,b)+DIST(b,d)
   end-if
 end-procedure

end-model

This implementation contains another example of the use of a subroutine in Mosel: the calculation of the distance data is carried out in the procedure calculate_dist. We thus use a subroutine to structure our model, removing secondary tasks from the main model formulation.

Results

The optimal solution to this problem has a total weighted distance of 2438. Since the region has a total of 185,000 inhabitants, the average distance per inhabitant is 2438/185 ≈ 13.178 km. The three offices are established at nodes 1, 6, and 11. The first serves cities 1, 2, 5, 7, the office in node 6 cities 3, 4, 6, 9, and the office in node 11 cities 8, 10, 11, 12.