110
Views
9
CrossRef citations to date
0
Altmetric
Section B

A mathematical programming approach to the construction of BIBDs

&
Pages 1067-1082 | Received 02 Aug 2009, Accepted 03 May 2010, Published online: 02 Feb 2011

Abstract

A balanced incomplete block design (BIBD) is instrumental in design of experiments. This is usually constructed by algebraic methods such as finite algebra or difference sets. However, in algebraic approaches, no unified method exists, and each BIBD has been constructed in some ad hoc ways. On the other hand, computer-based methods apply the same algorithm to all BIBDs; hence, these are unified approaches. Although various meta-heuristic algorithms have been tried, the success of these methods has been rather limited. This article presents an alternative approach to this problem that formulates the problem as a nonlinear mixed integer programming problem. We develop a branch-and-bound algorithm to solve this, and a tabu search algorithm to overcome some weakness in the former algorithm. We compare the performance of these algorithms against some previously developed algorithms, and demonstrate that our algorithms are competitive to these methods.

2000 AMS Subject Classifications :

1. Introduction

Let V be a set of v elements called points and B be a collection (i.e. multi-set) of b non-empty subsets of V called blocks. (V, B) is a balanced incomplete block design (BIBD Citation3 Citation8 Citation12 Citation22), if the following conditions are satisfied.

  • (i) Each point is contained in exactly r blocks.

  • (ii) Each block contains exactly k points.

  • (iii) Every pair of distinct points is contained in exactly λ blocks.

With positive integer parameters , this is denoted as BIBD. BIBD has its origin in statistical design of experiments Citation6, and recently it is also applied to the coding theory Citation2 Citation7, failure diagnosis and group testing Citation10, and scheduling sports leagues Citation1 among others.

In matrix notation, BIBD is a v×b binary matrix X=(x il ) such that

From these relations, the following is easily shown Citation12.

Proposition 1.1

In BIBD parameters satisfy

Then, b and r are determined by as

Thus, in literature, BIBD is often referred to as -BIBD, or alternatively 2- design Citation8 Citation25. Given a set of integer parameters satisfying EquationEquations (5) and Equation(6), we are concerned with the existence and construction of such a BIBD. For the limited case of k=3 or k=4, Hanani Citation13 Citation14 proved the following.

Theorem 1.2

Let v and λ be positive integers, and k=3 or k=4. Then, BIBD exists if and only if b and r determined by Equations (7) and (8) are both integers.

Triple system is the special case of BIBD with k=3, and Steiner triple system is with k=3 and λ=1 Citation8. These are denoted as TS(v, λ) and STS(v), respectively. For the Steiner triple system, the above result was shown by Kirkman Citation16 as early as in 1847, together with a construction method for STS(v). Hanani Citation14 proved similar results for k=5 and k=6 with some minor exceptions. However, except for these, neither necessary and sufficient conditions nor the construction method of BIBDs is known in general.

Various approaches have been tried for the construction of BIBDs. Among these, algebraic methods, such as finite projective (or affine) geometry and difference sets Citation12 Citation22, have been the main stream and most successful, with most of the BIBDs known today being discovered by these methods. Extensive list of BIBDs, as well as unknown BIBDs up to date, are available in the handbook by Colbourn and Dinitz Citation8. Unfortunately, these methods are quite complicated, and no unified approach exists for the construction of BIBDs by these methods. Thus, BIBDs have been constructed for each set of parameters by some ad hoc methods, in a case-by-case way.

Computational methods, such as constraint programming Citation20, neural networks Citation18 and meta-heuristic algorithms Citation5 Citation17, provide a unified framework. However, the ability of these methods for the construction of BIBDs has been rather limited, until recently some success stories have been reported. For example, by an exhaustive computer search over more than 90,000 days in total, Bilous et al. Citation4 proved that no (22,8,4)-BIBD exists. Similarly, Houghten et al. Citation15 declared the non-existence of (46, 6, 1)-BIBD. On the other hand, Morales Citation19 discovered six new BIBDs by an elaborate tabu search Citation11. In spite of these few successes, heuristic algorithms are usually quite time-consuming.

