NBS

REFERENCE

:::

' PUBLICATIONS

NBSiR 86-3472

AlllOE t31S13

NATL INST OF STANDARDS & TECH nal Revenue Service

)uty Location Modeling

System - Programmer’s Manual for

Pascal Solver

Paul D. Domich, Richard H. F. Jackson, Marjorie A. McClain

David M. Tate

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Applied Mathematics

Gaithersburg, MD 20899

July 1986

Issued February 1987

-OC y: 100 arch Division

.U56 evenue Service

86-3472 20224

1987
jII
NBSIR 86-3472

MBS

RESEARCH

INFORMATION CENTER THE INTERNAL REVENUE SERVICE

POST-OF-DUTY LOCATION MODELING

Paul D. Domich, Richard H. F. Jackson, Marjorie A. McClain

David M. Tate

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Applied Mathematics

Gaithersburg, MD 20899

July 1986

Issued February 1987

Sponsored by:

The Research Division

Internal Revenue Service

Washington, DC 20224

SYSTEM - PROGRAMMER’S MANUAL FOR

PASCAL SOLVER !ny

r

Q

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige, Secretary

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

ABSTRACT

This report is a programmer's manual for a microcomputer system designed

at the National Bureau of Standards for selecting optimal locations of IRS

Posts-of -Duty . The mathematical model is the uncapacitated, fixed charge,

facility location model which minimizes travel and facility costs. The

package consists of two sections of code, one in FORTRAN and the other in

PASCAL. The FORTRAN driver handles graphics displays and controls input and

output for the solution procedure. This report discusses the mathematical

techniques used to solve the mathematical model developed and includes a Greedy procedure, an Interchange procedure, and a Lagrangian approach to

the related linear program. A description of these PASCAL routines and

definitions of key data structures and variables are provided.

Key words: Uncapacitated fixed charge facility location problem. Greedy

heuristic. Interchange heuristic, Lagrangian Relaxation.

2
It-

, sn#,.

,!K:’%* -:ji!«. . -liiv' ..: ' li'' , ."; nj . kn.
TABLE OF CONTENTS

Section I. INTRODUCTION 4

Section II. METHODS 7

1. Greedy Heuristic 7

2. Teitz-Bart Interchange 8

3. Graph Coloring 9

4. Lower Bounds and Lagrangian Relaxation. 10

Section III. USER'S GUIDE 11

1. System Requirements, Hardware 11

2. Using the Package 12

Section IV. THE CODE 13

1. General Outline 13

2 . List of Functions and Procedures 14

3. Key Data Structures 16

4. Definition of Key Variables 16

5. Input/Output Processing 19

BIBLIOGRAPHY 21

APPENDIX: PROGRAM LISTING 23

3
I. INTRODUCTION

The Internal Revenue Service Post-of-Duty Location System is a microcomputer package designed to assist IRS district planners in selecting

locations for Post-of-Duty (POD's) that will minimize total costs. This

paper is part of a series of reports documenting the POD location system.

The reports in the series are as follows.

1) The Internal Revenue Service Post-of-Duty Location Modeling System:

Final Report.

This report describes the post-of-duty location problem and its

mathematical model. This report discusses the data used in

calculating costs, describes the solution procedures, and provides a brief introduction to the computer implementation of the model (NBS

Contact: Richard H. F. Jackson).

2) The Internal Revenue Serv^ice Post-of-Dutv Location Modeling System:

User's .Manual.

This report is a user's guide for the post-of-duty location computer

system. This report gives the hardware and software requirements, the

instructions for installing the system, a description of data files,

and detailed instructions for operating the system (NBS Contact:

Marjorie A. McClain).

3) The Internal Revenue Service Post-of-Dutv Location Modeling System:

Programmer's Manual for FORTRAN Driver.

This report describes the FORTRAN driver which handles graphics

displays and controls input and output for the solution procedure. An

alphabetical list of the FORTRAN routines includes a description of

purpose, a list of variables, and the calling sequence (NBS Contact:

Marjorie A. McClain).

4) The Internal Revenue Service Post-of-Dutv Location Modeling Svstem:

Programmer's Manual for PASCAL Solver.

This report is a programmer's manual for the PASCAL solver and

describes the mathematical techniques used to solve the facility

location problem. Included are a Greedy procedure, an Interchange

procedure, and a Lagrangian approach to the related linear program. A

description of these PASCAL routines and definitions of key data

structures and variables are provided (NBS Contact: Paul D. Domich)

.

4
For the Internal Revenue Service (IRS) , the facility location problem

involves the placement of Posts-of -Duty (POD's) for a given tax district,

according to the following model: locate k POD's so as to minimize the total

"cost" of the allocation. This cost is the sum of the fixed costs incurred

by opening or closing POD's, the operating costs for open POD sites, and the

travel costs incurred by taxpayers and IRS personnel. The interested reader

may also refer to any introductory textbook in Integer Programming (for

example, Garfinkel and Nemhauser (1972), Hu (1969)) or one of the many

papers on this subject (for example, Cornuejols, Fisher, and Nemhauser

(1977), Erlenkotter (1978)) for a general mathematical description of the

facility location model.

In the model, data is aggregated to a 5 -digit zip code level. The travel

cost of serving a given zip code is a function of the Euclidean distance to

the nearest POD in the solution, the workload for that zip code in the

period of interest (for example, one year), and the difficulty of travel

between that particular zip code and the zip code in which the POD is

located.

The disaggregated data for the problem comes in a variety of forms . The

first is map data which includes co-ordinates used as the centroid to the

zip code area, along with a list of zip code boundary points and of boundary

segments of adjacent zip code areas. This data is provided by contract to

IRS from Geographic Data Technology Inc. The centroid co-ordinates are used

in calculating distances between zip code areas and for displaying the map

of the state. Boundary points are used to draw the state map, while the

list of adjacent zip code areas can be used to display the POD service

regions for a given solution to the POD location problem.

A second source of data is workload data from the IRS Individual and

Business Master Files and includes Examination, Collection, Taxpayer Service

and Criminal Investigation workload data. Opening costs for new

"potential” POD sites or closing costs for currently "existing" POD sites,

and the cost of operating a POD facility in a particular zip code area are costs determined by the individual IRS District Offices. A more complete

description of these costs follows.

For each zipcode-POD pair, workload is combined with the distance and

travel difficulty factors between the two locations to produce a single

factor which represents the cost of servicing the zip code by that POD site

5
(see Report 1 for more information) . The distance from zip code to POD is

calculated using centroid co-ordinates from the geographic data mentioned

above, while difficulty factors are user-specified.

The fixed costs represent the cost of opening a potential site or closing

an existing site while operating costs are associated only with POD sites

determined by the solver routine to be open. These costs are included

directly in the zipcode-pair cost factors and are implicitly handled by the

SOLVER routine. The model correctly accommodates the interactive changes

made by the user without affecting the actual opening or closing costs for

the POD sites; i.e., a currently existing POD site changed to a potential

POD site by the user will not incur an opening cost, nor will a closing cost

exist for a potential POD site designated by the user as currently existing.

Thus , the initial POD file should reflect the current POD configuration and

accurately specify opening and closing costs. Opening costs for potential

POD sites not contained in the POD data file are interactively set by the

user

The operating cost for a POD site is computed in part from the zip code

areas it services. For each zip code area the nvunber of tax returns

received is translated into a floor space requirement at a particular POD

site. The cost of the floor space being different for each POD site

requires that this cost be included with the travel cost associated with

that zipcode-POD pair. Other related costs, for example overnight travel

costs and parking costs, may also be added to this factor. This cost data

can be used to determine the objective function coefficients for the

facility location problem.

Finally, a list of zip code areas designated as potential or currently

existing POD sites and a list of those zip code areas required to contain a POD site in any solution are required. The maximum distance allowed between

a POD site and the zip codes it serves determines which zipcode-POD pairs

are considered by the SOLVER routines. This distance represents the maximum

distance either an agent or taxpayer is expected to travel and varies from

region to region. This data along with the number of POD sites desired in

the final solution define the constraints to the facility location problem.

Note: Tradenames and products mentioned in this report are not

endorsed by the National Bureau of Standards nor does reference in this

report imply any such endorsement.

6
II. METHODS

The method for finding a "good" solution to the IRS POD location problem

is based on two well-known and dependable heuristic procedures. The first

is the Greedy heuristic (see, for example, Cornuejols, ^ , (1977)) and

the other is the Teitz-Bart Interchange heuristic (Teitz and Bart (1968)).

Also used by the procedure is a graph coloring algorithm, called the

Sequential Least-first Interchange Algorithm (see Matula et ^. , (1972)), to

display the final solution graphically. Each procedure is discussed below.

1. The Greedy Heuristic

In its simplest form, the Greedy heuristic for adding a POD to the current

configuration proceeds as follows (see procedure GreedyADD in the Appendix)

.

1) Choose the "cheapest" POD site and assign all workload to that site.

2) Choose k, the final niimber of open POD sites desired in the optimal

solution.

3) Among all allowable POD locations not currently in use, select that

site S which would most diminish the total assignment cost for the

problem, were it added to the current solution.

4) If this improvement is positive, and fewer than k sites are currently

active, add site S to the active POD set, let k=k+l, and go to 3).

5) Else, stop.

The GreedySUB routine for removing a POD from the current configuration

operates in a similar fashion.

The above procedure has been modified to accommodate the presence of

feasibility restrictions for the IRS model. Specifically, because of the

limit on the maximum travel distance from POD site to zip code area, an initial feasible solution must be provided by the user as input to the

solver routines. Without these travel distance restrictions, step 1 would

yield a feasible solution, though possibly at a large cost. Since the

Greedy heuristic restricts itself to feasible assignments, it assumes that a feasible solution exists prior to altering POD sites.

7
The current POD configuration is an adequate initial solution provided

the distance limit is properly specified. Unfortunately, the pre-specified

distance limit may be less than the actual distance traveled between a zip

code and a POD site in the current configuration. Feasibility can be

attained by increasing the distance limit to the maximum actual distance

traveled. Note that as a result of altering the distance limit, the number

of feasible zipcode/POD pairs changes, which consequently, affects the

complexity of the problem.

The second modification is that the target number of facilities, k,

supercedes objective function tests; the Greedy heuristic adds or subtracts

facilities from the current set as long as feasibility is maintained.

It is possible that an increase in cost may result after adding a facility.

This may be a result of a large fixed cost associated with a particular POD

site or a temporary aberration in the current assignment which will be

adjusted later in the algorithm. The procedure will add the site regardless

of the effect on the objective value. This provides the user with control

over the number of open POD sites in the final solution, against the chance

that the number of facilities desired may be influenced by factors not

incorporated in the mathematical model.

In the event of such objective value degeneration, a warning message will

be printed to the user. Note that such worsening does not necessarily

imply that fewer facilities will yield an eventual solution which is better

than that yielded by a larger number of facilities. Rather, the Greedy

heuristic has exhausted all other advantageous POD sites given its initial

allocation. The final application of the interchange heuristic will attempt

to correct this objective function value deterioration.

Should the Greedy heuristic fail to find a feasible solution at some

iteration, the program will advise the user and continue with the last known

feasible number of facilities as the target nvimber in all subsequent calculations .

2. Teitz-Bart Interchange.

Once the target number of facilities has been allocated by the Greedy

