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

Comput Optim Appl (2015) 62:565–588 DOI 10.1007/s10589-015-9739-3 An improved test set approach to nonlinear integer problems with applications to engineering design J. Gago-Vargas ·I. Hartillo ·J. Puerto · J. M. Ucha Received: 2 May 2014 / Published online: 14 March 2015 © Springer Science+Business Media New York 2015 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 exactmethodto obtainanoptimalpoint of an integerprogrammingproblemwith 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 J. Gago-Vargas ·J. M. Ucha Dpto. de Álgebra, Univ. de Sevilla, Apdo. 1160, 41080 Sevilla, Spain e-mail: [email protected] J. M. Ucha e-mail: [email protected] I. Hartillo (B) Dpto. de Matemática Aplicada I, E.T.S. de Ingeniería Informática, Univ. de Sevilla, Avda. de Reina Mercedes, s/n, 41012 Sevilla, Spain e-mail: [email protected] J. Puerto Dpto. de Estadística e Investigación Operativa, Facultad de Matemáticas, Univ. de Sevilla, C/ Tarfia, s/n, 41012 Sevilla, Spain e-mail: [email protected] 123 566 J. Gago-Vargas et al. applicability of this approach to deal with other problems with similar characteristics or to combine it with mainstream solvers to certify optimality. Keywords Non-linear integer programming ·Test set ·Gröbner basis · Chance constrained programming 1 Introduction The ever-increasing demand on operations researchers to lower production costs or to increasebenefitshaspromptedthespecializedcommunitytolookforrigorousmethods 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 occasionstheresultingmodels contain both nonlinearities and discrete variableswhich 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 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 [29] is the insertion in the pending nodes list by the ascending penalized cost. We will see in Sect.4the 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 schedulingjobs on parallel machines givenrestrictionsondemandsand 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 123 An improved test set approach to nonlinear integer problems 567 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 Sect.2the 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. Section4 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 [17,24,29] 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, Sect.5contains 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 [27,29]. 2.1 Test sets Consider an integer linear programming problem: (LIP(b)) min c(x)=ct·x s. t. A·x=b, x∈ZN +, where A∈Zd×N,b∈Zd,c∈RN. The notation (LIP(b)) denotes the integer linear programming problem with right-hand-side fixed to b.(LIP)denotes the set of all the integer linear programming problems obtained by varying the right-hand-side vector b,fixedAand the cost function c.Letπbe the map defined by π(x)=Ax.Givena vector b∈Zd,thesetπ−1(b)={u∈ZN +:π(u)=b}is the fiber of (LIP)over b. Fixed a term order <c, associated with the cost, defined as usual in [27], there exists a unique optimum βfor (LIP(b)).AsetG<c⊂ZNis a test set for the family of integer linear problems (LIP)with respect to the matrix Aand the order <cif – for each nonoptimal point αof (LIP(b)), there exists g∈G<csuch that α−gis a feasible solution for (LIP(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 (LIP)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: 123 568 J. Gago-Vargas et al. – 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=xα−xβ:Aα=Aβ,α,β∈ZN +⊂Q[x1,...,xN]. Given a term order <ccompatible with the cost function c(x), a reduced Gröbner basis G<cof IAyields a test set G<cfor (LP)(cf. [27]). If the reduced Gröbner 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∈Aare linear, and a test set can be obtained in practice with respect to the linear cost function and these restrictions; – the constraints x∈Bcan 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. 123 An improved test set approach to nonlinear integer problems 569 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 walk-back procedure, described in [29]. The above argument can be summarized in the following result: Theorem 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 [9,17,29]. 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∈Ato 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 [27], 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öbner 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 Sect. 3. 123 570 J. Gago-Vargas et al. 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. 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)=0forx∈A∩Band p(x)≥ γ>0ifx∈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 2 (cf. [20]Th. 2.11)Suppose that A\B=∅. Let c0be a lower bound of c(x)in Aand ρ>0be a lower bound of minx∈A\Bp(x). Then, for any μ>μ 0with μ0=opt(P)−c0 ρ we have opt(P)=opt(Ppen(μ)). 3.2 Ordering nodes by penalized cost Following the notation of [29] we describe in this subsection our algorithm. The walkbackprocedureprovidesalistoffeasiblepoints,notedP,fortherelaxedlinearproblem 123 An improved test set approach to nonlinear integer problems 571 (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 [17,29]. 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 [29]: 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 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: 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 [29] 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, 123 572 J. Gago-Vargas et al. –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 [29]is: (SP)minimize  i j (Kijzij +Lijyij) subject to m  j=1 yij =Mi,i=1,2,...,n,(1) Mizij ≥yij,i=1,2,...,n,j=1,2,...,m,(2) gj(z,y)= n  i=1 pij ˆ Di/Miyij + n  i=1 Sijzij ≤Cj,j=1,2,...,m,(3) g0(z,y) =Prob n  i=1 pij (Di/Mi)yij + n  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). The model (SP) admits a family of optimality cuts: zij ≤yij,i=1,...,n,j=1,...,m.(6) 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 Eqs.(3) and (4) are checked for every point and test feasibility for problem (ISP). Similarly to [29], the linear problem with restrictions (1), (2), (5) and (6) can be expressed with equalities by adding slack variables: 123 An improved test set approach to nonlinear integer problems 573 (LISP −II)minimize  i j (Kijzij +Lijyij) subject to m  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 3 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≤p= q≤m, bipyiqziqaMi−1 iq −zipyipbiqaMi−1 ip ,1≤p<q≤m, is the reduced Gröbner basis of the toric ideal corresponding to (LISP-II) with respect to the lexicographical order b >y>z>a>c. The letter y denotes the set of variables (y11,...,ymn)and similarly a,z,b and c denote the other sets of variables. Within a block, the variables are sorted lexicographically according to index. The proof is written in “Appendix” section. 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 [29]. 4.1.1 New penalized cost function Following the notation of Sect.3,letxbe 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}, 123 580 J. Gago-Vargas et al. Dn×N=⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ k1  1...1 k2  0...0... kn  0...0 0...01...1... 0...0 ... 0...00...0... 1...1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,Ck i×ki= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1−100··· 0 10−10··· 0 100−1··· 0 . . .. . ..... . . 1000···−1 01−10··· 0 010−1··· 0 . . .. . .. . ..... . . 0100··· −1 . . .. . .. . . 00001−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ where k i=ki 2. Under these conditions, we have the following result: Theorem 4 The rows of the matrix ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ IN(Dn×N)t−IN Ck 1×k1000 0−Ck 1×k1000 0Ck 2×k200 0 0 −Ck 2×k200 . . ..... . .... 00···Ck n×kn0000−Ck n×kn ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , form a test set of (LIRP)with respect to the monomial order defined by the cost function, where the leftmost block corresponds to the variables x11,...,x1k1,x21,...,xnkn, the center block to the variables d1,...,dNand the rightmost block to the variables t11,...,t1k1,t21,...,tnkn. For a complete proof of this result see [17]. 4.2.3 Penalized cost function The walk-back process for this problem is greatly improved by the knowledge of a feasible point y0of (RP), which provides a good initial upper bound. There are lots of heuristic methods to obtain such a point. In this case, a greedy algorithm similar to [15]or[24] is used, which obtains a feasible point y0, with a cost ij cijy0 ij =c0and a reliability R0. Obviously, the optimal solution of (RP) has a cost less than or equal to c0.Letβbe the optimum of (LIRP) with reliability Rβ. The penalized cost function can be defined in the following way, inspired by Theorem 2: T(x,μ)=ct·x+μ·max{0,(R0−R(x))},where μ=cY−ct·β R0−Rβ , 123 An improved test set approach to nonlinear integer problems 581 Table 4 No correlation, high values of rij ∈[0.99,0.998],cij ∈[10,20],uij =4,R0=0.90, average time in seconds and nodes Walk-back WB new ordering nkTime Nodes >LimitTimeNodes>Limit 17 4 0.2581 0.1231 19 2 40.1 3389 1 9.1 2255 1 20 2 173.7 7254 6 72.7 7177 5 with cYthe cost of the best point found for (RP). Initially, cY=c0. An improvement of the method consists in changing after a short period of time the penalized function. The list Pof pending nodes is ordered according to the penalized cost function for a fixed amount of nodes L. After that, when it is likely a good feasible point has been found, the ordering is changed to T2(x)=−ct·x, that is, a “Depth first search” [5]. This assures a fast search for the optimum and a set Pwith a medium size. The parameter used is described in the next subsection. 4.2.4 Computational results In all cases nand krefer to the number of subsystems, and the number of different components in each subsystem, respectively. For simplicity, the same number of components has been considered in each system ki=k. The complexity of the problem depends on iki, different values of kido not increase the problem size. The term “Walk-back” collects the results obtained with the algorithm published in [17]. The term “WB new ordering” summarizes the data for the new ordering given by the penalized cost function. For all the examples a summary table reports the average of a total of 30 instances run for each case. The columns show the average CPU time in seconds and nodes visited, both data computedfor thecaseswhen both algorithmscertifiedoptimality visiting80,000nodes or less. The parameter Ldescribed in the previous subsection for the improvement is set to 2000. The column “>Limit” counts how many cases did not reach optimality. 4.2.5 No correlation between cost and reliability, high reliability values Inthefirstcasethesampleisrandomlygeneratedfromauniformdistribution,asin[24], withuij =4,cij ∈[10,20],andvaluesofreliabilityverycloseto1:rij ∈[0.99,0.998]. The reliability objective is set to R0=0.90. There is no correlation between cost and reliability of components, so it is likely that certain components are easily discarded. In Table4the results of the experiment are summarized. 4.2.6 Correlation between cost and reliability, high reliability values In this case we use the hypothesis that the greater the reliability component, the greater cost. The size of the systems have changed, due to the increased difficulty of the cases. 123 582 J. Gago-Vargas et al. Table 5 Correlation, high values of rij ∈[0.99,0.998],cij ∈[10,20],uij =4,R0=0.90, average time in seconds and nodes Walk-back WB new ordering nkTime Nodes >Limit Time Nodes >Limit 15 3 322.8 6965 1 42.3 4795 0 15 4 571.117,629 13 432.621,843 8 17 2 92.7 6465 1 15.4 3813 1 Table 6 Correlation, low values rij ∈[0.90,0.99],cij ∈[10,20],uij =4,R0=0.90, average time in seconds and nodes Walk-back WB new ordering nkTime Nodes >Limit Time Nodes >Limit 75 23.6 6157 0 17.9 5301 0 8 4 593.537,728 7 373.434,754 1 The ranges of reliability and cost are similar to the first family of examples. Table5 contains the summarized data. 4.2.7 Correlation between cost and reliability, smaller reliability values The last family of examples has the reliability and cost also correlated. When the range of the reliability of components is changed, the size of the problem increases very quickly. However if the range of cost is changed, the difficulty of the problem is the same. Here we consider smaller reliability values of components rij in the range [0.90,0.99]. Table6summarizes the collected data. 4.2.8 Comments about the tables We highlight the following facts about the improvements: – In all the cases when the walk-back algorithm did not find any improvement to the initial estimation, the new ordering procedure gives a better point. – In almost the entire set of examples the optimum has been found with the new ordering algorithm in less than 1s. – 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 examplein Table 5theaveragetime needed in the walk-backalgorithm 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 [24]. 123 An improved test set approach to nonlinear integer problems 583 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 4.2.9 Comparison with other solvers In order to perform further comparisons with standard integer nonlinear solvers as Baron [28], Couenne [6]orBonmin [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,16,18] under these solvers. We have chosen the 30 instances for systems with n=8,k=4 from Table6, because it contains a good assortment of cases. A summary of the results are reported in Table7that shows average time to certify optimality and in Table 8that shows the time required to reach the optimum (Couenne is not considered in the last Table). 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 3s, as can be seen in Fig. 2. This chart compares the times used by walk-back, Bonmin and Baron to find the optimum. Couenne reports much larger times. 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 rightoptimalcost, and the CPU time is significantly reduced, as it can be seen in Fig.3. 123 584 J. Gago-Vargas et al. 0 200 400 600 Runs Time (s) WB new ordering Bonmin 0 20 40 60 80 Runs Time (s) WB new ordering Baron Fig. 2 Time to find an optimum: WB new ordering, Bonmin and Baron 0 50 100 Runs Time (s) Baron Baron+point Fig. 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 points to be inspected bya penalizedcostfunction. This newcost function tradesoffthe costofa pointagainst 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 123 An improved test set approach to nonlinear integer problems 585 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ón 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. Appendix: A Proof of Theorem 3 Let AbethematrixassociatedwiththerestrictionsforProblem(LSP-II).Then Acanbe 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 OOI 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.Similarfortheotherblocks. Consider the following cases: 1. K1={t∈S|ty=tb=0}.LetS=ker A, where A=⎛ ⎝ I−MI O OIO OII ⎞ ⎠→⎛ ⎝ I−MI O OIO OO 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}.LetS =ker A, where A =⎛ ⎜ ⎜ ⎝ DOOO II−MI O OO IO −IOII ⎞ ⎟ ⎟ ⎠. 123 586 J. Gago-Vargas et al. By the previous case, we can assume ty= 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 =⎛ ⎝ DOO IIO −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>psuch 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 yipaiqcip −yiqaipciq (12) divides the initial term of xt+−xt−. 3. Let S=ker A. By (a) and (b) we suppose tb= 0. Let Brs the leftmost nonzero component in tb. As in the previous case, take Brs >0. Because bis 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 brscrs −zrsaMr rs .(13) divides xt+. (b) Yrs <0. Then there exists Yrq,s= qwith Yrq =−Yrs >0. i. If Crq >0, then brsyrqcrq divides xt+. The initial term of brsyrqcrq −zrsyrsaMr−1 rs arq,s= q.(14) divides the initial term of xt+−xt−. 123 An improved test set approach to nonlinear integer problems 587 ii. If Crq ≤0 then 0 <Yrq ≤Zrq, which implies Brq <0sos<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 brsyrqzrqaMr−1 rq −zrsyrsbrqaMr−1 rs ,s<q(15) 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öbner basis. References 1. 4ti2 team. 4ti2—a software package for algebraic, geometric and combinatorial problems on linear spaces. www.4ti2.de 2. Achterberg, T.: SCIP: solving constraint integer programs. Math. Program. Comput. 1(1), 1–41 (2009) 3. Ahmadizar, F., Soltanpanah, H.: Reliability optimization of a series system with multiple-choice and budget constraints using an efficient ant colony approach. Expert Syst. Appl. 38(4), 3640–3646 (2011) 4. Ahmed, S., Shapiro, A.: Solving chance-constrained stochastic programs via sampling and integer programming. www2.isye.gatech.edu/people/faculty/Shabbir-Ahmed/cctutorial. Accessed 30 April 2014 5. Avis, D., Fukuda, K.: Reverse search for enumeration. Discret. Appl. Math. 65(1–3), 21–46 (1996) 6. Belotti, P. Couenne. www.coin-or.org/Couenne/ 7. Bonami, P., contributors.: Bonmin—basic open-source nonlinear mixed integer programming. https:// projects.coin-or.org/Bonmin 8. Cbc (Coin-or branch and cut). https://projects.coin-or.org/Cbc 9. Castro, F., Gago, J., Hartillo, I., Puerto, J., Ucha, J.M.: An algebraic approach to integer portfolio problems. Eur. J. Oper. Res. 210(3), 647–659 (2011) 10. Chern, M.S.: On the computational-complexity of reliability redundancy allocation in a series system. Oper. Res. Lett. 11(5), 309–315 (1992) 11. Coit, D.W., Smith, A.E., 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., Smith, A.E.: Reliability optimization of series-parallel systems using a genetic algorithm. IEEE Trans. Reliab. 45(2), 254–260 (1996) 13. Gurobi Optimization, Inc.: Gurobi optimizer reference manual. http://www.gurobi.com 14. Czyzyk, J., Mesnier, M.P., Moré, J.J.: The NEOS server. IEEE Comput. Sci. Eng. 5(3), 68–75 (1998) 15. Djerdjour, M., Rekab, K.: A branch and bound algorithm for designing reliable systems at a minimum cost. Appl. Math. Comput. 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., Ucha, J.M.: Exact cost minimization of a series-parallel system. Comput. Oper. Res. 40(11), 2752–2759 (2013) 18. Gropp, W., Moré, J.J.: Optimization Environments and the NEOS Server. Approximation theory and optimization. Cambridge University Press, Cambridge (1996) 19. Hemmecke, R., Schultz, R.: Decomposition of test sets in stochastic integer programming. Math. Program. 94(2–3), 323–341 (2003) 20. Li, D., Sun, X.: Nonlinear Integer Programming. Springer, New York (2006) 123 588 J. Gago-Vargas et al. 21. Luedtke, J., Ahmed, S., Nemhauser, G.: An Integer Programming Approach for Linear Programs with Probabilistic Constraints. Lecture notes in computer science 4513. Springer, Berlin (2007) 22. MOSEK ApS, Fruebjergvej 3, Boks 16, 2100 Copenhagen O, Denmark. The MOSEK optimization tools manual version 6.0, 2009. http://www.mosek.com/ 23. Ouzineb, M., Nourelfath, M., Gendreau, M.: Tabu search for the redundancy allocation problem of homogenous series-parallel multi-state systems. Reliab. Eng. Syst. Saf. 93(8), 1257–1272 (2008) 24. Ruan, N., Sun, X.: An exact algorithm for cost minimization in series reliability systems with multiple component choices. Appl. Math. Comput. 181(1), 732–741 (2006) 25. Ruszczy´nski, A.: Probabilistic programming with discrete distributions and precedence constrained knapsack polyhedra. Math. Program. 93(2), 195–215 (2002) 26. Schlueter, M., Gerdts, M.: The oracle penalty method. J. Global Optim. 47(2), 293–325 (2010) 27. Sturmfels, B.: Gröbner Bases and Convex Polytopes, vol. 8. American Mathematical Society, Providence (1996) 28. Tawarmalani, M., Sahinidis, N.V.: A polyhedral branch-and-cut approach to global optimization. Math. Program. 103(2), 225–249 (2005) 29. Tayur, S.R., Thomas, R.R., Natraj, N.R.: An algebraic-geometry algorithm for scheduling in presence of setups and correlated demands. Math. Program. 69(3), 369–401 (1995) 123 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.