In this article, we present an alternative approach for the construction of BIBDs that formulates the problem as a nonlinear mixed integer programming (MIP) problem. If this problem is solved to optimality, we either obtain a BIBD, or have the evidence that no BIBDs exist for the given set of parameters. In this sense, our approach is complete, as opposed to other heuristic methods where success of the algorithms is totally by chance.

The organization of this paper is as follows. In Section 2, after formulating the problem as a nonlinear MIP problem, we present a backtracking framework to solve this nonlinear optimization problem by solving a series of linear MIP problems. Based on this, in Section 3, we develop a branch-and-bound algorithm for this problem. We find that this method is successful for triple systems, but for some BIBDs, it takes very long computing time, or even unable to produce a solution. To overcome this weakness, we present a tabu search algorithm in Sections 4, and 5 to compare these algorithms against the local search method proposed by Prestwich Citation21 on the 86 benchmark instances given in his paper. From these 86 instances, we were able to solve 78 instances, while the local search algorithm of Prestwich Citation21 solved only 63 instances.

2. Backtracking algorithm

Finding a matrix satisfying EquationEquations (1)Equation(4) is a kind of constraint satisfaction problem which can be transformed into the following nonlinear 0–1 programming problem Citation23.

Note that equalities (1)–(3) are relaxed to inequalities (10)–(12), and the objective function is simply the sum of the left-hand sides of these inequalities. Then, the following is trivial.

Theorem

  • (i) P is always feasible.

  • (ii) For an arbitrary feasible solution X=(x ij ) to P, let z(X) denote its objective value. Then, we have

  • (iii) If Equation (14) is satisfied with equality, Equations (10)–(12) are also satisfied with equality, and thus X gives a BIBD

Therefore, if equality holds in EquationEquation (14) with respect to an optimal solution X* to P, we are done. Otherwise, if , no BIBD exists. However, since P is a nonlinear 0–1 programming problem, it is difficult to get an optimal solution. Thus, we propose an incremental approach that determines rows of X one by one.

Now, suppose that the first j rows of X is known as . From EquationEquations (1)Equation(4), this is a binary matrix that satisfies the followings.

Then, the (j+1)th row vector must satisfy

If such an is found, we can augment X j to a matrix

To obtain a binary vector satisfying EquationEquations (18)Equation(20), we formulate the following optimization problem.

Contrary to nonlinear P, P j (X j ) is a linear 0–1 programming problem which can be solved (in many cases) using free or commercial MIP solvers. Let the optimal objective value to P j (X j ) be . Then, since the objective function Equation(22) is the sum of the left-hand sides of EquationEquations (23) and Equation(25), any optimal solution to P j (X j ) gives a binary satisfying EquationEquations (18)Equation(20) if and only if

Without loss of generality, the first and second rows of X j can be assumed to be

and a BIBD is obtained if we can augment this through EquationEquation (21) to a v×b matrix X v . On the other hand, if
no BIBD exists as an extension of X j , and thus subproblem P j (X j ) is terminated.

Example 2.2 For BIBD (12,22,11,6,5), except for trivial inequalities and 0–1 conditions, P 2(X 2) is

Solving this, we have with , and
Continuing this for j=3, …, 8, we obtain
We stop here since , and restart with the eighth row of X 8 (denoted as ) changed to an optimal solution of P 7(X 7) other than . If no such solutions exist, we backtrack.

Here, we note that the optimal solution to P j (X j ) may not be unique even if EquationEquation (27) is satisfied. If this is the case, we do not know exactly which optimal solution to adopt in EquationEquation (21). In the branch-and-bound algorithm to be developed in the next section, we list up all the optimal solutions of P j (X j ) and define its children by appending each of these solution vectors to X j . Thus, we have a tree of subproblems rooted at P 2(X 2), and traversing all the subproblems obtain a BIBD, or prove that no BIBD exists.

To list up all the optimal solutions of P j (X j ), we introduce the following auxiliary problem.

Here, F and R are the disjoint subsets of variables which are fixed either at 1 and 0, respectively. We obtain an optimal solution to P j (X j ) by solving with an MIP solver. Let this be , with being the set of indexes such that . Then, by solving for s=1, 2, …, p, we obtain a set of distinct solutions to P j (X j ) other than , and repeating this recursively obtain the set of all optimal solutions to P j (X j ). This recursive procedure is included in the branch-and-bound algorithm to be given in the next section.

