A heuristic for finding embedded network structure in mathematical programmes

52

European Journalof Operational Research67 (1993)52-63 North-Holland

Theory and Methodology

A heuristic for finding embedded network structure in mathematical prograrnmes B a r r i e M. B a k e r a n d P.J. M a y e *

Coventry Polytechnic, Priory Street, Coventry CV1 5FB, England, UK ReceivedJuly 1990;revisedJanuary 1991 Abstract: An intuitive description is given of how to recognise when a mathematical programme can be converted to a network flow problem. In practice, it is more likely that a mathematical programme with a high degree of network structure will have additional complicating constraints. Thus, a heuristic is presented for identifying a maximal subset of the constraints which can be converted to a network flow problem. Computational results show that the heuristic is highly effective, and that much larger networks are identified by seeking a general network structure, rather than limiting the search to a subset of rows giving at most one + 1 and one - 1 in each column, as has been done previously. Keywords: Mathematical programming; Networks; Heuristics

L Introduction There is a growing interest in the problem of identifying network structure in mathematical programming problems. ?dgorithms which convert mathematical programmes into pure network flow problems, whenever such a conversion is possible, have been described by several authors, the most prominent being Bixby and Cunningham [2] and Iri [8-10]. More recent publications are due to Baston et al. [1] and B~by and Wagner [4]. The advantages of deriving such a network are: (i) greatly improved computational efficiency for finding an optimal solution to the problem, and (ii) a visual representation of the problem which may give improved insight. In practice, it is more likely that a mathematical programming problem will have a high degree of network structure, with additional complicating constraints which make a pure network representation impossible. Glover and Klingman [6] mention several publications concerned with solving linear programmes with special structure, including embedded network structure. However, there is little published work on the problem of finding embedded network structure in mathematical programmes. It is known that the problem of finding an embedded pure network of maximum row dimension in a mathematical programme is a member of the class of problems known as NP-hard (Wright [12]). There are many problems in this class, the travelling salesman problem probably being the most famous. It is * Now with the GovernmentStatistical Service. Correspondenceto: B.M. Baker, CoventryPolytechnic,PrioryStreet, CoventryCV1 5FB, England, UK. 0377-2217/93/$06.00 © 1993- ElsevierSciencePublishers B.V. All rights reserved

enerally believed that for all these xists for any problem in this class,

subject to

Cxjkni for i- 1,...,12, s, xjrO for j- l,...,?

(2) (3)

where IZ~represents the number of lecturers required during period i and Si represents the set of periods, j, s&h that a lecturer startin work at the beginning of period j would be lecturing durin period i. Table 2 shows how an initial simplex tableau could be constructed (temporarily i simple lower bound constraints x1 r 4 and x7 2 2). ss,. . . , sll are slack variables and the artificial variables CZ~, . . ., a,, comprise an initial basis. It is well known that in any minimum cost network flow problem, the arcs corresponding to a set of basic variables form a spanning tree of the network (ignoring direction). Any non-basic variable (arc) then forms a cycle with some basic arcs. A pivot would bring such an arc into the basis by sending a flow round this cycle which would alter the flow in one of the basic arcs to zero or to a lower or upper bound value (the variable leaving the basis). Thus, if a network representation of the formulation in Table 2 is possible, arcs ~~es~nding to us,, . . , ull must be drawn in the form of a tree.

Table 1 Time Minimum requirement

09-10

lo-11

11-12

12-13

13-14

14-15

15-16

16-17

17-18

18-19

19-20

20-Z

4

8

8

6

4

9

8

7

3

3

3

2

54

B.M. Baker, P.Z Maye / Embedded network structure heuristic in MP

Table 2 Simplex tableau for Braynepayne College XI

X2

!

1 1

1 1

1 1

X3

X4

X5

X6

X7

$2

$3

"'"

$10

$11

-1 I 1

1 1

a2

-1 1 1

1 1

1 1

! 1

1 1

a3

"'"

alO

all

1 1

1 1

1 1

-1 1

1 -1

1

8 8 6 4 9 8 7 3 3 3

