MIQL: A Fortran Subroutine for Convex Mixed
Contents
1. 1 see Leyffer and Fletcher 9 for details Thus there remain decisions by which the numerical performance of a branch and bound algorithm may become influenced 1 Select an integer variable with a non integer value in the actual continuous solution of a node for branching a maximal fractional branching Select an integer variable from solution values of a relaxed subproblem with largest distance from neighbored integer values b minimal fractional branching Select an integer variable from solution values of a relaxed subproblem with smallest distance from neighbored integer values 2 Search for a free node in the whole tree from where branching i e the generation of two new subproblems is initiated a best of all All free nodes are compared and a node with lowest objective function value is selected b 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 c depth first Select 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 very fast since the only difference between the both continuous programs is one2. R such that C RTR In case of numerical instabilities e g round off errors or a semi definite matrix C 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 12 15 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 N NMAX MNN NI IND C D A B XL XU X F ACC EPS IOPT ROPT LOPT IOUT IFAIL IPRINT RW LRW IW LIW LW LLW TTS TR ER ME TE Parameter Definitio
3. additional box constraint Usually the search tree is smaller i e there remain less unexplored nodes The disadvantage is that the child node is arbitrarily selected leading to a possibly larger number of subproblems to be solved After selecting an integer variable for branching and an unexplored node to continue the optimal solutions of the corresponding relaxed quadratic programs are computed This process is iterated until all nodes are either explored or fathomed The upper bound is infinity until the first integer feasible solution is found 3 Algorithmic Features of MIQL An essential component of a branch and bound is the ability to solve continuous quadratic programming problems very efficiently The primal dual method of Goldfarb and Idnani 5 is frequently used to solve convex quadratic programs To exploit as much information as possible from an existing solution within a branch and bound framework we extend the Fortran code QL of Schittkowski 13 This 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 A new node is a child of the 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 and Idnani successively adds cons
4. 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 execution no warm starts 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 the active set analogue to START 40 but continue afterwards with START 10 Table 1 Warm Start Modes The term dual warm start indicates a situation where a known solution is dual feasible
5. objective function of 1 is denoted by fean serwe e 2 2 y If n gt 0 the mixed integer quadratic program 1 is now solved by a branch and bound method This manual is organized as follows The subsequent section introduces the basic concept of a branch and bound method and describes the different options of the underlying Fortran subroutine Section 3 describes more details of the implementation especially the warm start features Section 4 summarizes some numerical results to compare different options Further implementation details and program documentation are found in Section 5 followed by an illustrative example in Section 6 2 The Branch and Bound Procedure A branch and bound algorithm starts at the solution of the relaxed problem in our case the optimal solution of a continuous quadratic program where the integrality condition y N is replaced by y R see also Dakin 3 for an early paper on mixed integer linear pro gramming An integer variable y I is selected and two different continuous subproblems are generated They are obtained by rounding the continuous value of yp k I to get two separate subproblems one with upper bound y 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
6. 2 1 IOPT 3 MAXNDS Cc C Set problem data G DO I 1 N DO J 1 N C I J 0 0D0 14 Q END DO c I I A 1 1 XL I XU I END DO DO I 1 NI IND I END DO A 1 1 A 1 5 B 1 D 1 D 2 D 3 D 4 D 5 Call MIQL CALL MIQL AS AN SK STOP END i H 7 56D0 0 5D0 39 1D0 21 98D0 1 26D0 61 39D0 5 3D0 101 3D0 M MNN A U ROPT RW LLW 1 0D0 0 0D0 100 0D0 100 0D0 ME MMAX NMAX NI IND C D zZ F ACC EPS IOPT LOPT IOUT IFAIL IPRINT LRW IW LIW LW The following output should appear on screen MIQL A branch and bound mixed integer quadratic solver kkk kkk kkk kkk kkk kkk INFO INFO INFO INFO INFO INFO Parameters Maximal Node selection strategy Branching variable selection Improved bounds calculation Lagrangian direction QL decomposition mode Maximal successive warmstarts MIQL Maximal number of successive warmstarts set to zero MIQL Improved bound computation deactivated MIQL Lagrangian direction deactivated MIQL Cholesky decomposition mode set to 1 MIQL Reduced QP dimension set to original dimension MIQL Initial number of branch and bound nodes set to zero branch and bound nodes 1000 orRFOOrFRF Start of the Mixed Integer Branch and Bound Code BFOUR Version 3 0 Mar 2013 15 Parameters Number of integer variables 5 Maximal number of nodes 1000 Integer tolerance 0 10
7. 4 Branch and bound root QP not solvable 3 Relaxed QP without feasible solution 2 Feasible solution could not be computed within maximal num ber of nodes 1 Feasible integer solution does not exist 0 Optimal solution found 1 Feasible solution found but the maximum number of nodes reached 2 Index in IND out of bounds 3 Termination due to internal inconsistency of branch and bound routine BFOUR 4 Length of a working array RW IW or LW too small 12 IPRINT RW LRW LRW IW LIW LIW 5 Parameters N M MNN ME or NI incorrect 6 Option of IOPT incorrectly set 7 Independent Lagrangian multipliers could not be calculated primal solution is valid Increase IOPT 3 for performing indepen dent Lagrangian multiplier calculation 8 IOUT or IPRINT incorrectly set 9 Lower variable bound greater than upper one 10 i Problem is continuous and could not be solved due to error i of QP solver QL see corresponding documentation 100 i Problem is continuous and infeasible constraint i could not be included in active set Input parameter defining the output level 0 No output to be generated 1 Root QP and final performance summary 2 Additional branch and bound iterations counter 3 Additional output from all generated subproblems fractional feasible node B new best integer feasible solution I infeasible node M marked node 3 Complete output for all nodes e g with
8. D 09 Branching strategy 1 Node selection 1 Convex problem No initial integer feasible solution known Output in the following order S status of current node B branched node I infeasible node M marked node other node IT iteration count ND index of current node F objective function value P index of parent node LB lower bound UB upper bound S IT ND F P LB UB 1 1 699651D 04 O 699651D 04 0 100000D 31 2 2 698324D 04 1 699651D 04 0 100000D 31 3 3 697573D 04 1 699651D 04 0 100000D 31 4 4 698317D 04 2 698324D 04 0 100000D 31 5 5 698306D 04 2 698324D 04 0 100000D 31 6 6 698312D 04 4 698317D 04 0 100000D 31 7 7 698292D 04 4 698317D 04 0 100000D 31 B 8 8 698285D 04 6 698312D 04 698285D 04 B 9 8 698309D 04 6 698312D 04 698309D 04 Optimal solution Number of explored subproblems 9 Highest index T F Y 6983 0900 Y 1 2 0000000 Y 2 1 0000000 Y 3 61 000000 Y 4 5 0000000 Y 5 100 00000 INFO MIQL QPs with numerical problems 0 Optimal solution of original problem Objective function value 6983 090000 X 1 2 0000000 X 2 1 0000000 X 3 61 000000 X 4 5 0000000 X 5 100 00000 Some further quadratic programs for testing a correct implementation are found in Hock and Schittkowski 7 16 References 1 Boland N L 1997 A dual active set algorithm for positive semi definite
9. MIQL A Fortran Subroutine for Convex Mixed Integer Quadratic Programming User s Guide T Lehmann K Schittkowski Address Prof K Schittkowski Siedlerstr 3 D 95488 Eckersdorf Germany Phone 49 921 32887 E mail klausQschittkowski de Web http www klaus schittkowski de Date October 2013 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 bound method The continuous subproblem 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 pre sented Keywords mixed integer quadratic programming quadratic optimization MIQP branch and bound 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 ming 9 7 a z BS Jepi 1 m j z 3 J e 9 x Ee R y e N P 1 a 7 a 20 j me 1 M TITL Tu with an n by n positive definite matrix C and n nj ne an n dimensional vector d an m by n matrix A a4 m and an m vector b Lower and upper bounds for the continuous and integer variables 1 u Y and y respectively are separately handled The
10. bounds 4 Full debug output 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 3 MMAX 8 N 3 MAXNDS 33 Integer working array of length LIW Input parameter for length of IW must be at least NMAX NI MMAX 6 MAXNDS 4 N 30 13 LW LLW Logical working array of length LLW LLW Input parameter for length of LW must be at least 3 NI 23 6 Example To give an example we consider the simple quadratic program min 45 gt a 21 987 1 262 61 3923 5 324 101 325 7 56a 0 525 39 1 gt 0 100 lt a lt 100 i 1 5 z1 T2 E R z3 4 5 EN The corresponding Fortran source code is listed below IMPLICIT NONE INTEGER NMAX MMAX MNNMAX LIW LRW LLW MAXNDS PARAMETER NMAX 5 MMAX 1 MAXNDS 1000 MNNMAX MMAX NMAX NMAX LIW 6 NMAX MMAX 6 MAXNDS 30 LRW 5 NMAX NMAX 2 25 NMAX 3 MMAX 3 MAXNDS 33 LLW 3 NMAX 23 INTEGER M ME N MNN NI IND NMAX IOUT IFAIL IPRINT IW LIW I J IOPT 20 DOUBLE PRECISION F ACC EPS A MMAX NMAX C NMAX NMAX B MMAX D NMAX XLCNMAX XUCNMAX X NMAX U MNNMAX RW LRW ROPT 20 LOGICAL LW LLW LOPT 20 Cc C Set parameters Cc N 5 M 1 ME 0 MNN MMAX N WN NI N ACC 1 0D 10 EPS 1 0D 12 IPRINT 2 IOUT 6 DO I 1 20 ROPT I 1 0D0 IOPT I 1 LOPT I TRUE END DO IOPT 1 1 IOPT
11. culation times over 0 1 seconds Table 2 shows some characteristic data of the corresponding nonlinear programs MITP57 MITP60 MITP62 MITP63 MITP66 MITP83 MITP115 MITP116 MIT PII Table 2 Problem Characteristics Most of those test examples are selected from the GAMS MINLP Library see Bussieck Drud and Meeraus 2 Some are simple case studies from petroleum industry Problems MITP115 MItp116 and MItp117 have multiple global optima 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 best of all depth first Table 3 Node Selection Strategy warm start no warm start warm start no warm start Table 4 Performance Improvements due to Warm Starts Table 3 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 4 For best of all warm starts can be neglected We observe that warm starts improve the perfo
12. for the subsequent continuous program START 10 is applicable to a child problem during a branching step Furthermore code QLX is able to handle the following special formulation of a quadratic program in an efficiently Suppose we want to solve a problem of the form min Lae gC d i x R a zEeER a where isann xn block diagonal matrix a o where C is a np X np matrix and Cp is a n np X n np diagonal matrix np lt n can be considered as reduced problem dimension is the j th column of a matrix Ap i e the constraint matrix A is decomposed in the form A A Ap 5 where A is a m x np matrix and Ap is an m x n np diagonal matrix This problem formulation 3 arises e g when a quadratic program is relaxed by the introduction of m additional slack variables to ensure feasibility or in special data mining problems related to support vector machines see Schittkowski 14 There are further features of our new version of our branch and bound code MIQL 1 The new depth first search benefits from warm starts 2 Since a node marked by a fathoming rule is not needed any longer the node data will become overwritten 3 In case of false termination when solving a quadratic program the whole process is not terminated and the error is treated in the same way as an infeasible leaf 4 Improved bounds according to Leyffer and Fletcher 9 are calculated if node
13. n M Input parameter for the total number of constraints including equal ity constraints 10 ME MMAX N NMAX MNN NI IND NI C NMAX N X N U MNN F ACC EPS IOPT 20 IOPT 1 IOPT 2 Input parameter for the number of equality constraints Input parameter defining the row dimension of array A containing linear constraints must be at least one and or equal to M 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 dimension 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 which should be symmetric and positive definite Double precision input vector for constant vector of objective func tion Double precision input matrix for linear constraints first ME rows for equality then M ME rows for inequality constraints Double precision input vector for constant values of linear con straints in same order as matrix A Double precision input 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 first lower then upper bounds Double precision out
14. nction matrix ibell3a 104 122 60 3 1 0 4 imas284 68 151 150 95 3 0 6 Table 7 Characteristics of Mittelmann Test Examples Numerical results i e number of nodes and calculation time are shown in Table 8 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 time 0 ibell3a maximal fraction best of all 61 588 ibell3a 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 maximal fraction depth first 142 434 Table 8 Results for Mittelmann Test Examples without Cuts 5 Program Documentation The Fortran subroutine MIQL solves the mixed integer quadratic program 1 with a pos itive definite objective function matrix by calling the general purpose branch and bound 9 code BFOUR see Lehmann and Schittkowski 8 Different rules for selecting variables and branching strategies are available The relaxed continuous quadratic program 1 is solved by a modification of the code QL of Schittkowski 13 which is based on the primal dual method of Goldfarb and Idnani 5 see also Powell 11 or Boland 1 for an extension to the semi definite case Initially a Cholesky decomposition of C is computed by an upper triangular matrix
15. nteger QC QP benchmark http plato asu edu ftp miqp html 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 NEPQL 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 17 15 Schittkowski K 2006 NZEPQLP 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 18
16. put parameter containing optimal objective function value at successful return see IFAIL Double precision input parameter for identifying integer values Double precision input parameter for termination accuracy of the QP solver QL should not be smaller than machine precision Integer option array Branching rule 1 Maximal fractional branching 2 Minimal fractional branching Node selection strategy 1 Best of all 11 IOPT 3 IOPT 4 IOPT 5 IOPT 6 IOPT 7 IOPT 12 IOPT 14 ROPT 20 LOPT 20 IOUT IFAIL 2 Best of two 3 Depth first Maximal number of nodes 100 000 Maximal number of successive warmstarts to avoid numerical in stabilities 100 Calculate improved bounds if best of all node 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 decomposition once and reuse it 1 Calculate new Cholesky decomposition whenever warmstart is not activate Input parameter for the reduced quadratic program dimension N Output parameter for total number of branch and bound nodes Double precision option array to remain compatible with previous versions Logical option array to remain compatible with previous versions Input parameter for desired output unit number i e all write statements start with WRITE IOUT Output parameter for termination reason
17. quadratic 10 11 12 13 14 programming Mathematical Programming Vol 78 1 27 Bussieck M R Drud A S Meeraus A 2007 MINLPLib A collection of test mod els for mixed integer nonlinear programming GAMS Development Corp Washington D C 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 Idnani A 1983 A numerically stable method for solving strictly convex quadratic programs Mathematical Programming Vol 27 1 33 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 2008 BFOUR A Fortran subroutine for integer op timization 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 Mittelmann H 2007 Mixed i
18. rmance significantly The branching direction is chosen by the value of the Lagrangian function in case of depth first to prevent blind diving Significant reduction of branch and bound nodes is achieved see Table 5 Table 6 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 9 are only calculated in case of the best of all strategy The results of Table 6 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 time sec time sec 122 055 02 15 Table 5 Branching Direction by Lagrangian Function improved bounds no improved bounds time see Table 6 Calculating Improved Bounds Information 4 2 Benchmark Test Cases of Mittelmann Next we investigate the performance of MIQL for some of the benchmark test examples of Mittelmann 10 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 7 for some data problem constraints number of number of non zero density non zero density variables integer of constraint of objective variables matrix fu
19. s are selected by best of all by using warm start options START 30 and START 41 of QLX 5 The branching direction for depth first search is determined by evaluating the La grangian function at y and where j is the selected branching variable and the solution of the previous branch and bound QP 6 Lagrangian multipliers are 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 multipliers would depend e g on the node selection strategy 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 4 since then Lagrangian multipliers are needed to calculate a BFGS update 4 Numerical Results 4 1 MISQP Subproblems All 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 4 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 cal
20. stage before trying all possible combinations by exploring the whole tree There are three e A subproblem is infeasible All subsequently generated subproblems obtained by branch ing from this node would also be infeasible e A subproblem has a feasible integer solution If the corresponding optimal objective function value is less than a known upper bound a better upper bound is provided A better solution cannot be generated starting from this node e Let the optimal value of a non integral solution of a subproblem be greater than the ac tual upper bound Then further branching from the node would only increase function values and there is no chance to find a better integer solution 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 6 An essential assumption is the convexity of the given nonlinear program or the possibility to compute global solutions at all nodes The integral node with minimal objective value represents the solution 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
21. traints 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 a complete new solution of the quadratic program 2 The new node has the same parent node as the actual one i e both problems differ only in one box constraint for an integer variable The algorithm identifies the solution of a quadratic program by an linearly 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 tries first to ensure dual feasibility This is achieved by saving the active 4 constraints of the common parent node They are compared to those of the child i e the previously 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 node 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
Download Pdf Manuals
Related Search