Example 2.3 Applying the above algorithm to P 7(X 7) with , we obtain , which is the eighth row of X 8 in EquationEquation (30). From the indexes of non-zero elements of , we have . Then, auxiliary problems are generated by adding constraints and to P 7(X 7) as {u s }), and the feasible region of is divided into mutually disjoint feasible regions of . Here for s=0, we mean , and from EquationEquation (23) is the unique optimal solution to this problem. Then, by solving recursively, we obtain all the optimal solutions to , together with .

3. A branch-and-bound method

We construct a branch-and-bound algorithm Citation23 Citation24 to find a BIBD in the following way. First of all, we define as the continuous relaxation of , and denotes an optimal solution to this problem with the corresponding objective value . Then, we can terminate subproblem if

We can check if EquationEquation (33) is satisfied by solving a linear programming (LP) problem. If EquationEquation (33) is not met, we solve the original integer programming (IP) problem to see if EquationEquation (27) is satisfied. Then, the algorithm is given as follows.

Algorithm

BAB_BIBD

Input: j, X j , F, R.

  • Step 1. If j=v, X j is a BIBD. Output this and stop.

  • Step 2.

    • (i) Solve an LP problem and obtain and

      • (a) If then terminate this subproblem and return.

      • (b) Otherwise, if the solution is a 0–1 vector, let , and go to Step 3.

    • Solve an IP problem and obtain and If terminate this subproblem and return.

  • Step 3. Append to X j and update this to X j+1 through EquationEquation (21). Call BAB_BIBD recursively.

  • Step 4. Let be the set of indexes such that and For do

    • (i) Let and

    • (ii) Call BAB_BIBD recursively.

  • Step 5. No BIBDs exist. Stop.

Some remarks are in order on the branch-and-bound strategies. First of all, linear relaxation is not mandatory, i.e. we may skip part (i) of Step 2. In the standard case, we solve an LP problem first, and then solve P j (X, F, R) if necessary, as described in the algorithm. This is referred to as the bounding strategy LP. If we skip part (i) in Step 2 and go directly to part (ii), we call this bounding strategy IP.

Next, recursive calls to subproblems can be carried out in the reverse order. For the problem in Example 2.3, the original algorithm calls subproblems in the order of . We call this the forward branching strategy (FWD). In the backward strategy (BWD), these subproblems are called reversely as . We denote these strategies explicitly as BAB_BIBD (bounding strategy, branching strategy), such as BAB_BIBD(LP, FWD).

Example 3.1 shows the behaviour of the algorithm under the bounding strategy IP and branching strategy FWD, i.e. BAB_BIBD(IP, FWD). As in Example 2.2, subproblem P 8 is terminated, and by solving we obtain an alternative solution that replaces the eighth row of X 8 in Example 2.2. In , this is P8. Restarting from this point, the process stops again at P 10. Here, we switch to P10, and finally obtain a BIBD as

after solving 13 IP problems.

Figure 1. Behaviour of BAB_BIBD(IP, FWD) for BIBD(12,22,11,6,5).

Figure 1. Behaviour of BAB_BIBD(IP, FWD) for BIBD(12,22,11,6,5).

Figure 2. Behaviour of BAB_BIBD(LP, FWD) for BIBD(12,22,11,6,5).

Figure 2. Behaviour of BAB_BIBD(LP, FWD) for BIBD(12,22,11,6,5).

Example 3.2

If we take bounding strategy LP instead of IP in the previous example, the algorithm BAB_BIBD(LP, FWD) produces the tree of subproblems as depicted in . Here, shadowed nodes show the subproblems where both of LP and IP problems are solved, while at the remaining nodes we solve only LP problems. Thus, we obtain a BIBD

after solving 12 LP and seven IP problems. Note that the obtained BIBD (EquationEquation (35)) is different from the one shown in EquationEquation (34), as early as at j=3.

Example 3.3

shows the behaviour of the backward branching strategy, i.e. BAB_BIBD(LP, BWD). The algorithm proceeds exactly the same way as in the previous example until again P 9 is terminated (compare and ). Then, auxiliary subproblems are generated and examined as in the previous example but in the reverse order. The first seven subproblems are terminated due to EquationEquation (33), we restart from and at P 11 obtain a BIBD

