SECTION 17 MATHEMATICAL PROGRAMMING
SECTION 17.1 LINEAR PROGRAMMING
The problem to be solved is the following:
Maximize c'x
subject to
Ax <= b
x >= 0
where A is an m x n array, c is a vector of length n, x is a vector of length n, and b is a vector of length m.
SECTION 17.1.1 WHEN A BASIC FEASIBLE SOLUTION IS AVAILABLE (LINPROA).
It is assumed in the program LINPROA that the starting vector b is nonnegative, i.e., b>=0. Then the constraints can be written as
Ax + Iu = b
where I is an identity matrix and u is a (nonnegative) vector of m components representing disposal activities. The initial basic feasible solution can then be written as
u = b - Ax with x=0.
Call the routine as
CALL LINPROA(A,IT,IL,IR,IB,LT,LL,LR,LB,M,N,ITER,IER,KPIV,
1 ASAF,A1,A1I,A2,CSAF,BSAF,CI,BI,TMP,A1IA2,IFR)
where
A = array dimensioned dimension A(0:M+1,0:N+1).
Input values:
A(0,0) should be initialized to zero. The elements in (A(i,j), i=1,m, j=1,n, should be the elements of -A above. The elements A(0,j), j=1,n, should the elements of the C-vector above. The elements of A(i,0), i=1,m should be the elements of the B-vector above.
Output values:
The element A(0,0) contains the value of the objective
function. The elements of A(i,0), i=1,m, are the solution values.
The negatives of the elements in A(0,j), j=1,n, are the dual
solution (shadow prices).
IT = integer array dimensioned IT(0:n) (output)
IL = integer array dimensioned IL(0:m) (output)
IR = integer array dimensioned IR(0:m) (output)
IB = integer array dimensioned IB(0:n) (output)
LT = character*2 array dimensioned LT(0:n) (output)
LL = character*2 array dimensioned LL(0:m) (output)
LR = character*2 array dimensioned LR(0:m) (output)
IB = character*2 array dimensioned LB(0:n) (output) In the
output, the primal productive variables are denoted by X, the
primal disposal (slack) activities by U, the dual variables by Y
and the dual slacks by V. Variables are numbered and the solution
is displayed in the output (if IPT and/or JPT are not zero; see Section 1.1) as in the following example:
Optimum reached Objective function = 0.28365529E+02 Primal solution
x- 2 0.52734045E+01 u- 2 0.20054093E+01 x- 5 0.12432843E+02 x- 4 0.24299007E+00 x- 3 0.10416291E+02
Dual solution
v- 1 0.12031298E+01 y- 1 0.60383894E-01 y- 5 0.34624878E+00 y- 4 0.16450454E+01 y- 3 0.78487476E+00
The above arrays beginning with "L" contain the appropriate letters (x,u,y,v) and the arrays beginning with "I" contain the number of the appropriate variable. The identifier for a variable is thus given by the combination of "L"-type and "I"-type arrays. The second letters of these arrays describe the positions of the identifiers as arranged along the borders of the condensed simplex tableau: T for top, L for left, R for right, and B for bottom. Thus, when the optimum has been reached, LL and IL contain the identifiers for the primal solution and LB and IB those for the dual solution. Any variable not listed among these has, of course, zero value in the optimal solution.
M = number of constraints (input)
N = number of activities in the primal problem (input)
ITER = number of iterations taken (output)
IER = error return (output) = 0 normal return,
IER =-98 problem is infeasible,
IER =-99 problem is unbounded
KPIV = pivot choice criterion (input) =1 steepest ascent,
KPIV =2 maximum improvement at each iteration
ASAF = M by N scratch array
A1 = M by M scratch array
A1I = M by M scratch array
A2 = M by N scratch array
CSAF = scratch vector of length N
BSAF = scratch vector of length M
CI = scratch vector of length M+N
BY = scratch vector of length M
TMP = scratch vector of length N
A1IA2 = M by N scratch array
IFR = inversion frequency (input) = 0 means do not reinvert
matrix = positive integer means reinvert matrix every IFR
iterations
The total number of variables including slacks is M+N. Now consider these M+N variables partitioned into subvectors x(1) (of length M) and x(2) (of length N) at some iteration, let the set x(1) represent the basic variables, and let the A matrix be partitioned conformably into A(1), A(2). Then the model can be written as
A(1)x(1) + A(2)x(2) = b
or
x(1) = A(1)-1 b - A(1)-1 A(2)x(2)
The current elements in the condensed simplex tableau are given by A(1) -1A(2) and if IFR is, for example, equal to 3, then at every third iteration the elements of the simplex tableau are replaced by the results of the matrix inversion and multiplication indicated above. If SUBROUTINE VERS determines that A(1) is singular or if the product of the inverse of A(1) and of b does not yield a vector with nonnegative elements, the simplext tableau is not replaced.
SECTION 17.1.2 WHEN A BASIC FEASIBLE SOLUTION IS NOT AVAILABLE (LINPROFS) The problem to be solved is the same as in Section 17.1, but unlike Section 17.1.1, it is now assumed that some of the elements of the b vector are negative. Partition the rows of the constraints, the A and the b vector so that the constraints can be written
A(1)x <= b(1)
A(2)x <=-b(2)
where A(1) is K by N, A(2) is M-K by N, b(1) is a K-vector, b(2) is an M-K vector, and both b(1) and b(2) are nonnegative. (Note that the A(1) and A(2) in this section are not the same as in Section 17.1.1). We introduce the usual slack variables, also partitioned into u(1) and u(2), and write
A(1)x + u(1) = b(1)
A(2)x + u(2) =-b(2)
It is obvious that a basic solution obtained by setting x = 0 is not feasible. Hence we introduce m artificial slack variables in a vector w, also partitioned into w(1) and w(2), and write the constraints as
A(1)x + u(1) + w(1) = b(1)
A(2)x + u(2) - w(2) =-b(2)
Our first feasible solution is then obtained by setting x = 0 and u = 0. We now maximize the (artificial) objective function -1'w, where 1' is a row vector of ones. (Since w(1) = b(1) - A(1)x - u(1) and w(2) = b(2) + A(2)x + u(2), we can replace the objective function by one written in terms of x's and u's.) If a maximum is reached with an objective function value of zero, the w's have been driven out of the solution and we have obtained a basic feasible solution to the original problem. Continued maximization will then lead to the optimal solution of the original problem, if one exists. If the optimal solution to the artificial problem has an objective function value that is positive, no feasible solution exists for the original problem. (In practice, because of rounding problems, this is tested for by checking whether all the w variables have become nonbasic, rather than by checking the value of the objective function.) The primal and dual variables in the original problem are designated as in Section 17.1.1; in addition, the artificial slacks are denoted by w and slack variables in the dual constraints corresponding to the inequalities adjoined to the inequalities of the original dual problem, which come into being in the artificial first phase problem, are denoted by z. Only steepest ascent is available as a pivot choice criterion, and no periodic reinversion of the matrix is provided for. The routine is called as
CALL LINPROFS(A,IT,IL,IR,IB,LT,LL,LR,LB,M,N,ITER,IER)
where A
A = is the coefficient matrix dimensioned A(-1:M+1,0:M+N+1)
with A(i,j), i=0,M, j=0,N being set the same as in Section 17.1.1
(input), with the other rows and columns being set by the program
IT = integer array dimensioned it(-1:M+N) (output)
IL = integer array dimensioned il(-1:1+M) (output)
IR = integer array dimensioned ir(-1:1+M) (output)
IB = integer array dimensioned ib(-1:M+N) (output)
LT = character*2 array dimensioned LT(-1:M+N) (output)
LL = character*2 array dimensioned LL(-1:1+M) (output)
LR = character*2 array dimensioned LR(-1:1+M) (output)
LB = character*2 array dimensioned LB(-1:M+N) (output)
M = number of constraints (input)
N = number of productive variables in original problem (input)
ITER = iteration count (output)
IER = error return (output)
IER = 0 normal return,
IER =-98 problem is infeasible,
IER =-99 problem is unbounded
SECTION 17.2 KNAPSACK PROBLEM OF INTEGER PROGRAMMING (KNAP)
The knapsack problem solves the following:
maximize c(1)x(1) + c(2)x(2) + ... + c(n)x(n)
subject to a(1)x(1) + a(2)x(2) + ... + a(n)x(n) <= b
x(j) >= 0 and integer, j = 1, 2, ...,n
See R. S. Garfinkel and G. L. Nemhauser, Integer Programming, 1972, pp. 214-248. KNAP is called from the MAIN program by
COMMON/BPRINT/IPT,NFILE,NDIG,JPT,MFILE
. . . . . . .
CALL DFLT
CALL KNAP(F,P,C,A,B,N,X)
where
F = is an integer array dimensioned F(0:N,0:B)
P = is an integer array dimensioned P(0:N,0:B)
C = an integer array dimensioned C(N)
A = an integer array dimensioned A(N)
X = an integer array dimensioned X(N)
B = a scalar integer
N = a scalar integer If an optimum solution is obtained, it will
be contained in the vector X and the optimal value of the
objective function will be given by F(N,B).