Suppose now that x I is chosen as entering basic variable. Examination of the xt-column in Table 2 shows that as xt is increased from zero, each of a2, a 5 and a 6 must decrease to preserve the constraint equations. Therefore, the arc corresponding to xt must be drawn so as to form a cycle with the arcs ~br a2, a s and a 6 and must be directed against these arcs in the cycle. Similarly, x z must form a cycle with a2, a s, a 6 and aT, s 2 forms a circuit with a z and similar requirements are identified from the remaining columns. It is easy to confirm that the network in Figure 1 satisfies these requirements. T h e initial basic solution sets a 2 = 8, a 3 = 8 . . . . . art ~ 3. If these values are inserted as flows in the network o f Figure 1, flows from the source or to the sink can be d e d u c e d for each vertex. These flows are

x.~Z a2

~3

s6

s7

Figure 1

3

B.M. Baker, PJ. Maye / Embedded r~twod¢ st.mc•re he'aristic in MP

-8

x2

55 +8

3

Figure 2

shown in Figure 2 as positive values or negative values, rather than drawing source and sink vertices and their connecting arcs. Note also that the artificial arcs can be dropped at this stage. All that remains is to associate the simple lower bounds with ares x~ and xT, and the costs per unit flow from the objective function (1) with ares xt . . . . . x7 and then a minimum cost network flow package can be used to obtain an optimal solution to the problem. Note that the existence of a minimum cost network flow representation of the LP together with the integer values n i in constraints (2) guarantees the existence of an integer-valued optimal solution of the LP, assuming that a feasible solution exists. In general, it is seen that if all coefficients of an LP formulation are + 1, - 1 or 0, the problem of identit~ing network structure can be approached by first writing the LP in standard form; that is with a set of columns which form an identity matrix to provide an initial basis. A tree must be drawn with one edge for each of these basic variables. The tree must be drawn so that for each non-basic column, an edge can be drawn to form a cycle with those basic variables which correspond to the non-zero entries in that column. Directions must then be inserted so that for each coefficient of + 1, the non-basic variable is directed against the basic variable in this cycle, and vice versa for each coefficient of - 1. If there are a few non-zero coefficients other than + 1 or - 1, then it may be possible to use row and column sealing to convert the coefficients to + 1 or - 1 . An important special case arises when each non-basic column contains at most a + 1 and a - 1 and all other values are zero. In this case a network representation is always possible, with the basis tree taking the form of a star (i.e. there is a vertex on which every are of the tree is incident). The problem addressed in this paper is that of finding the largest possible subset of rows of the LP formulation for which a network representation can be constructed as just described. The heuristic of

B.M. Baker,PJ. Maye / Embeddednetworkstructureheuristicin MP

56

Bixby and Fourer is based on finding a subset of rows for the case where every column contains at most a + 1 and a - 1. It is clear that the problem addressed here will be computationally much more demanding than the problem addressed by Bixby and Fourer. The additional computational effort must be weighed against the advantages of finding much larger networks, as shown in the Section 5, Computational Results.

3. Problem formulation Given an LP in standard form with constraint-set

[A, I ] X = b,

(4)

x~o,

(5)

it is assumed that the ( m x n)-matrix A satisfies ali E {0, 1, - 1} for all i and j, row and column scaling possibly having been carried out to achieve this. For simplicity, rows and columns of A containing only one non-zero entry can be omitted during application of the heuristic, since the former can be represented on a network as simple lower bound or simple upper bound constraints, and the latter merely mirror or duplicate arcs of the initial basis tree. Now consider the problem of converting the constraint matrix [A, 1] to a node-arc incidence matrix for a graph with m + n arcs, using simple row operations (i.e. converting to a matrix in which each column contains one + 1, one - 1 and all other entries are zero). If premultiplication by a matrix, T say, carries out such row operations, so that [TA, T] is a node-arc incidence matrix, then it is seen that T is the node-arc incidence matrix of the arcs corresponding to the variables in the initial basis, which must form a tree with m + 1 vertices. It is also straightforward to verify that the non-basic arcs, with the properties already deduced, have node-arc incidence matrix given by TA. Formally, the problem is to find an (m + 1) × m-matrix, 7, which is the node-arc incidence matrix of a tree, so that TA is also a node=arc incidence matrix. If such a T exists, then (4) can be re-wrRten as