after solving 18 LP and eight IP problems. Note that this BIBD is different again either from EquationEquation (34) or Equation(35), although the first eight rows of EquationEquations (35) and Equation(36) are identical.

Figure 3. The behaviour of BAB_BIBD(LP, BWD) for BIBD(12,22,11,6,5).

Figure 3. The behaviour of BAB_BIBD(LP, BWD) for BIBD(12,22,11,6,5).

Now, we evaluate the performance of the branch-and-bound algorithm for triple systems TS(v, λ). shows the result of BAB_BIBD(LP, FWD) on a DELL Precision 670 computer (CPU: Xeon ) with CPLEX 10.100 Citation9 to solve LP and IP problems. Here shown are the number of subproblems generated (#Subp), the number of LP and IP problems solved (#LP, #IP) and the CPU time in seconds. We solved three randomly selected instances for each value of , within v≤200 and b≤3000.

Table 1. Result of BAB_BIBD(LP, FWD) for triple systems TS(v, λ).

From this table, we see that in almost all instances , which means that backtracking seldom occurs in triple systems. Next, #IP is relatively small, i.e. in most cases we obtain BIBD by solving approximately v LP and some additional IP problems. This implies that we frequently obtain 0–1 solutions by solving only LP problems, and thus the bounding strategy LP is well suited to this case. Thus, the branch-and-bound method works effectively for triple systems.

BAB_BIBD is often inefficient for problems where backtracking occurs frequently. For example, in solving BIBD(10,15,6,4,2) by BAB_BIBD(LP, FWD), we had to solve 1,207,105 LP and 343,523 IP problems and the CPU time was 2068 s. We were unable to solve BIBD(16,16,6,6,2) within 3600 s, while the local search algorithm by Prestwich Citation21 solved this within 1 s on a slower computer.

The reason for this difficulty of BAB_BIBD for some BIBDs may be explained as follows. In the case of triple systems, usually there are many different (non-isomorphic) solutions for each set of parameters Citation8, and thus it does not much matter which optimal solution of P j (X j ) we choose in augmenting X j to X j+1 through EquationEquation (21). In most cases, any solution is satisfactory, and no backtracking is necessary. On the other hand, in general BIBDs, there are only a few number of solutions for each set of parameters. In such a case, we must be very careful in selecting a solution of P j (X j ). If we make a wrong choice and ‘bad’ rows are included in X j , it takes a long time until we backtrack and eliminate these bad rows.

4. Tabu search

To overcome the weakness of BAB_BIBD to some instances, we propose a tabu search method. The idea is that in X j we may have some bad rows, which means a 0–1 vector that is not included in any BIBDs of the given set of parameters. Then, after solving P j (X j ), we make a guess whether bad rows are included in X j . If we find a doubtful row, we eliminate it and restart with the reduced X j−1. More specifically, if we solve P j (X j ) and obtain an optimal solution with , we move forward by adding to X j as we did in the BAB_BIBD. On the other hand, if we have either

for some row i, or
In the case of EquationEquation (37), we consider as a bad row and eliminate it from X j . If no such row exists, we have EquationEquation (38). Then, we pick up a row vector in X j at random and remove it as well. In either of these cases, after eliminating such a row vector from X j , we restart with a reduced matrix with j−1 rows.

To prevent a removed vector from returning to X j in the subsequent few steps, we introduce a tabu list 𝒯 to keep the eliminated vectors ‘frozen’ for a certain period of steps. The size of 𝒯 is specified by a parameter TL (tabu length), and vectors in the tabu list are forbidden to be an optimal solution to P j (X j ). To realize this, let be the current tabu list, and consider the following 0–1 programming problem.

Let be the optimal objective value to this problem. Due to EquationEquation (39), no vectors in the tabu list satisfy EquationEquation (23) with equality. Then, if holds, the optimal is not tabooed, and thus we can expand X j by appending this vector to X j . Otherwise, we choose a row vector as stated before and eliminate it from X j . The tabu search algorithm Citation11, with the tabu length explicitly shown as TABU_BIBD(TL), is as follows.

Algorithm 1

TABU_BIBD

  • Step 1. Let j:=2, X 2 be given as Equation (28), and

  • Step 2. If j=v, then X j is a BIBD. Output this and stop.

  • Step 3. Solve an IP problem If , go to Step 4. Otherwise go to Step 5.

  • Step 4. Append the solution to X j and augment it to X j+1. Go back to Step 2 with

  • Step 5. If there exists a row i satisfying Equation (37), go to Step 6. Otherwise pick up a row vector from X j at random and go to Step 6.

  • Step 6. Eliminate the selected vector from X j , let and add the eliminated vector to the tabu list 𝒯. (If the oldest vector is removed from 𝒯.)

  • Step 7. If terminating condition is satisfied, then print ‘Tabu search failed’, and stop. Otherwise go back to Step 2.

Example 4.1 In applying the tabu search algorithm to BIBD(16,16,6,6,2) where BAB_BIBD failed to find a solution within 3600 CPU seconds, the tabu search produced a correct BIBD in 0.02 s as follows. Initially, the tabu list is initialized as . The algorithm behaves exactly the same way as BAB_BIBD(IP, FWD) until at j=8 we obtain

Here we have an optimal solution with . The first row conflicts against , since x . We remove from X 8 and solve with the reduced matrix X7. Again we have a conflicting row and eliminating this, obtain a matrix X6 with six rows. After this, 10 forward moves are repeated and at this point, we obtain a BIBD(16,16,6,6,2) as

Prestwich Citation21 proposed a local search algorithm, and tried 86 instances with v and b satisfying vb≤1000. compare his result against our branch-and-bound and tabu search methods with MIP solvers. In these tables, CPU time in seconds is shown for each method. The column of ‘Prestwich’ is the result of his method on a DEC Alphaserver 100A 5/300 computer (300 MHz), while ‘BAB’ is by the branch-and-bound algorithm BAB_BIBD(LP, FWD), and ‘TABU’ is by TABU_BIBDEquation(10) on a faster DELL Precision 670 machine . Computation was truncated (and shown in tables by –) at 2 million backtracks in ‘Prestwich’ and 3600 CPU seconds in our algorithms.

Table 2. Result of computation for Prestwich's 86 instances (Part 1).

Table 3. Result of computation for Prestwich's 86 instances (Part 2).

Table 4. Result of computation for Prestwich's 86 instances (Part 3).

From these tables, we see that the branch-and-bound method solves more instances than the Prestwich method. However, in some instances, it took much longer time, or even unable to solve, while the same instances were solved easily by the Prestwich local search algorithm. By the tabu search method of this paper, we were able to solve all the instances that were solved by Prestwich. The tabu search algorithm solved some additional instances that neither Prestwich's nor the branch-and-bound algorithms were able to solve.

gives a summary of the numbers of instances out of 86 solved by Prestwich's and our methods. The branch-and-bound method solved approximately 10 more instances than Prestwich's, irrespective of the bounding and branching strategies. Tabu search was able to solve 10 additional instances, irrespective of the tabu length.

Table 5. Summary of computation.

5. Conclusion

We have formulated the problem of constructing BIBDs as a nonlinear MIP problem, and presented branch-and-bound and tabu search algorithms that makes use of MIP solvers. The developed branch-and-bound method was successful for triple systems, but it sometimes encountered difficulty for general BIBDs. Tabu search improved this situation to some extent. In some instances, it succeeded in eliminating bad rows in earlier stages, and obtained BIBDs that the branch-and-bound method could not. Compared with heuristic algorithms proposed by previous researchers, our algorithms appear to be competitive. However, even with these algorithms we were unable to find new BIBDs. To solve problems with v≥30, some new bounding conditions are required, and thus make the algorithm more efficient.

Acknowledgements

The authors are grateful to two anonymous referees for their constructive criticisms and suggestions, which improved the exposition of the paper significantly.

References

  • Anderson , I. 2006 . Combinatorial Designs and Tournaments , New York : Oxford University Press .
  • Beth , T. , Jungnickel , D. and Lenz , H. 1999 . Design Theory , 2 , Cambridge : Cambridge University Press .
  • Biggs , N. L. , Lloyd , E. K. and Wilson , R. J. 1995 . The history of combinatorics, in Handbook of Combinatorics , Edited by: Graham , R. L. , Grötschel , M. and Lovász , L. Vol. II , 2163 – 2198 . Amsterdam : Elsevier .
  • Bilous , R. , Lam , C. W.H. , Thiel , L. H. , Li , B. P.C. , van Rees , G. H.J. , Radziszowski , S. P. , Holzmann , W. H. and Kharaghani , H. 2006 . There is no 2-(22, 8, 4) block design . J. Combin. Des. , 15 : 262 – 267 .
  • Bofill , P. , Guimera , R. and Torras , C. 2003 . Comparison of simulated annealing and mean field annealing as applied to the generation of block designs . Neural Netw. , 16 : 1421 – 1428 .
  • Box , G. E. , Hunter , W. G. and Hunter , J. S. 2005 . Statistics for Experimenters: Design, Innovation, and Discovering , 2 , Hoboken : John Wily & Sons .
  • Brouwer , E. 1995 . Block designs, in Handbook of Combinatorics , Edited by: Graham , R. , Gröetschel , M. and Lovász , L. Vol. I , 693 – 771 . Amsterdam : Elsevier .
  • Colbourn , C. J. and Dinitz , J. H. 2007 . Handbook of Combinatorial Designs , 2 , New York : Chapman & Hall/CRC .
  • CPLEX Ver. 11.0, ILOG, 2008. Available at http://www.ilog.com/products/cplex/
  • Du , D.-Z. and Hwang , F. K. 2000 . Combinatorial Group Testing and its Applications , 2 , Singapore : World Scientific Publication .
  • Glover , F. and Laguna , M. 1997 . “ Tabu Search ” . Norwell : Kluwer .
  • Hall , M. Jr. 1998 . Combinatorial Theory , 2 , New York : John Wiley & Sons .
  • Hanani , H. 1961 . The existence and construction of balanced incomplete block designs . Ann. Math. Stat. , 32 : 361 – 386 .
  • Hanani , H. 1975 . Balanced incomplete block designs and related designs . Discrete Math. , 11 : 255 – 369 .
  • Houghten , S. K. , Thiel , L. H. , Janssen , J. and Lamet , C. W.H. 2001 . There is no (46,6,1) block design . J. Combin. Des. , 9 : 60 – 71 .
  • Kirkman , T. 1847 . On a problem in combinations . Cambridge Dublin Math. J. , 2 : 191 – 204 .
  • Kreher , D. L. and Stinson , D. R. 1999 . Combinatorial Algorithms: Generation, Enumeration and Search , Boca Raton : CRC Press .
  • Kurokawa , T. and Takefuji , Y. 1992 . Neural network parallel computing for BIBD problems . IEEE Trans. Circuits Syst. II , 39 : 243 – 247 .
  • Morales , L. B. 2000 . Constructing difference families through an optimization approach: Six new BIBDs . J. Combin. Des. , 8 : 261 – 273 .
  • Prestwich , S. Balanced incomplete block design as satisfiability . Irish Conference on AI and Cognitive Science . Cork. Vol. 12 , pp. 189 – 198 .
  • Prestwich , S. 2002 . A local search algorithm for balanced incomplete block designs , Cork : University College Cork . Lecture Notes in Artificial Intelligence, Cork Constraint Computation Centre
  • Sachkov , V. N. 1977 . Combinatorial Methods in Discrete Mathematics , Cambridge : Cambridge University Press .
  • Schrijver , A. 1998 . Theory of Linear and Integer Programming , Amsterdam : John Wiley & Sons .
  • Sedegewick , R. 1998 . Algorithms in C , 3 , Boston : Addison-Wesley .
  • Stinson , D. R. 2004 . Combinatorial Designs: Construction and Analysis , New York : Springer-Verlag .

Reprints and Corporate Permissions

Please note: Selecting permissions does not provide access to the full text of the article, please see our help page How do I view content?

To request a reprint or corporate permissions for this article, please click on the relevant link below:

Academic Permissions

Please note: Selecting permissions does not provide access to the full text of the article, please see our help page How do I view content?

Obtain permissions instantly via Rightslink by clicking on the button below:

If you are unable to obtain permissions via Rightslink, please complete and submit this Permissions form. For more information, please visit our Permissions help page.