[37] | 1 | function output = mpcvx(p) |
---|
| 2 | %BMIBNB Branch-and-bound scheme for bilinear programs |
---|
| 3 | % |
---|
| 4 | % BMIBNB is never called by the user directly, but is called by |
---|
| 5 | % YALMIP from SOLVESDP, by choosing the solver tag 'bmibnb' in sdpsettings |
---|
| 6 | % |
---|
| 7 | % The behaviour of BMIBNB can be altered using the fields |
---|
| 8 | % in the field 'bmibnb' in SDPSETTINGS |
---|
| 9 | % |
---|
| 10 | % WARNING: THIS IS EXPERIMENTAL CODE |
---|
| 11 | % |
---|
| 12 | % bmibnb.lowersolver - Solver for lower bound [standard solver tag ('')] |
---|
| 13 | % bmibnb.uppersolver - Solver for upper bound [standard solver tag ('')] |
---|
| 14 | % bmibnb.lpsolver - Solver for LP bound tightening [standard solver tag ('')] |
---|
| 15 | % bmibnb.branchmethod - Branch strategy ['maxvol' | 'best' ('best')] |
---|
| 16 | % bmibnb.branchrule - Branch position ['omega' | 'bisect' ('omega')] |
---|
| 17 | % bmibnb.lpreduce - Improve variable bounds using LP [ real [0,1] (0 means no reduction, 1 means all variables) |
---|
| 18 | % bmibnb.lowrank - partition variables into two disjoint sets and branch on smallest [ 0|1 (0)] |
---|
| 19 | % bmibnb.target - Exit if upper found<target [double (-inf)] |
---|
| 20 | % bmibnb.roottight - Improve variable bounds in root using full problem [ 0|1 (1)] |
---|
| 21 | % bmibnb.vartol - Cut tree when x_U-x_L < vartol on all branching variables |
---|
| 22 | % bmibnb.relgaptol - Tolerance on relative objective error (UPPER-LOWER)/(1+|UPPER|) [real (0.01)] |
---|
| 23 | % bmibnb.absgaptol - Tolerance on objective error (UPPER-LOWER) [real (0.01)] |
---|
| 24 | % bmibnb.pdtol - A number is declared non-negative if larger than...[ double (-1e-6)] |
---|
| 25 | % bmibnb.eqtol - A number is declared zero if abs(x) smaller than...[ double (1e-6)] |
---|
| 26 | % bmibnb.maxiter - Maximum number of solved nodes [int (100)] |
---|
| 27 | % bmibnb.maxtime - Maximum CPU time (sec.) [int (3600)] |
---|
| 28 | |
---|
| 29 | % Author Johan Löfberg |
---|
| 30 | % $Id: callmpcvx.m,v 1.2 2005/05/07 13:53:20 joloef Exp $ |
---|
| 31 | |
---|
| 32 | % ******************************** |
---|
| 33 | % INITIALIZE DIAGNOSTICS IN YALMIP |
---|
| 34 | % ******************************** |
---|
| 35 | bnbsolvertime = clock; |
---|
| 36 | showprogress('Branch and bound started',p.options.showprogress); |
---|
| 37 | |
---|
| 38 | % ******************************* |
---|
| 39 | % Display-logics |
---|
| 40 | % ******************************* |
---|
| 41 | switch max(min(p.options.verbose,3),0) |
---|
| 42 | case 0 |
---|
| 43 | p.options.bmibnb.verbose = 0; |
---|
| 44 | case 1 |
---|
| 45 | p.options.bmibnb.verbose = 1; |
---|
| 46 | p.options.verbose = 0; |
---|
| 47 | case 2 |
---|
| 48 | p.options.bmibnb.verbose = 2; |
---|
| 49 | p.options.verbose = 0; |
---|
| 50 | case 3 |
---|
| 51 | p.options.bmibnb.verbose = 2; |
---|
| 52 | p.options.verbose = 1; |
---|
| 53 | otherwise |
---|
| 54 | p.options.bmibnb.verbose = 0; |
---|
| 55 | p.options.verbose = 0; |
---|
| 56 | end |
---|
| 57 | |
---|
| 58 | % ******************************* |
---|
| 59 | % Actual linear variables |
---|
| 60 | % ******************************* |
---|
| 61 | p.linears = find(sum(p.monomtable,2)==1); |
---|
| 62 | |
---|
| 63 | % ******************************* |
---|
| 64 | % PRE-CALCULATE INDICIES |
---|
| 65 | % FOR BILNEAR VARIABLES |
---|
| 66 | % ******************************* |
---|
| 67 | nonlinear = find(~(sum(p.monomtable~=0,2)==1 & sum(p.monomtable,2)==1)); |
---|
| 68 | nonlins = []; |
---|
| 69 | for i = 1:length(nonlinear) |
---|
| 70 | z = find(p.monomtable(nonlinear(i),:)); |
---|
| 71 | if length(z)==1 |
---|
| 72 | nonlins = [nonlins;nonlinear(i) z z]; |
---|
| 73 | else |
---|
| 74 | nonlins = [nonlins;nonlinear(i) z(1) z(2)]; |
---|
| 75 | end |
---|
| 76 | end |
---|
| 77 | p.nonlins = nonlins; |
---|
| 78 | |
---|
| 79 | p.options.saveduals = 0; |
---|
| 80 | |
---|
| 81 | % ******************************* |
---|
| 82 | % Select branch variables |
---|
| 83 | % ******************************* |
---|
| 84 | p.branch_variables = decide_branch_variables(p); |
---|
| 85 | |
---|
| 86 | % ******************************** |
---|
| 87 | % Extract bounds from model |
---|
| 88 | % ******************************** |
---|
| 89 | if isempty(p.ub) |
---|
| 90 | p.ub = repmat(inf,length(p.c),1); |
---|
| 91 | end |
---|
| 92 | if isempty(p.lb) |
---|
| 93 | p.lb = repmat(-inf,length(p.c),1); |
---|
| 94 | end |
---|
| 95 | if ~isempty(p.F_struc) |
---|
| 96 | [lb,ub,used_rows] = findulb(p.F_struc,p.K); |
---|
| 97 | if ~isempty(used_rows) |
---|
| 98 | lower_defined = find(~isinf(lb)); |
---|
| 99 | if ~isempty(lower_defined) |
---|
| 100 | p.lb(lower_defined) = max(p.lb(lower_defined),lb(lower_defined)); |
---|
| 101 | end |
---|
| 102 | upper_defined = find(~isinf(ub)); |
---|
| 103 | if ~isempty(upper_defined) |
---|
| 104 | p.ub(upper_defined) = min(p.ub(upper_defined),ub(upper_defined)); |
---|
| 105 | end |
---|
| 106 | % Remove linear bounds |
---|
| 107 | used_rows = used_rows(find(~any(p.F_struc(p.K.f+used_rows,1+nonlinear),2))); |
---|
| 108 | not_used_rows = setdiff(1:p.K.l,used_rows); |
---|
| 109 | for i = 1:length(p.KCut.l) |
---|
| 110 | p.KCut.l(i) = find(not_used_rows==p.KCut.l(i) ); |
---|
| 111 | end |
---|
| 112 | if ~isempty(used_rows) |
---|
| 113 | p.F_struc(p.K.f+used_rows,:)=[]; |
---|
| 114 | p.K.l = p.K.l - length(used_rows); |
---|
| 115 | end |
---|
| 116 | end |
---|
| 117 | end |
---|
| 118 | p = clean_bounds(p); |
---|
| 119 | p = updatenonlinearbounds(p); |
---|
| 120 | feasible = all(p.lb<=p.ub); |
---|
| 121 | |
---|
| 122 | % ******************************** |
---|
| 123 | % Remove empty linear rows |
---|
| 124 | % ******************************** |
---|
| 125 | if p.K.l > 0 |
---|
| 126 | empty_rows = find(~any(p.F_struc(p.K.f+1:p.K.f+p.K.l,2:end),2)); |
---|
| 127 | if ~isempty(empty_rows) |
---|
| 128 | if all(p.F_struc(p.K.f+empty_rows,1)>=0) |
---|
| 129 | p.F_struc(p.K.f+empty_rows,:)=[]; |
---|
| 130 | p.K.l = p.K.l - length(empty_rows); |
---|
| 131 | else |
---|
| 132 | feasible = 0; |
---|
| 133 | end |
---|
| 134 | end |
---|
| 135 | end |
---|
| 136 | |
---|
| 137 | % ******************************** |
---|
| 138 | % Tighten bounds at root |
---|
| 139 | % ******************************** |
---|
| 140 | if p.options.bmibnb.roottight & feasible |
---|
| 141 | lowersolver = eval(['@' p.solver.lowercall]); |
---|
| 142 | c = p.c; |
---|
| 143 | Q = p.Q; |
---|
| 144 | mt = p.monomtable; |
---|
| 145 | p.monomtable = eye(length(c)); |
---|
| 146 | i = 1; |
---|
| 147 | while i<=length(p.linears) & feasible |
---|
| 148 | j = p.linears(i); |
---|
| 149 | p.c = eyev(length(p.c),j); |
---|
| 150 | output = feval(lowersolver,p); |
---|
| 151 | if (output.problem == 0) & (output.Primal(j)>p.lb(j)) |
---|
| 152 | p.lb(j) = output.Primal(j); |
---|
| 153 | p = updateonenonlinearbound(p,j); |
---|
| 154 | p = clean_bounds(p); |
---|
| 155 | end |
---|
| 156 | if output.problem == 1 |
---|
| 157 | feasible = 0; |
---|
| 158 | else |
---|
| 159 | % p = updatenonlinearbounds(p,0,1); |
---|
| 160 | p.c = -eyev(length(p.c),j); |
---|
| 161 | output = feval(lowersolver,p); |
---|
| 162 | if (output.problem == 0) & (output.Primal(j) < p.ub(j)) |
---|
| 163 | p.ub(j) = output.Primal(j); |
---|
| 164 | p = updateonenonlinearbound(p,j); |
---|
| 165 | p = clean_bounds(p); |
---|
| 166 | end |
---|
| 167 | if output.problem == 1 |
---|
| 168 | feasible = 0; |
---|
| 169 | end |
---|
| 170 | i = i+1; |
---|
| 171 | end |
---|
| 172 | end |
---|
| 173 | p.lb(p.lb<-1e10) = -inf; |
---|
| 174 | p.ub(p.ub>1e10) = inf; |
---|
| 175 | p.c = c; |
---|
| 176 | p.Q = Q; |
---|
| 177 | p.monomtable = mt; |
---|
| 178 | end |
---|
| 179 | |
---|
| 180 | if feasible |
---|
| 181 | |
---|
| 182 | % ******************************* |
---|
| 183 | % Bounded domain? |
---|
| 184 | % ******************************* |
---|
| 185 | involved = unique(p.nonlins(:,2:3)); |
---|
| 186 | if isinf(p.lb(involved)) | isinf(p.ub(involved)) |
---|
| 187 | error('You have to bound all complicating variables explicitely (I cannot deduce bounds on all variables)') |
---|
| 188 | output.Primal = []; |
---|
| 189 | output.problem = -1; |
---|
| 190 | end |
---|
| 191 | |
---|
| 192 | % ******************************* |
---|
| 193 | % We don't need to save node data |
---|
| 194 | % ******************************* |
---|
| 195 | p.options.savesolverinput = 0; |
---|
| 196 | p.options.savesolveroutput = 0; |
---|
| 197 | |
---|
| 198 | % ******************************* |
---|
| 199 | % RUN BRANCH & BOUND |
---|
| 200 | % ******************************* |
---|
| 201 | [x_min,solved_nodes,lower,upper] = branch_and_bound(p); |
---|
| 202 | |
---|
| 203 | % ********************************** |
---|
| 204 | % CREATE SOLUTION |
---|
| 205 | % ********************************** |
---|
| 206 | output.problem = 0; |
---|
| 207 | if isinf(upper) |
---|
| 208 | output.problem = 1; |
---|
| 209 | end |
---|
| 210 | if isinf(-lower) |
---|
| 211 | output.problem = 2; |
---|
| 212 | end |
---|
| 213 | if solved_nodes == p.options.bnb.maxiter |
---|
| 214 | output.problem = 3; |
---|
| 215 | end |
---|
| 216 | else |
---|
| 217 | output.problem = 1; |
---|
| 218 | x_min = repmat(nan,length(p.c),1); |
---|
| 219 | solved_nodes = 0; |
---|
| 220 | end |
---|
| 221 | |
---|
| 222 | output.solved_nodes = solved_nodes; |
---|
| 223 | output.Primal = x_min; |
---|
| 224 | output.Dual = []; |
---|
| 225 | output.Slack = []; |
---|
| 226 | output.infostr = yalmiperror(output.problem,'BNB'); |
---|
| 227 | output.solverinput = 0; |
---|
| 228 | output.solveroutput =[]; |
---|
| 229 | output.solvertime = etime(clock,bnbsolvertime); |
---|
| 230 | |
---|
| 231 | function [x_min,solved_nodes,lower,upper] = branch_and_bound(p) |
---|
| 232 | |
---|
| 233 | % *************************************** |
---|
| 234 | % LPs ARE USED IN BOX-REDUCTION |
---|
| 235 | % (this is essentially a cutting plane pool) |
---|
| 236 | % *************************************** |
---|
| 237 | p.lpcuts = p.F_struc(1+p.K.f:1:p.K.l,:); |
---|
| 238 | |
---|
| 239 | % *************************************** |
---|
| 240 | % Create function handles to solvers |
---|
| 241 | % *************************************** |
---|
| 242 | try |
---|
| 243 | lowersolver = eval(['@' p.solver.lowercall]); % Local LMI solver |
---|
| 244 | uppersolver = eval(['@' p.solver.uppercall]); % Local BMI solver |
---|
| 245 | lpsolver = eval(['@' p.solver.lpcall]); % LP solver |
---|
| 246 | catch |
---|
| 247 | disp(' '); |
---|
| 248 | disp('The internal branch & bound solver requires MATLAB 6') |
---|
| 249 | disp('I am too lazy too do the changes to make it compatible') |
---|
| 250 | disp('with MATLAB 5. If you really need it, contact me...'); |
---|
| 251 | disp(' '); |
---|
| 252 | error(lasterr); |
---|
| 253 | end |
---|
| 254 | |
---|
| 255 | % ************************************************ |
---|
| 256 | % GLOBAL PROBLEM DATA |
---|
| 257 | % ************************************************ |
---|
| 258 | c = p.c; |
---|
| 259 | Q = p.Q; |
---|
| 260 | f = p.f; |
---|
| 261 | K = p.K; |
---|
| 262 | p.options.saveduals = 0; |
---|
| 263 | options = p.options; |
---|
| 264 | |
---|
| 265 | % ************************************************ |
---|
| 266 | % ORIGINAL PROBLEM (used in LOCAL BMI solver) |
---|
| 267 | % ************************************************ |
---|
| 268 | p_upper = p; |
---|
| 269 | |
---|
| 270 | % ************************************************ |
---|
| 271 | % Remove linear cutting planes from problem |
---|
| 272 | % ************************************************ |
---|
| 273 | p_upper.F_struc(p_upper.K.f+p_upper.KCut.l,:)=[]; |
---|
| 274 | p_upper.K.l = p_upper.K.l - length(p_upper.KCut.l); |
---|
| 275 | |
---|
| 276 | % ************************************************ |
---|
| 277 | % Remove sdp cutting planes from problem |
---|
| 278 | % ************************************************ |
---|
| 279 | if length(p_upper.KCut.s)>0 |
---|
| 280 | starts = p_upper.K.f+p_upper.K.l + [1 1+cumsum((p_upper.K.s).^2)]; |
---|
| 281 | remove_these = []; |
---|
| 282 | for i = 1:length(p_upper.KCut.s) |
---|
| 283 | j = p_upper.KCut.s(i); |
---|
| 284 | remove_these = [remove_these;(starts(j):starts(j+1)-1)']; |
---|
| 285 | end |
---|
| 286 | p_upper.F_struc(remove_these,:)=[]; |
---|
| 287 | p_upper.K.s(p_upper.KCut.s) = []; |
---|
| 288 | end |
---|
| 289 | |
---|
| 290 | % ************************************************ |
---|
| 291 | % INITILIZATION |
---|
| 292 | % ************************************************ |
---|
| 293 | p.depth = 0; |
---|
| 294 | p.dpos = 0; |
---|
| 295 | p.lower = NaN; |
---|
| 296 | upper = inf; |
---|
| 297 | lower = NaN; |
---|
| 298 | gap = inf; |
---|
| 299 | x_min = zeros(length(p.c),1); |
---|
| 300 | stack = []; |
---|
| 301 | solved_nodes = 0; |
---|
| 302 | solved_lower = 0; |
---|
| 303 | solved_upper = 0; |
---|
| 304 | solved_lp = 0; |
---|
| 305 | |
---|
| 306 | if isempty(p.x0) |
---|
| 307 | p.x0 = zeros(length(p.c),1); |
---|
| 308 | end |
---|
| 309 | |
---|
| 310 | x0 = evaluate_nonlinear(p,p.x0); |
---|
| 311 | upper_residual = resids(p,x0); |
---|
| 312 | x0_feasible = all(upper_residual(1:p.K.f)>=-options.bmibnb.eqtol) & all(upper_residual(1+p.K.f:end)>=options.bmibnb.pdtol); |
---|
| 313 | if p.options.usex0 & x0_feasible |
---|
| 314 | x_min = x0; |
---|
| 315 | upper = p.f+p.c'*x0+x0'*Q*x0; |
---|
| 316 | end |
---|
| 317 | |
---|
| 318 | % ************************************************ |
---|
| 319 | % Branch & bound loop |
---|
| 320 | % ************************************************ |
---|
| 321 | if options.bmibnb.verbose>0 |
---|
| 322 | fprintf('******************************************************************************************************************\n') |
---|
| 323 | fprintf('#node Was'' up gap upper node lower dpth stk Memory Vol-red\n') |
---|
| 324 | fprintf('******************************************************************************************************************\n') |
---|
| 325 | end |
---|
| 326 | |
---|
| 327 | doplot = 0; |
---|
| 328 | if doplot |
---|
| 329 | close all; |
---|
| 330 | hold on; |
---|
| 331 | end |
---|
| 332 | |
---|
| 333 | t_start = cputime; |
---|
| 334 | go_on = 1; |
---|
| 335 | while go_on |
---|
| 336 | |
---|
| 337 | if doplot;ellipplot(diag([200 50]),1,'y',[p.dpos;-p.depth]);drawnow;end; |
---|
| 338 | % ******************************************** |
---|
| 339 | % ASSUME THAT WE WON'T FATHOME |
---|
| 340 | % ******************************************** |
---|
| 341 | keep_digging = 1; |
---|
| 342 | % ******************************************** |
---|
| 343 | % REDUCE BOX |
---|
| 344 | % ******************************************** |
---|
| 345 | if ~options.bmibnb.lpreduce |
---|
| 346 | % [p.lb,p.ub] = tightenbounds(-p.F_struc(1+p.K.f:p.K.f+p.K.l,2:end),p.F_struc(1+p.K.f:p.K.f+p.K.l,1),p.lb,p.ub,[]); |
---|
| 347 | vol_reduction = 1; |
---|
| 348 | feasible = 1; |
---|
| 349 | else |
---|
| 350 | [p,feasible,vol_reduction] = boxreduce(p,upper,lower,lpsolver,options); |
---|
| 351 | end |
---|
| 352 | |
---|
| 353 | % ******************************************** |
---|
| 354 | % SOLVE LOWER AND UPPER |
---|
| 355 | % ******************************************** |
---|
| 356 | if feasible |
---|
| 357 | |
---|
| 358 | output = solvelower(p,options,lowersolver); |
---|
| 359 | |
---|
| 360 | info_text = ''; |
---|
| 361 | switch output.problem |
---|
| 362 | case 1 |
---|
| 363 | if doplot;ellipplot(diag([200 25]),1,'r',[p.dpos;-p.depth]);drawnow;end; |
---|
| 364 | info_text = 'Infeasible node'; |
---|
| 365 | keep_digging = 0; |
---|
| 366 | cost = inf; |
---|
| 367 | feasible = 0; |
---|
| 368 | |
---|
| 369 | case 2 |
---|
| 370 | cost = -inf; |
---|
| 371 | |
---|
| 372 | case {0,3,4} |
---|
| 373 | |
---|
| 374 | x = output.Primal; |
---|
| 375 | |
---|
| 376 | cost = f+c'*x+x'*Q*x; |
---|
| 377 | |
---|
| 378 | z = evaluate_nonlinear(p,x); |
---|
| 379 | p = addsdpcut(p,z); |
---|
| 380 | |
---|
| 381 | % Maybe the relaxed solution is feasible |
---|
| 382 | relaxed_residual = resids(p_upper,z); |
---|
| 383 | relaxed_feasible = all(relaxed_residual(1:p.K.f)>=-options.bmibnb.eqtol) & all(relaxed_residual(1+p.K.f:end)>=options.bmibnb.pdtol); |
---|
| 384 | if relaxed_feasible |
---|
| 385 | this_upper = f+c'*z+z'*Q*z; |
---|
| 386 | if (this_upper < (1-1e-5)*upper) & (this_upper < upper - 1e-5) |
---|
| 387 | x_min = z; |
---|
| 388 | upper = this_upper; |
---|
| 389 | info_text = 'Improved upper bound'; |
---|
| 390 | cost = cost-1e-10; % Otherwise we'll fathome! |
---|
| 391 | end |
---|
| 392 | end |
---|
| 393 | |
---|
| 394 | % UPDATE THE LOWER BOUND |
---|
| 395 | if isnan(lower) |
---|
| 396 | lower = cost; |
---|
| 397 | end |
---|
| 398 | if ~isempty(stack) |
---|
| 399 | lower =min(cost,min([stack.lower])); |
---|
| 400 | else |
---|
| 401 | lower = min(lower,cost); |
---|
| 402 | end |
---|
| 403 | |
---|
| 404 | if cost<upper |
---|
| 405 | output = solveupper(p,p_upper,x,options,uppersolver); |
---|
| 406 | xu = evaluate_nonlinear(p,output.Primal); |
---|
| 407 | upper_residual = resids(p_upper,xu); |
---|
| 408 | if output.problem == 0 | (all(upper_residual(1:p_upper.K.f)>=-options.bmibnb.eqtol) & all(upper_residual(1+p_upper.K.f:end)>=options.bmibnb.pdtol)) |
---|
| 409 | this_upper = f+c'*xu+xu'*Q*xu; |
---|
| 410 | if (this_upper < (1-1e-5)*upper) & (this_upper < upper - 1e-5) |
---|
| 411 | x_min = xu; |
---|
| 412 | upper = this_upper; |
---|
| 413 | info_text = 'Improved upper bound'; |
---|
| 414 | end |
---|
| 415 | end |
---|
| 416 | else |
---|
| 417 | if doplot;ellipplot(diag([200 25]),1,'b',[p.dpos;-p.depth]);drawnow;end |
---|
| 418 | info_text = 'Poor lower bound'; |
---|
| 419 | keep_digging = 0; |
---|
| 420 | end |
---|
| 421 | otherwise |
---|
| 422 | end |
---|
| 423 | else |
---|
| 424 | if doplot;ellipplot(diag([200 25]),1,'r',[p.dpos;-p.depth]);drawnow;end |
---|
| 425 | info_text = 'Infeasible during box-reduction'; |
---|
| 426 | keep_digging = 0; |
---|
| 427 | cost = inf; |
---|
| 428 | feasible = 0; |
---|
| 429 | end |
---|
| 430 | solved_nodes = solved_nodes+1; |
---|
| 431 | |
---|
| 432 | if ~isempty(stack) |
---|
| 433 | [stack,lower] = prune(stack,upper,options,solved_nodes); |
---|
| 434 | end |
---|
| 435 | lower = min(lower,cost); |
---|
| 436 | |
---|
| 437 | % ********************************** |
---|
| 438 | % CONTINUE SPLITTING? |
---|
| 439 | % ********************************** |
---|
| 440 | if keep_digging & max(p.ub(p.branch_variables)-p.lb(p.branch_variables))>options.bmibnb.vartol |
---|
| 441 | spliton = branchvariable(p,options,x); |
---|
| 442 | bounds = partition(p,options,spliton,x_min); |
---|
| 443 | |
---|
| 444 | node_1 = savetonode(p,spliton,bounds(1),bounds(2),-1,x,cost); |
---|
| 445 | node_2 = savetonode(p,spliton,bounds(2),bounds(3),1,x,cost); |
---|
| 446 | stack = push(stack,node_1); |
---|
| 447 | stack = push(stack,node_2); |
---|
| 448 | end |
---|
| 449 | |
---|
| 450 | % Pick and create a suitable node to continue on |
---|
| 451 | [p,stack] = selectbranch(p,options,stack,x_min,upper); |
---|
| 452 | |
---|
| 453 | if isempty(p) |
---|
| 454 | if ~isinf(upper) |
---|
| 455 | relgap = 0; |
---|
| 456 | end |
---|
| 457 | depth = 0; |
---|
| 458 | else |
---|
| 459 | relgap = 100*(upper-lower)/(1+abs(upper)); |
---|
| 460 | depth = p.depth; |
---|
| 461 | end |
---|
| 462 | if options.bmibnb.verbose>0 |
---|
| 463 | ws = whos; %Work-space |
---|
| 464 | Mb = sum([ws(:).bytes])/1024^2; %Megs |
---|
| 465 | showprogress(sprintf(['%3d ' info_text repmat(' ',1,35-length(info_text)) ' %8.2f%% %8.4f %8.4f %8.4f %3d %3d %5.2fMB %4.1f%% '],solved_nodes,relgap,upper,cost,lower,depth,length(stack),Mb,100-vol_reduction*100),options.bmibnb.verbose) |
---|
| 466 | end |
---|
| 467 | |
---|
| 468 | absgap = upper-lower; |
---|
| 469 | % ************************************** |
---|
| 470 | % Continue? |
---|
| 471 | % ************************************** |
---|
| 472 | time_ok = cputime-t_start < options.bmibnb.maxtime; |
---|
| 473 | iter_ok = solved_nodes < options.bmibnb.maxiter; |
---|
| 474 | any_nodes = ~isempty(p); |
---|
| 475 | relgap_too_big = (isinf(lower) | isnan(relgap) | relgap>100*options.bmibnb.relgaptol); |
---|
| 476 | absgap_too_big = (isinf(lower) | isnan(absgap) | absgap>options.bmibnb.absgaptol); |
---|
| 477 | taget_not_met = upper>options.bmibnb.target; |
---|
| 478 | go_on = taget_not_met & time_ok & any_nodes & iter_ok & relgap_too_big & absgap_too_big ; |
---|
| 479 | end |
---|
| 480 | if options.bmibnb.verbose>0 |
---|
| 481 | fprintf('******************************************************************************************************************\n') |
---|
| 482 | if options.bmibnb.verbose;showprogress([num2str2(solved_nodes,3) ' Finishing. Cost: ' num2str(upper) ' Gap: ' num2str(relgap) '%'],options.bnb.verbose);end |
---|
| 483 | fprintf('******************************************************************************************************************\n') |
---|
| 484 | end |
---|
| 485 | |
---|
| 486 | |
---|
| 487 | |
---|
| 488 | function stack = push(stackin,p) |
---|
| 489 | if ~isempty(stackin) |
---|
| 490 | stack = [p;stackin]; |
---|
| 491 | else |
---|
| 492 | stack(1)=p; |
---|
| 493 | end |
---|
| 494 | |
---|
| 495 | function [p,stack] = pull(stack,method,x_min,upper); |
---|
| 496 | if ~isempty(stack) |
---|
| 497 | switch method |
---|
| 498 | case 'maxvol' |
---|
| 499 | for i = 1:length(stack) |
---|
| 500 | vol(i) = sum(stack(i).ub(stack(i).branch_variables)-stack(i).lb(stack(i).branch_variables)); |
---|
| 501 | end |
---|
| 502 | [i,j] = max(vol); |
---|
| 503 | p=stack(j); |
---|
| 504 | stack = stack([1:1:j-1 j+1:1:end]); |
---|
| 505 | |
---|
| 506 | case 'best' |
---|
| 507 | [i,j]=min([stack.lower]); |
---|
| 508 | p=stack(j); |
---|
| 509 | stack = stack([1:1:j-1 j+1:1:end]); |
---|
| 510 | |
---|
| 511 | otherwise |
---|
| 512 | end |
---|
| 513 | else |
---|
| 514 | p = []; |
---|
| 515 | end |
---|
| 516 | |
---|
| 517 | |
---|
| 518 | function s = num2str2(x,d,c); |
---|
| 519 | if nargin==3 |
---|
| 520 | s = num2str(x,c); |
---|
| 521 | else |
---|
| 522 | s = num2str(x); |
---|
| 523 | end |
---|
| 524 | s = [repmat(' ',1,d-length(s)) s]; |
---|
| 525 | |
---|
| 526 | |
---|
| 527 | function res = resids(p,x) |
---|
| 528 | res= []; |
---|
| 529 | if p.K.f>0 |
---|
| 530 | res = -abs(p.F_struc(1:p.K.f,:)*[1;x]); |
---|
| 531 | end |
---|
| 532 | if p.K.l>0 |
---|
| 533 | res = [res;p.F_struc(p.K.f+1:p.K.f+p.K.l,:)*[1;x]]; |
---|
| 534 | end |
---|
| 535 | if (length(p.K.s)>1) | p.K.s>0 |
---|
| 536 | top = 1+p.K.f+p.K.l; |
---|
| 537 | for i = 1:length(p.K.s) |
---|
| 538 | n = p.K.s(i); |
---|
| 539 | X = p.F_struc(top:top+n^2-1,:)*[1;x];top = top+n^2; |
---|
| 540 | X = reshape(X,n,n); |
---|
| 541 | res = [res;min(eig(X))]; |
---|
| 542 | end |
---|
| 543 | end |
---|
| 544 | |
---|
| 545 | res = [res;min([p.ub-x;x-p.lb])]; |
---|
| 546 | |
---|
| 547 | |
---|
| 548 | function [stack,lower] = prune(stack,upper,options,solved_nodes) |
---|
| 549 | % ********************************* |
---|
| 550 | % PRUNE STACK W.R.T NEW UPPER BOUND |
---|
| 551 | % ********************************* |
---|
| 552 | if ~isempty(stack) |
---|
| 553 | toolarge = find([stack.lower]>upper*(1-1e-4)); |
---|
| 554 | if ~isempty(toolarge) |
---|
| 555 | if options.bnb.verbose;showprogress([num2str2(solved_nodes,3) ' Pruned ' num2str(length(toolarge)) ' nodes'],options.bnb.verbose-1);end |
---|
| 556 | stack(toolarge)=[]; |
---|
| 557 | end |
---|
| 558 | end |
---|
| 559 | if ~isempty(stack) |
---|
| 560 | lower = min([stack.lower]); |
---|
| 561 | else |
---|
| 562 | lower = upper; |
---|
| 563 | end |
---|
| 564 | |
---|
| 565 | function pcut = addmcgormick(p) |
---|
| 566 | |
---|
| 567 | pcut = p; |
---|
| 568 | top = 0; |
---|
| 569 | row = []; |
---|
| 570 | col = []; |
---|
| 571 | val = []; |
---|
| 572 | F_temp = []; |
---|
| 573 | for i = 1:size(p.nonlins,1) |
---|
| 574 | z = p.nonlins(i,1); |
---|
| 575 | x = p.nonlins(i,2); |
---|
| 576 | y = p.nonlins(i,3); |
---|
| 577 | x_lb = p.lb(x); |
---|
| 578 | x_ub = p.ub(x); |
---|
| 579 | y_lb = p.lb(y); |
---|
| 580 | y_ub = p.ub(y); |
---|
| 581 | top = 0; |
---|
| 582 | row = []; |
---|
| 583 | col = []; |
---|
| 584 | val = []; |
---|
| 585 | |
---|
| 586 | if x~=y |
---|
| 587 | row = [1;1;1;1;2;2;2;2;3;3;3;3;4;4;4;4]; |
---|
| 588 | col = [1 ; z+1 ; x+1 ; y+1 ; 1 ; z+1 ; x+1 ; y+1 ; 1 ; z+1 ; x+1 ; y+1 ; 1 ; z+1 ; x+1 ; y+1]; |
---|
| 589 | val = [x_lb*y_lb;1;-y_lb;-x_lb;x_ub*y_ub;1;-y_ub;-x_ub;-x_ub*y_lb;-1;y_lb;x_ub;-x_lb*y_ub;-1;y_ub;x_lb]; |
---|
| 590 | F_temp = [F_temp;sparse(row,col,val,4,size(pcut.F_struc,2))]; |
---|
| 591 | else |
---|
| 592 | |
---|
| 593 | nr = 3; |
---|
| 594 | row = [1;1;1;2;2 ;2; 3; 3; 3]; |
---|
| 595 | col = [1 ;z+1 ;x+1 ;1 ;z+1 ;x+1 ;1 ;z+1 ;x+1]; |
---|
| 596 | val = [-x_ub*x_lb;-1;x_lb+x_ub;x_lb*y_lb;1;-y_lb-x_lb;x_ub*y_ub;1;-y_ub-x_ub]; |
---|
| 597 | |
---|
| 598 | F_temp = [F_temp;sparse(row,col,val,nr,1+length(p.c))]; |
---|
| 599 | end |
---|
| 600 | bounds = [x_lb*y_lb x_lb*y_ub x_ub*y_lb x_ub*y_ub]; |
---|
| 601 | if x==y |
---|
| 602 | pcut.lb(pcut.nonlins(i,1)) = max(pcut.lb(pcut.nonlins(i,1)),max(0,min(bounds))); |
---|
| 603 | else |
---|
| 604 | pcut.lb(pcut.nonlins(i,1)) = max(pcut.lb(pcut.nonlins(i,1)),min(bounds)); |
---|
| 605 | end |
---|
| 606 | pcut.ub(pcut.nonlins(i,1)) = min(pcut.ub(pcut.nonlins(i,1)),max(bounds)); |
---|
| 607 | end |
---|
| 608 | |
---|
| 609 | keep = find(~isinf(F_temp(:,1))); |
---|
| 610 | F_temp = F_temp(keep,:); |
---|
| 611 | pcut.F_struc = [F_temp;pcut.F_struc]; |
---|
| 612 | pcut.K.l = pcut.K.l+size(F_temp,1); |
---|
| 613 | |
---|
| 614 | function [p,feasible,lower] = lpbmitighten(p,lower,upper,lpsolver) |
---|
| 615 | |
---|
| 616 | % Construct problem with only linear terms |
---|
| 617 | % and add cuts from lower/ upper bounds |
---|
| 618 | c = p.c; |
---|
| 619 | p_test = p; |
---|
| 620 | p_test.K.s = 0; |
---|
| 621 | p_test.F_struc = p_test.F_struc(1+p_test.K.f:1:p_test.K.l+p_test.K.f,:); |
---|
| 622 | |
---|
| 623 | if ~isnan(lower) |
---|
| 624 | p_test.F_struc = [-(p.lower-abs(p.lower)*0.01) p_test.c';p_test.F_struc]; |
---|
| 625 | end |
---|
| 626 | if upper < inf |
---|
| 627 | p_test.F_struc = [upper+abs(upper)*0.01 -p_test.c';p_test.F_struc]; |
---|
| 628 | end |
---|
| 629 | |
---|
| 630 | p_test.F_struc = [p_test.lpcuts;p_test.F_struc]; |
---|
| 631 | p_test.K.l = size(p_test.F_struc,1); |
---|
| 632 | |
---|
| 633 | % Add cuts for nonlinear terms |
---|
| 634 | p_test = addmcgormick(p_test); |
---|
| 635 | |
---|
| 636 | p_test.F_struc = [p.F_struc(1:1:p.K.f,:);p_test.F_struc]; |
---|
| 637 | |
---|
| 638 | |
---|
| 639 | feasible = 1; |
---|
| 640 | i = 1; |
---|
| 641 | |
---|
| 642 | p_test = clean_bounds(p_test); |
---|
| 643 | |
---|
| 644 | j = 1; |
---|
| 645 | n = ceil(max(p.options.bmibnb.lpreduce*length(p_test.linears),1)); |
---|
| 646 | res = zeros(length(p.lb),1); |
---|
| 647 | for i = 1:size(p.nonlins,1) |
---|
| 648 | res(p.nonlins(i,2)) = res(p.nonlins(i,2)) + abs( p.x0(p.nonlins(i,1))-p.x0(p.nonlins(i,2)).*p.x0(p.nonlins(i,3))); |
---|
| 649 | res(p.nonlins(i,3)) = res(p.nonlins(i,3)) + abs( p.x0(p.nonlins(i,1))-p.x0(p.nonlins(i,2)).*p.x0(p.nonlins(i,3))); |
---|
| 650 | end |
---|
| 651 | res = res(p.linears); |
---|
| 652 | [ii,jj] = sort(abs(res)); |
---|
| 653 | jj = jj(end-n+1:end); |
---|
| 654 | |
---|
| 655 | while feasible & j<=length(jj) |
---|
| 656 | i = p_test.linears(jj(j)); |
---|
| 657 | if abs(p.ub(i)-p.lb(i)>0.1) |
---|
| 658 | p_test.c = eyev(length(p_test.c),i); |
---|
| 659 | |
---|
| 660 | output = feval(lpsolver,p_test); |
---|
| 661 | if output.problem == 0 |
---|
| 662 | if p_test.lb(i) < output.Primal(i)-1e-5 |
---|
| 663 | p_test.lb(i) = output.Primal(i); |
---|
| 664 | p_test = updateonenonlinearbound(p_test,i); |
---|
| 665 | end |
---|
| 666 | p_test.c = -eyev(length(p_test.c),i); |
---|
| 667 | output = feval(lpsolver,p_test); |
---|
| 668 | if output.problem == 0 |
---|
| 669 | if p_test.ub(i) > output.Primal(i)+1e-5 |
---|
| 670 | p_test.ub(i) = output.Primal(i); |
---|
| 671 | p_test = updateonenonlinearbound(p_test,i); |
---|
| 672 | end |
---|
| 673 | end |
---|
| 674 | if output.problem == 1 |
---|
| 675 | feasible = 0; |
---|
| 676 | end |
---|
| 677 | end |
---|
| 678 | if output.problem == 1 |
---|
| 679 | feasible = 0; |
---|
| 680 | end |
---|
| 681 | p_test = clean_bounds(p_test); |
---|
| 682 | end |
---|
| 683 | j = j + 1; |
---|
| 684 | end |
---|
| 685 | p.lb = p_test.lb; |
---|
| 686 | p.ub = p_test.ub; |
---|
| 687 | |
---|
| 688 | |
---|
| 689 | function p = updateonenonlinearbound(p,changed_var); |
---|
| 690 | for i = 1:size(p.nonlins,1) |
---|
| 691 | x = p.nonlins(i,2); |
---|
| 692 | y = p.nonlins(i,3); |
---|
| 693 | if (x==changed_var) | (y==changed_var) |
---|
| 694 | z = p.nonlins(i,1); |
---|
| 695 | x_lb = p.lb(x); |
---|
| 696 | x_ub = p.ub(x); |
---|
| 697 | y_lb = p.lb(y); |
---|
| 698 | y_ub = p.ub(y); |
---|
| 699 | bounds = [x_lb*y_lb x_lb*y_ub x_ub*y_lb x_ub*y_ub]; |
---|
| 700 | if x==y |
---|
| 701 | p.lb(p.nonlins(i,1)) = max([p.lb(z) 0 min(bounds)]); |
---|
| 702 | p.ub(p.nonlins(i,1)) = min(p.ub(z),max(bounds)); |
---|
| 703 | else |
---|
| 704 | p.lb(p.nonlins(i,1)) = max(p.lb(z),min(bounds)); |
---|
| 705 | p.ub(p.nonlins(i,1)) = min(p.ub(z),max(bounds)); |
---|
| 706 | end |
---|
| 707 | end |
---|
| 708 | end |
---|
| 709 | |
---|
| 710 | |
---|
| 711 | function p = updatenonlinearbounds(p,changed_var,keepbest); |
---|
| 712 | % if nargin>1 |
---|
| 713 | % changed_var |
---|
| 714 | % else |
---|
| 715 | % i = 1:size(p.nonlins,1); |
---|
| 716 | % end |
---|
| 717 | for i = 1:size(p.nonlins,1) |
---|
| 718 | z = p.nonlins(i,1); |
---|
| 719 | x = p.nonlins(i,2); |
---|
| 720 | y = p.nonlins(i,3); |
---|
| 721 | x_lb = p.lb(x); |
---|
| 722 | x_ub = p.ub(x); |
---|
| 723 | y_lb = p.lb(y); |
---|
| 724 | y_ub = p.ub(y); |
---|
| 725 | bounds = [x_lb*y_lb x_lb*y_ub x_ub*y_lb x_ub*y_ub]; |
---|
| 726 | if x==y |
---|
| 727 | p.lb(p.nonlins(i,1)) = max([p.lb(z) 0 min(bounds)]); |
---|
| 728 | p.ub(p.nonlins(i,1)) = min(p.ub(z),max(bounds)); |
---|
| 729 | else |
---|
| 730 | p.lb(p.nonlins(i,1)) = max(p.lb(z),min(bounds)); |
---|
| 731 | p.ub(p.nonlins(i,1)) = min(p.ub(z),max(bounds)); |
---|
| 732 | end |
---|
| 733 | end |
---|
| 734 | return |
---|
| 735 | |
---|
| 736 | if nargin > 1 |
---|
| 737 | |
---|
| 738 | for i = 1:size(p.nonlins,1) |
---|
| 739 | z = p.nonlins(i,1); |
---|
| 740 | x = p.nonlins(i,2); |
---|
| 741 | y = p.nonlins(i,3); |
---|
| 742 | if isempty(changed_var) | (x==changed_var) | (y == changed_var) | nargin==3 |
---|
| 743 | bound_x1 = [p.lb(p.nonlins(i,2));p.ub(p.nonlins(i,2))]; |
---|
| 744 | bound_x2 = [p.lb(p.nonlins(i,3));p.ub(p.nonlins(i,3))]; |
---|
| 745 | bounds = [bound_x1(1)*bound_x2(1) bound_x1(1)*bound_x2(2) bound_x1(2)*bound_x2(1) bound_x1(2)*bound_x2(2)]; |
---|
| 746 | if nargin==3 |
---|
| 747 | if x==y |
---|
| 748 | p.lb(p.nonlins(i,1)) = max([p.lb(p.nonlins(i,1)) 0 min(bounds)]); |
---|
| 749 | p.ub(p.nonlins(i,1)) = min(p.ub(p.nonlins(i,1)),max(bounds)); |
---|
| 750 | else |
---|
| 751 | p.lb(p.nonlins(i,1)) = max(p.lb(p.nonlins(i,1)),min(bounds)); |
---|
| 752 | p.ub(p.nonlins(i,1)) = min(p.ub(p.nonlins(i,1)),max(bounds)); |
---|
| 753 | end |
---|
| 754 | else |
---|
| 755 | if x==y |
---|
| 756 | p.lb(p.nonlins(i,1)) = max(0,min(bounds)); |
---|
| 757 | p.ub(p.nonlins(i,1)) = max(bounds); |
---|
| 758 | else |
---|
| 759 | p.lb(p.nonlins(i,1)) = min(bounds); |
---|
| 760 | p.ub(p.nonlins(i,1)) = max(bounds); |
---|
| 761 | end |
---|
| 762 | end |
---|
| 763 | end |
---|
| 764 | end |
---|
| 765 | else |
---|
| 766 | for i = 1:size(p.nonlins,1) |
---|
| 767 | z = p.nonlins(i,1); |
---|
| 768 | x = p.nonlins(i,2); |
---|
| 769 | y = p.nonlins(i,3); |
---|
| 770 | bound_x1 = [p.lb(p.nonlins(i,2));p.ub(p.nonlins(i,2))]; |
---|
| 771 | bound_x2 = [p.lb(p.nonlins(i,3));p.ub(p.nonlins(i,3))]; |
---|
| 772 | bounds = [bound_x1(1)*bound_x2(1) bound_x1(1)*bound_x2(2) bound_x1(2)*bound_x2(1) bound_x1(2)*bound_x2(2)]; |
---|
| 773 | if x==y |
---|
| 774 | p.lb(p.nonlins(i,1)) = max( p.lb(p.nonlins(i,1)) ,max(0,min(bounds))); |
---|
| 775 | p.ub(p.nonlins(i,1)) = min( p.ub(p.nonlins(i,1)) ,max(bounds)); |
---|
| 776 | else |
---|
| 777 | p.lb(p.nonlins(i,1)) = max(p.lb(p.nonlins(i,1)),min(bounds)); |
---|
| 778 | p.ub(p.nonlins(i,1)) = min(p.ub(p.nonlins(i,1)),max(bounds)); |
---|
| 779 | end |
---|
| 780 | end |
---|
| 781 | end |
---|
| 782 | |
---|
| 783 | % ************************************* |
---|
| 784 | % DERIVE LINEAR CUTS FROM SDPs |
---|
| 785 | % THESE ARE ONLY USED IN BOXREDUCE |
---|
| 786 | % ************************************* |
---|
| 787 | function p = addsdpcut(p,x) |
---|
| 788 | if p.K.s > 0 |
---|
| 789 | top = p.K.f+p.K.l+1; |
---|
| 790 | newcuts = 1; |
---|
| 791 | newF = []; |
---|
| 792 | for i = 1:length(p.K.s) |
---|
| 793 | n = p.K.s(i); |
---|
| 794 | X = p.F_struc(top:top+n^2-1,:)*[1;x]; |
---|
| 795 | X = reshape(X,n,n); |
---|
| 796 | [d,v] = eig(X); |
---|
| 797 | for m = 1:length(v) |
---|
| 798 | if v(m,m)<0 |
---|
| 799 | for j = 1:length(x)+1; |
---|
| 800 | newF(newcuts,j)= d(:,m)'*reshape(p.F_struc(top:top+n^2-1,j),n,n)*d(:,m); |
---|
| 801 | end |
---|
| 802 | % max(abs(newF(:,2:end)),[],2) |
---|
| 803 | newF(newcuts,1)=newF(newcuts,1)+1e-6; |
---|
| 804 | newcuts = newcuts + 1; |
---|
| 805 | if size(p.lpcuts,1)>0 |
---|
| 806 | dist = p.lpcuts*newF(newcuts-1,:)'/(newF(newcuts-1,:)*newF(newcuts-1,:)'); |
---|
| 807 | if any(abs(dist-1)<1e-3) |
---|
| 808 | newF = newF(1:end-1,:); |
---|
| 809 | newcuts = newcuts - 1; |
---|
| 810 | end |
---|
| 811 | end |
---|
| 812 | end |
---|
| 813 | end |
---|
| 814 | top = top+n^2; |
---|
| 815 | end |
---|
| 816 | |
---|
| 817 | if ~isempty(newF) |
---|
| 818 | % Don't keep all |
---|
| 819 | m = size(newF,2); |
---|
| 820 | % size(p.lpcuts) |
---|
| 821 | p.lpcuts = [newF;p.lpcuts]; |
---|
| 822 | violations = p.lpcuts*[1;x]; |
---|
| 823 | p.lpcuts = p.lpcuts(violations<0.1,:); |
---|
| 824 | |
---|
| 825 | if size(p.lpcuts,1)>15*m |
---|
| 826 | violations = p.lpcuts*[1;x]; |
---|
| 827 | [i,j] = sort(violations); |
---|
| 828 | p.lpcuts = p.lpcuts(j(1:15*m),:); |
---|
| 829 | p.lpcuts = p.lpcuts(end-15*m+1:end,:); |
---|
| 830 | end |
---|
| 831 | end |
---|
| 832 | end |
---|
| 833 | |
---|
| 834 | |
---|
| 835 | function spliton = branchvariable(p,options,x) |
---|
| 836 | |
---|
| 837 | % Split if box is to narrow |
---|
| 838 | width = abs(p.ub(p.branch_variables)-p.lb(p.branch_variables)); |
---|
| 839 | if min(width)/max(width) < 0.1 |
---|
| 840 | [i,j] = max(width); |
---|
| 841 | spliton = p.branch_variables(j); |
---|
| 842 | else |
---|
| 843 | % res = zeros(length(p.lb),1); |
---|
| 844 | % for i = 1:size(p.nonlins,1) |
---|
| 845 | % res(p.nonlins(i,2)) = res(p.nonlins(i,2)) + abs( x(p.nonlins(i,1))-x(p.nonlins(i,2)).*x(p.nonlins(i,3))); |
---|
| 846 | % res(p.nonlins(i,3)) = res(p.nonlins(i,3)) + abs( x(p.nonlins(i,1))-x(p.nonlins(i,2)).*x(p.nonlins(i,3))); |
---|
| 847 | % end |
---|
| 848 | % |
---|
| 849 | % [ii,jj] = sort(abs(res)); |
---|
| 850 | % spliton = jj(end); |
---|
| 851 | |
---|
| 852 | res = x(p.nonlins(:,1))-x(p.nonlins(:,2)).*x(p.nonlins(:,3)); |
---|
| 853 | [ii,jj] = sort(abs(res)); |
---|
| 854 | v1 = p.nonlins(jj(end),2); |
---|
| 855 | v2 = p.nonlins(jj(end),3); |
---|
| 856 | |
---|
| 857 | acc_res1 = sum(abs(res(find((p.nonlins(:,2)==v1) | p.nonlins(:,3)==v1)))); |
---|
| 858 | acc_res2 = sum(abs(res(find((p.nonlins(:,2)==v2) | p.nonlins(:,3)==v2)))); |
---|
| 859 | |
---|
| 860 | if (acc_res1>acc_res2) & ismember(v1,p.branch_variables) |
---|
| 861 | spliton = v1; |
---|
| 862 | else |
---|
| 863 | spliton = v2; |
---|
| 864 | end |
---|
| 865 | end |
---|
| 866 | |
---|
| 867 | function bounds = partition(p,options,spliton,x_min) |
---|
| 868 | |
---|
| 869 | switch options.bmibnb.branchrule |
---|
| 870 | case 'omega' |
---|
| 871 | if ~isempty(x_min) |
---|
| 872 | bounds = [p.lb(spliton) 0.5*max(p.lb(spliton),min(x_min(spliton),p.ub(spliton)))+0.5*(p.lb(spliton)+p.ub(spliton))/2 p.ub(spliton)]; |
---|
| 873 | else |
---|
| 874 | bounds = [p.lb(spliton) (p.lb(spliton)+p.ub(spliton))/2 p.ub(spliton)]; |
---|
| 875 | end |
---|
| 876 | case 'bisect' |
---|
| 877 | bounds = [p.lb(spliton) (p.lb(spliton)+p.ub(spliton))/2 p.ub(spliton)]; |
---|
| 878 | otherwise |
---|
| 879 | bounds = [p.lb(spliton) (p.lb(spliton)+p.ub(spliton))/2 p.ub(spliton)]; |
---|
| 880 | end |
---|
| 881 | |
---|
| 882 | |
---|
| 883 | |
---|
| 884 | function [p,feasible,vol_reduction] = boxreduce(p,upper,lower,lpsolver,options); |
---|
| 885 | |
---|
| 886 | if options.bmibnb.lpreduce |
---|
| 887 | |
---|
| 888 | vol_start = prod(p.ub(p.branch_variables)-p.lb(p.branch_variables)); |
---|
| 889 | diag_before = sum(p.ub(p.branch_variables)-p.lb(p.branch_variables)); |
---|
| 890 | diag_before0 = diag_before; |
---|
| 891 | |
---|
| 892 | [pcut,feasible,lower] = lpbmitighten(p,lower,upper,lpsolver); |
---|
| 893 | diag_after = sum(pcut.ub(p.branch_variables)-pcut.lb(p.branch_variables)); |
---|
| 894 | iterations = 0; |
---|
| 895 | while (diag_after/(1e-18+diag_before) < 0.75 ) & feasible & iterations<4 |
---|
| 896 | [pcut,feasible,lower] = lpbmitighten(pcut,lower,upper,lpsolver); |
---|
| 897 | diag_before = diag_after; |
---|
| 898 | diag_after = sum(pcut.ub(p.branch_variables)-pcut.lb(p.branch_variables)); |
---|
| 899 | iterations = iterations + 1; |
---|
| 900 | end |
---|
| 901 | |
---|
| 902 | % Clean up... |
---|
| 903 | for i = 1:length(pcut.lb) |
---|
| 904 | if (pcut.lb(i)>pcut.ub(i)) & (pcut.lb-pcut.ub < 1e-3) |
---|
| 905 | pcut.lb(i)=pcut.ub(i); |
---|
| 906 | pcut = updatenonlinearbounds(pcut,i); |
---|
| 907 | end |
---|
| 908 | end |
---|
| 909 | p.lb = pcut.lb; |
---|
| 910 | p.ub = pcut.ub; |
---|
| 911 | |
---|
| 912 | % Metric = (V0/V)^(1/n) |
---|
| 913 | vol_reduction = max(0,min(1,(prod(p.ub(p.branch_variables)-p.lb(p.branch_variables))/(1e-12+vol_start))^(1/length(p.branch_variables)))); |
---|
| 914 | p.lb(p.lb<-1e12) = -inf; |
---|
| 915 | p.ub(p.ub>1e12) = inf; |
---|
| 916 | else |
---|
| 917 | vol_reduction = 1; |
---|
| 918 | feasible = 1; |
---|
| 919 | end |
---|
| 920 | |
---|
| 921 | function output = solvelower(p,options,lowersolver) |
---|
| 922 | |
---|
| 923 | % ******************************************** |
---|
| 924 | % Convex envelope |
---|
| 925 | % ******************************************** |
---|
| 926 | %p.binary_variables = []; |
---|
| 927 | p_with_bilinear_cuts = p; |
---|
| 928 | p_with_bilinear_cuts.F_struc(1:p.K.f,:)=[]; |
---|
| 929 | p_with_bilinear_cuts = addmcgormick(p_with_bilinear_cuts); |
---|
| 930 | p_with_bilinear_cuts.F_struc = [p.F_struc(1:p.K.f,:);p_with_bilinear_cuts.F_struc]; |
---|
| 931 | |
---|
| 932 | % ************************************** |
---|
| 933 | % SOLVE NODE PROBLEM |
---|
| 934 | % ************************************** |
---|
| 935 | if any(p_with_bilinear_cuts.ub<p_with_bilinear_cuts.lb) |
---|
| 936 | output.problem=1; |
---|
| 937 | else |
---|
| 938 | % We are solving relaxed problem (penbmi might be local solver) |
---|
| 939 | p_with_bilinear_cuts.monomtable = eye(length(p_with_bilinear_cuts.c)); |
---|
| 940 | |
---|
| 941 | % fix implied from mccormick |
---|
| 942 | % p_with_bilinear_cuts.lb(p.linears) = -inf; |
---|
| 943 | % p_with_bilinear_cuts.ub(p.linears) = inf; |
---|
| 944 | % p_with_bilinear_cuts.lb(p.nonlins(:,1)) = -inf; |
---|
| 945 | % p_with_bilinear_cuts.ub(p.nonlins(:,1)) = inf; |
---|
| 946 | |
---|
| 947 | output = feval(lowersolver,p_with_bilinear_cuts); |
---|
| 948 | |
---|
| 949 | relaxed_residual = resids(p_with_bilinear_cuts,output.Primal); |
---|
| 950 | % Minor clean-up |
---|
| 951 | output.Primal(output.Primal<p.lb) = p.lb(output.Primal<p.lb); |
---|
| 952 | output.Primal(output.Primal>p.ub) = p.ub(output.Primal>p.ub); |
---|
| 953 | end |
---|
| 954 | |
---|
| 955 | function [p,stack] = selectbranch(p,options,stack,x_min,upper); |
---|
| 956 | switch options.bmibnb.branchmethod |
---|
| 957 | case 'maxvol' |
---|
| 958 | [node,stack] = pull(stack,'maxvol',x_min,upper); |
---|
| 959 | case 'best' |
---|
| 960 | [node,stack] = pull(stack,'best',x_min,upper); |
---|
| 961 | otherwise |
---|
| 962 | [node,stack] = pull(stack,'best',x_min,upper); |
---|
| 963 | end |
---|
| 964 | % Copy node data to p |
---|
| 965 | if isempty(node) |
---|
| 966 | p = []; |
---|
| 967 | else |
---|
| 968 | p.depth = node.depth; |
---|
| 969 | p.dpos = node.dpos; |
---|
| 970 | p.lb = node.lb; |
---|
| 971 | p.ub = node.ub; |
---|
| 972 | p.lower = node.lower; |
---|
| 973 | p.lpcuts = node.lpcuts; |
---|
| 974 | p.x0 = node.x0; |
---|
| 975 | end |
---|
| 976 | |
---|
| 977 | |
---|
| 978 | |
---|
| 979 | function output = solveupper(p,p_original,x,options,uppersolver) |
---|
| 980 | |
---|
| 981 | p_upper = p_original; |
---|
| 982 | |
---|
| 983 | % Pick an initial point (this can be a bit tricky...) |
---|
| 984 | % Use relaxed point, shifted towards center of box |
---|
| 985 | if all(x<=p.ub) & all(x>=p.lb) |
---|
| 986 | p_upper.x0 = 0.1*x + 0.9*(p.lb+p.ub)/2; |
---|
| 987 | else |
---|
| 988 | p_upper.x0 = (p.lb+p.ub)/2; |
---|
| 989 | end |
---|
| 990 | % Shift towards interior for variables with unbounded lower or upper |
---|
| 991 | lbinfbounds = find(isinf(p.lb)); |
---|
| 992 | ubinfbounds = find(isinf(p.ub)); |
---|
| 993 | p_upper.x0(ubinfbounds) = x(ubinfbounds)+0.01; |
---|
| 994 | p_upper.x0(lbinfbounds) = x(lbinfbounds)-0.01; |
---|
| 995 | ublbinfbounds = find(isinf(p.lb) & isinf(p.ub)); |
---|
| 996 | p_upper.x0(ublbinfbounds) = x(ublbinfbounds); |
---|
| 997 | % ...expand the current node just slightly |
---|
| 998 | p_upper.lb = p.lb; |
---|
| 999 | p_upper.ub = p.ub; |
---|
| 1000 | p_upper.lb(~isinf(p_original.lb)) = 0.99*p.lb(~isinf(p_original.lb))+p_original.lb(~isinf(p_original.lb))*0.01; |
---|
| 1001 | p_upper.ub(~isinf(p_original.ub)) = 0.99*p.ub(~isinf(p_original.ub))+p_original.ub(~isinf(p_original.ub))*0.01; |
---|
| 1002 | p_upper.lb(isinf(p_original.lb)) = p_upper.lb(isinf(p_original.lb)) - 0.001; |
---|
| 1003 | p_upper.ub(isinf(p_original.ub)) = p_upper.ub(isinf(p_original.ub)) + 0.001; |
---|
| 1004 | p_upper.options.saveduals = 0; |
---|
| 1005 | |
---|
| 1006 | % Solve upper bounding problem |
---|
| 1007 | p_upper.options.usex0 = 1; |
---|
| 1008 | output = feval(uppersolver,p_upper); |
---|
| 1009 | % Project into the box (numerical issue) |
---|
| 1010 | output.Primal(output.Primal<p_upper.lb) = p_upper.lb(output.Primal<p_upper.lb); |
---|
| 1011 | output.Primal(output.Primal>p_upper.ub) = p_upper.ub(output.Primal>p_upper.ub); |
---|
| 1012 | |
---|
| 1013 | |
---|
| 1014 | % This one needs a lot of work |
---|
| 1015 | function p = nonlinear_constraint_propagation(p) |
---|
| 1016 | |
---|
| 1017 | for i = 1:size(p.nonlins,1) |
---|
| 1018 | x = p.nonlins(i,2); |
---|
| 1019 | y = p.nonlins(i,3); |
---|
| 1020 | z = p.nonlins(i,1); |
---|
| 1021 | |
---|
| 1022 | if y==x & p.ub(z)>0 |
---|
| 1023 | p.ub(x) = min(p.ub(x),sqrt(p.ub(z))); |
---|
| 1024 | p.lb(x) = max(p.lb(x),-sqrt(p.ub(z))); |
---|
| 1025 | end |
---|
| 1026 | |
---|
| 1027 | if p.lb(y)>0 & p.ub(z)>0 & p.ub(x)>0 |
---|
| 1028 | p.ub(x) = min(p.ub(x),p.ub(z)/p.lb(y)); |
---|
| 1029 | end |
---|
| 1030 | if p.lb(x)>0 & p.ub(z)>0 & p.ub(y)>0 |
---|
| 1031 | p.ub(y) = min(p.ub(y),p.ub(z)/p.lb(x)); |
---|
| 1032 | end |
---|
| 1033 | |
---|
| 1034 | if p.ub(y)>0 & p.lb(y)>0 & p.lb(z)>0 |
---|
| 1035 | p.lb(x) = max(p.lb(x),p.lb(z)/p.ub(y)); |
---|
| 1036 | end |
---|
| 1037 | if p.ub(x)>0 & p.lb(x)>0 & p.lb(z)>0 |
---|
| 1038 | p.lb(y) = max(p.lb(y),p.lb(z)/p.ub(x)); |
---|
| 1039 | end |
---|
| 1040 | end |
---|
| 1041 | |
---|
| 1042 | |
---|
| 1043 | function vars = decide_branch_variables(p) |
---|
| 1044 | |
---|
| 1045 | if p.options.bmibnb.lowrank==0 |
---|
| 1046 | nonlinear = find(~(sum(p.monomtable~=0,2)==1 & sum(p.monomtable,2)==1)); |
---|
| 1047 | vars = find(sum(abs(full(p.monomtable(nonlinear,:))),1)); |
---|
| 1048 | else |
---|
| 1049 | pool1 = p.nonlins(1,2); |
---|
| 1050 | pool2 = p.nonlins(1,3); |
---|
| 1051 | |
---|
| 1052 | for i = 2:size(p.nonlins,1) |
---|
| 1053 | v1 = p.nonlins(i,2); |
---|
| 1054 | v2 = p.nonlins(i,3); |
---|
| 1055 | if v1==v2 |
---|
| 1056 | % We are fucked |
---|
| 1057 | pool1 = [pool1 v1]; |
---|
| 1058 | pool2 = [pool2 v2]; |
---|
| 1059 | else |
---|
| 1060 | if ismember(v1,pool1) |
---|
| 1061 | pool2 = [pool2 v2]; |
---|
| 1062 | elseif ismember(v1,pool2) |
---|
| 1063 | pool1 = [pool1 v2]; |
---|
| 1064 | elseif ismember(v2,pool1) |
---|
| 1065 | pool2 = [pool2 v1]; |
---|
| 1066 | elseif ismember(v2,pool2) |
---|
| 1067 | pool1 = [pool1 v1]; |
---|
| 1068 | else |
---|
| 1069 | % No member yet |
---|
| 1070 | pool1 = [pool1 v1]; |
---|
| 1071 | pool2 = [pool2 v2]; |
---|
| 1072 | end |
---|
| 1073 | end |
---|
| 1074 | end |
---|
| 1075 | pool1 = unique(pool1); |
---|
| 1076 | pool2 = unique(pool2); |
---|
| 1077 | if isempty(intersect(pool1,pool2)) |
---|
| 1078 | if length(pool1)<=length(pool2) |
---|
| 1079 | vars = pool1; |
---|
| 1080 | else |
---|
| 1081 | vars = pool2; |
---|
| 1082 | end |
---|
| 1083 | else |
---|
| 1084 | nonlinear = find(~(sum(p.monomtable~=0,2)==1 & sum(p.monomtable,2)==1)); |
---|
| 1085 | vars = find(sum(abs(full(p.monomtable(nonlinear,:))),1)); |
---|
| 1086 | end |
---|
| 1087 | end |
---|
| 1088 | |
---|
| 1089 | |
---|
| 1090 | function x = evaluate_nonlinear(p,x); |
---|
| 1091 | x(p.nonlins(:,1)) = x(p.nonlins(:,2)).*x(p.nonlins(:,3)); |
---|
| 1092 | |
---|
| 1093 | function p = clean_bounds(p) |
---|
| 1094 | % Fix to improve numerica with integer bounds |
---|
| 1095 | %close = find(1e-6>abs(p.ub - round(p.ub))); |
---|
| 1096 | %p.ub(close) = round(p.ub(close)); |
---|
| 1097 | close = 1e-6>abs(p.ub - round(p.ub)); |
---|
| 1098 | p.ub(close) = round(p.ub(close)); |
---|
| 1099 | |
---|
| 1100 | close = 1e-6>abs(p.lb - round(p.lb)); |
---|
| 1101 | p.lb(close) = round(p.lb(close)); |
---|
| 1102 | |
---|
| 1103 | p.ub(p.binary_variables) = floor(p.ub(p.binary_variables) + 1e-2); |
---|
| 1104 | %p.lb(p.binary_variables) = ceil(p.lb(p.binary_variables) - 1e-2); |
---|
| 1105 | %p = updatenonlinearbounds(p); |
---|
| 1106 | |
---|
| 1107 | % Nothing coded to do non-linear propagation |
---|
| 1108 | %p = nonlinear_constraint_propagation(p); |
---|
| 1109 | p.lb(p.lb<-1e12) = -inf; |
---|
| 1110 | p.ub(p.ub>1e12) = inf; |
---|
| 1111 | |
---|
| 1112 | |
---|
| 1113 | |
---|
| 1114 | function node = savetonode(p,spliton,bounds1,bounds2,direction,x,cost); |
---|
| 1115 | |
---|
| 1116 | node.lb = p.lb; |
---|
| 1117 | node.ub = p.ub; |
---|
| 1118 | node.lb(spliton) = bounds1; |
---|
| 1119 | node.ub(spliton) = bounds2; |
---|
| 1120 | if direction == -1 |
---|
| 1121 | node.dpos = p.dpos-1/(2^sqrt(p.depth)); |
---|
| 1122 | else |
---|
| 1123 | node.dpos = p.dpos+1/(2^sqrt(p.depth)); |
---|
| 1124 | end |
---|
| 1125 | node.depth = p.depth+1; |
---|
| 1126 | node.x0 = x; |
---|
| 1127 | node.lpcuts = p.lpcuts; |
---|
| 1128 | node.lower = cost; |
---|
| 1129 | |
---|