heuristic, the solution procedure tries to determine a better solution with

the same specified nxamber of open POD sites. The procedure iteratively

8
locates pairs of POD sites, one which is presently selected and one not,

such that if the two are interchanged in the current configuration, the

overall cost is reduced. When no such pair exists, the routine terminates

with the last configuration. The following heuristic, which is a modified

version of that of Teitz and Bart (1968) is used;

1) Partition the set of allowable sites into two sets, A and B, where A

is the set of currently assigned sites and B is all other potential

POD sites.

2) Look for a pair of sites, a in A and b in B such that

(i) cost(A - {a} + (b)) is less than cost(A)

,

(ii) a is not required to be a POD site, and

(iii) A - (a) + {b} is feasible.

3) For all pairs satisfying 2, select that pair which produces the

largest improvement and exchange a for b in the set of active sites.

Go to 2

.

4) If no such pair exists, stop.

The modification of step 2 parts (ii) and (iii) are excluded in the

original reference which did not have the initial feasibility restrictions.

Because of the travel distance limit previously described, an initial

feasible solution is required. The combination of the Greedy heuristic

followed by the Interchange heuristic is well known to produce very good

solutions to the facility location problem (see, for example, Cornuejols,

Fisher, and Nemhauser (1977))

3. Graph Coloring.

To display in color the final assignment of zip code areas to POD

locations, it is necessary to ensure that no two adjacent POD service areas,

i.e., two areas sharing a common border, are colored with the same color.

This is a map coloring problem, where the regions involved are groups of

customers aggregated by their assigned POD facility. The problem is to

choose colors for the regions of a graph G, such that is not equal

to Cj if and are adjacent regions, and in such a fashion that a "small" nxomber of colors are used. Since all of the zip code maps can be

represented as planar graphs (i.e., graphs that can be drawn on a sheet of

paper so that no two edges cross)

,

theoretically all can be colored using

only four colors. In practice, to find a four-coloring is a very difficult

problem, so a five- or six-coloring is used. For a description of the

9
coloring algorithm, see Matula, ^ "Graph Coloring Algorithms", (1972).

The procedure used is called the Sequential Least- first Interchange

heuristic (SLI) and is presented in the Appendix.

4. Lower Bounds and Lagrangian Relaxation.

As previously stated, the Greedy heuristic and the Interchange heuristic

described above are well-known to produce good solutions to the facility

location problem. One drawback with these procedures involves determining

when the generated solution is in fact the optimal integral solution to the

described problem. One way to demonstrate the optimality of a solution

involves generating lower bounds to the optimal objective function value.

One bound can be obtained by solving the linear programming (LP) relaxation

of the original problem, i.e., the original problem without the integrality

constraints

.

In general the LP formulation of the facility location problem has a large number of constraints in the problem description and it, too, can be

difficult to solve. One Lagrangian relaxation of this LP problem removes

the requirement that a zip code is serviced by exactly one POD and adds a penalty to the objective value for any violation of these constraints. This

relaxation can produce the optimal LP objective function value in an

iterative manner, and provide at each step a lower bound to the optimal

solution to the original facility location problem. Further, by rounding

the possibly fractional real -valued solution produced by this method, an improved integral solution may be found as a by-product. The interested

reader may refer to the many articles in this subject (e.g., Cornuejols, ^

(1977), Fisher(1982))

.

The Lagrangian solution procedure will relax the constraints requiring

that a zip code area can be serviced by exactly one POD site while

penalizing the objective function for any violation in these constraints. A

feasible solution to the relaxed problem is found and the penalty factors

are modified in a manner which forces the relaxed constraints to be

satisfied. This iterative procedure generates a series of objective values

which are lower bounds to the optimal integral objective value to the

original problem. Further, by rounding the real -valued objective value

10
produced by this method, an improved integral objective value may be found

as a by-product of the procedure.

Often for the facility location problem, the optimal LP objective value

is equal to the optimal integral objective value (see, for example, Morris

(1978)) and therefore the optimality of the heuristic integral solution can

be demonstrated using the real-valued objective value. Otherwise, either

there exists a "gap" between the optimal LP objective value and the optimal

integral objective value, or the integral solution is nonoptimal. In the

latter case, the bound provides an estimate on the "goodness" of the

integral solution value.

III. USER'S GUIDE

1. System Requirements

The SOLVER package, and the graphics environment in which it runs, are

written specifically for systems running MS/DOS on an IBM PC compatible

(Intel 8088 -based) microcomputer with a math co-processor and a lOmb fixed

disk. When executing the FORTRAN driver routines, it is essential to have

the math co-processor to ensure correct type -matching in the input data

files produced by driver routines, as well as desirable speed of execution.

The graphics capability is provided via a number of different hardware and

software functions. Included are the following;

Graphics Display Monitor,

Graphics Expansion Card,

IBM Graphical Kernel System.

The SOLVER routines are written in TURBO PASCAL (Borland International

Inc., (1985)). There are several reasons for choosing Pascal as the

language for the SOLVER and these are summarized below.

1. Pascal has a dynamic storage capability, permitting a more efficient

use of core memory than is possible in static-allocation languages like

FORTRAN, BASIC, and APL. This is essential to solve large problems.

2. TURBO PASCAL compiles about one order of magnitude faster than other

available Pascal compilers, and several orders of magnitude faster than

available FORTRAN. As an example, the 1950 lines of code in the SOLVER

program compile and link in under 30 seconds, to a file of only 18K

11
bytes. A similar FORTRAN code requires over 6 minutes to compile and

link and has a much larger storage requirement.

3. Pascal supports pointer variables and structured data-t3rpes (user

defined records), making for much more legible, structured, and easily

altered code.

4. TURBO PASCAL is about one- fifth the price of most other Pascal or FORTRAN packages, and includes a number of graphics and utility programs

in this price. It runs its own developmental operating system, and traps

and locates run-time errors automatically, thus greatly enhancing program

development.

The flexibility provided by the Pascal programming language allows

development of a well -structured program which is easily understood. The

only limitation of the language in this application involved data transfers.

This problem was resolved using a FORTRAN unformatted write statement in the

preprocessor graphics routines which create the data files used by the

SOLVER routines. I/O issues are discussed in Section 4.

2. Using the Package

The SOLVER package is used as a subprogram to the IRS POD Location

Modeling System which performs all preprocessing of input data and

graphically displays workload data and SOLVER'S final solution. As input

data the SOLVER routine requires a single file (called STATE .DBL) that is

