[37] | 1 | function output = bmibnb(p) |
---|
| 2 | %BMIBNB Branch-and-bound scheme for bilinear programs |
---|
| 3 | % |
---|
| 4 | % BMIBNB is never called by the user directly, but is called by YALMIP from |
---|
| 5 | % SOLVESDP, by choosing the solver tag 'bmibnb' in sdpsettings |
---|
| 6 | % |
---|
| 7 | % The behaviour of BMIBNB can be altered using the fields in the field |
---|
| 8 | % 'bmibnb' in SDPSETTINGS |
---|
| 9 | % |
---|
| 10 | % bmibnb.lowersolver - Solver for lower bound [solver tag ('')] |
---|
| 11 | % bmibnb.uppersolver - Solver for upper bound [solver tag ('')] |
---|
| 12 | % bmibnb.lpsolver - Solver for LP bound tightening [solver tag ('')] |
---|
| 13 | % bmibnb.branchmethod - Branch strategy ['maxvol' | 'best' ('best')] |
---|
| 14 | % bmibnb.branchrule - Branch position ['omega' | 'bisect' ('omega')] |
---|
| 15 | % bmibnb.lpreduce - Improve variable bounds using LP [ real [0,1] (0 means no reduction, 1 means all variables) |
---|
| 16 | % bmibnb.lowrank - Partition variables into two disjoint sets and branch on smallest [ 0|1 (0)] |
---|
| 17 | % bmibnb.target - Exit if upper bound<target [double (-inf)] |
---|
| 18 | % bmibnb.roottight - Improve variable bounds in root using full problem [ 0|1 (1)] |
---|
| 19 | % bmibnb.vartol - Cut tree when x_U-x_L < vartol on all branching variables |
---|
| 20 | % bmibnb.relgaptol - Tolerance on relative objective error (UPPER-LOWER)/(1+|UPPER|) [real (0.01)] |
---|
| 21 | % bmibnb.absgaptol - Tolerance on objective error (UPPER-LOWER) [real (0.01)] |
---|
| 22 | % bmibnb.pdtol - A number is declared non-negative if larger than...[ double (-1e-6)] |
---|
| 23 | % bmibnb.eqtol - A number is declared zero if abs(x) smaller than...[ double (1e-6)] |
---|
| 24 | % bmibnb.maxiter - Maximum number nodes [int (100)] |
---|
| 25 | % bmibnb.maxtime - Maximum CPU time (sec.) [int (3600)] |
---|
| 26 | |
---|
| 27 | % Author Johan Löfberg |
---|
| 28 | % $Id: bmibnb.m,v 1.13 2006/09/20 12:43:45 joloef Exp $ |
---|
| 29 | |
---|
| 30 | bnbsolvertime = clock; |
---|
| 31 | showprogress('Branch and bound started',p.options.showprogress); |
---|
| 32 | |
---|
| 33 | % ************************************************************************* |
---|
| 34 | % INITIALIZE DIAGNOSTICS ETC |
---|
| 35 | % ************************************************************************* |
---|
| 36 | switch max(min(p.options.verbose,3),0) |
---|
| 37 | case 0 |
---|
| 38 | p.options.bmibnb.verbose = 0; |
---|
| 39 | case 1 |
---|
| 40 | p.options.bmibnb.verbose = 1; |
---|
| 41 | p.options.verbose = 0; |
---|
| 42 | case 2 |
---|
| 43 | p.options.bmibnb.verbose = 2; |
---|
| 44 | p.options.verbose = 0; |
---|
| 45 | case 3 |
---|
| 46 | p.options.bmibnb.verbose = 2; |
---|
| 47 | p.options.verbose = 1; |
---|
| 48 | otherwise |
---|
| 49 | p.options.bmibnb.verbose = 0; |
---|
| 50 | p.options.verbose = 0; |
---|
| 51 | end |
---|
| 52 | |
---|
| 53 | % ************************************************************************* |
---|
| 54 | % Assume feasible (this property can be changed in the presolve codes |
---|
| 55 | % ************************************************************************* |
---|
| 56 | p.feasible = 1; |
---|
| 57 | |
---|
| 58 | % ************************************************************************* |
---|
| 59 | % Extract bounds from model |
---|
| 60 | % ************************************************************************* |
---|
| 61 | [p.lb,p.ub,used_rows] = findulb(p.F_struc,p.K); |
---|
| 62 | |
---|
| 63 | % These two are needed for ex14_1_7 |
---|
| 64 | p = fixbounds(p); |
---|
| 65 | p = presolve_bounds_from_equalities(p); |
---|
| 66 | p = preprocess_eval_bounds(p); |
---|
| 67 | p = fixbounds(p); |
---|
| 68 | |
---|
| 69 | % ************************************************************************* |
---|
| 70 | % The bilinear solver does not support higher order polynomial (i.e. >2). |
---|
| 71 | % Hence, we need to bilinearize the problem. however, the upper bound |
---|
| 72 | % solver might support polynomial problems, so we should keep a copy of the |
---|
| 73 | % original problem also |
---|
| 74 | % ************************************************************************* |
---|
| 75 | n_in = length(p.c); |
---|
| 76 | [p_bilin,changed] = bilinearize(p); |
---|
| 77 | if (p.solver.uppersolver.constraint.equalities.polynomial & p.solver.uppersolver.objective.polynomial) |
---|
| 78 | p_bilin.originalModel = p; |
---|
| 79 | p = p_bilin; |
---|
| 80 | else |
---|
| 81 | p = p_bilin; |
---|
| 82 | p.originalModel = p_bilin; |
---|
| 83 | end |
---|
| 84 | |
---|
| 85 | % ************************************************************************* |
---|
| 86 | % Pre-calc lists of linear/bilinear/nonlinear variables |
---|
| 87 | % ************************************************************************* |
---|
| 88 | [p.linears,p.bilinears,p.nonlinears] = compile_nonlinear_table(p); |
---|
| 89 | |
---|
| 90 | % ************************************************************************* |
---|
| 91 | % Select branch variables. We should branch on all variables that are |
---|
| 92 | % involved in bilinear terms. |
---|
| 93 | % ************************************************************************* |
---|
| 94 | p.branch_variables = decide_branch_variables(p); |
---|
| 95 | |
---|
| 96 | % ************************************************************************* |
---|
| 97 | % Extract bounds from model |
---|
| 98 | % ************************************************************************* |
---|
| 99 | p = preprocess_bilinear_bounds(p); |
---|
| 100 | p = preprocess_eval_bounds(p); |
---|
| 101 | |
---|
| 102 | % ************************************************************************* |
---|
| 103 | % Now reduce the branch variables by removing bilinear terms that only have |
---|
| 104 | % been introduced due to a convex quadratic objective |
---|
| 105 | % ************************************************************************* |
---|
| 106 | p = reduce_bilinear_branching_variables(p); |
---|
| 107 | |
---|
| 108 | % ************************************************************************* |
---|
| 109 | % Is the given initial solution feasible, or is the zero-solution feasible? |
---|
| 110 | % ************************************************************************* |
---|
| 111 | [p,x_min,upper] = initializesolution(p); |
---|
| 112 | |
---|
| 113 | % ************************************************************************* |
---|
| 114 | % Simple pre-solve loop. The loop is needed for variables defined as w = |
---|
| 115 | % x*y, x = t*u,y=.. |
---|
| 116 | % ************************************************************************* |
---|
| 117 | i = 0; |
---|
| 118 | goon = 1; |
---|
| 119 | while goon |
---|
| 120 | start = [p.lb;p.ub]; |
---|
| 121 | i = i+1; |
---|
| 122 | p = updatenonlinearbounds(p); |
---|
| 123 | p = presolve_bounds_from_equalities(p); |
---|
| 124 | p = updatenonlinearbounds(p); |
---|
| 125 | p = updatenonlinearbounds(p); |
---|
| 126 | p = preprocess_eval_bounds(p); |
---|
| 127 | i = i+1; |
---|
| 128 | goon = (norm(start-[p.lb;p.ub],'inf') > 1e-1) & i < 50; |
---|
| 129 | start = [p.lb;p.ub]; |
---|
| 130 | end |
---|
| 131 | |
---|
| 132 | % ************************************************************************* |
---|
| 133 | % Default root bound improvement |
---|
| 134 | % ************************************************************************* |
---|
| 135 | %p.lb = (p.lb+1e5)-1e5; |
---|
| 136 | %p.ub = (p.ub+1e5)-1e5; |
---|
| 137 | close = find(abs(p.lb - p.ub) < 1e-12); |
---|
| 138 | p.lb(close) = (p.lb(close)+p.ub(close))/2; |
---|
| 139 | p.ub(close) = p.lb(close); |
---|
| 140 | p = root_node_tighten(p,upper); |
---|
| 141 | p = updatenonlinearbounds(p); |
---|
| 142 | p = preprocess_eval_bounds(p); |
---|
| 143 | p = presolve_bounds_from_equalities(p); |
---|
| 144 | p = updatenonlinearbounds(p); |
---|
| 145 | |
---|
| 146 | output = yalmip_default_output; |
---|
| 147 | if p.feasible |
---|
| 148 | |
---|
| 149 | % ******************************* |
---|
| 150 | % Bounded domain? |
---|
| 151 | % ******************************* |
---|
| 152 | if ~isempty(p.bilinears) |
---|
| 153 | involved = unique(p.bilinears(:,2:3)); |
---|
| 154 | if any(isinf(p.lb(involved))) | any(isinf(p.ub(involved))) |
---|
| 155 | output = yalmip_default_output; |
---|
| 156 | output.Primal = x_min; |
---|
| 157 | output.problem = -6; |
---|
| 158 | output.infostr = yalmiperror(-6); |
---|
| 159 | output.solved_nodes = 0; |
---|
| 160 | return |
---|
| 161 | end |
---|
| 162 | end |
---|
| 163 | |
---|
| 164 | % ******************************* |
---|
| 165 | % Save time & memory |
---|
| 166 | % ******************************* |
---|
| 167 | p.options.savesolverinput = 0; |
---|
| 168 | p.options.savesolveroutput = 0; |
---|
| 169 | p.options.saveduals = 0; |
---|
| 170 | p.options.dimacs = 0; |
---|
| 171 | |
---|
| 172 | % ******************************* |
---|
| 173 | % RUN BILINEAR BRANCH & BOUND |
---|
| 174 | % ******************************* |
---|
| 175 | [x_min,solved_nodes,lower,upper] = branch_and_bound(p,x_min,upper); |
---|
| 176 | |
---|
| 177 | % ******************************** |
---|
| 178 | % CREATE SOLUTION AND DIAGNOSTICS |
---|
| 179 | % ******************************** |
---|
| 180 | problem = 0; |
---|
| 181 | if isinf(upper) |
---|
| 182 | problem = 1; |
---|
| 183 | end |
---|
| 184 | if isinf(lower) & (lower<0) |
---|
| 185 | problem = 2; |
---|
| 186 | end |
---|
| 187 | if isinf(lower) & (lower>0) |
---|
| 188 | problem = 1; |
---|
| 189 | end |
---|
| 190 | if solved_nodes == p.options.bnb.maxiter |
---|
| 191 | problem = 3; |
---|
| 192 | end |
---|
| 193 | else |
---|
| 194 | problem = 1; |
---|
| 195 | x_min = repmat(nan,length(p.c),1); |
---|
| 196 | solved_nodes = 0; |
---|
| 197 | end |
---|
| 198 | output = yalmip_default_output; |
---|
| 199 | output.problem = problem; |
---|
| 200 | output.solved_nodes = solved_nodes; |
---|
| 201 | output.Primal = x_min(1:n_in); |
---|
| 202 | output.infostr = yalmiperror(output.problem,'BNB'); |
---|
| 203 | output.solvertime = etime(clock,bnbsolvertime); |
---|
| 204 | |
---|
| 205 | function model = fixbounds(model); |
---|
| 206 | polynomials = find(model.variabletype > 0); |
---|
| 207 | LU = max([abs(model.lb) abs(model.ub)],[],2); |
---|
| 208 | for i = 1:length(polynomials) |
---|
| 209 | j = polynomials(i); |
---|
| 210 | if j<=length(model.lb) |
---|
| 211 | vars = find(model.monomtable(j,:)); |
---|
| 212 | bound = 1; |
---|
| 213 | for k = 1:length(vars) |
---|
| 214 | bound = bound * LU(vars(k))^model.monomtable(j,vars(k)); |
---|
| 215 | end |
---|
| 216 | model.lb(j) = max(model.lb(j),-bound); |
---|
| 217 | model.ub(j) = min(model.ub(j),bound); |
---|
| 218 | end |
---|
| 219 | end |
---|
| 220 | |
---|
| 221 | |
---|