[37] | 1 | function output = bnb(p) |
---|
| 2 | %BNB General branch-and-bound scheme for conic programs |
---|
| 3 | % |
---|
| 4 | % BNB applies a branch-and-bound scheme to solve mixed integer |
---|
| 5 | % conic programs (LP, QP, SOCP, SDP) and mixed integer geometric programs. |
---|
| 6 | % |
---|
| 7 | % BNB is never called by the user directly, but is called by |
---|
| 8 | % YALMIP from SOLVESDP, by choosing the solver tag 'bnb' in sdpsettings. |
---|
| 9 | % |
---|
| 10 | % BNB is used if no other mixed integer solver is found, and |
---|
| 11 | % is only useful for very small problems, due to its simple |
---|
| 12 | % and naive implementation. |
---|
| 13 | % |
---|
| 14 | % The behaviour of BNB can be altered using the fields |
---|
| 15 | % in the field 'bnb' in SDPSETTINGS |
---|
| 16 | % |
---|
| 17 | % bnb.branchrule Deceides on what variable to branch |
---|
| 18 | % 'max' : Variable furthest away from being integer |
---|
| 19 | % 'min' : Variable closest to be being integer |
---|
| 20 | % 'first' : First variable (lowest variable index in YALMIP) |
---|
| 21 | % 'last' : Last variable (highest variable index in YALMIP) |
---|
| 22 | % 'weight' : See manual |
---|
| 23 | % |
---|
| 24 | % bnb.method Branching strategy |
---|
| 25 | % 'depth' : Depth first |
---|
| 26 | % 'breadth' : Breadth first |
---|
| 27 | % 'best' : Expand branch with lowest lower bound |
---|
| 28 | % 'depthX' : Depth until integer solution found, then X (e.g 'depthbest') |
---|
| 29 | % |
---|
| 30 | % solver Solver for the relaxed problems (standard solver tag, see SDPSETTINGS) |
---|
| 31 | % |
---|
| 32 | % inttol Tolerance for declaring a variable as integer |
---|
| 33 | % |
---|
| 34 | % feastol Tolerance for declaring constraints as feasible |
---|
| 35 | % |
---|
| 36 | % gaptol Exit when (upper bound-lower bound)/(1e-3+abs(lower bound)) < gaptol |
---|
| 37 | % |
---|
| 38 | % round Round variables smaller than bnb.inttol |
---|
| 39 | % |
---|
| 40 | % |
---|
| 41 | % See also SOLVESDP, BINVAR, INTVAR, BINARY, INTEGER |
---|
| 42 | |
---|
| 43 | % Author Johan Löfberg |
---|
| 44 | % $Id: bnb.m,v 1.12 2006/08/28 13:48:38 joloef Exp $ |
---|
| 45 | |
---|
| 46 | % ******************************** |
---|
| 47 | %% INITIALIZE DIAGNOSTICS IN YALMIP |
---|
| 48 | % ******************************** |
---|
| 49 | bnbsolvertime = clock; |
---|
| 50 | showprogress('Branch and bound started',p.options.showprogress); |
---|
| 51 | %p.options.verbose = 1; |
---|
| 52 | % ******************************** |
---|
| 53 | %% We might have a GP : pre-calc |
---|
| 54 | % ******************************** |
---|
| 55 | p.nonlinear = find(~(sum(p.monomtable~=0,2)==1 & sum(p.monomtable,2)==1)); |
---|
| 56 | p.nonlinear = union(p.nonlinear,p.evalVariables); |
---|
| 57 | |
---|
| 58 | % ******************************** |
---|
| 59 | %% Define infinite bounds |
---|
| 60 | % ******************************** |
---|
| 61 | if isempty(p.ub) |
---|
| 62 | p.ub = repmat(inf,length(p.c),1); |
---|
| 63 | end |
---|
| 64 | if isempty(p.lb) |
---|
| 65 | p.lb = repmat(-inf,length(p.c),1); |
---|
| 66 | end |
---|
| 67 | |
---|
| 68 | % ******************************** |
---|
| 69 | %% Extract bounds from model |
---|
| 70 | % ******************************** |
---|
| 71 | if ~isempty(p.F_struc) |
---|
| 72 | [lb,ub,used_rows] = findulb(p.F_struc,p.K); |
---|
| 73 | if ~isempty(used_rows) |
---|
| 74 | used_rows = used_rows(~any(full(p.F_struc(used_rows,1+p.nonlinear)),2)); |
---|
| 75 | if ~isempty(used_rows) |
---|
| 76 | lower_defined = find(~isinf(lb)); |
---|
| 77 | if ~isempty(lower_defined) |
---|
| 78 | p.lb(lower_defined) = max(p.lb(lower_defined),lb(lower_defined)); |
---|
| 79 | end |
---|
| 80 | upper_defined = find(~isinf(ub)); |
---|
| 81 | if ~isempty(upper_defined) |
---|
| 82 | p.ub(upper_defined) = min(p.ub(upper_defined),ub(upper_defined)); |
---|
| 83 | end |
---|
| 84 | p.F_struc(p.K.f+used_rows,:)=[]; |
---|
| 85 | p.K.l = p.K.l - length(used_rows); |
---|
| 86 | end |
---|
| 87 | end |
---|
| 88 | end |
---|
| 89 | |
---|
| 90 | % ******************************** |
---|
| 91 | %% ADD CONSTRAINTS 0<x<1 FOR BINARY |
---|
| 92 | % ******************************** |
---|
| 93 | if ~isempty(p.binary_variables) |
---|
| 94 | p.ub(p.binary_variables) = min(p.ub(p.binary_variables),1); |
---|
| 95 | p.lb(p.binary_variables) = max(p.lb(p.binary_variables),0); |
---|
| 96 | |
---|
| 97 | godown = find(p.ub(p.binary_variables) < 1);%-p.options.bnb.inttol); |
---|
| 98 | goup = find(p.lb(p.binary_variables) > 0);%p.options.bnb.inttol); |
---|
| 99 | p.ub(p.binary_variables(godown)) = 0; |
---|
| 100 | p.lb(p.binary_variables(goup)) = 1; |
---|
| 101 | end |
---|
| 102 | |
---|
| 103 | p.lb(p.integer_variables) = fix(p.lb(p.integer_variables)); |
---|
| 104 | p.ub(p.integer_variables) = fix(p.ub(p.integer_variables)); |
---|
| 105 | |
---|
| 106 | % ******************************* |
---|
| 107 | %% PRE-SOLVE (nothing fancy coded) |
---|
| 108 | % ******************************* |
---|
| 109 | pss=[]; |
---|
| 110 | if isempty(p.nonlinear) |
---|
| 111 | %if isempty(find(isinf([p.ub;p.lb]))) & p.K.l>0 & isempty(p.nonlinear) |
---|
| 112 | if p.K.f>0 |
---|
| 113 | Aeq = -p.F_struc(1:p.K.f,2:end); |
---|
| 114 | beq = p.F_struc(1:p.K.f,1); |
---|
| 115 | A = [Aeq;-Aeq]; |
---|
| 116 | b = [beq;-beq]; |
---|
| 117 | [p.lb,p.ub,redundant,pss] = tightenbounds(A,b,p.lb,p.ub,p.integer_variables,p.binary_variables,ones(length(p.lb),1)); |
---|
| 118 | end |
---|
| 119 | pss=[]; |
---|
| 120 | if p.K.l>0 |
---|
| 121 | A = -p.F_struc(1+p.K.f:p.K.f+p.K.l,2:end); |
---|
| 122 | b = p.F_struc(1+p.K.f:p.K.f+p.K.l,1); |
---|
| 123 | [p.lb,p.ub,redundant,pss] = tightenbounds(A,b,p.lb,p.ub,p.integer_variables,p.binary_variables,ones(length(p.lb),1)); |
---|
| 124 | if length(redundant)>0 |
---|
| 125 | pss.AL0A(redundant,:)=[]; |
---|
| 126 | pss.AG0A(redundant,:)=[]; |
---|
| 127 | p.F_struc(p.K.f+redundant,:)=[]; |
---|
| 128 | p.K.l = p.K.l - length(redundant); |
---|
| 129 | end |
---|
| 130 | end |
---|
| 131 | end |
---|
| 132 | |
---|
| 133 | % ******************************* |
---|
| 134 | %% PERTURBATION OF LINEAR COST |
---|
| 135 | % ******************************* |
---|
| 136 | p.corig = p.c; |
---|
| 137 | if nnz(p.Q)~=0 |
---|
| 138 | g = randn('seed'); |
---|
| 139 | randn('state',1253); %For my testing, I keep this the same... |
---|
| 140 | % This perturbation has to be better. Crucial for many real LP problems |
---|
| 141 | p.c = (p.c).*(1+randn(length(p.c),1)*1e-4); |
---|
| 142 | randn('seed',g); |
---|
| 143 | end |
---|
| 144 | |
---|
| 145 | % ******************************* |
---|
| 146 | %% Display logics |
---|
| 147 | % 0 : Silent |
---|
| 148 | % 1 : Display branching |
---|
| 149 | % 2 : Display node solver prints |
---|
| 150 | % ******************************* |
---|
| 151 | switch max(min(p.options.verbose,3),0) |
---|
| 152 | case 0 |
---|
| 153 | p.options.bnb.verbose = 0; |
---|
| 154 | case 1 |
---|
| 155 | p.options.bnb.verbose = 1; |
---|
| 156 | p.options.verbose = 0; |
---|
| 157 | case 2 |
---|
| 158 | p.options.bnb.verbose = 2; |
---|
| 159 | p.options.verbose = 0; |
---|
| 160 | case 3 |
---|
| 161 | p.options.bnb.verbose = 2; |
---|
| 162 | p.options.verbose = 1; |
---|
| 163 | otherwise |
---|
| 164 | p.options.bnb.verbose = 0; |
---|
| 165 | p.options.verbose = 0; |
---|
| 166 | end |
---|
| 167 | |
---|
| 168 | % ******************************* |
---|
| 169 | %% Figure out the weights if any |
---|
| 170 | % ******************************* |
---|
| 171 | try % Probably buggy first version... |
---|
| 172 | if ~isempty(p.options.bnb.weight) |
---|
| 173 | weightvar = p.options.bnb.weight; |
---|
| 174 | if isa(weightvar,'sdpvar') |
---|
| 175 | if (prod(size(weightvar)) == 1) |
---|
| 176 | weight = ones(length(p.c),1); |
---|
| 177 | for i = 1:length(p.c) |
---|
| 178 | weight(i,1) = full(getbasematrix(weightvar,p.used_variables(i))); |
---|
| 179 | end |
---|
| 180 | p.weight = weight; |
---|
| 181 | else |
---|
| 182 | error('Weight should be an SDPVAR scalar'); |
---|
| 183 | end |
---|
| 184 | else |
---|
| 185 | error('Weight should be an SDPVAR scalar'); |
---|
| 186 | end |
---|
| 187 | else |
---|
| 188 | p.weight = ones(length(p.c),1); |
---|
| 189 | p.weight(p.binary_variables) = (1./(1:length(p.binary_variables))); |
---|
| 190 | end |
---|
| 191 | catch |
---|
| 192 | disp('Something wrong with weights. Please report bug'); |
---|
| 193 | p.weight = ones(length(p.c),1); |
---|
| 194 | end |
---|
| 195 | |
---|
| 196 | % ******************************* |
---|
| 197 | %% START BRANCHING |
---|
| 198 | % ******************************* |
---|
| 199 | [x_min,solved_nodes,lower,upper] = branch_and_bound(p,pss); |
---|
| 200 | |
---|
| 201 | % ********************************** |
---|
| 202 | %% CREATE SOLUTION |
---|
| 203 | % ********************************** |
---|
| 204 | output.problem = 0; |
---|
| 205 | if isinf(upper) |
---|
| 206 | output.problem = 1; |
---|
| 207 | end |
---|
| 208 | if isinf(-lower) |
---|
| 209 | output.problem = 2; |
---|
| 210 | end |
---|
| 211 | if solved_nodes == p.options.bnb.maxiter |
---|
| 212 | output.problem = 3; |
---|
| 213 | end |
---|
| 214 | output.solved_nodes = solved_nodes; |
---|
| 215 | output.Primal = x_min; |
---|
| 216 | output.Dual = []; |
---|
| 217 | output.Slack = []; |
---|
| 218 | output.infostr = yalmiperror(output.problem,'BNB'); |
---|
| 219 | output.solverinput = 0; |
---|
| 220 | if p.options.savesolveroutput |
---|
| 221 | output.solveroutput.solved_nodes = solved_nodes; |
---|
| 222 | output.solveroutput.lower = lower; |
---|
| 223 | output.solveroutput.upper = upper; |
---|
| 224 | else |
---|
| 225 | output.solveroutput =[]; |
---|
| 226 | end |
---|
| 227 | output.solvertime = etime(clock,bnbsolvertime); |
---|
| 228 | %% -- |
---|
| 229 | |
---|
| 230 | function [x_min,solved_nodes,lower,upper] = branch_and_bound(p,pss) |
---|
| 231 | |
---|
| 232 | % ******************************* |
---|
| 233 | %% We don't need this |
---|
| 234 | % ******************************* |
---|
| 235 | p.options.savesolverinput = 0; |
---|
| 236 | p.options.savesolveroutput = 0; |
---|
| 237 | p.options.saveduals = 0; |
---|
| 238 | p.options.dimacs = 0; |
---|
| 239 | |
---|
| 240 | % ******************************* |
---|
| 241 | %% SET-UP ROOT PROBLEM |
---|
| 242 | % ******************************* |
---|
| 243 | p.depth = 0; |
---|
| 244 | p.lower = NaN; |
---|
| 245 | if p.options.usex0 |
---|
| 246 | [x_min,upper] = initializesolution(p); |
---|
| 247 | p.x0 = x_min; |
---|
| 248 | else |
---|
| 249 | upper = inf; |
---|
| 250 | x_min = zeros(length(p.c),1); |
---|
| 251 | p.x0 = zeros(length(p.c),1); |
---|
| 252 | end |
---|
| 253 | |
---|
| 254 | %if isinf(upper) & isequal(p.solver.lower.tag) |
---|
| 255 | %end |
---|
| 256 | |
---|
| 257 | % ******************************* |
---|
| 258 | %% Global stuff |
---|
| 259 | % ******************************* |
---|
| 260 | lower = NaN; |
---|
| 261 | stack = []; |
---|
| 262 | |
---|
| 263 | % ******************************* |
---|
| 264 | %% Create function handle to solver |
---|
| 265 | % ******************************* |
---|
| 266 | lowersolver = p.solver.lower.call; |
---|
| 267 | uppersolver = p.options.bnb.uppersolver; |
---|
| 268 | |
---|
| 269 | % ******************************* |
---|
| 270 | %% INVARIANT PROBLEM DATA |
---|
| 271 | % ******************************* |
---|
| 272 | c = p.corig; |
---|
| 273 | Q = p.Q; |
---|
| 274 | f = p.f; |
---|
| 275 | integer_variables = p.integer_variables; |
---|
| 276 | solved_nodes = 0; |
---|
| 277 | |
---|
| 278 | gap = inf; |
---|
| 279 | node = 1; |
---|
| 280 | |
---|
| 281 | if p.options.bnb.presolve |
---|
| 282 | savec = p.c; |
---|
| 283 | saveQ = p.Q; |
---|
| 284 | p.Q = p.Q*0; |
---|
| 285 | |
---|
| 286 | n = length(p.c); |
---|
| 287 | saveBinary = p.binary_variables; |
---|
| 288 | saveInteger = p.integer_variables; |
---|
| 289 | p.binary_variables = []; |
---|
| 290 | p.integer_variables = [];; |
---|
| 291 | |
---|
| 292 | for i = 1:length(c) |
---|
| 293 | p.c = eyev(n,i); |
---|
| 294 | output = feval(lowersolver,p); |
---|
| 295 | if output.problem == 0 |
---|
| 296 | p.lb(i) = max(p.lb(i),output.Primal(i)); |
---|
| 297 | end |
---|
| 298 | p.c = -eyev(n,i); |
---|
| 299 | output = feval(lowersolver,p); |
---|
| 300 | if output.problem == 0 |
---|
| 301 | p.ub(i) = min(p.ub(i),output.Primal(i)); |
---|
| 302 | end |
---|
| 303 | p.lb(saveBinary) = ceil(p.lb(saveBinary)-1e-3); |
---|
| 304 | p.ub(saveBinary) = floor(p.ub(saveBinary)+1e-3); |
---|
| 305 | end |
---|
| 306 | p.binary_variables = saveBinary; |
---|
| 307 | p.integer_variables = saveInteger; |
---|
| 308 | |
---|
| 309 | p.Q = saveQ; |
---|
| 310 | p.c = savec; |
---|
| 311 | end |
---|
| 312 | |
---|
| 313 | % ************************************************ |
---|
| 314 | % Some hacks to speed up solver calls |
---|
| 315 | % ************************************************ |
---|
| 316 | p.getsolvertime = 0; |
---|
| 317 | |
---|
| 318 | % ******************************* |
---|
| 319 | %% DISPLAY HEADER |
---|
| 320 | % ******************************* |
---|
| 321 | originalDiscrete = [p.integer_variables(:);p.binary_variables(:)]; |
---|
| 322 | originalBinary = p.binary_variables(:); |
---|
| 323 | |
---|
| 324 | if nnz(Q)==0 & (nnz(p.c-fix(p.c))==0) |
---|
| 325 | can_use_ceil_lower = all(ismember(find(p.c),originalDiscrete)); |
---|
| 326 | else |
---|
| 327 | can_use_ceil_lower = 0; |
---|
| 328 | end |
---|
| 329 | |
---|
| 330 | if p.options.bnb.verbose |
---|
| 331 | |
---|
| 332 | pc = p.problemclass; |
---|
| 333 | non_convex_obj = pc.objective.quadratic.nonconvex | pc.objective.polynomial; |
---|
| 334 | |
---|
| 335 | possiblynonconvex = non_convex_obj; |
---|
| 336 | if ~isequal(p.solver.lower.version,'') |
---|
| 337 | p.solver.lower.tag = [p.solver.lower.tag '-' p.solver.lower.version]; |
---|
| 338 | end |
---|
| 339 | |
---|
| 340 | disp('* Starting YALMIP integer branch & bound.'); |
---|
| 341 | disp(['* Lower solver : ' p.solver.lower.tag]); |
---|
| 342 | disp(['* Upper solver : ' p.options.bnb.uppersolver]); |
---|
| 343 | disp(['* Max iterations : ' num2str(p.options.bnb.maxiter)]); |
---|
| 344 | |
---|
| 345 | if possiblynonconvex & p.options.warning |
---|
| 346 | disp(' '); |
---|
| 347 | disp('Warning : The relaxed problem may be nonconvex. This means '); |
---|
| 348 | disp('that the branching process not is guaranteed to find a'); |
---|
| 349 | disp('globally optimal solution, since the lower bound can be'); |
---|
| 350 | disp('invalid. Hence, do not trust the bound or the gap...') |
---|
| 351 | end |
---|
| 352 | end |
---|
| 353 | if p.options.bnb.verbose; disp(' Node Upper Gap(%) Lower Open');end; |
---|
| 354 | |
---|
| 355 | if nnz(Q)==0 & nnz(c)==1 |
---|
| 356 | p.simplecost = 1; |
---|
| 357 | else |
---|
| 358 | p.simplecost = 0; |
---|
| 359 | end |
---|
| 360 | |
---|
| 361 | poriginal = p; |
---|
| 362 | p.cuts = []; |
---|
| 363 | |
---|
| 364 | %% MAIN LOOP |
---|
| 365 | p.options.rounding = [1 1 1 1]; |
---|
| 366 | |
---|
| 367 | % |
---|
| 368 | if p.options.bnb.nodefix & (p.K.s(1)>0) |
---|
| 369 | top=1+p.K.f+p.K.l+sum(p.K.q); |
---|
| 370 | for i=1:length(p.K.s) |
---|
| 371 | n=p.K.s(i); |
---|
| 372 | for j=1:size(p.F_struc,2)-1; |
---|
| 373 | X=full(reshape(p.F_struc(top:top+n^2-1,j+1),p.K.s(i),p.K.s(i))); |
---|
| 374 | X=(X+X')/2; |
---|
| 375 | v=real(eig(X+sqrt(eps)*eye(length(X)))); |
---|
| 376 | if all(v>=0) |
---|
| 377 | sdpmonotinicity(i,j)=-1; |
---|
| 378 | elseif all(v<=0) |
---|
| 379 | sdpmonotinicity(i,j)=1; |
---|
| 380 | else |
---|
| 381 | sdpmonotinicity(i,j)=nan; |
---|
| 382 | end |
---|
| 383 | end |
---|
| 384 | top=top+n^2; |
---|
| 385 | end |
---|
| 386 | else |
---|
| 387 | sdpmonotinicity=[]; |
---|
| 388 | end |
---|
| 389 | |
---|
| 390 | % Try to find sum(d_i) = 1 |
---|
| 391 | sosgroups = {}; |
---|
| 392 | sosvariables = []; |
---|
| 393 | if p.K.f > 0 & ~isempty(p.binary_variables) |
---|
| 394 | nbin = length(p.binary_variables); |
---|
| 395 | Aeq = -p.F_struc(1:p.K.f,2:end); |
---|
| 396 | beq = p.F_struc(1:p.K.f,1); |
---|
| 397 | notbinary_var_index = setdiff(1:length(p.lb),p.binary_variables); |
---|
| 398 | only_binary = ~any(Aeq(:,notbinary_var_index),2); |
---|
| 399 | Aeq_bin = Aeq(find(only_binary),p.binary_variables); |
---|
| 400 | beq_bin = beq(find(only_binary),:); |
---|
| 401 | % Detect groups with constraints sum(d_i) == 1 |
---|
| 402 | sosgroups = {}; |
---|
| 403 | for i = 1:size(Aeq_bin,1) |
---|
| 404 | if beq_bin(i) == 1 |
---|
| 405 | [ix,jx,sx] = find(Aeq_bin(i,:)); |
---|
| 406 | if all(sx == 1) |
---|
| 407 | sosgroups{end+1} = p.binary_variables(jx); |
---|
| 408 | sosvariables = [sosvariables p.binary_variables(jx)]; |
---|
| 409 | end |
---|
| 410 | end |
---|
| 411 | end |
---|
| 412 | end |
---|
| 413 | |
---|
| 414 | while ~isempty(node) & (solved_nodes < p.options.bnb.maxiter) & (isinf(lower) | gap>p.options.bnb.gaptol) |
---|
| 415 | |
---|
| 416 | % ******************************************** |
---|
| 417 | % Adjust variable bound based on upper bound |
---|
| 418 | % ******************************************** |
---|
| 419 | % p = Updatecostbound(p,upper); |
---|
| 420 | |
---|
| 421 | % This code typically never runs but can be turned on |
---|
| 422 | % using options.bnb.nodetight and bnb.nodefix. |
---|
| 423 | |
---|
| 424 | if ~isinf(upper) & ~isnan(lower) |
---|
| 425 | [p,poriginal,stack] = pruneglobally(p,poriginal,upper,lower,stack,x); |
---|
| 426 | [p,poriginal,stack] = fixvariables(p,poriginal,upper,lower,stack,x_min,sdpmonotinicity); |
---|
| 427 | end |
---|
| 428 | |
---|
| 429 | % ******************************************** |
---|
| 430 | % BINARY VARIABLES ARE FIXED ALONG THE PROCESS |
---|
| 431 | % ******************************************** |
---|
| 432 | binary_variables = p.binary_variables; |
---|
| 433 | |
---|
| 434 | % ******************************************** |
---|
| 435 | % ASSUME THAT WE WON'T FATHOME |
---|
| 436 | % ******************************************** |
---|
| 437 | keep_digging = 1; |
---|
| 438 | message = ''; |
---|
| 439 | |
---|
| 440 | % ************************************* |
---|
| 441 | % SOLVE NODE PROBLEM |
---|
| 442 | % ************************************* |
---|
| 443 | if any(p.ub<p.lb - 1e-12) |
---|
| 444 | x = zeros(length(p.c),1); |
---|
| 445 | output.Primal = x; |
---|
| 446 | output.problem=1; |
---|
| 447 | else |
---|
| 448 | |
---|
| 449 | relaxed_p = p; |
---|
| 450 | relaxed_p.integer_variables = []; |
---|
| 451 | relaxed_p.binary_variables = []; |
---|
| 452 | relaxed_p.ub(p.ub<p.lb) = relaxed_p.lb(p.ub<p.lb); |
---|
| 453 | output = solvelower(lowersolver,relaxed_p,upper+abs(upper)*1e-2+1e-4,lower); |
---|
| 454 | x = setnonlinearvariables(p,output.Primal); |
---|
| 455 | |
---|
| 456 | % ************************************** |
---|
| 457 | % Hmm, don't remember why this fix... |
---|
| 458 | % ************************************** |
---|
| 459 | if (p.K.l>0) & any(p.F_struc(p.K.f+1:p.K.f+p.K.l,:)*[1;x]<-1e-5) |
---|
| 460 | output.problem = 1; |
---|
| 461 | end |
---|
| 462 | end |
---|
| 463 | |
---|
| 464 | solved_nodes = solved_nodes+1; |
---|
| 465 | |
---|
| 466 | % ************************************** |
---|
| 467 | % THIS WILL BE INTIAL GUESS FOR CHILDREN |
---|
| 468 | % ************************************** |
---|
| 469 | p.x0 = x; |
---|
| 470 | |
---|
| 471 | % ************************************* |
---|
| 472 | % NODE HEURISTICS (NOTHING CODED) |
---|
| 473 | % ************************************* |
---|
| 474 | if output.problem==0 | output.problem==3 | output.problem==4 |
---|
| 475 | cost = f+c'*x+x'*Q*x ; |
---|
| 476 | if isnan(lower) |
---|
| 477 | lower = cost; |
---|
| 478 | end |
---|
| 479 | |
---|
| 480 | [upper1,x_min1] = feval(uppersolver,poriginal,output); |
---|
| 481 | |
---|
| 482 | if upper1 < upper |
---|
| 483 | x_min = x_min1; |
---|
| 484 | upper = upper1; |
---|
| 485 | [stack,stacklower] = prune(stack,upper,p.options,solved_nodes,p); |
---|
| 486 | lower = min(lower,stacklower); |
---|
| 487 | [p,poriginal,stack] = pruneglobally(p,poriginal,upper,lower,stack,x_min); |
---|
| 488 | [p,poriginal,stack] = fixvariables(p,poriginal,upper,lower,stack,x_min,sdpmonotinicity); |
---|
| 489 | end |
---|
| 490 | end |
---|
| 491 | p = adaptivestrategy(p,upper,solved_nodes); |
---|
| 492 | |
---|
| 493 | % ************************************* |
---|
| 494 | % ANY INTEGERS? ROUND? |
---|
| 495 | % ************************************* |
---|
| 496 | non_integer_binary = abs(x(binary_variables)-round(x(binary_variables)))>p.options.bnb.inttol; |
---|
| 497 | non_integer_integer = abs(x(integer_variables)-round(x(integer_variables)))>p.options.bnb.inttol; |
---|
| 498 | if p.options.bnb.round |
---|
| 499 | x(binary_variables(~non_integer_binary)) = round(x(binary_variables(~non_integer_binary))); |
---|
| 500 | x(integer_variables(~non_integer_integer)) = round(x(integer_variables(~non_integer_integer))); |
---|
| 501 | end |
---|
| 502 | non_integer_binary = find(non_integer_binary); |
---|
| 503 | non_integer_integer = find(non_integer_integer); |
---|
| 504 | |
---|
| 505 | % ************************************* |
---|
| 506 | % CHECK FATHOMING POSSIBILITIES |
---|
| 507 | % ************************************* |
---|
| 508 | feasible = 1; |
---|
| 509 | switch output.problem |
---|
| 510 | case 0 |
---|
| 511 | if can_use_ceil_lower |
---|
| 512 | lower = ceil(lower); |
---|
| 513 | end |
---|
| 514 | case {1,12} |
---|
| 515 | keep_digging = 0; |
---|
| 516 | cost = inf; |
---|
| 517 | feasible = 0; |
---|
| 518 | case 2 |
---|
| 519 | cost = -inf; |
---|
| 520 | otherwise |
---|
| 521 | % This part has to be much more robust |
---|
| 522 | cost = f+c'*x+x'*Q*x; |
---|
| 523 | end |
---|
| 524 | |
---|
| 525 | % ************************************** |
---|
| 526 | % YAHOO! INTEGER SOLUTION FOUND |
---|
| 527 | % ************************************** |
---|
| 528 | if isempty(non_integer_binary) & isempty(non_integer_integer) |
---|
| 529 | if (cost<upper) & feasible |
---|
| 530 | x_min = x; |
---|
| 531 | upper = cost; |
---|
| 532 | [stack,lower] = prune(stack,upper,p.options,solved_nodes,p); |
---|
| 533 | end |
---|
| 534 | p = adaptivestrategy(p,upper,solved_nodes); |
---|
| 535 | keep_digging = 0; |
---|
| 536 | end |
---|
| 537 | if cost>upper*(1-1e-6) |
---|
| 538 | keep_digging = 0; |
---|
| 539 | end |
---|
| 540 | |
---|
| 541 | % ********************************** |
---|
| 542 | % CONTINUE SPLITTING? |
---|
| 543 | % ********************************** |
---|
| 544 | if keep_digging & (cost<upper) |
---|
| 545 | |
---|
| 546 | % ********************************** |
---|
| 547 | % BRANCH VARIABLE |
---|
| 548 | % ********************************** |
---|
| 549 | [index,whatsplit] = branchvariable(x,integer_variables,binary_variables,p.options,x_min,[],p); |
---|
| 550 | |
---|
| 551 | % ********************************** |
---|
| 552 | % CREATE NEW PROBLEMS |
---|
| 553 | % ********************************** |
---|
| 554 | p0_feasible = 1; |
---|
| 555 | p1_feasible = 1; |
---|
| 556 | switch whatsplit |
---|
| 557 | case 'binary' |
---|
| 558 | |
---|
| 559 | [p0,p1,index] = binarysplit(p,x,index,cost,[],sosgroups,sosvariables); |
---|
| 560 | |
---|
| 561 | case 'integer' |
---|
| 562 | [p0,p1] = integersplit(p,x,index,cost,x_min); |
---|
| 563 | otherwise |
---|
| 564 | end |
---|
| 565 | |
---|
| 566 | % ********************************** |
---|
| 567 | % Only save varying data in the tree |
---|
| 568 | % ********************************** |
---|
| 569 | node1.lb = p1.lb; |
---|
| 570 | node1.ub = p1.ub; |
---|
| 571 | node1.depth = p1.depth; |
---|
| 572 | node1.lower = p1.lower; |
---|
| 573 | node1.x0 = p1.x0; |
---|
| 574 | node1.binary_variables = p1.binary_variables; |
---|
| 575 | |
---|
| 576 | node0.lb = p0.lb; |
---|
| 577 | node0.ub = p0.ub; |
---|
| 578 | node0.depth = p0.depth; |
---|
| 579 | node0.lower = p0.lower; |
---|
| 580 | node0.x0 = p0.x0; |
---|
| 581 | node0.binary_variables = p0.binary_variables; |
---|
| 582 | |
---|
| 583 | if p1_feasible |
---|
| 584 | stack = push(stack,node1); |
---|
| 585 | end |
---|
| 586 | if p0_feasible |
---|
| 587 | stack = push(stack,node0); |
---|
| 588 | end |
---|
| 589 | end |
---|
| 590 | |
---|
| 591 | % Lowest cost in any open node |
---|
| 592 | if ~isempty(stack) |
---|
| 593 | lower = min([stack.lower]); |
---|
| 594 | if can_use_ceil_lower |
---|
| 595 | lower = ceil(lower); |
---|
| 596 | end |
---|
| 597 | end |
---|
| 598 | |
---|
| 599 | % ********************************** |
---|
| 600 | % Get a new node to solve |
---|
| 601 | % ********************************** |
---|
| 602 | [node,stack] = pull(stack,p.options.bnb.method,x_min,upper); |
---|
| 603 | if ~isempty(node) |
---|
| 604 | p.lb = node.lb; |
---|
| 605 | p.ub = node.ub; |
---|
| 606 | p.depth = node.depth; |
---|
| 607 | p.lower = node.lower; |
---|
| 608 | p.x0 = node.x0; |
---|
| 609 | p.binary_variables = node.binary_variables; |
---|
| 610 | end |
---|
| 611 | gap = abs((upper-lower)/(1e-3+abs(upper)+abs(lower))); |
---|
| 612 | if isnan(gap) |
---|
| 613 | gap = inf; |
---|
| 614 | end |
---|
| 615 | |
---|
| 616 | %DEBUG if p.options.bnb.verbose;fprintf(' %4.0f : %12.3E %7.2f %12.3E %2.0f %2.0f %2.0f %2.0f %2.0f\n',solved_nodes,upper,100*gap,lower,length(stack)+length(node),sedd);end |
---|
| 617 | if p.options.bnb.verbose;fprintf(' %4.0f : %12.3E %7.2f %12.3E %2.0f \n',solved_nodes,upper,100*gap,lower,length(stack)+length(node));end |
---|
| 618 | |
---|
| 619 | end |
---|
| 620 | if p.options.bnb.verbose;showprogress([num2str2(solved_nodes,3) ' Finishing. Cost: ' num2str(upper) ],p.options.bnb.verbose);end |
---|
| 621 | |
---|
| 622 | |
---|
| 623 | function stack = push(stackin,p) |
---|
| 624 | if ~isempty(stackin) |
---|
| 625 | stack = [p;stackin]; |
---|
| 626 | else |
---|
| 627 | stack(1)=p; |
---|
| 628 | end |
---|
| 629 | |
---|
| 630 | %% |
---|
| 631 | function [p,stack] = pull(stack,method,x_min,upper); |
---|
| 632 | |
---|
| 633 | if ~isempty(stack) |
---|
| 634 | switch method |
---|
| 635 | case {'depth','depthfirst','depthbreadth','depthproject','depthbest'} |
---|
| 636 | [i,j]=max([stack.depth]); |
---|
| 637 | p=stack(j); |
---|
| 638 | stack = stack([1:1:j-1 j+1:1:end]); |
---|
| 639 | |
---|
| 640 | case 'breadth' |
---|
| 641 | [i,j]=min([stack.depth]); |
---|
| 642 | p=stack(j); |
---|
| 643 | stack = stack([1:1:j-1 j+1:1:end]); |
---|
| 644 | |
---|
| 645 | case 'best' |
---|
| 646 | [i,j]=min([stack.lower]); |
---|
| 647 | p=stack(j); |
---|
| 648 | stack = stack([1:1:j-1 j+1:1:end]); |
---|
| 649 | |
---|
| 650 | otherwise |
---|
| 651 | end |
---|
| 652 | else |
---|
| 653 | p = []; |
---|
| 654 | end |
---|
| 655 | |
---|
| 656 | % ********************************** |
---|
| 657 | %% BRANCH VARIABLE |
---|
| 658 | % ********************************** |
---|
| 659 | function [index,whatsplit] = branchvariable(x,integer_variables,binary_variables,options,x_min,Weight,p) |
---|
| 660 | all_variables = [integer_variables(:);binary_variables(:)]; |
---|
| 661 | |
---|
| 662 | switch options.bnb.branchrule |
---|
| 663 | case 'weight' |
---|
| 664 | interror = abs(x(all_variables)-round(x(all_variables))); |
---|
| 665 | [val,index] = max(abs(p.weight(all_variables)).*interror); |
---|
| 666 | case 'first' |
---|
| 667 | index = min(find(abs(x(all_variables)-round(x(all_variables)))>options.bnb.inttol)); |
---|
| 668 | case 'last' |
---|
| 669 | index = max(find(abs(x(all_variables)-round(x(all_variables)))>options.bnb.inttol)); |
---|
| 670 | case 'min' |
---|
| 671 | nint = find(abs(x(all_variables)-round(x(all_variables)))>options.bnb.inttol); |
---|
| 672 | [val,index] = min(abs(x(nint))); |
---|
| 673 | index = nint(index); |
---|
| 674 | case 'max' |
---|
| 675 | [val,index] = max(abs(x(all_variables)-round(x(all_variables)))); |
---|
| 676 | otherwise |
---|
| 677 | error |
---|
| 678 | end |
---|
| 679 | if index<=length(integer_variables) |
---|
| 680 | whatsplit = 'integer'; |
---|
| 681 | else |
---|
| 682 | index = index-length(integer_variables); |
---|
| 683 | whatsplit = 'binary'; |
---|
| 684 | end |
---|
| 685 | |
---|
| 686 | % ********************************** |
---|
| 687 | % SPLIT PROBLEM |
---|
| 688 | % ********************************** |
---|
| 689 | function [p0,p1,variable] = binarysplit(p,x,index,lower,options,sosgroups,sosvariables) |
---|
| 690 | p0 = p; |
---|
| 691 | p1 = p; |
---|
| 692 | |
---|
| 693 | variable = p.binary_variables(index); |
---|
| 694 | tf = ~(ismembc(p0.binary_variables,variable)); |
---|
| 695 | new_binary = p0.binary_variables(tf); |
---|
| 696 | |
---|
| 697 | % friends = []; |
---|
| 698 | % if ~isempty(sosvariables) |
---|
| 699 | % if ismember(variable,sosvariables) |
---|
| 700 | % i = 1; |
---|
| 701 | % while i<=length(sosgroups) |
---|
| 702 | % |
---|
| 703 | % if ismember(variable,sosgroups{i}) |
---|
| 704 | % friends = setdiff(sosgroups{i},variable); |
---|
| 705 | % break |
---|
| 706 | % else |
---|
| 707 | % i = i + 1; |
---|
| 708 | % end |
---|
| 709 | % end |
---|
| 710 | % end |
---|
| 711 | % end |
---|
| 712 | |
---|
| 713 | p0.ub(variable)=0; |
---|
| 714 | p0.lb(variable)=0; |
---|
| 715 | %if length(friends) == 1 |
---|
| 716 | % p0.ub(friends)=1; |
---|
| 717 | % p0.lb(friends)=1; |
---|
| 718 | %end |
---|
| 719 | |
---|
| 720 | p0.lower = lower; |
---|
| 721 | p0.depth = p.depth+1; |
---|
| 722 | p0.binary_variables = new_binary;%setdiff1D(p0.binary_variables,variable); |
---|
| 723 | %p0.binary_variables = setdiff(p0.binary_variables,friends); |
---|
| 724 | |
---|
| 725 | p1.ub(variable)=1; |
---|
| 726 | p1.lb(variable)=1; |
---|
| 727 | %p1.ub(friends)=0; |
---|
| 728 | %p1.lb(friends)=0; |
---|
| 729 | |
---|
| 730 | |
---|
| 731 | p1.binary_variables = new_binary;%p0.binary_variables;%setdiff1D(p1.binary_variables,variable); |
---|
| 732 | %p1.binary_variables = setdiff(p1.binary_variables,friends); |
---|
| 733 | p1.lower = lower; |
---|
| 734 | p1.depth = p.depth+1; |
---|
| 735 | |
---|
| 736 | % % ***************************** |
---|
| 737 | % % PROCESS MOST PROMISING FIRST |
---|
| 738 | % % (p0 in top of stack) |
---|
| 739 | % % ***************************** |
---|
| 740 | if x(variable)>0.5 |
---|
| 741 | pt=p1; |
---|
| 742 | p1=p0; |
---|
| 743 | p0=pt; |
---|
| 744 | end |
---|
| 745 | |
---|
| 746 | function [p0,p1] = integersplit(p,x,index,lower,options,x_min) |
---|
| 747 | |
---|
| 748 | variable = p.integer_variables(index); |
---|
| 749 | current = x(p.integer_variables(index)); |
---|
| 750 | lb = floor(current)+1; |
---|
| 751 | ub = floor(current); |
---|
| 752 | |
---|
| 753 | % xi<ub |
---|
| 754 | p0 = p; |
---|
| 755 | p0.lower = lower; |
---|
| 756 | p0.depth = p.depth+1; |
---|
| 757 | p0.x0(variable) = ub; |
---|
| 758 | p0.ub(variable)=min(p0.ub(variable),ub); |
---|
| 759 | |
---|
| 760 | % xi>lb |
---|
| 761 | p1 = p; |
---|
| 762 | p1.lower = lower; |
---|
| 763 | p1.depth = p.depth+1; |
---|
| 764 | p1.x0(variable) = lb; |
---|
| 765 | p1.lb(variable)=max(p1.lb(variable),lb); |
---|
| 766 | |
---|
| 767 | % ***************************** |
---|
| 768 | % PROCESS MOST PROMISING FIRST |
---|
| 769 | % ***************************** |
---|
| 770 | if lb-current<0.5 |
---|
| 771 | pt=p1; |
---|
| 772 | p1=p0; |
---|
| 773 | p0=pt; |
---|
| 774 | end |
---|
| 775 | |
---|
| 776 | |
---|
| 777 | function s = num2str2(x,d,c); |
---|
| 778 | if nargin==3 |
---|
| 779 | s = num2str(x,c); |
---|
| 780 | else |
---|
| 781 | s = num2str(x); |
---|
| 782 | end |
---|
| 783 | s = [repmat(' ',1,d-length(s)) s]; |
---|
| 784 | |
---|
| 785 | |
---|
| 786 | function [stack,lower] = prune(stack,upper,options,solved_nodes,p) |
---|
| 787 | % ********************************* |
---|
| 788 | % PRUNE STACK W.R.T NEW UPPER BOUND |
---|
| 789 | % ********************************* |
---|
| 790 | if ~isempty(stack) |
---|
| 791 | % toolarge = find([stack.lower]>upper*(1-1e-4)); |
---|
| 792 | toolarge = find([stack.lower]>upper*(1-options.bnb.prunetol)); |
---|
| 793 | if ~isempty(toolarge) |
---|
| 794 | stack(toolarge)=[]; |
---|
| 795 | end |
---|
| 796 | end |
---|
| 797 | |
---|
| 798 | if ~isempty(stack) |
---|
| 799 | lower = min([stack.lower]); |
---|
| 800 | else |
---|
| 801 | lower = upper; |
---|
| 802 | end |
---|
| 803 | |
---|
| 804 | function p = adaptivestrategy(p,upper,solved_nodes) |
---|
| 805 | % **********************************' |
---|
| 806 | % SWITCH NODE SELECTION STRATEGY? |
---|
| 807 | % **********************************' |
---|
| 808 | if strcmp(p.options.bnb.method,'depthproject') & (upper<inf) |
---|
| 809 | p.options.bnb.method = 'project'; |
---|
| 810 | end |
---|
| 811 | if strcmp(p.options.bnb.method,'depthbest') & (upper<inf) |
---|
| 812 | p.options.bnb.method = 'best'; |
---|
| 813 | end |
---|
| 814 | if strcmp(p.options.bnb.method,'depthbreadth') & (upper<inf) |
---|
| 815 | p.options.bnb.method = 'breadth'; |
---|
| 816 | end |
---|
| 817 | if strcmp(p.options.bnb.method,'depthest') & (upper<inf) |
---|
| 818 | p.options.bnb.method = 'est'; |
---|
| 819 | end |
---|
| 820 | |
---|
| 821 | function output = solvelower(lowersolver,relaxed_p,upper,lower) |
---|
| 822 | |
---|
| 823 | if all(relaxed_p.lb==relaxed_p.ub) |
---|
| 824 | x = relaxed_p.lb; |
---|
| 825 | res = resids(relaxed_p,x); |
---|
| 826 | if all(res>-1e-7) |
---|
| 827 | output.problem = 0; |
---|
| 828 | else |
---|
| 829 | output.problem = 1; |
---|
| 830 | end |
---|
| 831 | output.Primal = x; |
---|
| 832 | return |
---|
| 833 | end |
---|
| 834 | |
---|
| 835 | p = relaxed_p; |
---|
| 836 | removethese = p.lb==p.ub; |
---|
| 837 | if nnz(removethese)>0 & all(p.variabletype == 0) & isempty(p.evalMap)% ~isequal(lowersolver,'callfmincongp') & ~isequal(lowersolver,'callgpposy') |
---|
| 838 | |
---|
| 839 | |
---|
| 840 | if ~isinf(upper) & nnz(p.Q)==0 |
---|
| 841 | p.F_struc = [p.F_struc(1:p.K.f,:);upper-p.f -p.c';p.F_struc(1+p.K.f:end,:)]; |
---|
| 842 | p.K.l=p.K.l+1; |
---|
| 843 | end |
---|
| 844 | |
---|
| 845 | if ~isempty(p.F_struc) |
---|
| 846 | p.F_struc(:,1)=p.F_struc(:,1)+p.F_struc(:,1+find(removethese))*p.lb(removethese); |
---|
| 847 | p.F_struc(:,1+find(removethese))=[]; |
---|
| 848 | end |
---|
| 849 | |
---|
| 850 | p.c(removethese)=[]; |
---|
| 851 | if nnz(p.Q)>0 |
---|
| 852 | p.c = p.c + 2*p.Q(find(~removethese),find(removethese))*p.lb(removethese); |
---|
| 853 | p.Q(:,find(removethese))=[]; |
---|
| 854 | p.Q(find(removethese),:)=[]; |
---|
| 855 | else |
---|
| 856 | p.Q = spalloc(length(p.c),length(p.c),0); |
---|
| 857 | end |
---|
| 858 | p.lb(removethese)=[]; |
---|
| 859 | p.ub(removethese)=[]; |
---|
| 860 | p.x0(removethese)=[]; |
---|
| 861 | p.monomtable(:,find(removethese))=[]; |
---|
| 862 | p.monomtable(find(removethese),:)=[]; |
---|
| 863 | p.variabletype = []; % Reset, to messy to recompute |
---|
| 864 | |
---|
| 865 | output = feval(lowersolver,p); |
---|
| 866 | |
---|
| 867 | x=relaxed_p.c*0; |
---|
| 868 | x(removethese)=relaxed_p.lb(removethese); |
---|
| 869 | x(~removethese)=output.Primal; |
---|
| 870 | output.Primal=x; |
---|
| 871 | else |
---|
| 872 | output = feval(lowersolver,p); |
---|
| 873 | end |
---|
| 874 | |
---|
| 875 | |
---|
| 876 | function res = resids(p,x) |
---|
| 877 | res= []; |
---|
| 878 | if p.K.f>0 |
---|
| 879 | res = -abs(p.F_struc(1:p.K.f,:)*[1;x]); |
---|
| 880 | end |
---|
| 881 | if p.K.l>0 |
---|
| 882 | res = [res;p.F_struc(p.K.f+1:p.K.f+p.K.l,:)*[1;x]]; |
---|
| 883 | end |
---|
| 884 | if (length(p.K.s)>1) | p.K.s>0 |
---|
| 885 | top = 1+p.K.f+p.K.l; |
---|
| 886 | for i = 1:length(p.K.s) |
---|
| 887 | n = p.K.s(i); |
---|
| 888 | X = p.F_struc(top:top+n^2-1,:)*[1;x];top = top+n^2; |
---|
| 889 | X = reshape(X,n,n); |
---|
| 890 | res = [res;min(eig(X))]; |
---|
| 891 | end |
---|
| 892 | end |
---|
| 893 | res = [res;min([p.ub-x;x-p.lb])]; |
---|
| 894 | |
---|
| 895 | function p = Updatecostbound(p,upper); |
---|
| 896 | if p.simplecost & ~isinf(upper) |
---|
| 897 | ind = find(p.c); |
---|
| 898 | if p.c(ind)>0 |
---|
| 899 | p.ub(ind) = min(p.ub(ind),(upper-p.f)/p.c(ind)); |
---|
| 900 | else |
---|
| 901 | p.lb(ind) = max(p.lb(ind),(p.f-upper)/abs(p.c(ind))); |
---|
| 902 | end |
---|
| 903 | end |
---|
| 904 | |
---|
| 905 | function [x_min,upper] = initializesolution(p); |
---|
| 906 | |
---|
| 907 | x_min = zeros(length(p.c),1); |
---|
| 908 | upper = inf; |
---|
| 909 | if p.options.usex0 |
---|
| 910 | z = p.x0; |
---|
| 911 | residual = resids(p,z); |
---|
| 912 | relaxed_feasible = all(residual(1:p.K.f)>=-1e-12) & all(residual(1+p.K.f:end)>=-1e-6); |
---|
| 913 | if relaxed_feasible |
---|
| 914 | upper = p.f+p.c'*z+z'*p.Q*z; |
---|
| 915 | x_min = z; |
---|
| 916 | end |
---|
| 917 | else |
---|
| 918 | p.x0 = zeros(length(p.c),1); |
---|
| 919 | x = p.x0; |
---|
| 920 | z = evaluate_nonlinear(p,x); |
---|
| 921 | residual = resids(p,z); |
---|
| 922 | relaxed_feasible = all(residual(1:p.K.f)>=-p.options.bmibnb.eqtol) & all(residual(1+p.K.f:end)>=p.options.bmibnb.pdtol); |
---|
| 923 | if relaxed_feasible |
---|
| 924 | upper = p.f+p.c'*z+z'*p.Q*z; |
---|
| 925 | x_min = x; |
---|
| 926 | end |
---|
| 927 | end |
---|
| 928 | |
---|
| 929 | |
---|
| 930 | |
---|
| 931 | function [p,poriginal,stack] = pruneglobally(p,poriginal,upper,lower,stack,x); |
---|
| 932 | |
---|
| 933 | if isempty(p.nonlinear) & (nnz(p.Q)==0) & p.options.bnb.nodetight |
---|
| 934 | pp = poriginal; |
---|
| 935 | |
---|
| 936 | if p.K.l > 0 |
---|
| 937 | A = -pp.F_struc(1+pp.K.f:pp.K.f+pp.K.l,2:end); |
---|
| 938 | b = pp.F_struc(1+p.K.f:p.K.f+p.K.l,1); |
---|
| 939 | else |
---|
| 940 | A = []; |
---|
| 941 | b = []; |
---|
| 942 | end |
---|
| 943 | |
---|
| 944 | if (nnz(p.Q)==0) & ~isinf(upper) |
---|
| 945 | A = [pp.c';-pp.c';A]; |
---|
| 946 | b = [upper;-(lower-0.0001);b]; |
---|
| 947 | else |
---|
| 948 | % c = p.c; |
---|
| 949 | % Q = p.Q; |
---|
| 950 | % A = [c'+2*x'*Q;A]; |
---|
| 951 | % b = [2*x'*Q*x+c'*x;b]; |
---|
| 952 | end |
---|
| 953 | |
---|
| 954 | [lb,ub,redundant,pss] = milppresolve(A,b,pp.lb,pp.ub,pp.integer_variables,pp.binary_variables,ones(length(pp.lb),1)); |
---|
| 955 | |
---|
| 956 | if ~isempty(redundant) |
---|
| 957 | if (nnz(p.Q)==0) & ~isinf(upper) |
---|
| 958 | redundant = redundant(redundant>2)-2; |
---|
| 959 | else |
---|
| 960 | % redundant = redundant(redundant>1)-1; |
---|
| 961 | end |
---|
| 962 | if length(redundant)>0 |
---|
| 963 | poriginal.K.l=poriginal.K.l-length(redundant); |
---|
| 964 | poriginal.F_struc(poriginal.K.f+redundant,:)=[]; |
---|
| 965 | p.K.l=p.K.l-length(redundant); |
---|
| 966 | p.F_struc(p.K.f+redundant,:)=[]; |
---|
| 967 | end |
---|
| 968 | end |
---|
| 969 | if ~isempty(stack) |
---|
| 970 | keep = ones(length(stack),1); |
---|
| 971 | for i = 1:length(stack) |
---|
| 972 | stack(i).lb = max([stack(i).lb lb]')'; |
---|
| 973 | stack(i).ub = min([stack(i).ub ub]')'; |
---|
| 974 | if any(stack(i).lb>stack(i).ub) |
---|
| 975 | keep(i) = 0; |
---|
| 976 | end |
---|
| 977 | end |
---|
| 978 | stack = stack(find(keep)); |
---|
| 979 | end |
---|
| 980 | poriginal.lb = max([poriginal.lb lb]')'; |
---|
| 981 | poriginal.ub = min([poriginal.ub ub]')'; |
---|
| 982 | p.lb = max([p.lb lb]')'; |
---|
| 983 | p.ub = min([p.ub ub]')'; |
---|
| 984 | end |
---|
| 985 | |
---|
| 986 | |
---|
| 987 | function [p,poriginal,stack] = fixvariables(p,poriginal,upper,lower,stack,x_min,monotinicity) |
---|
| 988 | % Fix variables |
---|
| 989 | |
---|
| 990 | if p.options.bnb.nodefix & (p.K.f == 0) & (nnz(p.Q)==0) & isempty(p.nonlinear) |
---|
| 991 | |
---|
| 992 | A = -poriginal.F_struc(poriginal.K.f + (1:poriginal.K.l),2:end); |
---|
| 993 | b = poriginal.F_struc(poriginal.K.f + (1:poriginal.K.l),1); |
---|
| 994 | c = poriginal.c; |
---|
| 995 | [fix_up,fix_down] = presolve_fixvariables(A,b,c,poriginal.lb,poriginal.ub,monotinicity); |
---|
| 996 | % |
---|
| 997 | poriginal.lb(fix_up) = 1; |
---|
| 998 | p.lb(fix_up) = 1; |
---|
| 999 | |
---|
| 1000 | % not_in_obj = find(p.c==0); |
---|
| 1001 | % constrained_blow = all(poriginal.F_struc(1:poriginal.K.l,1+not_in_obj)>=0,1); |
---|
| 1002 | % sdp_positive = sdpmonotinicity(not_in_obj) == -1; |
---|
| 1003 | % can_fix = not_in_obj(find(constrained_blow & sdp_positive)); |
---|
| 1004 | % |
---|
| 1005 | % still_on = find(p.lb==0 & p.ub==1); |
---|
| 1006 | % p.lb(intersect(can_fix,still_on)) = 1; |
---|
| 1007 | % still_on = find(poriginal.lb==0 & poriginal.ub==1); |
---|
| 1008 | % poriginal.lb(intersect(can_fix,still_on)) = 1; |
---|
| 1009 | |
---|
| 1010 | if ~isempty(stack) & ~(isempty(fix_up) & isempty(fix_down)) |
---|
| 1011 | keep = ones(length(stack),1); |
---|
| 1012 | for i = 1:length(stack) |
---|
| 1013 | stack(i).lb = max([stack(i).lb poriginal.lb]')'; |
---|
| 1014 | stack(i).ub = min([stack(i).ub poriginal.ub]')'; |
---|
| 1015 | if any(stack(i).lb>stack(i).ub) |
---|
| 1016 | keep(i) = 0; |
---|
| 1017 | end |
---|
| 1018 | end |
---|
| 1019 | stack = stack(find(keep)); |
---|
| 1020 | end |
---|
| 1021 | end |
---|
| 1022 | |
---|
| 1023 | |
---|
| 1024 | |
---|
| 1025 | |
---|
| 1026 | |
---|
| 1027 | |
---|