automatically generated by the graphics package (see the IRS Post-Of-Duty

Location Modeling System: User's Manual). This file defines the facility

location problem and contains information about the individual zip code

areas and also specifies assignment costs from zip code areas to the

feasible POD sites. For computational efficiency, this file is written in

binary format. The exact commands needed to call SOLVER from the main

program are discussed in detail in the report mentioned above.

Once the driver routine generates the input files for the solver

routines, the user is provided with a summary of the problem characteristics, followed by a query to the user for additional information on the

number of POD's desired in the final solution. Once the current problem is

fully described, control is passed to the SOLVER routines and the following

steps occur. The following text illustrates this phase of the program:

12
Total number of zip codes is

Nvnnber of possible POD's is

Number of existing POD's is

non.

111 .

Enter the desired number of POD’s in the final solution: kkk

where nnn, mmm, 111, and kkk are integer values. After the last prompt has

been answered, the solver proceeds to solve the POD location problem. An

in-depth examination of the solver routine is given in the next section.

IV . THE CODE

1. General Outline

The structure of the solver routines involves four basic program units.

The first performs the input of the facility location problem as defined in

the pre-processor graphics package (see the IRS Post-Of-Duty Location

Modeling System: User's Manual). The problem file is read and entered into

the data structures and, from the existing configuration of POD's, an initial interchange is performed so as to locate the best possible solution

given the original number of POD facilities. Next, the number of POD's is

altered by adding or deleting POD's as required via the Greedy heuristic.

Upon termination of the Greedy heuristic, a final interchange is performed

which seeks the best possible solution of the given size. Finally, a graph

coloring is performed so as to display the POD service areas in the final

solution.

The code for the solver portion of the package is found in the files;

SOLVER. PAS,

INIT.PAS

,

DSTRUCT . PAS

,

GREEDY . PAS

INTCHG . PAS

,

FIVCLR.PAS,

PODCLR.PAS,

LGRN . PAS

.

13
The SOLVER file contains the driver program as well as routines to compute

the cost of an allocation (i.e., assign customers to their nearest facility)

and to output the current solution. The routine INIT performs array, set,

and pointer initializations. The sparse-matrix data structures determined

from the input data are set up by DSTRUCT. The procedure GREEDY performs

the Greedy heuristic calculations, while INTCHG is the interchange heuristic

algorithm. General utility routines, denoted as InsertPOD and DeletePOD,

add and delete POD's from the data structure, re-establishing the data

structure for the new set of POD's. The LGRN routine determines a lower

bound on the optimal integer-valued solution to the problem and can be used

to verify the solution found by the Greedy and interchange heuristics. The

FIVCLR and PODCLR routines are used to determine a coloring of the final

solution map for displaying the POD service regions.

2. List of Functions and Procedures;

The following is a list of procedures and functions, and their purposes:

FUNCTIONS

SwapVal(old,new) Returns the change in objective function value

associated with an exchange of facility "new"

for facility "old" in the current set of

facilities

Exist(filename) Boolean function returning true if the string

"filename" is the name of a current disk file,

false otherwise.

PROCEDURES

Match

ComputeCost

ListCurrent

Associates with each customer area (zip code)

the nearest facility in the current set of

facilities. Values are set for the arrays

BestPOD[], NextBestPOD[ ] , CurrentCost [ ]

NextCost [

]

Adds all assignment costs to find the current

objective function value.

Sends a list of current POD assignments to the

default list device.

Initialize Zeroes arrays, empties sets, and NIL's

pointers prior to program execution.

14
CreateDataStructures

Greedy

GreedyADD

GreedySUB

Interchange

InsertPOD

Remove POD

GraphColor

Lagrangian_

Quick_Sort

Reads zip code data from the special file

STATE .DBL then establishes the sparse array

data structure which has cost and feasibility

data for specific POD allocations. Rows of

the array are pointed to by the vector of

pointers Map[], and columns are pointed to by

the pointer fields of the vector of records

CanBe[]. As each record of data is read for a feasible zipcode/POD pair, an entry in the

sparse array is created, specifying the zip

code index and the zip code index of the POD

site involved, the cost, and a pointer to the

next zip code entry for that POD site and a pointer to the next POD site for that zip

code. This record is inserted in the data

structure ordered by increasing cost. Rows

correspond to all POD sites which may feasibly

serve a given zip; columns correspond to all

zips which may be feasibly served by a given

POD site.

Performs the Greedy heuristic as described

in Section II. 1.

Increases the number of POD's by one,

according to the Greedy heuristic.

Decreases number of POD's by one, according to

the Greedy heuristic (if feasible).

Performs interchange heuristic on problem, as described above

.

Performs the insertion of a POD to the current

set and updates the BestPOD[], NextBes tPOD [ ]

,

CurrentCost [ ] , NextCost[] arrays.

Performs the removal of a POD from the current

set and updates the BestPOD[], NextBestPOD [ ]

,

CurrentCost [] , NextCost[] arrays.

Performs the sequential least- first

interchange coloring algorithm on the graph of

the final solution, coloring POD "spheres of

influence" to avoid having identical adjacent

colors

.

dual Computes a lower bound on the best possible

solution to the problem. Can be used to

verify the optimality of the heuristic

solutions

.

Performs a sort of a vector of real numbers.

15
3. Key Data Structures:

The key data structure in the solver program is a doubly linked- list for

maintaining the zipcode/POD pair data. The basic element of this structure

is a five-field record, defined as follows:

(1) node,

(2) target,

(3) cost,

(4) nextZip, and

(5) nextPOD.

The "node" field is the index of the zip code for this record. The

"target" field is the index of the potential POD site to which this node

refers. "Cost" is the cost of assigning zip "node" to POD "target" (if

node=*target , then this also includes the fixed operating cost of having a POD at target) . The entry "nextZip" is a pointer to the next record which

refers to POD site "target", and "nextPOD" is a pointer to the next record

which refers to zip code index "node".

Map [ 1 . .MaxZips ] is an array whose entries for any given zip, are pointers

to the linked records by POD, and CanBe [ 1 . .MaxPossible ] is an array whose

entries are records, one field of which is, for any given allowable POD

site, a pointer to the linked records by zip code index. Thus, starting

with Map [27] and following the "nextPOD" links results in a linked list of

records corresponding to all possible POD's which can serve zip code index

#27 with their associated costs. This linked list of potential POD sites is

sorted in order of increasing cost.

Similarly, starting with CanBe [11] .next (the pointer field of the 27th

entry of array CanBe) and following the NextZip links produces a linked list

of records corresponding to all zip code indices which can be served from

the 11th allowable POD site. Both of these data structures are static, in

the sense that once they are created (by procedure CreateDataS true Cures)

they will never change.

4. Definition of Key Variables:

There are certain global variables in the program that the programmer

should be familiar with before attempting to modify the code. This section

16
will list the most important variables and their definitions and structures

(if any) . First the various Pascal constants and variable t}rpes are introduced.

CONSTANTS

:

MaxZips = 1000; Maximum number of zip code areas

allowed. This constant may be

changed.

MaxPossible = 85; Maximum number of possible POD sites

allowed. This constant may be

changed but can not exceed 256.

These two constants determine the size of the various storage arrays used in

the SOLVER routines. Consequently, limiting the size of these constants

will lower the storage requirements for the system.

VARIABLE TYPES:

Zcode = 0.. MaxZips; Integer type in the range [0,Maxzips]

ZipSet = set of 1 . .Maxpossible

;

NOTE: The "set" data type is an

implementation dependent type. TURBO

PASCAL allows set types up to 255

distinct possible elements. This

means that in no case can MaxPossible

be set to a value of more than 256.

Link » ^Neighbor; A Pascal pointer type for record of

type Neighbor.

Neighbor

Site

Target

Record type which includes:

zip code index of type Zcode,

zip code index of POD's of type

Zcode

,

Cost

NextZip

NextPOD

cost of Site-Target assignment,

pointer to the next zip code,

pointer to the next POD zip code.

PODsite

Where

Must

Next

Record type which includes:

the zip code index of type Zcode,

boolean flag for a required POD site,

a Pascal link to the first of its

neighbors

.

SingleZip

zip code

SType

Data record type which includes:

actual zip code number,

site type;

17
FixCost

SType=0 => never a POD site,

SType=l => can be a POD site,

SType=2 => must be a POD site,

FixCost Opening/Closing cost for a POD site.

PairOfZips Data record type which includes:

Number zip code index of zip code area,

PODnum zip code index of POD site,

Cij cost of area to site assignment.

Co lumnPo intArray Array type of length MaxPossible of

PODsite

.

RowPo intArray Array type of length MaxZips of Link.

IndexArray Array type of length MaxZips of

Zcode

.

ValueArray Array type of length MaxZips of real.

FileString Character string of length 15.5.

VARIABLE DEFINITIONS

CurrentPODs The set of POD's in the current

assignment

.

PossiblePODs The set of all possible POD sites.

BestPOD[l. .MaxZips] The POD index which is the nearest

POD in the current solution.

NextBestPOD [ 1 . .MaxZips] The POD index which is the second

best POD in the current solution.

CurrentCost [ 1 . .MaxZips] Value of the cost of the BestPOD.

NextCost[l. .MaxZips] Value of the cost of the NextbestPOD.

CanBe[l. .MaxPossible] This is an array of Co lumnpo intArray

storing the index of the POD site and

a pointer to the linked list of zip

code areas reached from that POD

site

.

Index[l. .MaxZips] Pointer for all possible POD sites to

records in array CanBe

.

Map [ 1 . .MaxZips

]

This is an array of RowPo intArray

pointing to the start of the linked

list of feasible POD sites for a zip

code area.

ZCreal[l. .MaxZips] The actual zip code number.

18
CurrentNumber

EndNumber

Nzips

Nposs

Switch

TotalCost

Error

Change

ErrLoc

The current number of POD sites.

The desired number of POD sites in

the final solution.

The total number of zip code areas.

The total number of possible POD

sites

.

0 if graph- coloring is used, and

1 otherwise.

The current objective function value.

A flag to warn that the solver has

run into a situation where the user's

wishes cannot be satisfied; e.g. no feasible solution exists using only

EndNumber POD sites.

Flag indicating whether any swapping

was performed by the interchange

heuristic

.

The site of ERROR if true.

MinCode The smallest zip code number in the

state

.

StateNumber

StateNameFile

The two-digit state code number.

The name of the state.

6. Input/Output Processing:

Input to the SOLVER routines comes from the STATE .DEL file where

refers the index of the tax district (1 to 76) . The DBL file is written by

the driver routines using an unformatted FORTRAN write statement. This file

consists of sets of records, each set preceded by, and followed by, a two

byte word indicating the total number of bytes used in that set (see the IBM

Professional FORTRAN Reference Guide) . The following is a representation of one such set:

Wordl , i

,

ZIP. , r

j, , C .

. , type . , Wordl

.

19
The first parameter, Wordl is used by the CreateDataStructures routine to

determine the number of elements in the set. The set involves index i

having zip code ZIP^ which is of type type, and has feasible POD assignments

to j,, j, , at a cost of C , C.. , C.. , respectively. The

costs are in decreasing sorted order except possibly for the last record

which, if the index is also a POD site, contains the operating cost for that

site

This type of data transfer is very efficient. Alternative methods of

transfering large amounts of data from a FORTRAN to a Pascal program

consumed nearly twice as much time. Further, all of the problem information

for the SOLVER routines is contained in a single file. This includes travel

costs, floor space rental costs, operating costs and fixed opening and

closing costs. The latter two costs are included into the C. . factors above

ij

before the data transfer is performed.

The Pascal input is performed in a pairwise form. Each pair consists of

a two byte integer followed by an eight byte real number. The

CreateDataStructures procedure reads 116 pairs at a time and processes the

vector of information sequentially. The length of the vector is arbitrary.

To ensure proper sequencing of the Pascal read statements with the DEL file,

additional zero entries are inserted during the FORTRAN write statement.

Output from the SOLVER routines is stored in the STATE .SOL file.

Included in this output is the index of the zip code, its assigned POD, and

a mamber indicating the color determined by the graph coloring algorithm for

this zip code area. The STATE .SOL file is used by the driver package to

display the final solution.

20
BIBLIOGRAPHY

Borland International Inc., TurboPascal version 3.0, Reference Manual,

Scotts Valley Drive, Scotts Valley, CA.

,

1985.

Bradeau, M.L.

,

and Chiu, S.S., "Sequential Location and Allocation;

Characterization and Estimation of Globally Optimal Solutions", Department

of Engineering-Economic Systems, working paper, Stanford University, 1984.

Cornuejols, G. , Fisher, M.L.

,

and Nemhauser, G.L.

,

"Location of Bank

Accounts to Optimize Float; An Analytical Study of Exact and Approximate

Algorithms", Management Science . vol. 23, 789-810, 1977.

Domich, P.D.

,

Hoffman, K.L.

,

Jackson, R.H.F., and McClain, M.A.

,

"The

Internal Revenue Service Post-of Duty Location Modeling System;

Final Report", National Bureau of Standards Technical Report NBSIR 86-3482,

Gaithersburg, MD, July, 1986a.

Domich, P.D.

,

Jackson, R.H.F., and McClain, M.A.

,

"The Internal Revenue

Service Post-of Duty Location Modeling System; Programmer's Manual for the

FORTRAN Driver", National Bureau of Standards Technical Report NBSIR 86-

3473, Gaithersburg, MD, July, 1986b.

Domich, P.D., Jackson, R.H.F., and McClain, M.A.

,

"The Internal Revenue

Service Post-of Duty Location Modeling System; User's Manual", National

Bureau of Standards Technical Report NBSIR 86-3471, Gaithersburg, MD, July,

1986c.

Erlenkotter, D. , "A Dual-Based Procedure for Uncapacitated Facility

Location", Operations Research , vol. 26, no. 6, 992-1009, 1978.

Fisher, M.L.

,

"The Lagrangian Relaxation Method for Solving Integer

Programming Problems", Management Science , vol. 27, no. 1, 1982.

Francis, R.L.

,

and White, J.A., Facility Lavout and Location . Prentice-Hall,

Englewood Cliffs, NJ , 1974.

Geographic Data Technology, Inc., 13 Dartmouth College Highway, Lyme, NH.

Garfinkel, R.S. and Nemhauser, G.L.

,

Integer Programming . John Wiley & Sons,

New York, NY, 1972.

Hu, T.C., Integer Programming and Network Flows . Addison-Wesley , Reading MA,

1969.

Morris, J.G., "On the Extent to which Certain Fixed-Charge Depot Location

Problem? Can Be Solved by LP"

,

Journal of the Operational Research Society . vol. 29, no. 1, 71-6, 1978.

Matula, D.W.

,

Marble, G. , and Isaacson, J.D., "Graph Coloring Algorithms,"

in R.C. Read, Graph Theory and Computing . Academic Press, New York, NY,

1972.

21
Ryan-McFarland Corporation, IBM Personal Computer Professional FORTRAN:

Reference, for IBM Corp . Boca Raton, Florida, 1984.

Teitz, M.B. and Bart, P. , "Heuristic Methods for Estimating the Generalized

Vertex Median of a Weighted Graph," Operations Research . 16, 955-61, 1968.

Wagner, H.M. , Principles of Managment Science, Prentice -Hall , Englewood

Cliffs, NJ, 1975.

22
APPENDIX: Program Listing

Program Solver;

{ This is the main driver program for the package that finds good

heuristic solutions to IRS Post-of-Duty (POD) location problem.

The program takes data from specially formatted data files, which

have been created by a separate pre-processing package. The final

solution is colored, and the result may be saved for graphic display on a map of the region in question.

}

Const MaxZips = 800;

MaxPossible = 40;

{ Maximum number of zip -code areas allowed and the maximum number of

possible POD sites allowed. These numbers are somewhat flexible,

although 1000 may not be large enough for some districts and 75 is

probably more than we need for any district )

Type

Zcode =» 0.. MaxZips;

ZipSet =- set of 1 . .Maxpossible

;

( ! !N0TE! ! : The "set" data type is an implementation-dependent type.

TURBO Pascal allows set types up to 255 distinct possible elements.

This means that in no case can MaxPossible be set to a value of more

than 255 .

]

Link = ^Neighbor; { a pointer to a record of type Neighbor ) Neighbor = record

Site, Target

Cost

NextZip, NextPOD

end;

{Each "Neighbor" record is one entry in the sparse matrix of

information relating zip codes to POD sites. The field SITE indicates

which zip-code area the information in the record applies to. The

field TARGET tells with reference to which POD site. COST gives the

cost of travel between SITE and TARGET (which will always be less than

the user- supplied upper limit on travel distance for any customer.

NEXTZIP is a pointer to the record which holds the next-nearest zip

code to TARGET (after SITE), if there is one. NEXTPOD is a pointer to

the record which holds the POD possible location which is the nextnearest (after TARGET) to SITE, if one exists.)

PODsite - record

Where ; Zcode; { Which site is this? ) Must : boolean; { Is it a required site? ) Next : Link; ( a pointer to the first of its neighbors)

end;

SingleZip - record

Zcode

;

Real

;

Link;

23
ZipCode : real;

Slype : integer;

end;

{actual zipcode niomber, a real var because integers in TURBO are not

large enough. Site type:

SType=l => not a POD site, ever. SType=2 => can be a POD site.

SType=3 => is now a POD site.

SType=4 => must be a POD site. )

PairOfZips * record

Number

,

end;

PODnum

Cij:

integer; (site, target *POD, edge)

real; { type: Cij weights]

( i,j assignment.)

Co lumnPo intArray

RowPo intArray

IndexArray

ValueArray

FileString

array [ 1 . .MaxPossible ] of PODsite;

array[l. .MaxZips] of Link;

array[l. .MaxZips] of Zcode;

array [1. .MaxZips] of real;

string[14]

;

var

CurrentPODs , PossiblePODs : Zipset;

(CurrentPODs is the set of all POD's assigned in the current solution.

PossiblePODs is the set of all possible POD sites.)

BestPOD, NextBestPOD : IndexArray;

{BestPOD holds the zip which is the nearest POD in the current

solution. NextBestPOD holds the secondbest current POD for each zip.)

CurrentCost, NextCost : ValueArray;

(CurrentCost [zip] holds the cost from zip to BestPOD [zip] in the

current solution. Similarly, NextCost [zip ] is the cost from zip to

NextBestPOD [zip] .

}

CanBe : ColumnPointArray

;

{CanBe is an array which allows us to find all the pertinent data

concerning the Jth potential POD site: which site it is, which zips

can be served from it, and how much that would cost. Its field NEXT

points to a column of Neighbor records, along the NEXTZIP links.)

Index : Array [1 . .MaxZips ] of Zcode;

{Index[i] tells which entry in CanBe refers to POD i )

Map : RowPo intArray;

(Map is an array which lets us find, for any zip area, which POD sites

24
can serve it and how much that costs. Each entry of Map is a pointer

to a row of Neighbor records, along the NEXTPOD links.)

ZCreal : ValueArray;

{Actual zip code number in real format)

RawDat : SingleZip;

{Will hold individual zip area data for the construction of the data

structures .

)

CurrentNumber , Switch, Nzips, Nposs, EndNumber : integer;

{ CurrentNumber is the number of POD sites assigned in the current

solution. Switch is 0 if graph coloring is to be used, 1 otherwise.

Nzips is the niomber of zip code areas ( <= MaxZips ) . NPOSS is the

number of possible POD sites ( <= MaxPossible ) . EndNumber is the

number of POD sites the user has requested be in the final solution.

)

TotalCost, Limit : real;

{TotalCost is the current objective function value. Limit is the usersupplied upper bound on travel distance.)

Error, Changes, Stuck : boolean;

{ERROR is a flag to warn that the solver has run into a situation where

the user's wishes cannot be done; e.g., no feasible solution exists

using only EndNtunber POD sites. Changes indicates whether any

swapping has been done in the interchange heuristic.)

ErrLoc : Zcode; { Site of ERROR if true ) MinCode : real; { Smallest zip-code in state ) StateNumber : string[2] ; { Two-digit state code number ) StateNameFile : text;

function exist(fn: FileString) :boolean;

{TURBO PASCAL FUNCTION returns true if file fn already exists } var fil:file;

begin

assign( f il , fn)

;

{$!-}

reset(fil)

;

($i+):

exist (lOresult - 0);

end

;

25
procedure Match;

( Given the contents of CurrentPODs and the arrays of neighbor data, this

procedure determines the nearest and next-nearest currently assigned POD for

each individual zip-code area, and the associated costs.)

var

base

zip , pod

empty, done

ipod, izip

begin

TotalCost := 0.0;

error := false;

for zip :=• 1 to Nzips do

{ find the first current POD in zip's list of possible POD's, and assign

zip to it. } begin

done :=• false;

base := map[zip]

;

if base=nil then

done : = true

;

while not done do

if base=nil then { no POD is close enough, so this is illegal: ) begin

done := true;

error := true;

writeln( ' feasiblity error at ',zip:5);

end

else

begin

pod :=« base^ . target

;

ipod := Index [pod]

;

if ipod in CurrentPODs then { pod is the best choice: } begin

done : = true

;

BestPOD[zip] :=- pod;

CurrentCost[zip] :=• base"'. cost;

NextBestPOD [ zip ] := 0;

base := base'' .nextpod;

empty false;

while not empty do { see if there's a next-best POD: ) if base=»nil then { there isn’t a next-best: ) empty := true

else if Index[base'' . target ] in CurrentPODs then

begin {this is next best ) NextBestPOD [zip] :=» base^ . target

;

NextCost [zip] :=• base^.cost;

emp ty : =» true

;

end

else

base := base'' . nextpod; { keep looking for a next-best ) end ( if POD in CurrentPODs . . .

)

link;

zcode

;

boolean;

integer

;

26
else

{writeln(' pod not in CurrentPODS ',pod:5);)

base := base^ .nextpod; { keep looking for a best POD }

end; {while not done...}

end; { for zip := 1 to . . .

}

end; { Procedure Match )

( }

procedure ConiputeCost

;

( Just add up all CurrentCost values to find the new total cost. Recall

that fixed costs are included in the i,j assignment costs)

var

zip : Zcode;

begin

TotalCost := 0.0;

for zip := 1 to Nzips do

TotalCost TotalCost + CurrentCost [zip] ; end;

{ ^

) procedure Lis tCurrent;

( For larger problems, modify this to only print out POD sites}

var i : Zcode

;

begin

{ writeln( ' Current zip-code assignments:');

for i:= 1 to NZips do

writeln(i: 5

,

' at ' ,BestP0D[ i]

, ' : cost - CurrentCost [ i] : 3 : 2) ;

}

CoraputeCost

;

writeln( {LST ,

}

' Total cost of this allocation is ' .TotalCost : 12 : 2)

;

end;

( }

procedure dumps true t;

{ This is a diagnostic procedure which prints out the contents of the sparse

matrix structure set up in procedure CreateDataStructures

}

var zip.i ; zcode;

ptr : link;

begin

27
{ for zip :=» 1 to Nzips do

begin

ptr:- Map [zip]

;

while ptrOnil do

begin

write(LST,ptr'' . target : 5)

;

ptr := ptr^ . nextpod;

end;

writeln(LST)

;

writeln(LST)

end; ) for zip := 1 to Nposs do

begin

ptr :=» CanBe [ zip ]. next

;

write (LST,CanBe[ zip] .where: 5)

;

while ptrOnil do

begin

write(LST,ptr^ .site: 5)

;

ptr := ptr'' .nextzip;

end;

writeln(LST)

;

writeln(LST)

end;

end;

{ k-k'k'k-kic'k'k'k-k-Jrk-ic'k MAIN PROGRAM '**'*’*'*'*'*'**'**'*'*'*'*'**'*'*'*'5t*'*-!t'*-5!r'5t'5Sr-*-A-

]

{$i init.pas } {$! dstruct.pas)

[$I greedy. pas } {$! intchg.pas } {$! fivclr.pas } {$! podclr.pas ) {$1 Igrn.pas }

{ Include array initializations )

{ Include data -structure initialization package }

{ Include greedy heuristic routines }

{ Include interchange routines }

{ Include graph- coloring algorithm }

{ Include POD-coloring algorithm )

{ Include Lagrangian Dual algorithm)

{ Read the number of the state under consideration }

begin

Assign(StateNameFile , 'NAMES ' )

;

reset(StateNameFile)

;

read(StateNameFile , StateNumber)

;

close(StateNameFile)

;

Initialize

;

CreateDataStructures

;

writeln;

write(' Enter the desired number of POD’s in the final solution:');

readln(EndNumber)

;

Match;

ClrScr

;

if error then

writeln( { 1st ,)' Initial allocation is not feasible- -program aborted.')

else

begin

28
writeln( ' **********•***** INITIAL ASSIGNMENT ' ) ; ComputeCost

;

ListCurrent

writeln;

writeln(' it'k'k-k-k-k-k-k'k-k-k-kic^'k INITIAL INTERCHANGE k-k-k-k-k^-k-k-k-k-k-k-k'k'k-k writeln;

Interchange

;

if not changes then

writeln( {LST ,

)

' No interchanges were necessary.');

ComputeCost;

ListCurrent

;

if EndNumber O CurrentNumber then

begin'

writeln;

write In (

' GREEDY HEURISTIC k-k-k-k-k-k-k-k-k-k'k'k'k'k'k-k ' ) j writeln;

Greedy;

if changes then

begin

if error then

writeln( ' Greedy heuristic solution is not feasible')

else

begin

writeln;

writeln( '  •** ***•*•***•*•* FINAL INTERCHANGE •*•***•*'*******•****• ' ) ; writeln;

Interchange

;

if not changes then

writeln( {LST , )

' No interchanges were necessary.');

end;

end

end;

Match;

ComputeCost

;

ListCurrent;

{perform the optional routine to find a lower bound on the optimal solution)

{ writeln( ' start lagrangian dual ' )

;

lagrangian_dual{ (NextCost , NextBestPOD) ) { ; writeln( ' end lagrangian dual ' ) ; }

if Switch-0 then

begin

writeln;

writeln( ' Calculating colors for solution map - Please wait');

GraphColor

end

else

PODColor;

end;

end.

{ )

29
procedure initialize;

{This procedure initializes various data arrays, pointers and sets

used by the solver package.)

var

i,j : integer;

begin

for i;=l to MaxZips do

begin

BestPOD[i] ;= 0;

CurrentCost [ i] := 0.0;

NextBestPOD[ i] :=• 0;

NextCost[i] := 0.0;

Index [ i ]

; = 0

Map[i] :=• nil;

end;

for i:= 1 to MaxPossible do

CanBe[i] .next :=• nil;

PossiblePODs := [];

CurrentPODs := [];

end;

{ ) procedure CreateDataStructures

( This procedure creates the sparse matrix structure which holds the

information concerning which zip code area can be served from which POD

sites, and at what cost. The data structure is a cross-linked array, with

row links joining all PODs which can serve a given zip, and column links

joining all zips which can be served by a given POD site. The entries are ordered along both row and column lists in order of increasing cost. )

type

pair =• record

iteml: integer;

item2: real;

end;

pair_vec = array [ 1. . 116] of pair;

pair_rec = record

pair_vector; pair_vec;

end;

var pairs

i, j ,k, pointer

pair_file

hold

C_ij

filename

TempPair

DoubleFile

count , POD , t

pair_rec

integer;

file of pair_rec;

Array [1. .Maxpossible] of integer;

Array [1. .Maxpossible] of real;

Files tring;

PairOfZips

file of PairOfZips;

30
pt,pl,p2

scanning

PODnum , Number

ipod, izip integer;

integer;

link;

boolean;

begin {MAIN PROCEDURE)

count := 0;

Nposs := 0;

Nzips :=> 0;

mincode := 99999.0;

CurrentNumber := 0;

filename ;=' STATE' +Statenumber+' .ADJ'

;

if exist (filename) then

Switch := 0

else

Switch ;= 1;

Assign(Pair_file ,' STATE '+StateNumber+' .DBL' )

;

reset (pair_file)

;

read (pair_file , pairs )

;

pointer := 1;

while pairs .pair_vector [pointer] . iteml O 0 do

k ;= pairs .pair_vector[ pointer] . iteml;

k : k div 10 - 2

;

for j 1 to k do

begin

pointer := pointer+1;

if pointer < 117 then

begin

hold[ j ] pairs .pair_vector [pointer] . iteml

;

[ j ] pairs .pair_vector [pointer ]. item2

;

end

else

begin

pointer 1;

read( pair_file, pairs);

hold[ j ] pairs .pair_vector[ pointer ], iteml

;

C_ij [ j

]

pairs .pair_vector [pointer ]. item2

;

end;

end;

for j 1 to k do writeln(hold[ 1] : 3

, '

' ,hold[ j ] ; 3 , C_ij [ j ] ) ;

)

if pointer - 115 then

begin

pointer 1;

T pairs .pair_vector[ 116 ]. iteml

;

read( pair_file, pairs);

end

else if pointer - 116 then

begin

pointer 2;

read( pair_file, pairs);

T pairs . pair_vector [ 1 ]. iteml

;

begin

31
end

else

begin

T pairs .pair_vec tor [pointer+1] . iteml

;

pointer := pointer+2;

end;

{ start entering the i,j data )

{ i,j entries are assumed to be in sorted order, except for i,i)

number;— hold[l];

count := succ(count);

ZCreal [number ] := C_ij[l];

if C_ij[l] < mincode then mincode := C_ij[l];

if Nzips < number then Nzips ;= number;

if T > 1 then { this is a possible POD site } begin

if Index [number] = 0 then

begin

Nposs succ(Nposs)

;

Index [number] := Nposs;

POD ;= Nposs

end

else

POD := Index [number]

;

CanBe[POD] .where ;= number;

if T=4 then

CanBe [POD] .must := true

else

CanBe [ POD] . must :=» false;

PossiblePODs :=• PossiblePODs + [POD] ; case T of

3,4 ; begin

CurrentPODs := CurrentPODs + [POD]

;

CurrentNumber := succ(CurrentNumber)

end;

end;

end;

for j :=2 to k do

begin

new(pt)^

with pt'' do

begin

site :=• number;

target ;= hold[ j ] ; cost := C_ij [j ] ; nextzip:- nil;

nextpod;- nil;

end;

pi := Map[number];

if pi = nil then

Map [number] := pt

32
else

begin

{ if pi''. cost < C_ij [ j ] then

writeln( ' scanning ' ,pl^ .cost ,C_ij [ j ] ,

' index [n\imber] ) ; } Map[ntimber] :=• pt;

pt''.nextpod :=» pi;

end; { if pl=nil . . . else . . .

}

if Index[hold[ j ] ] = 0 then

begin

Nposs :=• succ(Nposs);

Index[hold[ j ] ] :=• Nposs;

POD :=- Nposs

end

else

POD :=• Index[hold[ j ] ] ;

pi ;=» CanBe[POD] .next;

if pi =• nil then

CanBe [ POD ]. next : =»pt

else

begin

p2 :=«pl^ .nextzip

;

if pi'', site O pi'', target then

begin

CanBe [POD] .next:=»pt;

pt'' .nextzip :=pl; • end

else

begin

pt'" . nextzip :=p2 ; pl^ . nextzip : =pt

;

end;

end;

end;

end;

close(pair_file)

;

{ LST sends output to default list device (printer) )

{ Remove braces to change from console to printer }

vrriteln( { 1st , )

'

' ) ; writeln( ( LST ,

}

' Total number of zipcodes is ' , count : 5

, '

. ' )

;

writeln( { LST ,

)

'

Number of possible POD's is '

, Nposs : 5

, '. ' )

;

writeln( [LST, )

' Number of existing POD’s is '

, CurrentNximber : 5 , '

. ' ) ;

{ dumpstruct; ) { diagnostic only--prints out entire data structure ) end;

{

) procedure InsertPOD ( new_POD: Zcode );

{ this procedure performs the actual addition of pod site 'new_POD'

33
and updates all pertinent data in arrays such as GurrentCost [ ]

BestPOD[], NextBestPOD[ ] , NextCost[], and the set CurrentPODs .

}

var

base

opt, old

new_cost

begin

opt := Index [new_POD] ; base := CanBe [opt] .next;

new_cost :=• base^.cost; { include operating cost } CurrentPODs :=• CurrentPODs + [opt] ; base :=• base^ . nextzip

while base O nil do

begin

with base'' do

begin

zip site;

old_POD := BestPOD[zip]

if ( cost < CurrentCost[zip] ) and ( old_POD O zip ) then

begin

NextCost [zip] := CurrentCost [ zip ]

CurrentCost [zip] :=» cost;

NextBestPOD[zip] := old_POD;

BestPOD[zip] :=• new_POD;

end

else if ( cost < NextCost [zip] ) then

begin

NextCost [zip] ;= cost;

NextBestPOD[zip] := new_POD;

end;

end; { with base'' do. . . } base; “base'' .nextzip;

end;

NextBestPOD [new_POD] := BestPOD[new_POD]

NextCost [new_POD] := CurrentCost [new_POD]

CurrentCost [new_POD] := new_cost;

BestPOD[new_POD] := new_POD;

end; { procedure InsertPOD )

{

procedure RemovePOD( old_POD: integer);

{ this procedure performs the actual removal of pod site 'old'

and updates all pertinent data in arrays such as CurrentCost []

BestPOD[], NextBestPOD [] , NextCost[], and the set CurrentPODs.)

var

base, base2 ; link;

id, idold : Zcode

: link;

pod, zip : Zcode;

: real;

34
looking : boolean;

Begin

idold := Index [old_POD]

;

CurrentPODs ;= CurrentPODs - [idold];

base := CanBs [ idold] .next

;

while baseOnil do

begin

with base'' do

begin

if NextBestPOD [ site ] = old_POD then { NextBestPOD corrected here } begin

base2 :=• map[site];

looking := true;

while looking do

begin

id := base2' . target

;

if (Index[id] in CurrentPODs) and (id O BestPOD [ site ]) then

begin

looking :=• false;

NextBestPOD [ site ] := base2'' . target

;

NextCost [site] ;= base2^.cost;

end

else if base2=nil then

begin

looking :=» false;

NextBestPOD [site] :=• 0;

end

else

base2 := base2^ .nextpod;

end;

end;

if BestPOD [site] - old_P0D then

begin

BestPOD [ site] NextBestPOD] site]

;

CurrentCost [site] := NextCost[site] ; base2 niap[site];

looking!” true;

while looking do

begin

id base2^ . target

;

if (Index [id] in CurrentPODs) and (id O NextBestPOD [ site ] ) then

begin

looking false;

NextBestPOD [ site ] base2^ . target

;

NextCost [ site ] base2'' . cost

;

end

else if base2-nil then

begin

looking false;

NextBestPOD [ site ] 0;

end

else

base2 base2'' , nextpod;

35
}

end;

end;

end;

base := base^ .nextzip

end;

end; ( procedure RemovePOD() }

{

procedure GreedyADD;

{ Locates one new POD by greedy heuristic } var

Zip_Code

opt, zip, ind,

new_POD, i, old_POD, z impr, addval

base

done

string[5]

Zcode

real

link;

boolean;

{ Procedure GreedyADD )

{ find the index of the best choice )

begin

opt := 1;

impr := lE+38;

for ind := opt to Nposs do

begin

zip :=• CanBe [ ind] .where

if BestPOD[zip] O zip then ( zip isn't a POD site at the moment } begin

base :=• CanBe [ ind] . next

addval ;= base''. cost - CurrentCost[zip] ; { Operating Cost ) base := base^ .nextzip

while base O nil do

with base'' do

begin

z := site;

if (cost < CurrentCost [ z ] ) and (z O BestPOD[z]) then

addval := addval - currentcost [ z ] + cost;

base ;=• nextzip;

end;

if addval < impr then { this is the best choice so far ) begin

impr :=• addval;

opt := zip;

end;

end; { if BestPOD, . . } end; { for ind ... } if impr < lE+38 then { add the new POD to the current set ) begin

InsertPOD(opt)

str (ZCreal [ opt ] : 5 : 0 , Zip_Code)

writeln( {LST, }' Adding ',Zip_Code,' to active set of PODs .

' )

end

I

II

I

I

I

!l

!l

I

I

I

I

I

I

36
else

begin

stuck := true;

writeln(' No POD can be added to the current allocation.' );

end;

end; { procedure GreedyADD }

procedure GreedyDEL;

{ This procedure subtracts a POD from the currently assigned set

according to the greedy heuristic.)

var

Zip_Code

base

tset

minval , change

i, opt, ipod

String[5]

Link

ZipSet

Real

Zcode

begin

opt := 0;

minval:- lE+32;

for i := 1 to Nposs do

begin

if ( i in CurrentPODs ) and not ( Canbe [ i ] .must ) then

begin {ith site is a candidate for deletion } base Canbe [ i ]. next

ipod base^.site;

change : - 0

while base O nil do

with base'' do

if BestPOD[site] - ipod then

if NextBestPOD[site] O 0 then

begin

change change + NextCost [ site ] - CurrentCost [ site

base base^ .nextzip

end

else

begin

base nil;

change lE+32;

end

else

base base^ . nextzip

if ( change < minval ) then

begin

opt ipod;

minval : - change

end;

end; { if i in CurrentPODs ... ) end; { for .... ) if opt O 0 then { we're not stuck: delete opt from CurrentPODs )

37
begin

str(ZCreal[opt] : 5 : 0 , Zip_Code)

writeln( ' Deleting ',Zip_Code,' from active set of PODs');

Remove POD (opt)

end

else

begin

stuck := true;

writeln(' No POD can legally be deleted from current allocation.' );

end;

end; {procedure GreedyDEL}

procedure greedy;

( This procedure, given an initial and a final number of POD sites, adds

or subtracts sites using the greedy heuristic until the desired number

remain. Procedures GreedyADD and GreedyDEL are called.

{

) begin { MAIN PROCEDURE }

{ writeln(' Entering Greedy heuristic..,');}

stuck := false;

changes := false;

while (CurrentNumber O EndNumber) and not stuck do

if (CurrentNumber < EndNiomber) then

begin

GreedyADD;

ComputeCost

ListCurrent

changes := true;

if not stuck then

CurrentNumber := succ( CurrentNumber );

end

else

begin

GreedyDEL;

ComputeCost

ListCurrent;

changes :=• true;

if not stuck then

CurrentNumber :=• pred( Currentnumber );

end;

if not changes then

writeln(' No changes in greedy heuristic.');

end; { Procedure Greedy ]

procedure interchange;

{This procedure performs the Teitz/Bart interchange heuristic for the

38
Simple Plant Location Problem. Given an initial allocation of

customers to service sites, the heuristic checks to see if it would be

advantageous to exchange one currently assigned service site for one potential service site not currently assigned. The best such exchange

is performed, and the heuristic repeats 'ontil no advantageous

exchanges exist. }

done

POD , TestOut , Swapin, SwapOut ,

mincost,val

Zip_Code_out , Zip_Code_in

•.boolean;

: integer

: real

:string[5]

( )

function SwapVal(old_POD,new_POD; integer) :real;

{This function computes the value of a potential site-exchange of site

'new' for site 'old'. )

const

failure-100 . 0

var

contr

illegal , looking

idold, idnew

base, base2

real

boolean;

Zcode

link;

begin

contr := 0.0;

illegal false;

idold Index [old_P0D]

idnew Index [ new_P0D ]

if CanBe [ idold] .must then ( permanent POD Site ) illegal true

else

begin

base : -canbe [ idnew] . next

contr :-base^ . cost-CurrentCost [new_P0D] ; { Operating Cost for new_P0D 1 base ; -base'' . nextzip ; while baseOnil do

begin

with base'' do

if (cost<currentcost[site] ) then ( make new assignment ] if (siteOBestPOD [ site ] ) and (old_PODOBestPOD [ site ] ) then

contr contr + coct - CurrentCost [ site]

base base^ .nextzip

end;

base ; - CanBe [ idold] . next

while base O nil do

begin

with base^ do

if BestPOD [ site ] - old_P0D then

39
begin

base2 := Map[site];

if site =• new_pod then

looking :=» false

else

looking :=• true;

while looking do

begin

looking := false;

if base2=nil then

illegal ;=« true

else if ( Index [base2^ . target ] in CurrentPODs ) and

(base2^ . target O old_POD) then

contr ;=• contr + base2''.cost - CurrentCost [ site

else if (base2'' . target = new_POD) then

contr := contr + base2''.cost - CurrentCost [site]

else

begin

looking : = true

base2 :=• base2^ . nextpod;

end;

end;

end;

base := base'' .nextzip';

end;

end; [else clause from top)

if illegal then

Swapval := failure

else

Swapval ;= contr;

end; { function swapval () }

I }

{ MAIN PROCEDURE } begin

{ writeln (

Entering interchange heuristic.

changes := false;

repeat

done ; = true

for POD := 1 to nposs do

begin

if not (POD in CurrentPODs) then { POD is a candidate: } begin

Swapin := CanBe [ POD] .where

SwapOut : =“ 0

mincost := 0.0;

for i:=» 1 to Nposs do

begin

if i in CurrentPODs then

begin

TestOut ;= canbe [ i] .where

il

ill

I

i

I

I

i

I

I

I

I

I

I

I

P

40
val := swapval (TestOut , Swapin)

if vaKmincost then { this is the best swap so far } begin

mincost val;

SwapOut :=" TestOut;

end;

end;

end;

if mincost<0.0 then { go ahead and make the best swap: } begin

str( ZCreal [ SwapOut] : 5 ; 0 , Zip_Code_out)

str( ZCreal [ Swapin ]:5:0, Zip_Code_in );

writeln( {LST, }' Swapping ', Zip_Gode_out ,' out'

,

Zip_Code_in, ' in. ' )

InsertPOD( Swapin)

Remove?OD( SwapOut)

changes : = true;

done :=false;

end;

end;

end;

until done

{ writeln( ' Leaving interchange heuristic:'); ) end; { interchange heuristic procedure )

{-- - -

) procedure PODColor;

{ This procedure creates a solution file which does not use the

graph-coloring algorithm. )

var

zip, clr, pod : integer;

WriteFile : text;

base : link;

begin

Assign (WriteFile , ' STATE' +StateNumber+' . SOL' )

Rewrite(WriteFile)

write (Write file , totalcost : 12 : 2)

writeln(Writefile)

clr : =0

for pod:-l to Nzips do

if Index [pod] O 0 then

begin

base :-CanBe [ Index [pod] ] .next;

while baseOnil do

begin

zip : - base^ . site

if BestPOD [ zip ] -pod then

begin

write (WriteFile , zip : 5 , ZcReal [zip] :6:0,clr:3, pod: 5 )

writeln(WriteFile)

end

41
base : =*base^ . nextzip

end;

end;

Close (Writefile)

end; { PODColor )

{ -

procedure GraphColor;

{ This procedure computes the adjacency of POD service -regions in

the current solution to the POD location problem, and colors the

zips in these regions such that no two adjacent regions use the same color. At most six (five?) colors will be used. For a good

description of the coloring algorithm, see David W. Matula et al,

"Graph Coloring Algorithms", in Ronald C. Read, "Graph Theory and

Computing", 1972 Academic Press, N.Y. The idea for the algorithm

is based on the 'two-color chain' proof of the five-color theorem.

The solution may be saved to a file, if desired. }

type

ptr_type = ''adj_list_el

adj_list_el = record

V : integer;

next : ptr_type

end;

graph_type - array [1. .MaxPossible] of adj_list_el;

Node_array - array [ 1 . .MaxPossible] of integer;

set_type = set cf 1 . .Maxpossible

var

graph : graph_type

{ contains adjacency list representation of the graph }

Node_num, color, ordering : Node_array;

{ color [vi] is the color assigned to vertex vi,

ordering[] stores the order in which vertices should be

colored. Node_num[i] tells which vertex in Graph corresponds to

zip site number i }

last_color, num_nodes , count, pod, spot ; integer;

{ total # of colors used, number of nodes in graph

count , pod, spot are temporary variables ]

answer : string[3];

filename : FileString;

WriteFile : text;

GettingName : boolean;

base : link;

{

)

42
procedure init_graph;

var

nodes, count, w : integer;

begin

for count 1 to CurrentNumber do

graph [count] .next := NIL;

nodes := 0;

for count := 1 to Nposs do

if count in CurrentPODs then

begin

nodes := succ (nodes);

w := CanBe[ count] .where;

graph [nodes] .v :=» w; Node_num[w] ;= nodes;

end;

end; { init_graph }

procedure write_out_graph (var graph : graph_type; num_nodes : var

count : integer;

temp : ptr_type;

begin { write_out_graph ) writeln; writeln;

for count := 1 to n’am_nodes do begin

write (' adj acency list for node ', count,' is : ');

temp graph [count] .next;

while (temp O nil) do begin

if temp O nil then write ( temp'' .v) ; temp temp ''.next;

if temp O nil then write(',');

end; { while ) writeln;

end; ( for ) end; { write_out_graph )

{

) procedure read_in_graph

var adj num, zipnum, neighbor , count : integer;

adj_file : text;

{ ) procedure Add_to_list(z ,n: integer)

var podl , pod2 : integer;

ptr, p ; ptr_type;

begin

podl Node_num[ BestPOD [ z ]]

pod2 Node_num[ Bes tPOD [ n] ]

ptr nil;

integer)

43
new(ptr)

ptr^ .V := podl;

ptr''.next :=• graph [pod2] .next;

graph [pod2] .next := ptr;

ptr := nil;

new (ptr)

ptr'' .V := pod2

ptr^.next := graph [podl] .next;

graph [podl] .next :=• ptr;

end; { Add to_list )

{ -

)

begin

assign(adj_file ,

' STATE '+StateNumber+' .ADJ ' ) ; reset (adj_file)

while not EOF(adj_file) do

begin

read(adj_file , zipnum, adjnum)

readln(adj_file)

for count := 1 to adjnxim do

begin

read(adj_file .neighbor)

if BestPOD[zipnum]OBestPOD[neighbor] then

Add_to_list( zipnum, neighbor)

end; { for ) readln(adj_file)

end; ( while ) end; { read_in_graph )

{ -

} procedure delete_node ( node : ptr_type; var list_ptr : ptr_type)

var

temp ; ptr_type

begin

temp ;=• list_ptr;

if (node - list_ptr) then begin

writeln(' error');

list_ptr node^.next;

temp : = node

dispose (temp)

end [if]

else begin

while (temp''. next O node) do

temp ;= temp ''.next;

temp next := node ''.next;

temp : - node

dispose (temp)

end; [else]

end; { delete node )

{

)

44
procedure clean_up (""^ar graph : graph_type ; num_nodes ; integer) ;

{ eliminates duplications from the adjacency list of each vertex } var

node , temp_node : ptr_type;

index , current : integer;

adjacent : set_type;

begin

for index := 1 to num__nodes do begin

adj acent ; = [ ] ; node := graph [ index] .next

while (node O nil) do begin

current := node'^.v;

if (current IN adjacent ) then begin

temp_node ;= node;

node := node ''.next;

delete_node ( temp_node , graph [ index] . next)

end ( if ) else begin

adjacent := adjacent + [current]

node ;= node ''.next;

end; { else } end; { while } end; { for } end; { clean_up )

{ } procedure find_min_degree (var vertex : integer; var graph : graph_type

num_nodes : integer; var deleted : set_type)

var

v_count, degree ,min_degree : integer;

temp : ptr_t}rpe;

begin

min_degree MaxPossible;

for v_count := 1 to ni.im_nodes do

if not (v_count IN deleted) then begin

temp graph [v_count] .next;

degree := 0;

while (temp O NIL) do begin

if not (temp''.v IN deleted) then

degree degree + 1;

temp temp''. next;

end; { while ) if (degree < min_degree) then begin

vertex v_count;

min_degree ; - degree;

end; { if } end; { if ) end; { find min degree )

{ )

45
procedure order_graph (var graph : graph_type; num_nodes : integer;

var ordering : Node_array)

var

deleted : set_type;

count, vertex : integer;

begin

deleted : [ ] ; for count ;=• num_nodes downto 1 do begin

f ind_min_degree (vertex, graph, num_nodes , deleted)

ordering [ count ] := vertex;

deleted := deleted + [vertex]

end; { for ) end; { order_graph )

{ -

} procedure write_colors (var color : Node_array; last_color , num_nodes

: integer)

var count ; integer;

begin

writeln; writeln;

for count :=» 1 to num_nodes do

writeln( ' POD '

, graph[count] . v, ' is colored in color

' ,color[count] )

writeln(' This coloring used '

, last_color ,' colors.');

end; { write_out_graph }

{ } procedure find_available( var color : Node_array; v_point : ptr_type;

vertex : integer; var first_not_used : integer);

var

temp : ptr_type;

begin

temp := v_point;

f irst_not_used := 1;

while (temp O NIL) do

if ( color [ temp^ .v] =• first_not_used)and( temp^ , v O vertex)

begin

first_not_used first_not_used + 1;

temp :=• v_point;

end { if . . ) else temp := temp^.next;

end; { find_available }

{ }

procedure determine_colors_used(point : ptr_type;

last_color : integer; var color, used : Node_array)

var

temp : ptr_type;

current : integer;

begin

for current := 1 to last color do

then

46
used [current] := 0;

temp :=» point;

while (temp O NIL) do begin

current ;= color [temp'' .v] ; if (used [ current ] = 0) then

used[current] ;= temp^.v

else if (used [current] > 0) then

used [current] := -1;

temp := temp^.next;

end; { while }

{ writeln( 'determining the colors used out of '

, last_color)

for current := 1 to last_color do

writeln( current , ' : ' ,used [current] ) ; ) end; { determine_colors_used }

{ }

procedure change_colors ( var mark, color : Node_array;

colorl , color2 , num_marked : integer);

var

index , countl : integer;

begin

for countl :=• 1 to num_marked-l do begin

index := mark [countl]

if (color [ index] = colorl) then

color [index] := color2

else if (color [ index] =• color2) then

color[index] := colorl

else writeln( ' ERROR in change_colors , ignore the results')

end; { for } end; ( change_colors }

{ ) procedure try_swap (a,b : Integer; var graph : graph_type ; var color

Node_array; var success ; boolean)

var

visited, sub_colors ; set_type;

mark : Node_array

num_marked , current , colorl , color2 ; integer;

{

) procedure find_component (var current : integer)

var

temp : ptr_type;

begin

{ writeln( 'visited vertex '.current,' '.success);}

if (current - b) then success false

else begin

visited visited + [current]

mark [ num_marked] current;

num_marked num_marked + I; temp graph [ current ]. next

47
while (temp O NIL) and (success) do begin

current temp^.v;

{ writeln( ' should we visit '.current,' ?'); ) if (not (current IN visited) )and

(color [current] IN sub_colors) then begin

find_component(current)

{ writeln( ' popped up'); } end;

temp := temp^.next;

end; { while ) end; { else } end; { find_component

{ } begin { try_swap } success := true;

num_marked : = 1

current a; visited : =>

[ ] ; colorl color[a]

color2 color[bj;

sub_colors :» [colorl] + [color2] ; find_component( current)

if (success) then change_colors (mark, color , colorl , color2 , num_marked)

end; { try_swap }

{

procedure try_interchange (var graph : graph_type; vertex ; integer;

var color : Node_array; var first , last_color : integer);

var

used : Node_array;

countl , count2 , trial_color ; integer;

change_successful : boolean;

begin

( writeln ('trying interchange for #' .vertex, ' ..');}

determine_colors_used(graph[vertex] .next, last_color , color .used)

change_successful false;

countl := 1;

while (countl <- last_color) and(not change_successful) do begin

if (used [ countl ] >0) then begin

count2 := countl + 1;

while (count2 <= last_color) and (not change_successful) do

begin

if (used[count2j > 0) then begin

{ writeln( ' considering '.countl,', ',count2);}

trial_color := color [countl] ; try_swap(used[ countl] ,used[count2] , graph , color

change_successful)

if (change_successful) then first := trial_color;

end; { if . . } count2 ;= count2 + 1;

end; { while count2 ... )

48
end; { if . . countl } countl := countl + 1;

end; { while ] end; { try_interchange }

{ )

procedure color_graph (var graph ; graph_type; nura_nodes : integer;

var color : Node_array; var last_colcr : integer)

var

vertex, first_not_used, counter : integer;

begin

for counter ;= 1 to n\im_nodes do

color [counter ] :=» 0;

last_color := 1;

for counter ;= 1 to num_nodes do begin

vertex := ordering [ counter ]

{ writeln( ' now coloring vertex number vertex) ; } find_avai lab le (color , graph [vertex] . next .vertex , first_not_used)

if (f irst_not_used > last_color) then

try_interchange(graph,vertex, color, first_not_used, last_color)

color [vertex] ;= f irst_not_used;

if (f irst_not_used > last_color) then

last_color := last_color + 1;

end; { for } end; [ color_graph }

{ - MMN PROGRAM }

begin { main ) init_graph;

n;am_nodes CurrentNumber ; ( CurrentNumber is a global variable

telling how many PODs are assigned ] read_in_graph

c 1 e an_up ( graph , num_node s )

( write_out_graph(graph,num_nodes) ; } {diagnostic only)

order_graph( graph , num_nodes , ordering)

CO lor_graph( graph , n\om_nodes , color , last_color) ;

{ write_colors (color , last_color , num_nodes ) ;

)

{ writeln(' Do you wish to save this solution and coloring on the disk');

writeln( ' for later graphic output ? (Yes or No)');

readln(answer) ; ] answer [1] :=-'y' ; if (answer [l]-'y' ) or (answer [ 1]->'Y' ) then

begin

{ writeln( ' Enter the filename under which you wish to save the data:')

readln( filename) ; ) filename ' STATE '+S tat eNumber+' , SOL'

GettingName exist(filename)

while GettingName do

begin

writeln( ' NOTE: file '.filename,' already exists:');

writeln( ' Write over this file ? ');

readln(answer)

49
if (answer [l]0'y' ) and (answer [1]0'Y' ) then

begin

writeln( ' Enter new filename:');

readln( filename)

GettingName ;= exist (filename)

end

else

GettingName := false;

end; { while } Assign (WriteFile , filename)

Rewrite (WriteFile)

write (Writefile , totalcost : 12 : 2)

writeln(Writefile)

for pod:=l to Nzips do

if Index [pod] O 0 then

begin

base :=»CanBe [ Index [pod] ] .next;

while baseOnil do

begin

count:— base^.site;

if BestPOD[ count ]=pod then

begin

spot := Node_num[pod]

write (Writefile , count : 5 , ZcReal [count ] : 6 : 0 , Color [ spot] : 3 , pod: 5)

writeln(WriteFile)

end;

base: -base .nextzip;

end;

end;

Close(Writefile)

{ writeln( ' Data have been saved in file '.filename);}

end;

end; { GraphColor )

{ )

Procedure Lagrangian_dual ; { (var Dual_var : ValueArray ; mono_Xij : Indexarray) ; )

{ This procedure attempts to locate a lower bound on the optimal

allocation. )

Type

RealArray = array[l. .Maxpossible] of real;

IntegerArray = array [ 1 . .Maxpossible ] of integer;

Const

w_eps =• 0.001;

var

Hold_POD_set

Running_Average

red_cost , work

pod_id

Dual var

zipset

Array [1.. 5] of real;

RealArray;

IntegerArray;

ValueArray;

50
Sum_Xij : ValueArray ; {Indexarray;

}

mono_Xij: Indexarray

multiplier, mean, delta_s, delta_u,

scale_factor , norm_factor, min_dif,

w_target, w_previous , w_new, w, s : real;

Sum_Xj j , monotone, nits, iter,

pod, ipod, zip, i, tick

error. Existing

base

dual_file

filename

Out file

integer

boolean;

1 ink

file of real

Filestring;

Text

procedure Assignment;

{Given the contents of CurrentPODs and the arrays of neighbor data,

this procedure determines the nearest currently assigned POD for

each individual zip-code area, and the associated costs. } var

base

zip , pod

empty, done

ipod, izip

link;

zcode

boolean;

integer

begin

error false;

TotalCost :=* 0.0;

for zip :=“ 1 to Nzips do { find the first current POD in zip's list of

possible pod's, and assign zip to it.

begin

done false;

base map [zip]

if base-nil then

done true

while not done do

if base-nil then { this zipcode will be skipped } begin

done : - true

error true;

writeln(' Feasibility error at ',zip:5);

end

else

begin

pod base'' . target

ipod Index [pod]

if ipod in CurrentPODs then { pod is the best choice: ) begin

done : - true

BestPOD[zip] pod;

CurrentCost[zip] base^.cost;

end ( if POD in CurrentPODs . . .

)

else

base base^ . nextpod ; ( keep looking for a best POD )

51
end; {while not done. . .

end; { for zip 1 to . . .

}

end; { Procedure Assignment )

Procedure Quicksort (Var value :RealArray;Var index: IntegerArray;N: Integer)

Procedure Exchange(I , J : Integer);

{ Change records I and J )

var

temp : real

indx: integer;

Begin

temp : =value [ i]

indx;=index[ i]

value [ i ]

: =value [ j ]

index [ i ]

: =« index [ j ] ; value [ j ]

; =t emp

index [j ]

: =indx;

End;

Const

MaxStack =20; { Log2(N) - MaxStack, i. e. for MaxStack = 20

it is possible to sort 1 million records ) Var

t The stacks } LStack : Array [ 1 . .MaxStack] Of Integer; ( Stack of left index } RStack : Array [1 . .MaxStack] Of Integer; { Stack of right index ) Sp : Integer; { Stack SortPointer ) M,L,R,I,J : Integer;

X : Real

Begin

{ The quicksort algorithm } If N>0 Then

Begin

LStack[l] :=1;

RStack[l] :=N;

Sp:=l

End

Else

Sp:=0;

While Sp>0 do

Begin

{ Pop(L,R) } L :=LStack[Sp]

R :=RStack[Sp]

Sp : =Sp- 1

Repeat

I:=L; J:=R;

M:=(I+J) shr 1;

X:=Value[M]

{writeln('l r m x '

, 1 : 5 , r : 5 ,m: 5 , x) ;

52
R.epeat

while (K=J) and (Value [I] < x) do

I: =1+1;

while (K=J) and (Value[J] > x) do

J:=J-1;

If K=J Then

begin

{writeln('i j v[i] v[j] ' , i : 5 , j : 5 , value [ i] , value [ j ] ) ;

}

If iOj then Exchange(I , J) ; i : =i+l

j

end

Until I>J;

{ Push longest interval on stack ) If J-L < R-I Then

Begin

If I<R Then

Begin

{ Push(I,R) ) Sp:=Sp+l;

LStack[Sp] ;=I

RStack[Sp] :=R;

{writeln('sp i r ' , sp: 5 , i: 5 , r ; 5) ;

End;

R:=I-1

End

Else

Begin

If L<J Then

Begin

{ Push(L,J) ] Sp;=Sp+l;

LStack[Sp] :=L;

RStack[Sp] ;-J

{writeln('sp Ij ',sp:5,l:5,j;5);}

End;

L;-J+l

End;

Until L>-R

End;

End { Quicksort } ;

begin

( Initialize Lagrangian Solution ) tick 0;

nits ;= 50;

Error false;

scale_factor : — 1;

w_previous 0;

w_new : - 2+w_eps

w_target totalCost+10

hold_POD_set : - CurrentPODs;

53
{ Compute a dual value for the number of open POD sites )

s : -totalcost

GreedyADD

ComputeCos t

s : —to talcos t- s

CurrentPODs ; =hold_POD_set

match;

s:--1035;

{ Perform File Initialization ) assign(dual_file , ' dual . var ' )

reset (dual_file)

assign(out_file ,' out . fil ' )

rewrite (out_file)

writeln(Out_file , scale_factor , s)

{ Initialize dual variables using best and nextbest costs }

for i;-l to nzips do

if map[i]Onil then

begin { determine an interval for the dual variable ) if nextbestPOD[ i] = 0 then

nextcos t [ i ]

: =2'*currentcost [ i ] else if currentcost [ i] > nextcost[i] then

currentcos t [ i ]

; -nextcos t [ i ] /2

{ Estimate an initial dual value ) w := CurrentCost[i]/nextcost[i]

nextcost[i] nextcost[i] -currentcost[i]

if nextcost[i] > abs(s) then nextcost[i] :=abs(s)

dual_var [ i] :=CurrentCost [ i] + w*NextCost [ i]

mono_Xij [ i ]

: =0

end

else

dual_var [ i ]

; =0 . 0

{ read(dual_file , scale_factor , s) ;

for i;=l to nzips do

begin

{ read(dual_file , dual_var [ i] ) ;

{ if index[i] O 0 then]

writeln(Out_file , '

'

, dual_var [ i] : 7 : 1

, ' NextCost [ i ] ; 7 : 1

,

CurrentCost [ i] : 7 : 1

, ' ' , index[ i] : 3

, '

'

, i ; 4 , BestPOD [ i ] : 5)

end;

for i;=l to 5 do

Running_average [ i ]

: -1 . 0

{ Begin the main loop }

while abs (w_previous-w_new) > w_eps do

begin

clrscr

monotone ; =0

for iter:—1 to nits do

54
Compute Reduced Costs ] be sin

tick := tick + 1;

if iter = 1 then

for i:=l to Nposs do

begin { Compvite the reduced costs from scratch } bas e : = CanBe [ i ] . next

ipod := base^.site;

work[i] ;= base'' . cost-dual_var[ ipod] -s ; base := base'' . nextzip ; while base O nil do

with base'' do

begin

zip : =site

if cost-dual_var [ zip ] < 0 then

work[i] :=>work[ i]+cost-dual_var [zip] ; base ; =base'' . nextzip

end;

dual_var

end

else { Compute the same thing only faster }

i]

begin

for i:= 1 to nposs do

work[i] :=work[i]+delta_s*(Sum_Xj j -Endnumber)

for i:= 1 to nzips do

if Sum_Xij[i] O 0 then { Examine only those which changed ) begin

dual_var [ i ]

; =• dual_var [ i ] +Sum_Xij [ i ] ; base : =* map [ i ] ; while base O nil do

with base'' do

begin

ipod:- Index [ target ]

if site-target then

work[ipod] work[ ipod] -Sum_Xij [ i]

else if cost-dual_var [ i] < 0 then

begin

if cost-dual_var [ i

] +Sum_Xij [ i ] <0 then

work[ipod] work[ipod] - Sum_Xij[i]

else

work[ipod] work[ipod] + cost -

end

else if cos t-dual_var [ i

] +Sum_Xij [ i

] <0 then

work[ipod] work[ipod] - cost + dual_var[i]

- Sum_Xij [i]

base nextpod;

end;

end;

end;

{ find the best k POD set

for i:-l to Nposs do

55
begin

pod_id[ i] :=i

red_cost[i] :=work[i] ; end;

Quicksort(red_cost ,pod_id,nposs)

{ find a feasible solution }

Sum_Xj j ;=0;

w_previous := w_new;

w_new := s"*EndNumber

CurrentPODs : = [ ]

for i:= 1 to nposs do

if ( red_cost[i] < 0 ) then

begin

Sum_Xjj :=Suin_Xj j+1

w_new : =w_new+red_cos t [ i ]

CurrentPODs :=CurrentPODs+[pod_id[ i] ]

end;

{ Compute the Xij

'and Objective Value }

for i;=l to Nzips do

if map[i]=nil then

Sum_Xij [ i ]

: 0

else

Sum_Xij[i] :=•-!;

for i:=l to Nposs do

begin

base := CanBe [ i ] . next

pod := base^.site;

ipod := index[pod]

if ipod in currentPODs then

while base O nil do

with base^ do

begin

if (site-pod) or (cos t- dual_var [site] <0) then

Sum_Xij [site] Sum_Xij [ site ] +1

base :-nextzip

end;

end;

norm_factor : =0 . 0

for i:=l to nzips do

begin

w_new ; =w_new+dual_var [ i ]

{ norm_factor ; =norm_factor+Sum_Xij [ i

] *Sum_Xij [ i ] ;

}

no rm_factor : -norm_factor+abs ( Sum_Xij [ i ] )

end;

if norm_factor = 0 then norm_factor ;=0 . 9

Running_average [ iter mod 3+1] := norm_factor;

56
ni6S.ri I = (Running_average [ 1 ] +Rurming_average [ 2 ] +Running_average [ 3 ] ) /3

{ compute a new scale factor and compute new dual variables }

if w_new-v;_previous < -w_eps then

monotone : =0

else

monotone ; =monotone+l

if (monotone >= 5) and (scale_factor < 0.5) then

begin

scale_factor := 2 . 0*scale_factor

writeln(LST, ' 2*scale_factor ', scale_factor : 7:6);

monotone : =0

end;

multiplier : = scale_factor''f (w_target-w_new)/mean + v_eps;

{ compute the min_dif needed to change Sum_Xij by 1

} min_dif : =»lE+38

for i:=l to Nzips do

begin

if Sum_Xij[i] < 0 then

begin

if mono_Xij[ij > 0 then

mono_Xij [ i ]

: -0

else

mono_Xij [ i ]

: =mono_Xij

[ i

] - 1

base :=> map [ i] ; while base O nil do

with base^ do

begin

pod:- target;

if Index [pod] in currentPODs then

if (cost-dual_var [ i] >0) then

if (cost-dual_var[ i]<min_dif ) then

min_dif : -cost-dual_var [ i]

base : -nextpod

end;

end

else if Sum_Xij[i] > 0 then

begin

if mono_Xij

[ i ] <0 then

mono_Xij

[ i ]

: -0

else

mono_Xij [ i ]

: -mono_Xij [ i

] +1

base : - map [ i ]

while base O nil do

with base^ do

begin

pod:- target;

if Index[pod] in currentPODs then

if (cost-dual_var [ i ] < 0) then

57
if ( dual_var [ i

] - cost<min_dif) then

min_dif :-»dual_var [ i] - cost;

base :=»nextpod;

end;

end

else

mono_Xij [ i ]

: =0

end;

multiplier ;= scale_factor*(w_target-w_new)/mean + w_eps;

if (min_dif<lE+38) and (min_dif>multiplier) then

begin

multiplier := min_dif;

writeln(LST, ' min_dif applied ' ,min_dif : 10 ; 5)

end;

mean : =w_new/w_targe t

if (tick mod 5 = 0) then

writeln(LST , tick: 3

, '

' ,w_target : 12 , ' ',w_new:12,' ' ,norm_factor : 9

' multiplier : 9

, ' '.mean: 12,' ',s:6:l,' '

, Sum_Xj j : 3)

GotoXY(l.l)

writeln( ' Upgraded lower bound on the optimal solution value is '

w_new:7:0,' '.iter: 3);

writeln( ' iter w_target w_new norm_fact multiplier',

' min_dif s Sum_Xj j

' )

writeln(tick: 3

, '

' ,w_target : 12 , ' ',w_new:12,' ' ,norm_factor : 9

'

' .multiplier : 9

, ' ' ,min_dif : 9

, ' ',s:6:l,' '

, Sum_Xj j : 3)

for i:=»l to nzips do

begin

if abs(mono_Xij [ i] ) >- 5 then

begin

Sum_Xij [ i ]

: =»Sum_Xij [ i ] *2

writeln(' 2*Sum_Xij on ',i:4);

end;

if min_dif =» multiplier then

Sum_Xij [i] :=“-multiplier'*Sum_Xij [i]

else if multiplier < nextcost[i] then

Sum_Xij [i] :=-multiplier*Sum_Xij [i]

else

Sum_Xij [i] :=»-Sxam_Xij [ i]*scale_factor*Nextcost [ i] ; end;

if multiplier > 10 then

delta_s : =-scale_factor

else if (Sum_Xj j -EndNumber>0) and (multiplier<-red_cost[Sum_Xj j ]

then

delta_s:=abs(red_cost[Sum_Xj j ])

else if (Sum_Xj j -EndNumber<0) and (multiplier<red_cost [Sum_Xj j+1]

then

delta_s :=abs(red_cost[Sum_Xj j+1]

)

else

delta_s : “multiplier

s :=s-delta_s*(Sum_Xj j -EndNumber)

58
end;

writeln('i pod sumXij red_cost dual_var ' )

for i:=l to nposs do

begin

ipod:=pod_id[i]

pod :=canbe [ ipod] .where

writeln(i: 5

, ' ',ipod:5,' '

, sum_Xij [pod] : 5

,

red_cost[i] : 5:0, '

' , dual_var [ ipod] : 5:0)

end;

Assignment

if error then

begin

writeln(' current lagrangian solution not feasible ');

end’

else

begin

listcurrent

if (totalCost < w_target) and (Sum_Xjj = EndNumber) then

begin

mean:=(w_target-w_new) /( totalCost -w_new)

if mean > 2 then

begin

writeln(lst, ' scale_factor adjusted'

scale_factor : 5 : 4)

scale_factor : =scale_factor'’=mean;

if scale_factor > 1 then scale_factor :=1;

end;

w_target : =totalCost

hold_POD_set :=current?ODs

writeln( 'upgraded w_target' ,w_target)

end;

end;

icer:"=nits shr 1;

nits :=”iter+10;

scale_factor : -scale_factor/2 , 0+w_eps

{ writeln( ' enter nits scalefactor ' ,nits : 6

, '

'

, scale_factor : 8)

readln(nits , scale_factor )

writeln( ' w_new ' ,w_new, ' norm_factor '

, norm_factor )

writeln(' nits ',nits,' scale_factor ' , scale_factor ) ; ) reset (dual_file)

write (dual_file , scale_factor , s) ; for i:-l to nzips do

write (dual_file , dual_var [ i] )

end;

close(dual_file) ; close(out_file)

end;

59
NBS*n4A (REV. 2-dC) U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instruction s)

1. PUBLICATION OR

REPORT NO.

NBSIR 86-3472

2. Performing Organ. Report No. 3. Publication Date FEBRUARY 1987

4. TITLE and subtitle

The Internal Revenue Service Post-of-Duty Location Modeling System - Programmer's Manual for Pascal SOLVER

5. AUTHOR(S)

Paul D. Domich, Richard H. F. Jackson, Marjorie A. McClain, David M. Tate

6. performing organization (If joint or other than NBS, see instructions)

national bureau of standards

DEPARTMENT OF COMMERCE

WASHINGTON, O.C. 20234

7. Contract/Grant No.

S. Type of Report & Period Covered

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street, City. State. ZIP)

The Research Division

U. S. Internal Revenue Service

1201 E Street

Washington, DC 20224

10. SUPPLEMENTARY NOTES

I I Document describes a computer program; SF-185, FIPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bibliography or literature survey, mention it here)

This reports documents a project undertaken by the National Bureau of Standards to

develop a mathematical model which identifies optimal locations of Internal Revenue

Service Posts-of-Duty . The mathematical model used for this problem is the uncapacita

ted, fixed charge, facility location model which minimizes travel and facility costs,

given a specified level of activity. This reports discusses the mathematical techniques used to solve the mathematical model developed and includes a Greedy procedure,

an Interchange procedure, and a Lagrangian Relaxation of the related linear program.

A description of the Pascal routines, definitions of key data structures and

variables is provided. Data sources identified and used are also described.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolon S) uncapacitated fixed charge facility location problem. Greedy heuristic.

Interchange heuristic, Lagrangian Relaxation

13. AVAILABILITY

iX I Unlimited

I I For Official Distribution. Do Not Release to NTIS

I I Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.

20402.

(2^ Order From National Technical information Service (NTlS). Springfield. VA. 22161

14. NO. OF

PRINTED PAGES

62

15. Price

$13.95

jSCOmm-OC «0.J-*t0
IIIII1

fl