[TA, T I X = Tb,

(6)

so that the LP can be solved as a minimum cost network flow problem with node-arc incidence matrix [TA, T] and with source and sink flows given by Tb. If no such T exists, then we wish to partition the constraint set into 2 subsets, with coefficient matrices A l ( p X n) and A2 ((m - p ) × n), so that such a T, with conformable dimensions, does exist with respect to A~, and so that p is made as large as possible. Without loss of generality, assume that A 1 consists of the first p rows of A and the constraint set will then have the form [A,, Ip, OIX= bl,

(7)

[As, o, ~m_.]x= t,~,

(8)

X>__O,

(9)

which may then be viewed as a minimum cost network flow problem with node-arc incidence matrix [TA 1, T] and source and sink flows given by Tb l, which must be solved subject to the additional complicating constraints (8).

B.M. Baker, P.Z Maye / Embedded network str~cturebeurislic in MP

4. Outline of the heuristic Consider the ~bUowing example: 1 1

2

-1

-1

2

1

3 4 A-~

3 4 5 1 -1

6

7 -1

9 -1

1

1

1

-1

1 -1

1

s

1

-1

1

6

-1

7 8

8

1

-1

9

-11

10

1

1

1

1 -1

If a pure network representation for A is possible, then from column 1, arcs 1, 4 and 8 of the basis tree must form a chain, the order of appearance of these arcs in the chain being unknown at this stage. From column 9, arcs 1, 4, 9 and 10 must form a chain. It follows that arcs 1 and 4 must be adjacent in the basis tree. This observation is extended in the first step of the heuristic as follows. For every pair of rows in A, i and k say (i < k), check whether there is at least one column, j, with a o ~ 0 and akj ~ O. If such a column can be found, then form a list Lik ~- {r : arj ~ 0 and i ~ r ~ k}. Then, if any further column(s), j ' say, can be found with a o, ~ 0 and a k f ~ 0 and then if there exists r ~ Lik for which ar~, = 0, delete r from L,.k. Thus, in the case of rows 1 and 4 of A, initially Lt4 ~- {8}, and then consideration of column 9 reduces this to L~4 = ~. On completion of this process, arcs i and k in the basis tree can be separated only by arcs listed in L~k, and further, the arcs of L~k together with arcs i and k must form a chain. The lists obtained by this procedure are shown in Table 3. On the assumption that a pure network representation is possible, a variety of further steps can be used to determine more adjaceneies (rows i and k for which Lik ~ ~J) (Maye [11]). For example, Lla = ¢ and Ll4 ~ ~ together imply that arcs 3 and 4 can be separated at most by are 1. But La4 -- {6}, and so this could be replaced by L34 = ~. However, for non-pure network problems, such results could prove misleading and, in tests, the additional computation time for techniques such as this has not proved worthwhile.

Table 3 Adjacent arcs and arc chains in network representation

L.:

1 2 3 4 5 6 7 8 9

2 ~

3 ~

4 ~[ 6

5 ~ 1 1

6

7

4 3

8 4 1

9 8,9

7,9 6,9

9 ~ 1 I, 10 1,10 7 7,8

10 9 1, 9 1,9

6, 7 1

58

B.M. Baker, PJ. Maye / Embedded network structure heuristic in MP

t

9

2

5

(t~)

(i~i)

(i,,)

~

9

t

(t)

9

(v)

(v~)

b,

9

2

4

Figure3

The heuristic then chooses 2 arcs which are adjacent to initialise the basis tree, say arcs 1 and 2 (Figure 30)). Relative directions are inserted (say from column 2 of A, where a~2 = -a22, indicating opposing directions). If none of the defined lists satisfies Lik = ~, then a list of minimum cardinality is chosen to initialise the tree. Now L13, L~4, Lls and L~9 are all empty, indicating 2 possible positions for each of arcs 3, 4, 5 and 9 to be appended. Arc 5 can be placed uniquely, either by observing L25 or by referring to column 2 of A (Figure 30i)). The direction can be derived from column 2. Arc 3 can now be placed using L35 in conjunction with Ll3 (Figure 3(iii)), followed by arc 4 using L34 in conjunction with Ll4 (Figure 3(iv)). L1~ in combination with L29 now locates arc 9, L79 in conjunction with L57 locates arc 7 and LIa in conjunction with LTs or La9 locates arc 8 (Figures 3(v)-3(vii)). L36 in conjunction with L67 now shows that arc 6 cannot be appended consistently and so must be omitted. Finally, Ll,10 in conjunction with L3,1o locates arc 10 (Figure 3(viii)).

59

B.M. Baker, P.& Maye / Embedded network strucm~ beur~tic in MP

In general, there is no guarantee that any arc can be placed uniquely at each incremental stage. In such cases, the heuristic can proceed by choosing at random a position for an are which cannot be placed uniquely. Clearly, the number of arcs in the final network may depend on the pair of arcs chosen i n i t i a l . Hence the heuristic uses different starting combinations to generate several networks, from which the best result is chosen. No backtracking mechanism is incort~rated, as this would require an e x c e ~ c e amount of computer time. A full statement of the heuristic is given in Appendix A. In Appendix B, it is shown that for an m x n coefficient matrix, the heuristic has a worst case performance of O ( n m 4 ) . T h e heuristic of B£xby and Fourer, with its much simpler objective, has a worst ease performance of O(mn). As seen in the above example, the hearistie works mainly by identifying pairs of arcs which must be adjacent. Alternative forms of hidden information can be deduced and used to construct a network. For example, Gould [7] has produced an algorithm for pure network problems, that finds sets of ares which must form chains within the network, and from these sets, the order of arcs in these chains is deduced.

5. Computational results To test the heuristic a procedure was devised for generating linear programmes cvnvertible to pure network flow problems using a random number generator. This procedure is described in Appendix C. A number of rows in each problem were then corrupted so that a pure network structure would not be obtainable. However, it is then not possible to state the maximum number of rows for which conversion to a network is possible. Hence, the heuristic was applied to each problem both in its pure network form and its corrupted form. Table 4 shows the results for 10 such problems. For comparison, the heuristic of Bixby and Fourer was applied to each problem. The purpose of this comparison is to see whether there is any significant advantage to seeking the general network structure as opposed to seeking a subset of rows containing at most one + 1 and one - 1 in each column. The two heuristics are not comparable in any other sense. The results in Table 4 show that the heuristic presented here is very effective for finding embedded network rows. Further, there is a considerable advantage due to seeking the general network structure.

Table 4 Computational results for heuristic Problem

Rows

Columnsa

Rows/Corrupted

Non-zeros

1 2 3 4 5 6 7 8 9 10

110 120 130 140 150 160 170 180 190 200

500 500 600 600 700 700 800 800 900 990

11 15 24 21 31 40 51 55 67 77

2782 2671 2970 2994 3828 4036 5008 4804 4861 5592

Arcs in network b. Baker & Maye c 97 (110) 109 (120) 105 (130) 127 (140) 118 (150) 133 (160) 124 (170) 136 (180) 128 (190) 133 (200)

a Number of columns includingthe identity matrix. b Number of network rows found for the corrupted LP, and for the uncorrupted LP in brackets. c Best result using 5 different starting points. d Best result using 4 different starting points.

Bixby& Fourer a 51 (52) 63 (64) 65 (69) 73 (73) 68 (74) 77 (76) 84 (90) 97 (97) 99 (99) 96 (95)

60

B.M. Baker, P.Z Maye / Embedded network structure heuristic in MP

6. Conclusion For linear programmes with a high degree of network structure, but with a small number of complicating (non-network) constraints present, a relatively simple heuristic has been found to be highly effective at identifying maximal subsets of constraints which can be converted to pure network problems. It has also been shown that much larger networks are identified by seeking a general network structure, rather than limiting the search to a subset of rows giving at most one + 1 and one - 1 in each column, as has been done previously. For linear programmes which can be converted to pure network problems, there is a high probability that the heuristic will find such a network. Appendix A. Formal statement of the heuristic

Step 1. Initialise dik *-- O, i = I . . . . . m - I; k ~- i + 1. . . . . m. Step 2. For j ~- I , . . . , n, for i ~ 1. . . . . m - 1, if a~j ~ 0 for k ~ - i + 1. . . . . m, if akj -~ O, if d~k ~ O, L~k ~- {r: i ~ r ~ k and a~j ~ 0}, 1

dik <--

Step 3.

Step 4. Step 5.

Step 6.

-- 1

if a~j-~ak~, if ajj ~ --ak~,

else if dik ~ 2, if aij ~ dik " ~ikj, dlk ~" 2, undcfine Lik, e[SC Lik ~- L~k n {r: a,j ~ 0}; next k; next i; next j. Choose any i I and i2 such that Lili2 ~- fJ or has minimum cardinality in the event of there being 110 Li~i2 ~ ~. Initialise the tree T with arc i, and add arc i 2 taking into account directions using dili 2. S ! ~ ~, S 2 ~ ~. For each i ~ T, check the lists Lik (k > i) and Lki (k < i) for k ~ T U S ! u S 2. If no such list exists, then: if S ! ~ ~, k ~ k s for any k s ~ S l and go to Step 7; else if $2 ~ ~, k ~ k~ for any k~ ~ S 2 and go to Step 7; else if there exists any k ~ T, append arc k to T at the root; go to Step 8. If such list(s) exist, then: If possible, choose k such that L~k c T and go to Step 6; else if S ! ~ {~, k ~- k~ for any k~ E S~ and go to Step 7; else choose k such that Lik - T has minimum cardinality. Try adding arc k to each end of the chain formed by the arcs of (Lik U {i}) n T. For each end check (i) every d~k (or dkr) for r ~ T to see that there is a direction for arc k which is consistent with every are in T, and (ii) each column j of A for which akj ~ 0, for any contradiction.

RM. Baker,P.Z Maye / Embeddednetworkstntcturehair, tic in MP

61

If arc k cannot be appended at either end, then k is omitted from further consideration for appending to T and Step 5 is repeated. If arc k can be appended at just one end, then go to Step 7. If are k can be appended at both ends, then if L~ c T, S t ~ S t U {k}; if L~k ¢ T, S z *'- S z U {k}; go to Step 5. Step Z Append arc k to T, choosing its position at random if there are 2 positions in which it can be placed. Step 8. If all arcs have been appended, stop. Else, if all arcs not appended have been omitted at Step 6, then repeat from Step 3 with different i t and i 2 or choose the best result found so far and stop. Else go to Step 4.

Appendix B. Complexity analysis For an m × tz coefficient matrix A, the heuristic presented in Appendix A has a worst case performance of O(nm4). Step 2 has O(nm4): For any pair of arcs i, k in a column, the greatest effort is checking whether arcs already in Lik are also in the column. This has O(m2). In eaehcolnmn there are at most (m - 1) + (m 2) + . . . + 2 + 1 = -[m(m - I) pairs to check. As there are n columns, Step 2 has O(nm4). Step 3 has O(m f) as there are ½m(m - 1) possible pairs of arcs to consider. Step 4 has O(m) as the maximum possible number of elements in S l or S 2 is m - 2. Steps 5 to 8 have O(nm4): These steps are repeated several times according to how many arcs are in the tree being constructed. The worst case would be that at each stage of construction of the tree, every arc not in the tree is considered for appending to the tree. For a particular arc to be considered for appending to the tree it may be necessary to examine every cell in the adjacency tableau which relates to that arc and arcs already in the tree. To examine a single cell takes O(m), to see if all elements in the cell are in the tree. There will be at most t cells to examine, where t is the number of arcs in the tree at the stage. Although t < m, it makes no difference to the complexity analysis if we assume that t ~ m. Hence, checking one arc at Step 5 as a candidate for appending to the tree has O(m2). At Step 6, to check a candidate are for consistent direction has O(m). To check a column of the matrix to verify that the corresponding chain is consistent with the tree takes at most O(m 2) since to cheek such a chain is equivalent to ordering m numbers. Thus, to check n columns has O(nm2). Step 7 has O(1) and Step 8 has O(m) as there are at most m arcs to considered. Therefore, the dominant effort in Steps 5 to 8 is checking an arc for appending to the tree, which has O(nm2). The number of arcs which must be checked is (m-2)+(m-1)+'..

+2+l----½(m-1)(m-2),

giving a complexity of O(nm4). Thus, the overall complexity of the heuristic is O(nm4).

Appendix C. Generating test problems Suppose we wish to generate an LP with 5 rows an~ 12 columns, which is convertible to a pure network.

B.M. Baker, P.J. Maye / Embedded network structure heuristic in MP

1

8

(I)

(11)

(al)

(iv)

Figure 4

1. A tree with 5 arcs (6 nodes) is generated as follows. Node 1 is drawn as a root node and node 2 is drawn with node 1 as predecessor. For i = 3 . . . . . 6, generate a random number, 1 g RN A i - 1, as predecessor to node i. Directions are inserted randomly. E.g. Figure 4(i) shows the tree corresponding to the random number sequence 1, 3, 3, 2. 2. If possible, for each node i an arc is generated with tail vertex i and head vertex not equal to i, so as not to duplicate or mirror any arc already generated. (This may not be possible, but for large problems it usually is.) E.g. Figure g(ii) shows how the network may now appear. 3. Any further arcs (one in this case) are generated randomly, with only Simple loops excluded (as in Figure 4(iii)). 4. Node numbers 2 . . . . . 6 are randomly permuted and arcs 1. . . . . 5 are re-named correspondingly (as in Figure 4(iv)). 5. An incidence matrix is constructed: Arc: 1

1

2

2

3

4

1

5

6

1

1

7

8

9

10

-1

1 -1

3

-1

-1

11

-1

-1

1

1

Node:

4

1

5 6

-1 1

1

-1

1

-1 -1

12

-1

1

-1

1

B.M. Baker, EL Maye / Embedded networkstructuret~uristicin MP

63

6. F r o m the incidence matrix, the coefficient matrix o f the L P for the network o f Figure 4(iv), with basis tree as shown, is obtained as follows: 12345 2i l l me: 3 4 5

6

1

7 -1

I

8

9 l0

1 -I

1 -- | 1

--I 1

-1 1

11

1 -1

-1

--

(Add rows 1, 2, 5, 6) (Row 4) (Row 5)

1 1

12

|

(Add rows I, 3, 4)

F o r the larger problems used in testing the heuristic, to corrupt a row, a random number, 1 _< R N _< 10, is generated and if the number of non-zero entries in the row exceeds RN, then R N o f these entries are chosen at random and re-set to zero.

References [1] Boston, V..LD,, Rahmouni, M.K, and Williams, H.P,, "The practical conversion of linear programmes to network flaw models", in: Universityof Swattmmpton Preprint Series, 1987. [2] BLxby,R.E., and Cunningham, W.H., "Converting linear programs to network problems", Mathematics of OperationsResearch 5 (1980) 321~357. [3] Bixby, R.E., and Fourer, R., "Finding embedded network rows in linear programs", Management Science 34 (1988) 342-376. [4] Bixby,R.E., and Wagner, D.K, "An almost linear time algorithm for graph realization", Mathematics of Operatious Research 13 (1988) 99-123. [5] Fugishige, S., "An efficient PQ-graph algorithm for solving the graph realization problem", Journalof Computerand System Sciences 21/1 (1980) 63-86. [6] Glover, F., and Klingman, D., "The Simplex SON Algorithm for LP/embedded network problems", MathematicalPregramming Study 15 (1981) 148-176. [7] Gould, R., "Graphs and vector spaces", Journalof Mathematics and Physics 37 (1958) 193-214. [8] Iri, M., "A necessary and sufficient condition for a matrix to be the loop or cut-set matrix of a graph and a practical method for the topological synthesis of a network", in: RAAG Research Notes, Third Series, No. 50, 1962. [9] lri, M., "A criterion for the reducibility of a linear programming problem to a linear network flow problem", RAAG Research Notes, Third Series, No. 98, 1966. [10] Iri, M., "On the synthesis of loop and cutset matrices and the related problems", RAAG Memeirs 4/A-XIII (1968) 4-38. [11] Maye,P..I., "Exact and heuristic methods of finding embedded network structure in mathematical programmes", Ph.D. Thesis, Dept. of Stats & OR, Coventry Polytechnic, England, 1990. [12] Wright, W., "Automatic identification of network rows in large-scale optimization models", M.S. Thesis, Naval Postgraduate School, USA, 1980.