An improved test set approach to nonlinear integer problems with applications to engineering design

AN IMPROVED TEST SET APPROACH TO NONLINEAR INTEGER PROBLEMS WITH APPLICATIONS TO ENGINEERING DESIGN J. GAGO-VARGAS, M.I. HARTILLO-HERMOSO, J. PUERTO-ALBANDOZ, AND J.M. UCHA-ENR´ IQUEZ Abstract. Many problems in engineering design involve the use of nonlinearities and some integer variables. Methods based on test sets have been proposed to solve some particular problems with integer variables, but they have not been frequently applied because of computation costs. The walk-back procedure based on a test set gives an exact method to obtain an optimal point of an integer programming problem with linear and nonlinear constraints, but the calculation of this test set and the identification of an optimal solution using the test set directions are usually computationally intensive. In problems for which obtaining the test set is reasonably fast, we show how the effectiveness can still be substantially improved. This methodology is presented in its full generality and illustrated on two specific problems: (1) minimizing cost in the problem of scheduling jobs on parallel machines given restrictions on demands and capacity, and (2) minimizing cost in the series parallel redundancy allocation problem, given a target reliability. Our computational results are promising and suggest the applicability of this approach to deal with other problems with similar characteristics or to combine it with mainstream solvers to certify optimality. Non-linear Integer Programming and test set and Gr¨obner basis and chance constrained programming 1. Introduction The ever-increasing demand on operations researchers to lower production costs or to increase benefits has prompted the specialized community to look for rigorous methods of decision making, such as optimization methods, to design and process both economically and efficiently most engineering systems. Optimization techniques, having reached a degree of maturity over the last years, are being used in a wide spectrum of engineering applications, including reliability of systems, finance, scheduling, as well as various contributions in aerospace, automotive, chemical, electrical or manufacturing industries. For some of those applications Operations Research has developed efficient algorithms, especially when the problems can be modeled either as continuous (linear or nonlinear) programs or discrete linear programs. However, in many occasions the resulting models contain both nonlinearities and discrete variables which make their resolution challenging. The purpose of this paper is to improve a general technique to handle nonlinear integer problems and to show some applications to engineering optimization problems in a simple manner. This technique combines three essential elements: 1) it borrows some tools from Computational Algebra which are applied to generate Date: Received: date / Accepted: date. 1 2J. GAGO-VARGAS, M.I. HARTILLO-HERMOSO, J. PUERTO-ALBANDOZ, AND J.M. UCHA-ENR´ IQUEZ test sets for linear integer programs, 2) it uses relaxed problems to obtain seeds for a dual search procedure, and 3) it defines an order on the search tree of feasible solutions that uses the classical idea of a penalty strategy that should take into account both the objective function value and the distance to the feasible region. This third ingredient is specifically the new one that makes the general dual search (or walk-back procedure) much more efficient in the experiments. The general theory of penalty functions usually considers the transformation of a constrained problem into unconstrained one(s) and afterward, starting from some convenient seed, uses suitable directions in order to reach an optimum. Our approach is a natural generalization in this context: once the optimum for a relaxed integer linear problem of the original problem is reached, the process starts from this optimum and chooses the most promising nodes balancing the cost and the distance to the whole feasible region. The main advance of our algorithm with respect to the one described in [30] is the insertion in the pending nodes list by the ascending penalized cost. We will see in Section 4 the effectiveness of this change in some examples. The method is described in its full generality and then it is tested in two engineering design problems taken from the literature, namely, the minimization of the cost of scheduling jobs on parallel machines given restrictions on demands and capacity and in the minimization of the cost in the series parallel redundancy allocation problem given a target reliability. For the first problem, with stochastic restrictions, it is remarkable that our method can manage the examples tested whereas other mainstream nonlinear solvers cannot deal with the nonlinear formulation. The results obtained in the second problem are competitive with the nonlinear solvers, with a remarkable performance in the time needed to reach the optimum. The paper is organized as follows. In Section 2 the basics about the walk-back procedure using a test set associated with an integer linear program are presented. Section 3 contains the main technique followed to improve the pure walk-back procedure. In this section it is also introduced the penalized cost function that can be employed to order, in an efficient way, the feasible points of the relaxed integer linear problem. Section 4 is devoted to the examples where the method has been tested. In all the examples tables with execution data (CPU time, number of processed nodes) and comparison with previous works ([30], [25], [17]) and other nonlinear integer programming solvers are included. Additionally, we show how to exploit a combined strategy with other solvers to speed up certification of optimality. Finally, Section 5 contains the conclusions. 2. Preliminaries This section contains a brief summary of the concepts and algorithms used to solve integer linear programming problems from an algebraic point of view. In addition it is recalled the walk back procedure for nonlinear integer programming problems based on test sets. To this end, we have followed [28] and [30]. 2.1. Test sets. Consider an integer linear programming problem: (LP(b)) min c(x) = ct·x s. t. A·x=b, x∈ZN +, AN IMPROVED TEST SET APPROACH TO NONLINEAR INTEGER PROBLEMS 3 where A∈Zd×N,b∈Zd,c∈RN. The notation (LP(b)) denotes the integer linear programming problem with right-hand-side fixed to b. (LP) denotes the set of all the integer linear programming problems obtained by varying the right-hand-side vector b, fixed Aand the cost function c. Let πbe the map defined by π(x) = Ax. Given a vector b∈Zd, the set π−1(b) = {u∈ZN +:π(u) = b}is the fiber of (LP) over b. Fixed a term order <c, associated with the cost, defined as usual in [28], there exists a unique optimum βfor (LP(b)). A set G<c⊂ZNis a test set for the family of integer linear problems (LP) with respect to the matrix Aand the order <cif •for each nonoptimal point αof (LP(b)), there exists g∈G<csuch that α−gis a feasible solution for (LP(b)) and α−g<cα, •for the optimal point β,β−gis not a feasible point for any g∈G<c. In this way, a test set for (LP) connects the fiber as a directed graph with a unique sink. It gives an obvious algorithm to solve an integer program, starting from a feasible solution to this problem. At every step of this algorithm, we have two different cases: •There exists an element in the test set which, when subtracted from the current point, yields an improved point. •There does not exist such an element in the set, so that the point is the optimum of the fiber. 2.2. Toric ideal. The toric ideal associated with A, denoted as IAis defined as IA=hxα−xβ:Aα=Aβ,α,β∈ZN +i ⊂ Q[x1,...,xN]. Given a term order <ccompatible with the cost function c(x), a reduced Gr¨obner basis G<cof IAyields a test set G<cfor (LP) (cf. [28]). If the reduced Gr¨obner basis is formed by binomials G<c={xαi−xβi, i = 1,2,...,r},with in<c(xαi−xβi) = xαi, where in<cstands for the greatest monomial with respect to <c, then the test set is expressed as G<c={αi−βi, i = 1,2,...,r}. 2.3. Walk back procedure. The walk back procedure is an algorithm which computes the optimum for a nonlinear integer programming problem under some conditions. Our problem is of the form: (P) min c(x) = ct·x x∈ A ⊂ Zn +, x∈ B ⊂ Zn +, where •the constraints x∈ A are linear, and a test set can be obtained in practice with respect to the linear cost function and these restrictions; •the constraints x∈ B can be linear and nonlinear. We can assume without loss of generality that c(x) is a linear cost function because if it were nonlinear we would transform the problem introducing the new constraint c(x)≤z, that would be included in the set x∈ B, and the objective function min z. The problem would now fit within the original hypothesis, but the test set will change because the cost function that defines the term order is different. 4J. GAGO-VARGAS, M.I. HARTILLO-HERMOSO, J. PUERTO-ALBANDOZ, AND J.M. UCHA-ENR´ IQUEZ Let (LIP) denotes the relaxed linear integer problem, namely: (LIP) min c(x) x∈ A ⊂ Zn +. Let Gbe the test set associated with the relaxed linear problem (LIP). Starting from some solution of (LIP) and adding vectors from G, every feasible point of the relaxed problem (LIP) can be reached. So, in particular, every feasible point of problem (P) can be reached as well. But we do not need to complete an exhaustive enumeration of the feasible points of (P): if p′is a feasible point for (P), with cost c(p′), all the paths starting from points whose cost is greater than c(p′) can be pruned. Remember that subtracting the elements of the test set always reduces the cost, so adding them (or subtracting the elements of the reverse test set) provides points that do not improve the value c(p′). When a path is pruned, all its points are discarded. Points with negative components are deleted too. Along this process new points are generated, called pending nodes, and the paths starting from them must be analyzed. At the end all the possibilities will have been processed, and the final result will be an optimum of (P), or a certificate that problem (P) has an empty feasible region if no feasible point could be reached. This is the walkback procedure, described in [30]. The above argument can be summarized in the following result: Theorem 2.1. The walk-back procedure solves the problem (P), that is, it obtains an optimum or shows that the problem has no feasible solution. To the best of our knowledge, there have been few practical examples where the walk-back procedure has been successfully applied, as in [30], [17], [9]. This lack of applicability may be due to its two main drawbacks: the computation of the test set and the time required for visiting the points x∈ A to eventually obtain an optimum for (P). •The computation of the test set has, in general, a high cost. We have used the toric ideal approach [28], which is implemented in a very efficient way in the software 4ti2 [1]. In some cases, there is no need of such computation, because it is possible to give a closed form of a Gr¨obner basis of the problem, as it is made for example in [17]. •The points are ordered taking into account their cost. As soon as a feasible point for (P) is reached, new points are found quickly, and one of them is often an optimum. However, the processing time to reach such feasible points is usually long and the list of pending nodes is huge. The question is whether it is possible to order the points that have to be visited in an alternative way, so that new feasible better points are obtained earlier. This approach is explained in Section 3. 3. Combining test sets and a penalized strategy 3.1. Penalty functions. Traditionally, methods using penalty functions transform a constrained problem into a single unconstrained problem or into a sequence of unconstrained problems. Unconstrained problems are easier to manage in general. The general idea is to place the constraints into the objective function with a penalty parameter that penalizes properly any violation of these constraints. AN IMPROVED TEST SET APPROACH TO NONLINEAR INTEGER PROBLEMS 5 Given a problem (P), (P) min c(x) x∈ A ⊂ Zn +, x∈ B ⊂ Zn +, apenalty function p(x) is a function such that p(x) = 0 for x∈ A ∩ B and p(x)≥γ > 0 if x∈ A \ B. Moreover, the statement γ > 0 comes from the discrete nature of (P). For a problem (P) the penalty problem (Ppen(µ)) is defined as minx∈AT(x, µ) = c(x) + µp(x),for a certain µ > 0. Since T(x, µ) = c(x) if x∈ A ∩ B, it is clear that opt(P) ≥opt(Ppen(µ)). In the problems that we deal with in this work the region Bis possibly an enormous, but finite, subset of Zn. It can be proved that, in this situation the problem is equivalent to a penalized one. Theorem 3.1 (cf. [20] Th. 2.11).Suppose that A \B 6=∅. Let c0be a lower bound of c(x)in Aand ρ > 0be a lower bound of minx∈A\B p(x). Then, for any µ > µ0 with µ0=opt(P) −c0 ρ we have opt(P) = opt(Ppen(µ)). 3.2. Ordering nodes by penalized cost. Following the notation of [30] we describe in this subsection our algorithm. The walk-back procedure provides a list of feasible points, noted P, for the relaxed linear problem (LIP) and a guided search procedure that allows us to get an exact optimum for the original problem (P). However, the way in which the points are scanned is fundamental to create an algorithm that performs the computations in a reasonable time. A first approach is to order the pending nodes list according to their cost ([30], [17]). A second approach, presented here, consists of ordering the pending nodes by a penalized cost function of the form T(x, µ) = c(x) + µp(x), where the cost of the point is added to a value associated with the distance to the region defined by the constraints of Bin problem (P). Instead of ordering the points to be visited by the original cost function c, the penalized cost function T(x, µ) is used, expecting to reach faster some points that are feasible for the whole problem (P). This is the main advance of our algorithm compared to the one described in [30]: the insertion in the pending nodes list ordered by the new penalized cost function. It is worthwhile to emphasize that this procedure is exact and certifies the optimality as proved in Theorem 2.1, because it is only a different way of visiting the points required to solve the problem. It starts from the optimum βof (LIP) (that we assume infeasible for (P) since otherwise we are done), and produces recursively new points xthat are in two possible situations: 6J. GAGO-VARGAS, M.I. HARTILLO-HERMOSO, J. PUERTO-ALBANDOZ, AND J.M. UCHA-ENR´ IQUEZ (1) xis feasible for (P), so we include it in the set Yof feasible points for (P) obtained so far. It is not necessary to follow paths from xsince the cost would be worse; (2) xis not feasible for (P) and we do not have points in Ywith better cost, then xis inserted according to the function T(x, µ) in the set Pof nodes to examine. 4. Applications 4.1. Correlated setup. The walk back procedure was introduced in [30] and applied in a problem of assignment of jobs to machines, with given production and correlated setup costs, capacity constraints and probability to reach a given demand. We deal with the same problem by adding the improvement of the penalized cost and a set of new cuts that does not alter the optimum. The notation for the model is the following: •nnumber of job types, indexed by i, •mnumber of machines, indexed by j, •(D1,...,Dn) random vector of demands, •Nsize of the sample set used to estimate the probability, •Cjcapacity (time) for each machine, •(ˆ D1,..., ˆ Dn) means vector of the probability distribution of demand, •Sij setup time for job type ion machine j, •Kij setup cost for job type ion machine j, •Milot splitting, •L′ ij the cost of producing a unit of product type ion machine j, •Lij = ( ˆ Di/Mi)L′ ij, •pij processing time for a unit of job type ion machine j, •γprobability of no shortfall. •zij equals 1 if job type iis scheduled on machine j, 0 otherwise, •yij multiples of 1/Miof demand of product iare scheduled on machine j. The model presented in [30] is: (SP) minimize X iX j (Kijzij +Lijyij) subject to m X j=1 yij =Mi, i = 1,2,...,n,(1) Mizij ≥yij, i = 1,2,...,n,j = 1,2,...,m,(2) gj(z,y) = n X i=1 pij ˆ Di/Miyij + n X i=1 Sijzij ≤Cj, j = 1,2,...,m,(3) g0(z,y) = Prob (n X i=1 pij (Di/Mi)yij + n X i=1 Sijzij ≤Cj, j = 1,2,...,m)≥γ,(4) zij ∈ {0,1}, yij ∈ {0,1,...,Mi}(5) The condition (4) will be computed in a sample dataset of size N. As usual we note (LSP) the linear relaxed problem defined by constraints (1), (2) and (5). AN IMPROVED TEST SET APPROACH TO NONLINEAR INTEGER PROBLEMS 7 The model (SP) admits a family of optimality cuts: (6) zij ≤yij, i = 1,...,n,j = 1,...,m. It is clear that if the job iis assigned to the machine j, then something is produced in it. These constraints reduce the feasible region, but do not alter the optimum of problem (SP). We call problem (ISP) the resulting model obtained from (SP) after adding the constraints (6). The test set is associated with the linear constraints (1), (2), (5) and (6), which defines a relaxed problem (LISP). The equations (3) and (4) are checked for every point and test feasibility for Problem (ISP). Similarly to [30], the linear problem with restrictions (1), (2), (5) and (6) can be expressed with equalities by adding slack variables: (LISP −II) minimize X iX j (Kijzij +Lijyij) subject to m X j=1 yij =Mi, i = 1,2,...,n,(7) yij +aij −Mizij = 0, i = 1,2,...,n,j = 1,2,...,m,(8) zij +bij = 1, i = 1,2,...,n,j = 1,2,...,m,(9) −yij +zij +cij , i = 1,2,...,n,j = 1,2,...,m,(10) yij, aij, bij, cij ∈Z+, i = 1,2,...,n,j = 1,2,...,m.(11) Theorem 4.1. The set of binomials yipaiqcip −yiqaipciq, i = 1,2,...,n,1≤p < q ≤m, bipcip −zipaMi ip , i = 1,2,...,n,p= 1,2,...,m, bipyiqciq −zipyipaMi−1 ip aiq, i = 1,2,...,n,1≤p6=q≤m, bipyiqziqaMi−1 iq −zipyipbiqaMi−1 ip ,1≤p < q ≤m, is the reduced Gr¨obner basis of the toric ideal corresponding to (LISP-II) with respect to the lexicographical order b > y > z > a > c. The letter ydenotes the set of variables (y11,...,ymn)and similarly a, z, b and cdenote the other sets of variables. Within a block, the variables are sorted lexicographically according to index. The proof is written in Appendix A. However, this basis is not the one needed to build the test set with respect to the cost function. We have used 4ti2 to get the basis in a straightforward way, and the computation cost has been always under one second for all the considered instances, which are of the same size of those proposed in [30]. 4.1.1. New penalized cost function. Following the notation of Section 3, let xbe the set of variables yij, zij. Multiplying by a suitable constant, we can assume that all the coefficients in gj(x) for j= 1,2,...,n are integers. We consider the functions G0(x) = γ−g0(x), Gj(x) = gj(x)−Cj. Let Abe the finite set in Zn+mdefined by the linear constraints (1), (2), (5) and (6). We call B={x∈Zn+m|Gj(x)≤0, j = 0,1,...,m}, 8J. GAGO-VARGAS, M.I. HARTILLO-HERMOSO, J. PUERTO-ALBANDOZ, AND J.M. UCHA-ENR´ IQUEZ and define the penalty function p(x) = m X j=0 max(Gj(x),0). Theorem 3.1 gives a constant µ0which frames a parameter region where the same objective value for the nonlinear problem and the penalized version is obtained. Let c(x) = PiPj(Kij zij +Lijyij) be the cost function and ˆ xthe optimal solution of the relaxed problem (LISP). Proposition 4.1. Let c1:= PiPj(Kij +Lij),Nbe the size of the sample set and ρ0=1 2N. Let (ISP)′be the problem that results from (ISP) after replacing constraint (4) by g0(x)≥γ′, with γ′=⌈Nγ⌉ N−1 2N. The objective values of problems (ISP) and (ISP)′coincide. Moreover, for any µ≥µ0=c1−c(ˆx) ρ0, Problem (ISP) and the penalized version of (ISP)′with objective function c(x) + µp(x)have the same optimal value. Proof. First of all, we observe that g0(x) is a stepwise function that can only assume the values k Nfor k= 1,...,N. Therefore, the set of feasible points of (ISP) that satisfy g0(x)≥⌈Nγ⌉ Nfor any γ∈(0,1) is the same that those that satisfy g0(x)≥⌈Nγ⌉ N−1 2N. This proves the first statement of the proposition. From now on we only consider (ISP)′and its penalized version. In order to apply Theorem 2, we need to prove that ρ0is a valid lower bound of p(x) on A \ B. Indeed, we note that if x6∈ B then either Gj(x)>1 for some j= 1,...,m or g0(x)<⌈Nγ⌉ N−1 2N. In the former case, we have p(x)≥Pm j=1 Gj(x)≥1> ρ0. In the latter case, due to ⌈Nγ⌉−1 N<⌈Nγ⌉ N−1 2N<⌈Nγ⌉ Nand taking into account that g0(x) only assumes values k N, then the constraint g0(x)<⌈Nγ⌉ N−1 2Nis equivalent to g0(x)≤⌈Nγ⌉−1 N. So γ′−g0(x)≥1 2Nand hence p(x)≥G0(x)≥γ′−g0(x)≥1 2N. Theorem 2 provides a constant µ0which depends on the unknown optimal value of (ISP)′. To avoid this inconvenience, we shall replace µ0by something greater. We substitute the unknown value of (ISP)′by c1, a valid upper bound of that value. Hence, we can apply Theorem 3.1 to (ISP)′and its penalized version to conclude the second statement of the proposition.   In spite of the above result, our algorithm does not solve the penalized version. Rather than that we use the order induced by its values to guide the search strategy in our branching tree. However, the above estimation gives a good estimate of µ in our penalty term because a big value of µyields a poor performance of our algorithm. It is clear that c0=c(ˆ x) is a lower bound of minx∈A c(x). In our computational results, after some experimentation, we have used a different estimate of the value of minx∈A\B p(x) than the one given in Proposition 4.1 because of its better performance. Indeed, the values of the functions Gj(x) for j= 1,2,...,m are greater than one, while the values of G0(x) are less than 1. We let ¯ρ=G0(ˆ x), and take ¯µ0=c1−c0 ¯ρ, α =⌊log( c0 ¯µ0 )⌋,and µ=α¯µ0. AN IMPROVED TEST SET APPROACH TO NONLINEAR INTEGER PROBLEMS 9 This quotient is intended to avoid that the coefficient ¯µ0takes on values very large with respect to c0. The penalized cost is then T(x, µ) = c(x) + µp(x). We have considered other penalized cost functions, as the adaptive function described in [11], or the oracle function of [27], but the computational results have been slightly worse than the experiments with the previously defined function T(x, µ). 4.1.2. Computational results. In order to evaluate the performance of the new ordering in the list Pof pending nodes, several families of instances have been considered and compared with the previous walk-back procedure. Both algorithms have been coded in Matlab and run on a Intel Xeon X5660 (2.8 GHz) with 32 GB RAM. The data have been randomly generated with the conditions described in [30, Example 4], that is, n= 7 jobs and m= 4 machines, with capacities (20,20,40,40) (there is a misprint in the original article), mean demands equal to µ= (20,16,12,8,8,4,4), lots M= 2 and covariance matrix V=           36 0 −10.8 0 0 0 7.34 0 23 0 5.9 11.500 −10.8 0 13 0 0 −4.4−4.4 0 5.9 0 6 0 0 0 0 11.5 0 0 23 0 0 0 0 −4.4 0 0 6 0 7.34 0 −4.4 0 0 0 6           . All setup times Sij and production times pij are equal to 1. The setup and production costs (Kij , L′ ij) are provided in Table 1, similarly to [30, Fig. 3], but with the errors corrected (the values were swapped). Table 1. Setup and production costs for n= 7 jobs, m= 4 machines job 1 job2 job 3 job 4 job 5 job 6 job 7 machine 1 5, 3 3, 4 3, 3 1, 2 1, 3 3, 3 5, 4 machine 2 1, 4 1, 4 4, 2 5, 2 3, 4 1, 3 2, 2 machine 3 4, 3 2, 3 5, 2 4, 2 4, 4 2, 2 4, 4 machine 4 5, 3 4, 4 4, 2 4, 3 1, 3 5, 4 2, 2 Different sample sets have been considered, with 250 records, mean µand covariance matrix V. There are no differences among the different samples, we have chosen the more descriptive one. The results appear in Table 2 for different values of the probability γ. We compare the runs for the original walk-back procedure of [30], ordered only by the cost (column “Walk-back” in tables), and the runs with the previously defined penalized cost function (column “WB new ordering” for function T(x, µ)). The column “Total” groups the number of visited nodes and CPU time in seconds. The maximum number of points to be considered has been set to 120000, (a proxy for a CPU time approximately equivalent to 90 minutes of computation), and the label “Max” denotes that the run has reached this limit and it has been terminated. The column “Optimum” groups the number of visited nodes and CPU time to get an optimum. When the process has reached the maximum number of points, the value of the best point found is shown, but it is not possible to guarantee that it is an optimum. If the process has not even found a feasible point we write “NNP”. The computation of the Gr¨obner basis is done 16J. GAGO-VARGAS, M.I. HARTILLO-HERMOSO, J. PUERTO-ALBANDOZ, AND J.M. UCHA-ENR´ IQUEZ •When the walk-back procedure certifies optimality, the optimum is usually found after a large number of nodes explored during the search. The new ordering method finds it very quickly and most part of the time is used to certify optimality. •The required time to certify optimality has highly decreased in all the cases. This fact is due to the decrease in the time needed to reach an initial feasible point. It contributes to speeding up the process through pruning paths. As an illustrative example in Table 5 the average time needed in the walkback algorithm for n= 15, k= 4 is 571.1 s, and the new ordering method requires 432.6 s, an improvement of 25%. •The performance is better than that appeared in [25]. 4.2.9. Comparison with other solvers. In order to perform further comparisons with standard integer nonlinear solvers as Baron ([29]), Couenne ([6]) or Bonmin ([7]), the problem (RP) has been modeled in GAMS format. This encoding has allowed us to launch our data sets in the Neos Server ([14], [18], [16]) under these solvers. We have chosen the 30 instances for systems with n= 8, k = 4 from Table 6, because it contains a good assortment of cases. A summary of the results are reported in Table 7 that shows average time to certify optimality and in Table 8 that shows the time required to reach the optimum (Couenne is not considered in the last Table). Table 7. Comparison of solvers for the case n= 8, k = 4, time to certify optimality WB new ordering Baron Couenne Bonmin Avg. Time 373.4 53.9 1020.9 272.8 N/F 0 0 9 0 Table 8. Comparison of solvers for the case n= 8, k = 4, time to rearch the optimum WB new ordering Baron Bonmin Avg. Time 2.8 25.7 209.4 The columns identify the solvers. The row “Avg. Time” contains the average time of CPU reported by the outputs for the 30 instances. The row “N/F” indicates the number of examples that the system has not been able to finish. We note the following facts: •The packages Bonmin and Baron have, in general, a better performance to certify optimality. •In some cases we have observed that Baron returns a nonoptimal point, although it is very close to the minimum cost. The walk-back procedure and Bonmin always return an actual optimal point. •The walk-back with new ordering found the optimum very quickly, usually under 3 seconds, as can be seen in Figure 2. This chart compares the times used by walk-back, Bonmin and Baron to find the optimum. Couenne reports much larger times. AN IMPROVED TEST SET APPROACH TO NONLINEAR INTEGER PROBLEMS 17 0 200 400 600 Runs Time (s) WB new ordering Bonmin 0 20 40 60 80 Runs Time (s) WB new ordering Baron Figure 2. Time to find an optimum: WB new ordering, Bonmin and Baron A natural question is how to combine the speed of the walk-back procedure with new ordering to find the optimum, and the good approach of a nonlinear solver like Baron to certify optimality of the point. This has been done including the point found by the walk-back procedure in the GAMS code as an initial feasible point. In all the cases, the point is certified with the right optimal cost, and the CPU time is significantly reduced, as it can be seen in Figure 3. 0 50 100 Runs Time (s) Baron Baron+point Figure 3. Time to certify: Baron and Baron with initial point given by WB 5. Conclusions The walk-back procedure is based on extracting certain linear constraints from a nonlinear integer problem for which a test set can be easily computed. Starting from an optimum for the relaxed problem, and using the reverse test set, all the feasible points for the original problem can be reached, and an optimum can be found. The method has been improved by giving a criterion to sort the set of 18J. GAGO-VARGAS, M.I. HARTILLO-HERMOSO, J. PUERTO-ALBANDOZ, AND J.M. UCHA-ENR´ IQUEZ points to be inspected by a penalized cost function. This new cost function trades off the cost of a point against its distance to the feasible region. As a result, new better points are obtained faster than with the order only by cost. Two real design engineering applications have been presented: the scheduling of jobs on parallel machines given restrictions on demands and capacity to minimize costs, and the series parallel redundancy allocation problem given a target reliability. The first problem includes a probabilistic constraint with random variables in the technology matrix, that can not be treated with nonlinear solvers. Our approach deals with any computable constraint. This problem can be reformulated as a linear integer program, and still in certain cases our approach is competitive with the best MILP solvers. The second application considered contains a nonlinear restriction given by the reliability function. This nonlinear model can be compared with other solvers with remarkable results. We conclude, through extensive computations, that both applications show that the combination of the walk-back procedure and an improved ordering is a promising tool as an exact method to solve nonlinear integer programs. Acknowledgments This paper has been partially supported by Junta de Andaluc´ıa under grant FQM-5849, and Ministerio de Ciencia e Innovaci´on MTM2010-19336, MTM2010- 19576, MTM2013-46962-C2-1-P and FEDER. The authors would like to thank the anonymous reviewers for their valuable comments and suggestions to improve the quality of the paper. References [1] 4ti2 team. 4ti2—a software package for algebraic, geometric and combinatorial problems on linear spaces. Available at www.4ti2.de. [2] Achterberg, T. SCIP: Solving constraint integer programs. Mathematical Programming Computation, 1(1):1–41, (2009). [3] Ahmadizar, F. and Soltanpanah, H. Reliability optimization of a series system with multiplechoice and budget constraints using an efficient ant colony approach. Expert Systems with Applications, 38(4):3640–3646, (2011). [4] Ahmed, S. and Shapiro, A. Solving Chance-Constrained Stochastic Programs via Sampling and Integer Programming. Available at www2.isye.gatech.edu/people/faculty/Shabbir- Ahmed/cctutorial.pdf. Accessed 30 April 2014. [5] Avis, D. and Fukuda, K. Reverse search for enumeration. Discrete Applied Mathematics, 65(1-3):21–46, (1996). [6] Belotti, P. Couenne. Available at www.coin-or.org/Couenne/. [7] Bonami, P. and contributors. Bonmin — basic open-source nonlinear mixed integer programming. Available at https://projects.coin-or.org/Bonmin. [8] Cbc (Coin-or branch and cut). Available at https://projects.coin-or.org/Cbc. [9] Castro, F., Gago, J., Hartillo, I., Puerto, J., and Ucha, J.M. An algebraic approach to integer portfolio problems. European Journal of Operational Research, 210(3):647–659, (2011). [10] Chern, M.S. On the computational-complexity of reliability redundancy allocation in a series system. Operations Research Letters, 11(5):309–315, (1992). [11] Coit, D. W., Smith, A.E. and Tate, D.M. Adaptive penalty methods for genetic optimization of constrained combinatorial problems. INFORMS J.Comput., 8(2):173–182, (1996). [12] Coit, D.W. and Smith, A.E. Reliability optimization of series-parallel systems using a genetic algorithm. IEEE Transactions on Reliability, 45(2):254–260, (1996). [13] Gurobi Optimization, Inc. Gurobi Optimizer Reference Manual. Available at http://www.gurobi.com. [14] Czyzyk, J., Mesnier, M.P. and Mor´e, J.J. The NEOS server. IEEE Computational Science and Engineering, 5(3):68–75, (1998). AN IMPROVED TEST SET APPROACH TO NONLINEAR INTEGER PROBLEMS 19 [15] Djerdjour, M. and Rekab, K. A branch and bound algorithm for designing reliable systems at a minimum cost. Applied Mathematics and Computation, 118(2-3):247–259, (2001). [16] Dolan, E. The NEOS server 4.0 administrative guide. Technical Report Technical Memorandum ANL/MCS-TM-250, Mathematics and Computer Science Division, Argonne National Laboratory, (2001). [17] Gago, J., Hartillo, I., Puerto, J., and Ucha, J.M. Exact cost minimization of a series-parallel system. Computers and Operations Research, 40(11):2752–2759, (2013). [18] Gropp, W. and Mor´e, J.J. Optimization environments and the NEOS server, pages 167– 182. Approximation theory and optimization (Cambridge, 1996). Cambridge Univ. Press, Cambridge, (1997). [19] Hemmecke, R. and Schultz, R. Decomposition of test sets in stochastic integer programming. Mathematical Programming, 94(2-3):323–341, (2003). [20] Li, D. and Sun, X. Nonlinear integer programming. Springer, New York, (2006). [21] Li, Q., Guo, Y.K., Darlington, J. and Ida, T. Minimised geometric Buchberger algorithm for integer programming. Annals of Operations Research, 108(1-4), (2001). [22] Luedtke, J., Ahmed, S., and Nemhauser, G. An integer programming approach for linear programs with probabilistic constraints. Lecture Notes in Computer Science 4513, pages 410–423, Berlin, 2007. Springer-Verlag. [23] MOSEK ApS, Fruebjergvej 3, Boks 16, 2100 Copenhagen O, Denmark. The MOSEK optimization tools manual version 6.0, 2009. Available at http://www.mosek.com/. [24] Ouzineb, M., Nourelfath, M., and Gendreau, M. Tabu search for the redundancy allocation problem of homogenous series-parallel multi-state systems. Reliability Engineering and System Safety, 93(8):1257–1272, (2008). [25] Ruan, N. and Sun, X. An exact algorithm for cost minimization in series reliability systems with multiple component choices. Applied Mathematics and Computation, 181(1):732–741, (2006). [26] Ruszczy´nski, A. Probabilistic programming with discrete distributions and precedence constrained knapsack polyhedra. Mathematical Programming, 93 (2), 195–215 (2002). [27] Schlueter, M. and Gerdts, M. The oracle penalty method. Journal of Global Optimization, 47(2):293–325, (2010). [28] Sturmfels, B. Gr¨obner bases and convex polytopes, volume 8. American Mathematical Society, Providence, RI, (1996). [29] Tawarmalani, M. and Sahinidis, N.V. A polyhedral branch-and-cut approach to global optimization. Mathematical Programming, 103(2):225–249, (2005). [30] Tayur, S.R., Thomas, R.R., and Natraj, N.R. An algebraic-geometry algorithm for scheduling in presence of setups and correlated demands. Mathematical Programming, 69(3):369–401, (1995). Appendix A. Proof of Theorem 4.1 Let Abe the matrix associated with the restrictions for Problem (LSP-II). Then Acan be divided in five blocks according to the five sets of variables yij , aij , zij, bij and cij. Let Idenote the identity matrix of size mn and Odenote a matrix with all entries zero. Let Dbe the block diagonal (0,1)-coefficient matrix associated with constraints (7). Let −MI denote the coefficient matrix of the variables zij, i = 1,...,n,j = 1, . . . , m in constraints (8). We can then write Aas the (n+3mn)×5mn integer matrix A=    DO O O O II−MI O O O O I I O −IOIOI     We note t= (ty, ta, tz, tb, tc), S ={t|At = 0}and the block tyrepresents the variables (Y11,...,Ymn), noted in capital letters. Similar for the other blocks. Consider the following cases: 20J. GAGO-VARGAS, M.I. HARTILLO-HERMOSO, J. PUERTO-ALBANDOZ, AND J.M. UCHA-ENR´ IQUEZ (1) K1={t∈S|ty=tb= 0}. Let S′= ker A′, where A′=  I−MI O OIO OI I  →  I−MI O OIO O O I , which is not singular. Then K1= 0 and there is no binomial xα−xβsuch that contains only the variables aij , zij , cij . (2) K2={t∈S|tb= 0}. Let S′′ = ker A′′, where A′′ =    DO O O I I −MI O O O IO −IOI I    . By the previous case, we can assume ty6= 0. The rows in the third block imply that tz= 0. Then t∈K2if and only if (ty, ta, tc)∈S′′′ = ker A′′′, for A′′′ =  DO O I I O −IOI . Let Yip be the leftmost nonzero component in ty. We can assume Yip >0, because S′′′ is a vector space. The row block corresponding to Dimplies that there exists q > p such that Yiq <0. Then yip divides xt+, yiq divides xt−. The second block of rows in A′′′ implies Aip =−Yip <0, Aiq =−Yiq >0. From the third block of rows, Cip =Yip >0, Ciq =Aiq <0. Then yipaiqcip divides xt+, yiqaipciq divides xt−. The initial term in xt+−xt−with respect to >is xt+because yip divides xt+and yip it the greatest monomial in the binomial. Then the initial term of (12) yipaiqcip −yiqaipciq divides the initial term of xt+−xt−. (3) Let S= ker A. By (a) and (b) we suppose tb6= 0. Let Brs the leftmost nonzero component in tb. As in the previous case, take Brs >0. Because b is the set of greatest variables for our term order, xt+is the leading term. We also have Zrs =−Brs <0. (a) If Yrs ≥0 then −Yrs +Zrs +Crs = 0 then brscrs divides xt+ and the initial term of (13) brscrs −zrsaMr rs . divides xt+. (b) Yrs <0. Then there exists Yrq, s 6=qwith Yrq =−Yrs >0 . AN IMPROVED TEST SET APPROACH TO NONLINEAR INTEGER PROBLEMS 21 (i) If Crq >0, then brsyrqcrq divides xt+. The initial term of (14) brsyrqcrq −zrsyrsaMr−1 rs arq, s 6=q. divides the initial term of xt+−xt−. (ii) If Crq ≤0 then 0 < Yrq ≤Zrq, which implies Brq <0 so s < q, because Brs was the first nonzero element. The rows in the second block imply that Yrq +Arq −MrZrq = 0,so Arq =MrZrq −Yrq ≥MrYrq −Yrq = (Mr−1)Yrq. Then the initial term of (15) brsyrqzrqaMr−1 rq −zrsyrsbrqaMr−1 rs , s < q divides the initial term of xt+−xt−. We have reviewed all the cases, so the sets (12), (13), (14) and (15) form the reduced Gr¨obner basis. Depto. de ´ Algebra, Universidad de Sevilla. Apdo. 1160, E-41080 Sevilla (Spain) E-mail address:[email protected] Dpto. de Matem´ atica Aplicada I, E.T.S. de Ingenier´ ıa Inform´ atica, Av. Reina Mercedes, s/n, 41012 Sevilla, Spain E-mail address:[email protected] Dpto. de Estad´ ıstica e I.O., Facultad de Matem´ aticas, apdo. 1160, 41080 Sevilla, Spain E-mail address:puer[email protected] Depto. de ´ Algebra, Universidad de Sevilla. Apdo. 1160, E-41080 Sevilla (Spain) E-mail address:[email protected]