MIQL: A Fortran Subroutine for Convex Mixed
Contents
1. Computational study of a family of mixed integer quadratic pro gramming problems Columbia University New York USA Bixby R E Fenelon M Gu Z Rothberg E Wunderling R 2004 Mixed integer programming A progress report in Martin Gr tschel ed The sharpest cut The impact of Manfred Padberg and his work MPS SIAM Series on Optimization Vol 4 Boland N L 1997 A dual active set algorithm for positive semi definite quadratic programming Mathematical Programming Vol 78 1 27 26 14 15 16 18 19 20 21 22 23 24 25 26 21 Bussieck M R Drud 5 Meeraus A 2007 MINLPLib A collection of test mod els for mixed integer nonlinear programming GAMS Development Corp Washington D C USA Ceria S Pataki G 1998 Solving integer and disjunctive programs by lift and project Lecutre Notes in Computer Science Vol 1412 271 283 Ceria S Soares J 1997 Disjunctive cut generation for mized 0 1 programs duality and lifting GSB Columbia University New York USA Dakin R J 1965 A tree search algorithm for mixed integer programming problems Computer Journal Vol 8 250 255 Exler O Schittkowski K 2006 MISQP A Fortran implementation of a trust region SQP algorithm for mixed integer nonlinear programming User s guide version 2 0 Report Department of Computer Science University of Bayreuth Germany Goldfarb D I2. NI IND C D A B XL XU X ACC EPS ROPT LOPT IOUT IFAIL IPRINT RW LRW IW LIW LW LLW SS ALS Hj Parameter Definition Input parameter for the total number of constraints including equal ity constraints ME Input parameter for the number of equality constraints MMAX Input parameter defining the row dimension of array A containing linear constraints must be at least M CMAX CMAX maximum number of cuts 16 NMAX MNN IND NI D N A MMAX N B MMAX XL N XU N X N U MNN EPS IOPT 20 IOPT 1 IOPT 2 Input parameter for the number of optimization variables Input parameter defining the row dimension of C at least one and greater or equal to N Input parameter equal to MMAX N N when calling QL dimen sion of U Input parameter for the number of integer variables Integer input vector for indices of integer variables Double precision input matrix for objective function matrix Double precision input vector for constant vector of objective func tion Double precision input and output matrix for linear constraints first ME rows for equality then M ME rows for inequality con straints last MMAX M rows set to cut constraints on return Double precision input and output vector for constant values of linear constraints and the cut constraints on return Double precision input
3. whereas the variables u Uo v Vo ensure the validity of this constructed cut The cutting plane obtained by solving CGLP only regards the integer condition on variable yj Therefore this cut can be strengthened by a procedure of Balas and Jeroslow 7 that also considers the integrality condition on other variables The result of solving CGLP is one single violated cut There exists various suggestions how to generate more than just one cut out of a disjunctive relaxation corresponding to disjunction 6 see Perregaard 28 for details Furthermore it is possible to restrict the disjunctive relaxation to a subspace of the original one neglecting all variables that are fixed at a bound This reduces the size of CGLP significantly The resulting cutting plane of the reduced cut generating linear program can be lifted to be valid in the full space see Perregaard 28 In our current implementation of MIQL we replace the linear program 8 by a quadratic min VR a ATu uper gt 0 8 E R xR a ATu gt 0 R x B b u ugly 0 v vo E R x R 2 bv vof 0 A 10 oy us Uo Du 0 1 ug 9 00 gt 0 where e gt 0 R is a suitable small value Program 10 is solved by the quadratic solver QL of Schittkowski 31 Furthermore our implementation of disjunct
4. Feasible solution found but the maximum number of nodes reached Termination tolerance ACC negative Length of a working array RW IW or LW too small Parameters N M MNN ME or NI incorrect Option of IOPT incorrectly set Lagrangian multipliers could not be calculated IOUT or IPRINT incorrectly set Lower variable bound greater than upper one ON no w Input parameter defining the output level 0 No output to be generated 1 Output of relaxed solution and final performance summary 2 Intermediate reduced i e one line output for all generated subproblems i e all nodes including node status fractional feasible node B new best integer feasible solution I infeasible node M marked node 19 RW LRW LRW IW LIW LIW LW LLW LLW 3 Complete output for all nodes e g with bounds Double precision working array of length LRW Input parameter for length of RW must be at least 5 NMAX NMAX 2 7 NI 10 NMAX 16 MMAX 15 N 3 MAXNDS 14 Integer working array of length LIW Input parameter for length of IW must be at least 4 MMAX NMAX 5 NI 6 MAXNDS 6 N 11 Logical working array of length LLW Input parameter for length of LW must be at least MMAX N 3 NI 3 The continuous quadratic programming solver QL see Schittkowski 31 is extended to satisfy the requirements of Section 4 Usage CALL QLX M ME MMAX NMAX MNN C D
5. Vi eIg 4 i j i icr jeJ The inequality constraints of the original problem 1 are transformed into equalities by introducing slack variables within the CMIR cut generation procedure The cut generation phase tries to find violated complemented mixed integer rounding inequalities Fp Eui Fp 9 rl 5 2 mg 2 EE 1 f for XMK and to improve them and U are a partition of I and 6 Q 0 Furthermore we define es 04 fa c F d d 2 d d 20 which is called the MIR function with 0 lt 1 and max d 0 Marchand and Wolsey 26 show that the complemented mixed integer rounding inequal ity 5 is always valid for We also observe that X is a relaxation of the original mixed integer feasible region X and that every cut valid for is therefore valid for X If any complemented mixed integer rounding inequality 5 is violated by the relaxed solution of MIQP 1 then we have fund a cutting plane By varying the partition T and U and by trying out different values for 6 the quality of CMIR cuts can be improved Different aggregated constraints can lead to different cuts The drawback is that there is no guarantee to find violated cuts by this procedure But computational experience shows that it is often possible to find deep violated cuts Wolters 35 reports that complemented mixed integer rounding cuts are the most
6. actual one i e the new quadratic program is identical to the previous one apart from one additional box constraint subject to an integer variable Since the primal dual method of Goldfarb Idnani successively adds constraints to the active set we need in most cases only one additional iteration to get a new optimal solution Even if more iterations are needed almost always we save some calculation time compared to solving the complete quadratic program 9 2 The new node has the same parent node as the actual one 1 both problems differ only in one box constraint for an integer variable The algorithm identifies the solution of a quadratic program by an linear independent active set A problem is solved by adding constraints to and eliminating constraints from the active set until all constraints are satisfied and each iterate is dual feasible To enable a warm start in the present situation one has to ensure dual feasibility first This is achieved by saving the active constraints of the common parent node They are compared to the ones of the child ie the last solved quadratic program All other active constraints are eliminated Then a dual feasible point is obtained and a new iteration of the algorithm is started by searching for a violated constraint Time savings are achieved since in most cases only few active constraints have to be deleted and added afterwards 3 The new node is the child of a node which has the same parent nod
7. integer variable from the solution of the relaxed subproblem with smallest distance from next integer value 2 Searching the tree for a free node from where branching i e the generation of two new subproblems is started a best of two The optimal objective function values of the two child nodes are compared and the node with a lower value is chosen If both are leafs i e if the corresponding solution is integral or if the corresponding problem is infeasible or if there is already a better integral solution strategy best of all is started b best of all All free nodes are compared and a node with lowest objective function value is selected c depth first This search strategy selects a child node whenever possible If the node is a leaf the best of all strategy is applied The depth first strategy has the advantage that the quadratic problem of a child node can be solved faster than any other one since the only difference is one additional box constraint Moreover the search tree usually is smaller i e there are less unexplored nodes The disadvantage is that the child node is arbitrarily selected leading to a larger number of subproblems to be solved Until the first integer feasible solution is found the upper bound is infinity After selecting an integer variable for branching and an unexplored not to continue the optimal solutions of the corresponding relaxed quadratic programs are computed This process is iterated
8. matrix C n an n dimensional vector d an by n matrix a1 am and an m vector b Lower and upper bounds for the continuous and integer variables r u and y respectively are separately handled To ensure a consistant notation we introduce the two index sets J and J with UJ 1 n and IN J 0 where y N for k I I n and x R for 1 J J ne The objective function of 1 is denoted by fee setane 5 er 5 2 y If n gt 0 the mixed integer quadratic program 1 is solved by a branch and cut method Branch and cut combines two solution strategies a method for generating cutting planes and a branch and bound procedure We apply cut generation only at the root node of the branch and bound search tree i e before starting the enumeration process For the present program version we have implemented two different cut generation procedures known as complemented mixed integer rounding cutting planes and disjunctive cutting planes This manual is organized as follows The subsequent section introduces the basic con cept of a branch and bound method and describes the different options of the underlying Fortran subroutine Section 3 deals with the generation of cutting planes whereas Section 4 describes more details of the implementation especially the warm start features Section 5 summarizes some numerical results to compare various options of the code Further im plem
9. parameters N 5 M E ME 0 MNN MMAX N N NI N ACC 1 0D 10 EPS 1 0D 15 IOPT 1 1 IOPT 2 3 IOPT 3 MAXNDS IOPT 4 1000 IOPT 5 1 IOPT 6 1 IOPT 7 1 IOPT 8 3 IOPT 9 1 IOPT 10 1 11 2 IOPT 12 N IPRINT 2 IOUT 6 Set problem data DO I 1 N DO J 1 N C I J 0 000 ENDDO 1 1 0DO 1 1 0 000 1 100 0D0 XU I 100 000 ENDDO DO I 1 NI 1 I ENDDO A 1 1 7 56D0 A 1 5 0 50D0 B 1 39 1D0 D 1 21 9800 D 2 1 26 0 D 3 61 39D0 4 5 3 D 5 101 300 Calling MIQL CALL MIQL M ME MMAX N NMAX MNN NI IM C D A B XL XU X U FMIN ACC EPS IOPT ROPT LOPT IOUT IFAIL IPRINT W LWMAX IW LIWMAX LW LLWMAX 23 STOP END The following output should appear on screen MIQL A branch and cut mixed integer quadratic solver INFO MIQL Improved bounds calculation deactivated Parameters Maximal branch and bound nodes 1000 Node selection strategy 3 Branching variable selection 1 Improved bounds calculation 0 Lagrangian direction 1 QL decomposition mode 1 Cut generation mode 3 Primal heuristic mode 2 Maximal successive warmstarts 1000 Maximal number of cutting planes 9 Disjunctive cutting rounds 1 CMIR cutting rounds 1 Problem transformed to positive orthant BFOUR A branch and bound framework Parameters Number of integer variables 5 Maximal number of nodes
10. suc cessful cutting planes of the MILP solver SCIP of Achterberg 1 For the code MIQL we try to construct as many violated cuts as possible Thus we allow a large number of different starting constraints for the aggregation process Together with some further modifications we are able to construct significantly more cuts compared to the original aggregation scheme number of rounds is controlled by the user large number of cutting rounds is only beneficial for MIQP problems which are difficult to solve 3 3 Disjunctive Cuts Disjunctive cuts are cutting planes that are derived from a disjunction of the feasible region of a linearly constrained mixed integer program A disjunction in general decomposes the feasible region of the original optimization problem into disjoint polyhedra One natural disjunction in the context of mixed integer programming is a single disjunction for an integer variable where k I is a fractional component of the relaxed solution Z 7 yk lt Yel V ye I 6 disjunction gives rise to a disjunctive relaxation of the original problem feasible region of this disjunctive relaxation is the union of the disjoint polyhedra corresponding to the considered disjunction If the feasible region of the continuous relaxation of the original mixed integer problem is given by E ann i the disjunctive relaxation corresponding to the standard disjunction 6 is the union of the polyhedra P
11. until all nodes are either explored or fathomed 3 Cut Generation 3 1 Motivation of Mixed Integer Cutting Planes Cutting plane methods for mixed integer linear programming have been studied since a famous paper of Gomory 20 in 1958 The purpose of a cutting plane is to cut off the solution of a continuous relaxation of the original problem while retaining all feasible integer solutions The combination of a branch and bound enumeration procedure and a cutting plane method is called branch and cut Most cutting plane techniques relay either on the fact that the current iteration point to be cut off is a feasible basis solution i e a corner of the feasible polyhedral set or on special structures of some constraints or the whole feasible region The tremendous possible speed up for MILP problems is estimated by Bixby 12 see Table 1 In general the solution of a relaxed quadratic program does not need to be a basis solution at all nor is special information about the constraints available in most cases This speedup factor cutting planes presolve branching rules heuristics node presolve probing on dives Table 1 Speed up by Cutting in Linear Programming is the main reason why there exists only very few results on cutting planes for solving mixed integer quadratic programming MIQP problems of the form 1 Bienstock 11 suggests a branch and cut method based on disjunctive cuts for special MIQP problems arising in port
12. 1000 Integer tolerance 0 100000D 09 Branching variable selection 1 Node selection strategy 3 Convex problem Initial integer feasible solution 0 1750D 05 Output in the following order 5 status of current node iteration number ND index of current node objective function value P index of parent node LB lower bound UB upper bound S IT ND F FA LB UB 1 1 175108D 05 0 175108D 05 175035D 05 2 2 175107D 05 1 175108D 05 175035D 05 3 175082D 05 2 175108D 05 175035D 05 4 4 175081D 05 3 175108D 05 175035D 05 B 5 175081D 05 4 175108D 05 175081D 05 6 5 175106D 05 1 175108D 05 175081D 05 M 7 6 175081D 05 175108D 05 175081D 05 M 8 6 175007D 05 2 1751070 05 175081 05 9 6 175005D 05 5 175106D 05 175081D 05 24 10 6 175079D 05 3 175082D 05 175081D 05 M 11 6 175079D 05 4 175081D 05 175081D 05 Number of explored subproblems 11 Highest used index 5 F Y 17508 0900000 YC 1 98 000000 2 101 00000 3 39 000000 4 95 000000 5 0 0000000 QPs exposed to numerical problems 0 Generated cutting planes 2 Reduced gap 0 797323D 00 Optimal solution of original problem Objective value 6983 09000000 XC 1 2 0000000 XC 2 1 0000000 XC 3 61 000000 XC 4 5 0000000 XC 5 100 00000 Some further quadratic programs for testing a correct implementation are found in Hock a
13. 3 Generating disjunctive cuts for mixed integer programs Doc toral Dissertation Carnegie Mellon University Pittsburgh USA Powell M J D 1983 ZQPCVX A Fortran subroutine for convex quadratic program ming Report DAMTP 1983 NA17 University of Cambridge England Schittkowski K 1985 86 NLPQL A Fortran subroutine solving constrained nonlin ear programming problems Annals of Operations Research Vol 5 485 500 Schittkowski K 2003 QL A Fortran code for convex quadratic programming User s guide Report Department of Mathematics University of Bayreuth Germany Schittkowski K 2005 Optimal parameter selection in support vector machines Jour nal of Industrial and Management Optimization Vol 1 465 476 Schittkowski K 2006 NLPQLP A Fortran implementation of a sequential quadratic programming algorithm with distributed and non monotone line search user s guide version 2 2 Report Department of Computer Science University of Bayreuth Ger many Westerlund T Petterson F 1995 An extended cutting plane mathod for solving convex MINLP problems Computers and Chemical Engineering Vol 21 131 136 Wolters K 2006 Implementation of cutting plane separators for mixed integer pro grams Diploma Thesis Zuse Institute Berlin Germany 28
14. B XL XU X U EPS MODE IOUT IPRINT RW LRW IW LIW NP ADIAG j CDIAG START DELC NRDC STEPS Definition of the parameters ME MMAX N NMAX Number of constraints Number of equality constraints Row dimension of array A containing linear constraints MMAX must be at least one and greater or equal to M Number of optimization variables Row dimension of C NMAX must be at least one and greater or equal to N 20 MNN C NMAX N EPS MODE IOUT IFAIL IPRINT RW LRW Must be equal to M N WN when calling QL dimension of U Objective function matrix which should be symmetric and positive definite If MODE 0 C is supposed to be the upper triangular factor of a Cholesky decomposition of the objective function matrix Contains the constant vector of the quadratic objective function Matrix of the linear constraints first ME rows for equality then M ME rows for inequality constraints Constant values of linear constraints in the same order On input the one dimensional arrays XL and XU must contain the upper and lower bounds of the variables On return X contains the optimal solution On return U contains the multipliers subject to the linear con straints first M locations and bounds At the optimal solution all multipliers of inequality constraints are nonnegative Final termination accuracy e g 1 0D 12 The parameter value should not b
15. MIQL Fortran Subroutine for Convex Mixed Integer Quadratic Programming User s Guide T Lehmann K Schittkowski T Spickenreuther Address Department of Computer Science University of Bayreuth D 95440 Bayreuth Phone 921 553278 E mail klaus schittkowski uni bayreuth de Web http www klaus schittkowski de Date May 2010 Abstract The Fortran subroutine MIQL solves strictly convex mixed integer quadratic pro gramming problems subject to linear equality and inequality constraints by a branch and cut method At the root node of the branch and bound search tree disjunctive and complemented mixed integer rounding cuts are generated The continuous sub problem solutions are obtained by the primal dual method of Goldfarb and Idnani The code is designed for solving small to medium size mixed integer programs Its usage is outlined and an illustrative example is presented Keywords mixed integer quadratic programming quadratic optimization MIQP branch and bound branch and cut disjunctive cutting planes complemented mixed integer round ing cutting planes numerical algorithms 1 Introduction The Fortran subroutine MIQL solves strictly convex mixed integer quadratic programming problems subject to linear equality and inequality constraints of the form nins 7 ee 5 TI lude 1 m 5 J 3 1 5 1 with an n by n positive definite
16. dnani 1983 A numerically stable method for solving strictly convex quadratic programs Mathematical Programming Vol 27 1 33 Gomory R E 1958 An algorithm for integer solutions to linear programms Prince ton Mathematics Research Project Technical Report No 1 Princeton Univer sity USA Goncalves J P M Ladanyi L 2005 An implementation of a separation procedure for mixed integer rounding inequalities IBM Research Report RC23686 W0508 022 Gupta O K Ravindran A 1985 Branch and bound experiments in convex nonlinear integer programming Management Science Vol 31 1533 1546 Hock W Schittkowski K 1981 Test Examples for Nonlinear Programming Codes Lecture Notes in Economics and Mathematical Systems Vol 187 Springer Lehmann T Schittkowski K Spickenreuther T 2008 BFOUR A Fortran subrou tine for integer optimization by branch and bound user s guide Report Department of Computer Science University of Bayreuth Germany Leyffer S Fletcher R 1998 Numerical experiments with lower bounds for MIQP branch and bound SIAM Journal on Optimization Vol 8 604 616 Marchand H Wolsey L A 2001 Aggregation and mixed integer rounding to solve mixed integer programs Operations Research Vol 49 363 371 Mittelmann H 2007 Mixed integer QC QP benchmark http plato asu edu ftp miqp html 27 28 29 30 3l 32 33 34 35 oo Zi Perregaard M 200
17. e as the actual one nephew The warm start procedure is similar to the previous one The only difference is that active constraints in the two child nodes are compared This option is particularly efficient in case of the best of two strategy It is possible that warm starts lead to numerical instabilities based on round off errors To avoid this situation it is possible to limit the number of successive warm starts In case of a numerical error the problem is automatically resolved Warm starts are indicated by parameter START of QLX that can obtain the following values START Warm start performed by QLX 0 Normal QL run warm start performed 10 Dual warm start where only bounds are changed 11 Dual warm start where additional constraints are included 30 Warm start and early termination only one QL iteration to be performed 31 Warm start and early termination more than one QL iteration to be performed as provided by the parameter STEPS 40 Remove active constraints from the active set where the number of these constraints is stored in the calling parameter NRDC and the indices of the constraints to be removed are specified in the array DELC 41 Remove active constraints from the active set analogue to START 40 but continue afterwards with START 30 42 Remove active constraints from the active set analogue to START 40 but continue afterwards with START 31 43 Remove active constraints from t
18. e smaller than the underlying machine precision 0 Initial Cholesky factorization of C provided by the calling pro gram stored in the upper triangular part of array C 1 Cholesky decomposition internally computed start with Output unit number i e all write statements WBITEGIOUT Termination reason Optimality conditions satisfied Termination after too many iterations 40 N M Termination accuracy insufficient Inconsistency division by zero own Length of a working array too short 100 Inconsistent constraints and IFAIL 100 ICON where ICON denotes a constraint causing the conflict Specification of the desired output level 0 No output of the program 1 Only a final error message is given Double precision working array of length LRW 21 LRW IW LIW LIW NP ADIAG CDIAG START NACT DELC NRDC NRDC STEPS 7 Example Input parameter for length of RW must be at least 3 NMAX NMAX 2 10 NMAX 2 MMAX 1 Integer working array of length LIW Input parameter for length of IW must be at least N Input parameter defining the dimension of the reduced problem 11 Input array defining the diagonal entries of the constraint matrix of the reduced problem formulation 11 Input array defining the diagonal entries of the objective matrix of the reduced problem formulation 11 Input parameter specifying the type of warmstart to be performed more deta
19. entation details and program documentation are found in Section 6 followed by an illustrative example in Section 7 2 The Branch and Bound Procedure A branch and bound algorithm starts at the solution of the relaxed problem In our case this is the corresponding continuous problem where the integrality condition Vk J is removed See also Dakin 17 for an early paper on mixed integer linear programming An integer variable I is selected and two different continuous subproblems are generated They are obtained by rounding the continuous value of k I to get two separate subproblems one with upper bound another one with lower bound y 1 Each subproblem determines a node in a binary search tree which is internally created step by step A particular advantage of branch and bound methods is that certain subtrees can be eliminated or fathomed at an early stage before trying all possible combinations by exploring the whole tree There are three situations where fathoming is possible e A subproblem is infeasible All further subproblems obtained by branching from this node would also be infeasible i e the node can be eliminated e The subproblem has a feasible integer solution The corresponding optimal objective function value if less than the known upper bound provides a better upper bound for the optimal solution Further branching from this node is not required e The continuous solution of a s
20. ficient implementation or Boland 13 for an extension to the semi definite case Initially a Cholesky decomposition of is computed by an upper triangular matrix R such that C RTR In case of numerical instabilities e g round off errors or a semi definite matrix a certain multiple of the unit matrix is added to C to get a positive definite matrix for which a Cholesky decomposition is computed Subsequently the known triangular factor is saved and used again whenever solving another relaxed subproblem Successively violated constraints are added to an active set until a solution is obtained In each step the minimizer of the objective function subject to the new active set is computed If an iterate satisfies all linear constraints and bounds the optimal solution is obtained and the algorithm terminates A particular advantage of a dual method is that phase I of a primal algorithm i e the computation of an initial feasible point satisfying all linear constraints and bounds can be avoided The algorithm proves its robustness and efficiency in the framework of sequential quadratic programming algorithm see Schittkowski 30 33 A further advantage is the possibility of handling warm starts efficiently as pointed out in the previous section The calling sequence and the meaning of the parameter settings are described below where default values as far as applicable are set in brackets Usage CALL MIQL M ME MMAX NMAX MNN
21. folio optimization where the structure of one single constraint is exploited to construct cuts In general a cutting plane method tries to tighten the feasible region of a mixed integer program Although they are applied to solve convex nonlinear programs see e g Wester lund 34 their main application area are problems whose feasible regions are described by a polyhedron Cutting plane methods successively add linear inequalities to the constraint set They are called cutting planes and have to be valid for all feasible integral points 1 for every integral point satisfying the constraints of 1 On the other hand they have to cut off the relaxed feasible solution 2 7 of a continuous relaxation Cutting plane restrictions are added to the originally given set of constraints In the sequel we generate cuts only at the root node of the branch and bound search tree which implies that they are valid for each subproblem within the enumeration process One reason is that the construction of these cutting planes is quite expansive and pays out only if they are globally valid Furthermore the generation of cutting planes within the search process requires an efficient cut management and selection procedure which is a challenging task In our implementation cuts are constructed in rounds We solve the relaxed original problem then construct the cut planes and solve the tightened problem again This process is continued until no more cuts are fo
22. he active set analogue to START 40 but continue afterwards with START 10 Table 2 Warm Start Modes 10 The term dual warm start indicates situation where known solution is dual feasible for the subsequent continuous QP START 10 is applicable to a child problem during a branching step while START 11 is suitable for integrating cutting planes By now we do not handle the situation where a previous solution is primal feasible but not optimal This situation can be bypassed by first removing all active constraints preventing dual feasibility Furthermore QLX can handle the following special QP formulation in an efficient way min s x e ee i 2 iz b 0 z eR a x Tz b gt 0 4 11 C is an n x n matrix possessing the following structure C ET 12 where C is a n x n matrix and CD is a n np x n np diagonal matrix lt n can be considered as reduced problem dimension is the vector AD j where e is the j th unit vector i e the constraint matrix A with rows 1 01 exhibits the subsequent structure A A AD 13 where A is a m x n matrix and AD is a m x n n diagonal matrix The diagonal elements of the matrices CD and AD are stored in the input arrays CDIAG and ADIAG and n has to be given to QLX by input parameter NP This problem formulation 11 arises e g when a quadratic program is relaxed b
23. ils are given in Section 4 Input and output parameter indicating the number of active con straints For a correct warm start it has to possess the value corre sponding to the solution of the previous quadratic program Input array containing the indices of the active constraints to be removed by a corresponding warm start Input parameter defining the number of active constraints to be deleted by a corresponding warm start Therefore it has to be smaller than or equal to NACT Input parameter defining the number of QL iterations to be per formed before an early termination To give an example we consider the simple quadratic program 21 22 E R 23 24 25 min 1352 22 21 9811 1 26 61 3923 5 324 101 325 7 56 1 0 525 39 1 gt 0 100 lt z lt 100 i 1 5 The corresponding Fortran source code is listed below IMPLICIT NONE INTEGER NMAX MMAX MNNMAX LIWMAX LWMAX LLWMAX MAXNDS PARAMETER NMAX 5 MMAX 10 MAXNDS 1000 MNNMAX MMAX NMAX NMAX 4 USUS LIWMAX LWMAX 5 NMAX NMAX 2 32 NMAX 16 MMAX 16 NMAX 6 MAXNDS 4 MMAX 11 3 MAXNDS 14 22 LLWMAX 4 NMAX 3 INTEGER I J IPRINT IOUT N M ME NI MNN IOPT 20 IFAIL IM NMAX IW LIWMAX DOUBLE PRECISION FMIN ACC EPS A MMAX NMAX C NMAX NMAX B MMAX D NMAX XL NMAX XU NMAX X NMAX UCMNNMAX W LWMAX ROPT 20 LOGICAL LW LLWMAX LOPT 20 Set
24. ive cuts in MIQL creates and solves one cut generating quadratic program for each integer variable that is fractional at the relaxed solution of the mixed integer quadratic program 1 All generated cuts are strengthened if possible by taking into account the integrality condition of all other integer variables At present we only generate one single cut per disjunctive relaxation We restrict each CGQP to the subspace of those variables that are not fixed at their lower bound with value zero Afterwards we lift the generated disjunctive cut to the full space which significantly reduces the effort for generating disjunctive cuts 4 Algorithmic Features of MIQL As outlined in the previous sections the two main components of a branch and cut approach are the branch and bound enumeration method and the cut generation procedures In ad dition we have to take into account the algorithm by which the relaxed subproblems i e continuous convex quadratic programs are solved The primal dual method of Goldfarb and Idnani 19 is perfectly suited for this task For a better integration into the branch and cut framework we have extended the Fortran code QL of Schittkowski 31 This extended version is called QLX and its calling sequence is described in the next section When called from MIQL QLX performs warm starts depending on the relationship between successive subproblems e g in a branch and bound procedure 1 The new node is a child of the
25. k A gt b gt 0 y lower y y Uk Lyx Doe dem fa z gt 2 20 n gt In In general any optimization program with a constraint set consisting of the union of two or more polyhedra is called a disjunctive program z Essential for the theory behind disjunctive programs is the existence of a compact repre sentation of the convex hull of the union of polyhedra in a higher dimensional space It can be projected back to the original space This projection can be used to generate cuts which are called disjunctive cutting planes see Balas 2 3 4 5 Balas Ceria and Cornuejols 6 Balas and Perregaard 8 9 10 Ceria and Pataki 15 Ceria and Soares 16 for details In practice this is achieved by solving the so called cut generating linear program CGLP corresponding to disjunction 6 a ATu ugeg gt 0 a B R xR u uo R x R a voer gt 0 8 vw e R xR 6 bu uolgr 0 B 00 vof 0 uU ug 9 Vo gt 0 together with a normalization constraint of the form X utut Sut 1 9 The linear program 8 constructs the cutting plane a 5 gt p that is violated most at the current relaxed solution z y The constraints of problem 8 restrict the choice to those cutting planes that are valid for the union of the polyhedra P and P Therefore the optimization variables a 8 correspond to the resulting cutting plane
26. ling the general purpose branch and bound code BFOUR see Lehmann Schittkowski and Spickenreuther 24 Different rules for selecting variables and branching strategies are available Prior to the call of BFOUR the cut genera tion routines are executed complemented mixed integer rounding cuts are constructed by subroutine CMIR based is the corresponding cut generation routine of Achterberg 1 They are adapted for example to construct more cuts than the original implementation routine CMIR is called in combination with the routine DJCUT which generates cuts for each fractional variable by calling QLX to solve the corresponding cut generation quadratic program 10 T he necessary storage for generating the quadratic program is taken from the one that is reserved for the branch and bound process since the cuts are generated prior to the enumeration process The number of 2n 2m 6 maximum branch and bound nodes guarantees the generation of disjunctive cuts If sufficient storage is not available the generation of disjunctive cuts is omitted CMIR and disjunctive cutting plane generation are both executed within a given maximum number of rounds A suitable value depends heavily on the complexity of the MIQP problem The relaxed continuous quadratic program 1 is solved by an extention of the code QL of Schittkowski 31 which is based on the primal dual method of Goldfarb and Idnani 19 see 15 also Powell 29 for an ef
27. me data problem constraints number of number of non zero density non zero density variables integer of constraint of objective variables matrix function matrix Table 8 Characteristics of Mittelmann Test Examples Numerical results 1 number of nodes and calculation time are shown in Table 9 Cuts are not generated The option depth first is superior to the other node selection strategies since it benefits from warm starts Furthermore an integer feasible solution is often found quickly enabling fathoming of branch and bound branches 14 time sed ibell3a mazimal fraction best of all 61 588 ibella maximal fraction best of two 63 777 ibell3a maximal fraction depth first 74 720 imas284 maximal fraction best of all 142 241 imas284 maximal fraction best of two 186 412 imas284 mazimal fraction depth first 142 434 Table 9 Results for Mittelmann Test Examples without Cuts Table 10 includes the number of cuts and the reduced gap Calculation time is signifi cantly reduced by about 6096 for problem ibell3a and 2396 for imas284 problem branching strategy reduced gap ibell3a maximal fraction depth first 29 808 E z 51 5 imas284 maximal fraction depth first 127 849 1 494 11 6 Table 10 Results for Mittelmann Test Examples with Cuts 6 Program Documentation The Fortran subroutine MIQL solves the mixed integer quadratic program 1 with a posi tive definite objective function matrix by cal
28. nd Schittkowski 23 25 References 1 Achterberg T 2004 SCIP a framework to integrate constraint and mixed integer programming Technical Report 04 19 Zuse Institute Berlin Germany Balas E 1971 Intersection cuts a new type of cutting planes for integer program ming Operations Research Vol 19 19 39 Balas E 1975 Disjunctive programming cutting planes from logical conditions Nonlinear Programming Vol 2 330 382 Balas E 1979 Disjunctive programming Annals of Discrete Mathematics Vol 5 3 51 Balas E 1998 Disjunctive programming Properties of the convex hull of feasible points Discrete Applied Mathematics Vol 89 1 44 Balas E Ceria 5 Cornuejols 1993 A lift and project cutting plane algorithm for mixed 0 1 programs Mathematical Programming Vol 58 295 324 Balas E Jeroslow R 1980 Strengthening cuts for mixed integer programs European Journal of Operations Research Vol 4 224 234 Balas E Perregaard M 2001 Generating cuts from multiple term disjunctions Lecture Notes in Computer Science Vol 2081 348 360 Balas E Perregaard M 2002 Lift and project for mixed 0 1 programming recent progress Discrete Applied Mathematics Vol 123 129 154 Balas E Perregaard M 2003 A precise correspondence between lift and project cuts simple disjunctive cuts and mixed integer gomory cuts Mathematical Programming Vol 94 221 2245 Bienstock D 1994
29. omposition once and reuse it 1 Calculate new Cholesky decomposition whenever warmstart is not activate Control the cutting process 0 0 No cuts 1 Disjunctive cuts only 2 Complemented mixed integer rounding CMIR cuts only 3 Both disjunctive and CMIR cuts Maximal number of cycles for disjunctive cuts 1 Maximal number of cycles for CMIR cuts 1 Primal heuristic mode 0 0 No primal heuristics 1 Nearest integer 2 Feasibility pump Input parameter for the reduced quadratic program dimension N Output parameter for total number of generated cuts Output parameter for total number of branch and bound nodes Double precision option array Output parameter for time spent for generating cuts Output parameter for time spent for branch and bound 18 ROPT 3 ROPT 4 LOPT 20 LOPT 1 IOUT IFAIL IPRINT Output parameter for gap reduced by cutting planes compared to original relaxed solution Output parameter for relaxed optimal objective value Logical option array Output parameter indicating cut generation Input parameter for desired output unit number i e all write statements start with WRITE IOUT Output parameter for termination reason 4 Branch and bound root QP not solvable 3 Relaxed QP without feasible solution 2 Feasible solution not within maximal number of nodes 1 Feasible integer solution does not exist 0 Optimal solution found 1
30. pliers depend on the solution strategy ie the node selection strategy cut generation and the branching rule Independent multipliers are especially useful if MIQL is used as subproblem solver for the mixed integer nonlinear solver MISQP of Exler and Schittkowski 18 since then Lagrangian multipliers are needed to calculate a BFGS update 7 Primal heuristics are now used for finding a first integer feasible solution before the start of the enumeration process diving and shifting procedure is implemented which often finds a integer feasible solution close to the optimal one by successively fixing one integer variable Afterwards the values of some variables are shifted according to the corresponding value of the partial derivative of the Lagrangian function Furthermore we have implemented a cheap routine to fix all integer variables by rounding them up or down and trying to retain feasibility 5 Numerical Results 5 1 MISQP Subproblems numerical results are obtained on a Intel E6600 Dual Core CPU with 2 4 GHz under Windows XP 64 using Intel Fortran Compiler version 9 1 For testing the nonlinear mixed integer code MISQP a large number of test problems have been implemented by Exler and Schittkowski 18 Since in each iteration of the SQP based method a mixed integer quadratic programming problem must be solved by MIQL a large number of test examples is available A subset of 35 MIQP subproblems is selected with calculation
31. times over 0 1 seconds Table 3 shows some characteristic data of the corresponding nonlinear programs Most of those test examples are selected from the GAMS MINLP Library see Bussieck Drud and Meeraus 14 Some are simple case studies from petroleum industry Prob 12 problem continuous integer inequality equality Peale bes constraints oma MISQP57 2 4 0 1 MISQP60 0 MISQP62 MISQP63 MISQP66 MISQP83 MISQP115 MISQP116 MISQP117 Table 3 Problem Characteristics Table 4 Node Selection Strategy lems MISQP115 MISQP116 and MISQP117 have multiple global optima and the numerical results might become incompatible In the subsequent tables the total number of branch and bound nodes and the total calculation times are reported For comparing all other options we always choose the fastest alternative Table 4 shows the influence of different node selection strategies Although best of all leads to the smallest amount of continuous QP problems that have to be solved during the branch and bound process calculation time is considerably higher than for the other two strategies because of very few warm starts Strategies best of two and depth first guided by the Lagrangian function value exhibit nearly identical performance Next we consider the influence of warm starts subject to the node selection strategies best of two and depth first see Table 5 For best of all warm starts can be neglected We observe that
32. ubproblem has a minimal solution value greater than the actual upper bound Further branching from the node would only increase function values and there is no chance to find a better integer solution The node is eliminated A branch and bound method represents a general framework which is easily adapted to other classes of mathematical optimization problems say linear or nonlinear programming see for example Gupta and Ravindran 22 Essential is convexity of the original problem or the possibility to compute global solutions at all nodes The integral node with minimal objective value represents the solution of 1 or we get the information that no feasible mixed integer solution exists Note that although the relaxed version of 1 is assumed to be a strictly convex optimization problem the solution of the mixed integer problem may not be unique There is also a possibility to exploit lower bounds within a branch and bound algorithm by investigation the dual of 1 see Leyffer and Fletcher 25 for details Thus there remain two open questions leading to different options by which a branch and bound algorithm can be tuned 1 Selecting an integer variable with a non integer value in the continuous solution of a given node for branching a maximal fractional branching We select an integer variable from the solution of the relaxed subproblem with largest distance from next integer value b minimal fractional branching We select
33. und or the maximal number of rounds is reached At the moment two different cut generation procedures are implemented which will be described in detail in the remainder of this section First we will describe complemented mixed integer rounding cutting planes which go back to Marchand and Wolsey 26 Sub sequently we introduce disjunctive programming and the resulting disjunctive cuts both heavily influenced by Balas 3 3 2 Complemented Mixed Integer Rounding Cuts Complemented mixed integer rounding CMIR cuts have been developed by Marchand and Wolsey 26 see also Goncalves and Ladanyi 21 and Wolters 35 To simplify the notation the m constraint indices are denoted by P 1 m To generate valid cuts one has to transform 1 such that all variables are positive The generation of CMIR cutting planes consists of three steps aggregation bound sub stitution and cut generation Within the aggregation step a selection of original constraints is aggregated and if possible some continuous variables are eliminated During bound sub stitution continuous variables are substituted by 2 u 2 or j 2 The resulting constraint is relaxed to yield a mixed knapsack set X defined by _ ts ER x N Yoy lt 8 8 y Xy Vic r 3 icI where f y ges j J 0 original feasible set X is defined by X 4 x y E Re x N X aly X VLE P Yi yu
34. vectors for upper and lower bounds of the variables Double precision output vector for the optimal solution Double precision output vector for the multipliers subject to linear constraints at the first M positions and for bounds but without cut constraints Double precision output parameter containing optimal objective function value Double precision input parameter for identifying integer values set to machine precision in case of ACC 0 Double precision input parameter for termination accuracy of the QP solver QL should not be smaller than machine precision set to machine precision in case of EPS 0 Integer option array Branching rule 1 Maximal fractional branching 2 Minimal fractional branching Node selection strategy 1 Best of all 17 IOPT 3 IOPT 4 IOPT 5 IOPT 6 IOPT T IOPT 8 IOPT 9 IOPT 10 IOPT 11 IOPT 12 IOPT 13 IOPT 14 ROPT 20 1 ROPT 2 2 Best of two 3 Depth first Maximal number of nodes should be at least 2 N 2 M 6 2 to guarantee generation of disjunctive cuts and must be at least 2 4 2 for recomputing Lagrange multipliers 100 000 Maximal number of successive warmstarts to avoid numerical in stabilities 100 Calculate improved bounds if best of all selection strategy is used 0 Select direction for depth first according to value of Lagrange func tion 0 Cholesky decomposition mode 1 0 Calculate Cholesky dec
35. warm starts improve the performance significantly The branching direction is chosen by the value of the Lagrangian function in case of depth best of two depth first Table 5 Performance Improvements due to Warm Starts 13 time time 12255 302 115 Table 6 Branching Direction by Lagrangian Function improved bounds no improved bounds nodes time sec time sec 115 588 102 Table 7 Calculating Improved Bounds Information first to prevent blind diving Significant reduction of branch and bound nodes is achieved see Table 6 Table 7 shows the effect of calculating improved bounds leading to faster fathoming The improved bounds exploit dual information and are calculated by an additional QL iteration Improved bounds see Leyffer and Fletcher 25 are only calculated in case of the best of all strategy The results of Table 7 show that although the number of branch and bound nodes is reduced by early fathoming the calculation of the improved bound information slows down the overall performance 5 2 Benchmark Test Cases of Mittelmann Next we investigate the performance of MIQL for some of the benchmark test examples of Mittelmann 27 Since MIQL is a dense solver and cannot exploit sparsity patterns some of the problems are not suitable In some other situations the continuous relaxation is not solvable without special adjustments Thus we only consider examples ibell3a and imas284 see Table 8 for so
36. y the introduction of m additional slack variables to ensure feasibility or in special data mining problems related to support vector machines see Schittkowski 32 In addition to the extended version of QL there are some further improvements of our new version of MIQL 1 The new depth first search benefits from warm starts 2 storing procedure for nodes in the branch and bound process is improved Since a node marked by a fathoming rule is not needed any longer its memory is used again for storing a new node 11 3 Error handling is improved In case of false termination when solving a quadratic program the whole process is not terminated as in the previous version of MIQL Instead the error is treated in the same way as an infeasible leaf 4 Improved bounds according to Leyffer and Fletcher 25 are calculated if nodes are selected by best of all using warm start options START 30 and START 41 5 The branching direction for depth first search is determined by evaluating the La grangian function at j and where j is the selected branching variable and 2 7 the solution of the previous branch and bound QP 6 Lagrangian multipliers can be calculated independently of the solution method by a modified least squares approach which uses an extended QP formulation to guarantee feasibility This calculation is executed if storage for 2n m 4 branch and bound nodes is provided Otherwise the returned multi
Download Pdf Manuals
Related Search