[37] | 1 | function [x_min,solved_nodes,lower,upper] = branch_and_bound(p,x_min,upper) |
---|
| 2 | |
---|
| 3 | % ************************************************************************* |
---|
| 4 | % Create handles to solvers |
---|
| 5 | % ************************************************************************* |
---|
| 6 | lowersolver = p.solver.lowersolver.call; % For relaxed lower bound problem |
---|
| 7 | uppersolver = p.solver.uppersolver.call; % Local nonlinear upper bound |
---|
| 8 | lpsolver = p.solver.lpsolver.call; % LP solver for bound propagation |
---|
| 9 | |
---|
| 10 | % ************************************************************************* |
---|
| 11 | % GLOBAL PROBLEM DATA (these variables are the same in all nodes) |
---|
| 12 | % ************************************************************************* |
---|
| 13 | c = p.c; |
---|
| 14 | Q = p.Q; |
---|
| 15 | f = p.f; |
---|
| 16 | K = p.K; |
---|
| 17 | options = p.options; |
---|
| 18 | |
---|
| 19 | % ************************************************************************* |
---|
| 20 | % DEFINE UPPER BOUND PROBLEM. Basically just remove the cuts |
---|
| 21 | % ************************************************************************* |
---|
| 22 | p_upper = cleanuppermodel(p); |
---|
| 23 | |
---|
| 24 | % ************************************************************************* |
---|
| 25 | % Active constraints in main model |
---|
| 26 | % 0 : Inactive constraint (i.e. a cut which unused) |
---|
| 27 | % 1 : Active constraint |
---|
| 28 | % inf : Removed constraint (found to be redundant) |
---|
| 29 | % ************************************************************************* |
---|
| 30 | p.InequalityConstraintState = ones(p.K.l,1); |
---|
| 31 | p.InequalityConstraintState(p.KCut.l,1) = 0; |
---|
| 32 | p.EqualityConstraintState = ones(p.K.f,1); |
---|
| 33 | |
---|
| 34 | % ************************************************************************* |
---|
| 35 | % LPs ARE USED IN BOX-REDUCTION |
---|
| 36 | % ************************************************************************* |
---|
| 37 | p.lpcuts = p.F_struc(1+p.K.f:1:p.K.l+p.K.f,:); |
---|
| 38 | p.cutState = ones(p.K.l,1); |
---|
| 39 | p.cutState(p.KCut.l,1) = 0; % Don't use to begin with |
---|
| 40 | |
---|
| 41 | % ************************************************************************* |
---|
| 42 | % INITIALITAZION |
---|
| 43 | % ************************************************************************* |
---|
| 44 | p.depth = 0; % depth in search tree |
---|
| 45 | p.dpos = 0; % used for debugging |
---|
| 46 | p.lower = NaN; |
---|
| 47 | lower = NaN; |
---|
| 48 | gap = inf; |
---|
| 49 | stack = []; |
---|
| 50 | solved_nodes = 0; |
---|
| 51 | numGlobalSolutions = 0; |
---|
| 52 | |
---|
| 53 | % ************************************************************************* |
---|
| 54 | % Silly hack to speed up solver calls |
---|
| 55 | % ************************************************************************* |
---|
| 56 | p.getsolvertime = 0; |
---|
| 57 | |
---|
| 58 | if options.bmibnb.verbose>0 |
---|
| 59 | disp('* Starting YALMIP bilinear branch & bound.'); |
---|
| 60 | disp(['* Upper solver : ' p.solver.uppersolver.tag]); |
---|
| 61 | disp(['* Lower solver : ' p.solver.lowersolver.tag]); |
---|
| 62 | if p.options.bmibnb.lpreduce |
---|
| 63 | disp(['* LP solver : ' p.solver.lpsolver.tag]); |
---|
| 64 | end |
---|
| 65 | disp(' Node Upper Gap(%) Lower Open'); |
---|
| 66 | end |
---|
| 67 | |
---|
| 68 | t_start = cputime; |
---|
| 69 | go_on = 1; |
---|
| 70 | |
---|
| 71 | while go_on |
---|
| 72 | |
---|
| 73 | % ******************************************** |
---|
| 74 | % ASSUME THAT WE WON'T FATHOME |
---|
| 75 | % ******************************************** |
---|
| 76 | keep_digging = 1; |
---|
| 77 | |
---|
| 78 | % ******************************************** |
---|
| 79 | % Reduce size of current box (bound tightening) |
---|
| 80 | % ******************************************** |
---|
| 81 | p = propagatequadratics(p,upper,lower); |
---|
| 82 | p = domain_reduction(p,upper,lower,lpsolver); |
---|
| 83 | p = presolve_bounds_from_equalities(p); |
---|
| 84 | p = preprocess_eval_bounds(p); |
---|
| 85 | |
---|
| 86 | % ******************************************** |
---|
| 87 | % Detect redundant constraints |
---|
| 88 | % ******************************************** |
---|
| 89 | p = remove_redundant(p); |
---|
| 90 | |
---|
| 91 | % ******************************************** |
---|
| 92 | % SOLVE LOWER AND UPPER |
---|
| 93 | % ******************************************** |
---|
| 94 | if p.feasible |
---|
| 95 | [output,cost] = solvelower(p,options,lowersolver); |
---|
| 96 | |
---|
| 97 | % Cplex sucks... |
---|
| 98 | if output.problem == 12 |
---|
| 99 | pp = p; |
---|
| 100 | pp.c = pp.c*0; |
---|
| 101 | [output2,cost2] = solvelower(pp,options,lowersolver); |
---|
| 102 | if output2.problem == 0 |
---|
| 103 | output.problem = 2; |
---|
| 104 | else |
---|
| 105 | output.problem = 1; |
---|
| 106 | end |
---|
| 107 | end |
---|
| 108 | |
---|
| 109 | |
---|
| 110 | |
---|
| 111 | info_text = ''; |
---|
| 112 | switch output.problem |
---|
| 113 | case {1,12} % Infeasible |
---|
| 114 | info_text = 'Infeasible'; |
---|
| 115 | keep_digging = 0; |
---|
| 116 | cost = inf; |
---|
| 117 | feasible = 0; |
---|
| 118 | |
---|
| 119 | case 2 % Unbounded (should not happen!) |
---|
| 120 | cost = -inf; |
---|
| 121 | x = output.Primal; |
---|
| 122 | |
---|
| 123 | case {0,3,4} % No problems (disregard numerical problems) |
---|
| 124 | |
---|
| 125 | x = output.Primal; |
---|
| 126 | |
---|
| 127 | % UPDATE THE LOWER BOUND |
---|
| 128 | if isnan(lower) |
---|
| 129 | lower = cost; |
---|
| 130 | end |
---|
| 131 | if ~isempty(stack) |
---|
| 132 | lower = min(cost,min([stack.lower])); |
---|
| 133 | else |
---|
| 134 | lower = min(lower,cost); |
---|
| 135 | end |
---|
| 136 | |
---|
| 137 | if cost<upper-1e-5 |
---|
| 138 | |
---|
| 139 | z = evaluate_nonlinear(p,x); |
---|
| 140 | |
---|
| 141 | % Manage cuts etc |
---|
| 142 | p = addsdpcut(p,z); |
---|
| 143 | p = addlpcuts(p,x); |
---|
| 144 | |
---|
| 145 | if numGlobalSolutions < p.options.bmibnb.numglobal |
---|
| 146 | [upper,x_min,cost,info_text,numGlobalSolutions] = heuristics_from_relaxed(p_upper,x,upper,x_min,cost,numGlobalSolutions); |
---|
| 147 | [upper,x_min,info_text,numGlobalSolutions] = solve_upper_in_node(p,p_upper,x,upper,x_min,uppersolver,info_text,numGlobalSolutions); |
---|
| 148 | end |
---|
| 149 | else |
---|
| 150 | keep_digging = 0; |
---|
| 151 | info_text = 'Poor bound'; |
---|
| 152 | end |
---|
| 153 | otherwise |
---|
| 154 | cost = -inf; |
---|
| 155 | x = (p.lb+p.ub)/2; |
---|
| 156 | end |
---|
| 157 | else |
---|
| 158 | info_text = 'Infeasible'; |
---|
| 159 | keep_digging = 0; |
---|
| 160 | cost = inf; |
---|
| 161 | feasible = 0; |
---|
| 162 | end |
---|
| 163 | solved_nodes = solved_nodes+1; |
---|
| 164 | |
---|
| 165 | % ************************************************ |
---|
| 166 | % PRUNE SUBOPTIMAL REGIONS BASED ON UPPER BOUND |
---|
| 167 | % ************************************************ |
---|
| 168 | if ~isempty(stack) |
---|
| 169 | [stack,lower] = prune(stack,upper,options,solved_nodes,p); |
---|
| 170 | end |
---|
| 171 | if isempty(stack) |
---|
| 172 | if isinf(cost) |
---|
| 173 | lower = upper; |
---|
| 174 | else |
---|
| 175 | lower = cost; |
---|
| 176 | end |
---|
| 177 | else |
---|
| 178 | lower = min(lower,cost); |
---|
| 179 | end |
---|
| 180 | |
---|
| 181 | % ************************************************ |
---|
| 182 | % CONTINUE SPLITTING? |
---|
| 183 | % ************************************************ |
---|
| 184 | if keep_digging & max(p.ub(p.branch_variables)-p.lb(p.branch_variables))>options.bmibnb.vartol |
---|
| 185 | spliton = branchvariable(p,options,x); |
---|
| 186 | bounds = partition(p,options,spliton,x); |
---|
| 187 | for i = 1:length(bounds)-1 |
---|
| 188 | node = savetonode(p,spliton,bounds(i),bounds(i+1),-1,x,cost,p.EqualityConstraintState,p.InequalityConstraintState,p.cutState); |
---|
| 189 | node.bilinears = p.bilinears; |
---|
| 190 | node = updateonenonlinearbound(node,spliton); |
---|
| 191 | stack = push(stack,node); |
---|
| 192 | end |
---|
| 193 | lower = min([stack.lower]); |
---|
| 194 | end |
---|
| 195 | |
---|
| 196 | % ************************************************ |
---|
| 197 | % Pick and create a suitable node |
---|
| 198 | % ************************************************ |
---|
| 199 | [p,stack] = selectbranch(p,options,stack,x_min,upper); |
---|
| 200 | |
---|
| 201 | if isempty(p) |
---|
| 202 | if ~isinf(upper) |
---|
| 203 | relgap = 0; |
---|
| 204 | end |
---|
| 205 | if isinf(upper) & isinf(lower) |
---|
| 206 | relgap = inf; |
---|
| 207 | end |
---|
| 208 | depth = 0; |
---|
| 209 | else |
---|
| 210 | relgap = 100*(upper-lower)/(1+abs(upper)); |
---|
| 211 | depth = p.depth; |
---|
| 212 | end |
---|
| 213 | if options.bmibnb.verbose>0 |
---|
| 214 | fprintf(' %4.0f : %12.3E %7.2f %12.3E %2.0f %s \n',solved_nodes,upper,relgap,lower,length(stack)+length(p),info_text); |
---|
| 215 | end |
---|
| 216 | |
---|
| 217 | absgap = upper-lower; |
---|
| 218 | % ************************************************ |
---|
| 219 | % Continue? |
---|
| 220 | % ************************************************ |
---|
| 221 | time_ok = cputime-t_start < options.bmibnb.maxtime; |
---|
| 222 | iter_ok = solved_nodes < options.bmibnb.maxiter; |
---|
| 223 | any_nodes = ~isempty(p); |
---|
| 224 | relgap_too_big = (isinf(lower) | isnan(relgap) | relgap>100*options.bmibnb.relgaptol); |
---|
| 225 | absgap_too_big = (isinf(lower) | isnan(absgap) | absgap>options.bmibnb.absgaptol); |
---|
| 226 | taget_not_met = upper>options.bmibnb.target; |
---|
| 227 | go_on = taget_not_met & time_ok & any_nodes & iter_ok & relgap_too_big & absgap_too_big ; |
---|
| 228 | end |
---|
| 229 | if options.bmibnb.verbose>0 |
---|
| 230 | if options.bmibnb.verbose;showprogress([num2str(solved_nodes) ' Finishing. Cost: ' num2str(upper) ' Gap: ' num2str(relgap) '%'],options.bnb.verbose);end |
---|
| 231 | end |
---|
| 232 | |
---|
| 233 | % ************************************************************************* |
---|
| 234 | % Stack functionality |
---|
| 235 | % ************************************************************************* |
---|
| 236 | function stack = push(stackin,p) |
---|
| 237 | if ~isempty(stackin) |
---|
| 238 | stack = [p;stackin]; |
---|
| 239 | else |
---|
| 240 | stack(1)=p; |
---|
| 241 | end |
---|
| 242 | |
---|
| 243 | function [p,stack] = pull(stack,method,x_min,upper,branch_variables); |
---|
| 244 | if ~isempty(stack) |
---|
| 245 | switch method |
---|
| 246 | case 'maxvol' |
---|
| 247 | for i = 1:length(stack) |
---|
| 248 | vol(i) = sum(stack(i).ub(branch_variables)-stack(i).lb(branch_variables)); |
---|
| 249 | end |
---|
| 250 | [i,j] = max(vol); |
---|
| 251 | p=stack(j); |
---|
| 252 | stack = stack([1:1:j-1 j+1:1:end]); |
---|
| 253 | |
---|
| 254 | case 'best' |
---|
| 255 | [i,j]=min([stack.lower]); |
---|
| 256 | p=stack(j); |
---|
| 257 | stack = stack([1:1:j-1 j+1:1:end]); |
---|
| 258 | |
---|
| 259 | otherwise |
---|
| 260 | end |
---|
| 261 | else |
---|
| 262 | p =[]; |
---|
| 263 | end |
---|
| 264 | |
---|
| 265 | function [stack,lower] = prune(stack,upper,options,solved_nodes,p) |
---|
| 266 | if ~isempty(stack) |
---|
| 267 | toolarge = find([stack.lower]>upper*(1+1e-4)); |
---|
| 268 | if ~isempty(toolarge) |
---|
| 269 | stack(toolarge)=[]; |
---|
| 270 | end |
---|
| 271 | if ~isempty(stack) |
---|
| 272 | indPOS = find(p.c>0); |
---|
| 273 | indNEG = find(p.c<0); |
---|
| 274 | LB = [stack.lb]; |
---|
| 275 | UB = [stack.ub]; |
---|
| 276 | LOWER = p.c([indPOS(:);indNEG(:)])'*[LB(indPOS,:);UB(indNEG,:)]; |
---|
| 277 | toolarge = find(LOWER > upper*(1+1e-4)); |
---|
| 278 | stack(toolarge)=[]; |
---|
| 279 | end |
---|
| 280 | end |
---|
| 281 | if ~isempty(stack) |
---|
| 282 | lower = min([stack.lower]); |
---|
| 283 | else |
---|
| 284 | lower = upper; |
---|
| 285 | end |
---|
| 286 | |
---|
| 287 | function node = savetonode(p,spliton,bounds1,bounds2,direction,x,cost,EqualityConstraintState,InequalityConstraintState,cutState); |
---|
| 288 | node.lb = p.lb; |
---|
| 289 | node.ub = p.ub; |
---|
| 290 | node.lb(spliton) = bounds1; |
---|
| 291 | node.ub(spliton) = bounds2; |
---|
| 292 | node.lb(p.integer_variables) = ceil(node.lb(p.integer_variables)); |
---|
| 293 | node.ub(p.integer_variables) = floor(node.ub(p.integer_variables)); |
---|
| 294 | node.lb(p.binary_variables) = ceil(node.lb(p.binary_variables)); |
---|
| 295 | node.ub(p.binary_variables) = floor(node.ub(p.binary_variables)); |
---|
| 296 | |
---|
| 297 | if direction == -1 |
---|
| 298 | node.dpos = p.dpos-1/(2^sqrt(p.depth)); |
---|
| 299 | else |
---|
| 300 | node.dpos = p.dpos+1/(2^sqrt(p.depth)); |
---|
| 301 | end |
---|
| 302 | node.depth = p.depth+1; |
---|
| 303 | node.x0 = x; |
---|
| 304 | node.lpcuts = p.lpcuts; |
---|
| 305 | node.lower = cost; |
---|
| 306 | node.InequalityConstraintState = InequalityConstraintState; |
---|
| 307 | node.EqualityConstraintState = EqualityConstraintState; |
---|
| 308 | node.cutState = cutState; |
---|
| 309 | |
---|
| 310 | % ************************************* |
---|
| 311 | % DERIVE LINEAR CUTS FROM SDPs |
---|
| 312 | % ************************************* |
---|
| 313 | function p = addsdpcut(p,x) |
---|
| 314 | if p.K.s > 0 |
---|
| 315 | top = p.K.f+p.K.l+1; |
---|
| 316 | newcuts = 1; |
---|
| 317 | newF = []; |
---|
| 318 | for i = 1:length(p.K.s) |
---|
| 319 | n = p.K.s(i); |
---|
| 320 | X = p.F_struc(top:top+n^2-1,:)*[1;x]; |
---|
| 321 | X = reshape(X,n,n); |
---|
| 322 | [d,v] = eig(X); |
---|
| 323 | for m = 1:length(v) |
---|
| 324 | if v(m,m)<0 |
---|
| 325 | for j = 1:length(x)+1; |
---|
| 326 | newF(newcuts,j)= d(:,m)'*reshape(p.F_struc(top:top+n^2-1,j),n,n)*d(:,m); |
---|
| 327 | end |
---|
| 328 | % max(abs(newF(:,2:end)),[],2) |
---|
| 329 | newF(newcuts,1)=newF(newcuts,1)+1e-6; |
---|
| 330 | newcuts = newcuts + 1; |
---|
| 331 | if size(p.lpcuts,1)>0 |
---|
| 332 | dist = p.lpcuts*newF(newcuts-1,:)'/(newF(newcuts-1,:)*newF(newcuts-1,:)'); |
---|
| 333 | if any(abs(dist-1)<1e-3) |
---|
| 334 | newF = newF(1:end-1,:); |
---|
| 335 | newcuts = newcuts - 1; |
---|
| 336 | end |
---|
| 337 | end |
---|
| 338 | end |
---|
| 339 | end |
---|
| 340 | top = top+n^2; |
---|
| 341 | end |
---|
| 342 | |
---|
| 343 | if ~isempty(newF) |
---|
| 344 | % Don't keep all |
---|
| 345 | m = size(newF,2); |
---|
| 346 | % size(p.lpcuts) |
---|
| 347 | p.lpcuts = [newF;p.lpcuts]; |
---|
| 348 | p.cutState = [ones(size(newF,1),1);p.cutState]; |
---|
| 349 | violations = p.lpcuts*[1;x]; |
---|
| 350 | p.lpcuts = p.lpcuts(violations<0.1,:); |
---|
| 351 | |
---|
| 352 | if size(p.lpcuts,1)>15*m |
---|
| 353 | disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'); |
---|
| 354 | violations = p.lpcuts*[1;x]; |
---|
| 355 | [i,j] = sort(violations); |
---|
| 356 | %p.lpcuts = p.lpcuts(j(1:15*m),:); |
---|
| 357 | %p.cutState = lpcuts = p.lpcuts(j(1:15*m),:); |
---|
| 358 | %p.lpcuts = p.lpcuts(end-15*m+1:end,:); |
---|
| 359 | end |
---|
| 360 | end |
---|
| 361 | end |
---|
| 362 | |
---|
| 363 | function p = addlpcuts(p,z) |
---|
| 364 | inactiveCuts = find(~p.cutState); |
---|
| 365 | violation = p.lpcuts(inactiveCuts,:)*[1;z]; |
---|
| 366 | need_to_add = find(violation < -1e-4); |
---|
| 367 | if ~isempty(need_to_add) |
---|
| 368 | p.cutState(inactiveCuts(need_to_add)) = 1; |
---|
| 369 | end |
---|
| 370 | inactiveCuts = find(p.InequalityConstraintState == 0 ); |
---|
| 371 | violation = p.F_struc(p.K.f+inactiveCuts,:)*[1;z]; |
---|
| 372 | need_to_add = find(violation < -1e-4); |
---|
| 373 | if ~isempty(need_to_add) |
---|
| 374 | p.InequalityConstraintState(inactiveCuts(need_to_add)) = 1; |
---|
| 375 | end |
---|
| 376 | |
---|
| 377 | |
---|
| 378 | % ************************************************************************* |
---|
| 379 | % Strategy for deciding which variable to branch on |
---|
| 380 | % ************************************************************************* |
---|
| 381 | function spliton = branchvariable(p,options,x) |
---|
| 382 | % Split if box is to narrow |
---|
| 383 | width = abs(p.ub(p.branch_variables)-p.lb(p.branch_variables)); |
---|
| 384 | if (min(width)/max(width) < 0.1) | (size(p.bilinears,1)==0) % |
---|
| 385 | [i,j] = max(width);%.*p.weight(p.branch_variables)); |
---|
| 386 | spliton = p.branch_variables(j); |
---|
| 387 | else |
---|
| 388 | res = x(p.bilinears(:,1))-x(p.bilinears(:,2)).*x(p.bilinears(:,3)); |
---|
| 389 | [ii,jj] = sort(abs(res)); |
---|
| 390 | v1 = p.bilinears(jj(end),2); |
---|
| 391 | v2 = p.bilinears(jj(end),3); |
---|
| 392 | |
---|
| 393 | acc_res1 = sum(abs(res(find((p.bilinears(:,2)==v1) | p.bilinears(:,3)==v1)))); |
---|
| 394 | acc_res2 = sum(abs(res(find((p.bilinears(:,2)==v2) | p.bilinears(:,3)==v2)))); |
---|
| 395 | |
---|
| 396 | if ~ismember(v2,p.branch_variables) | (acc_res1>acc_res2) |
---|
| 397 | spliton = v1; |
---|
| 398 | else |
---|
| 399 | spliton = v2; |
---|
| 400 | end |
---|
| 401 | end |
---|
| 402 | |
---|
| 403 | % ************************************************************************* |
---|
| 404 | % Strategy for diving the search space |
---|
| 405 | % ************************************************************************* |
---|
| 406 | function bounds = partition(p,options,spliton,x_min) |
---|
| 407 | x = x_min; |
---|
| 408 | if isinf(p.lb(spliton)) |
---|
| 409 | p.lb(spliton) = -1e6; |
---|
| 410 | end |
---|
| 411 | if isinf(p.ub(spliton)) |
---|
| 412 | p.ub(spliton) = 1e6; |
---|
| 413 | end |
---|
| 414 | |
---|
| 415 | switch options.bmibnb.branchrule |
---|
| 416 | case 'omega' |
---|
| 417 | if ~isempty(x_min) |
---|
| 418 | U = p.ub(spliton); |
---|
| 419 | L = p.lb(spliton); |
---|
| 420 | x = x(spliton); |
---|
| 421 | 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)]; |
---|
| 422 | else |
---|
| 423 | bounds = [p.lb(spliton) (p.lb(spliton)+p.ub(spliton))/2 p.ub(spliton)]; |
---|
| 424 | end |
---|
| 425 | case 'bisect' |
---|
| 426 | bounds = [p.lb(spliton) (p.lb(spliton)+p.ub(spliton))/2 p.ub(spliton)]; |
---|
| 427 | otherwise |
---|
| 428 | bounds = [p.lb(spliton) (p.lb(spliton)+p.ub(spliton))/2 p.ub(spliton)]; |
---|
| 429 | end |
---|
| 430 | if isnan(bounds(2)) %FIX |
---|
| 431 | if isinf(p.lb(spliton)) |
---|
| 432 | p.lb(spliton) = -1e6; |
---|
| 433 | end |
---|
| 434 | if isinf(p.ub(spliton)) |
---|
| 435 | p.ub(spliton) = 1e6; |
---|
| 436 | end |
---|
| 437 | bounds(2) = (p.lb(spliton)+p.ub(spliton))/2; |
---|
| 438 | end |
---|
| 439 | |
---|
| 440 | function [p,stack] = selectbranch(p,options,stack,x_min,upper) |
---|
| 441 | switch options.bmibnb.branchmethod |
---|
| 442 | case 'maxvol' |
---|
| 443 | [node,stack] = pull(stack,'maxvol',x_min,upper,p.branch_variables); |
---|
| 444 | case 'best' |
---|
| 445 | [node,stack] = pull(stack,'best',x_min,upper); |
---|
| 446 | otherwise |
---|
| 447 | [node,stack] = pull(stack,'best',x_min,upper); |
---|
| 448 | end |
---|
| 449 | % Copy node data to p |
---|
| 450 | if isempty(node) |
---|
| 451 | p = []; |
---|
| 452 | else |
---|
| 453 | p.depth = node.depth; |
---|
| 454 | p.dpos = node.dpos; |
---|
| 455 | p.lb = node.lb; |
---|
| 456 | p.ub = node.ub; |
---|
| 457 | p.lower = node.lower; |
---|
| 458 | p.lpcuts = node.lpcuts; |
---|
| 459 | p.x0 = node.x0; |
---|
| 460 | p.InequalityConstraintState = node.InequalityConstraintState; |
---|
| 461 | p.EqualityConstraintState = node.EqualityConstraintState; |
---|
| 462 | p.cutState = node.cutState; |
---|
| 463 | end |
---|
| 464 | |
---|
| 465 | |
---|
| 466 | |
---|
| 467 | function [output,cost] = solvelower(p,options,lowersolver) |
---|
| 468 | |
---|
| 469 | removeThese = find(p.InequalityConstraintState==inf); |
---|
| 470 | p.F_struc(p.K.f + removeThese,:) = []; |
---|
| 471 | p.K.l = p.K.l - length(removeThese); |
---|
| 472 | |
---|
| 473 | removeThese = find(p.EqualityConstraintState==inf); |
---|
| 474 | p.F_struc(removeThese,:) = []; |
---|
| 475 | p.K.f = p.K.f - length(removeThese); |
---|
| 476 | |
---|
| 477 | p_with_bilinear_cuts = p; |
---|
| 478 | |
---|
| 479 | if ~isempty(p.bilinears) |
---|
| 480 | p_with_bilinear_cuts.F_struc(1:p.K.f,:)=[]; |
---|
| 481 | p_with_bilinear_cuts = addmcgormick(p_with_bilinear_cuts); |
---|
| 482 | p_with_bilinear_cuts.F_struc = [p.F_struc(1:p.K.f,:);p_with_bilinear_cuts.F_struc]; |
---|
| 483 | end |
---|
| 484 | |
---|
| 485 | if ~isempty(p.evalMap) |
---|
| 486 | p_with_bilinear_cuts = addEvalVariableCuts(p_with_bilinear_cuts); |
---|
| 487 | end |
---|
| 488 | |
---|
| 489 | % ************************************** |
---|
| 490 | % SOLVE NODE PROBLEM |
---|
| 491 | % ************************************** |
---|
| 492 | if any(p_with_bilinear_cuts.ub+1e-8<p_with_bilinear_cuts.lb) |
---|
| 493 | output.problem=1; |
---|
| 494 | cost = inf; |
---|
| 495 | else |
---|
| 496 | % We are solving relaxed problem (penbmi might be local solver) |
---|
| 497 | p_with_bilinear_cuts.monomtable = eye(length(p_with_bilinear_cuts.c)); |
---|
| 498 | |
---|
| 499 | if p.solver.lowersolver.objective.quadratic.convex |
---|
| 500 | % Setup quadratic |
---|
| 501 | for i = 1:size(p.bilinears,1) |
---|
| 502 | if p_with_bilinear_cuts.c(p.bilinears(i,1)) |
---|
| 503 | p_with_bilinear_cuts.Q(p.bilinears(i,2),p.bilinears(i,3)) = p_with_bilinear_cuts.c(p.bilinears(i,1))/2; |
---|
| 504 | p_with_bilinear_cuts.Q(p.bilinears(i,3),p.bilinears(i,2)) = p_with_bilinear_cuts.Q(p.bilinears(i,3),p.bilinears(i,2))+p_with_bilinear_cuts.c(p.bilinears(i,1))/2; |
---|
| 505 | p_with_bilinear_cuts.c(p.bilinears(i,1)) = 0; |
---|
| 506 | end |
---|
| 507 | end |
---|
| 508 | if ~all(eig(full(p_with_bilinear_cuts.Q))>-1e-12) |
---|
| 509 | p_with_bilinear_cuts.Q = p.Q; |
---|
| 510 | p_with_bilinear_cuts.c = p.c; |
---|
| 511 | end |
---|
| 512 | end |
---|
| 513 | |
---|
| 514 | fixed = p_with_bilinear_cuts.lb >= p_with_bilinear_cuts.ub; |
---|
| 515 | if isempty(fixed) |
---|
| 516 | output = feval(lowersolver,p_with_bilinear_cuts); |
---|
| 517 | cost = output.Primal'*p_with_bilinear_cuts.Q*output.Primal + p_with_bilinear_cuts.c'*output.Primal + p.f; |
---|
| 518 | % Minor clean-up |
---|
| 519 | output.Primal(output.Primal<p.lb) = p.lb(output.Primal<p.lb); |
---|
| 520 | output.Primal(output.Primal>p.ub) = p.ub(output.Primal>p.ub); |
---|
| 521 | else |
---|
| 522 | pp = p_with_bilinear_cuts; |
---|
| 523 | removethese = fixed; |
---|
| 524 | if ~isempty(p_with_bilinear_cuts.F_struc) |
---|
| 525 | p_with_bilinear_cuts.F_struc(:,1)=p_with_bilinear_cuts.F_struc(:,1)+p_with_bilinear_cuts.F_struc(:,1+find(fixed))*p_with_bilinear_cuts.lb(fixed); |
---|
| 526 | p_with_bilinear_cuts.F_struc(:,1+find(fixed))=[]; |
---|
| 527 | |
---|
| 528 | rf = find(~any(p_with_bilinear_cuts.F_struc,2)); |
---|
| 529 | rf = rf(rf<=(p_with_bilinear_cuts.K.f + p_with_bilinear_cuts.K.l)); |
---|
| 530 | p_with_bilinear_cuts.F_struc(rf,:) = []; |
---|
| 531 | p_with_bilinear_cuts.K.f = p_with_bilinear_cuts.K.f - nnz(rf<=p_with_bilinear_cuts.K.f); |
---|
| 532 | p_with_bilinear_cuts.K.l = p_with_bilinear_cuts.K.l - nnz(rf>p_with_bilinear_cuts.K.f); |
---|
| 533 | end |
---|
| 534 | p_with_bilinear_cuts.c(removethese)=[]; |
---|
| 535 | if nnz(p_with_bilinear_cuts.Q)>0 |
---|
| 536 | p_with_bilinear_cuts.c = p_with_bilinear_cuts.c + 2*p_with_bilinear_cuts.Q(find(~removethese),find(removethese))*p_with_bilinear_cuts.lb(removethese); |
---|
| 537 | p_with_bilinear_cuts.Q(:,find(removethese))=[]; |
---|
| 538 | p_with_bilinear_cuts.Q(find(removethese),:)=[]; |
---|
| 539 | else |
---|
| 540 | p_with_bilinear_cuts.Q = spalloc(length(p_with_bilinear_cuts.c),length(p_with_bilinear_cuts.c),0); |
---|
| 541 | end |
---|
| 542 | |
---|
| 543 | if ~isempty(p_with_bilinear_cuts.binary_variables) |
---|
| 544 | new_bin = []; |
---|
| 545 | new_var = find(~fixed); |
---|
| 546 | for i = 1:length(p_with_bilinear_cuts.binary_variables) |
---|
| 547 | temp = find(p_with_bilinear_cuts.binary_variables(i) == new_var); |
---|
| 548 | new_bin = [new_bin temp(:)']; |
---|
| 549 | end |
---|
| 550 | p_with_bilinear_cuts.binary_variables = new_bin; |
---|
| 551 | end |
---|
| 552 | if ~isempty(p_with_bilinear_cuts.integer_variables) |
---|
| 553 | new_bin = []; |
---|
| 554 | new_var = find(~fixed); |
---|
| 555 | for i = 1:length(p_with_bilinear_cuts.integer_variables) |
---|
| 556 | temp = find(p_with_bilinear_cuts.integer_variables(i) == new_var); |
---|
| 557 | new_bin = [new_bin temp(:)']; |
---|
| 558 | end |
---|
| 559 | p_with_bilinear_cuts.integer_variables = new_bin; |
---|
| 560 | end |
---|
| 561 | |
---|
| 562 | p_with_bilinear_cuts.lb(removethese)=[]; |
---|
| 563 | p_with_bilinear_cuts.ub(removethese)=[]; |
---|
| 564 | p_with_bilinear_cuts.x0(removethese)=[]; |
---|
| 565 | p_with_bilinear_cuts.monomtable(:,find(removethese))=[]; |
---|
| 566 | p_with_bilinear_cuts.monomtable(find(removethese),:)=[]; |
---|
| 567 | output = feval(lowersolver,p_with_bilinear_cuts); |
---|
| 568 | x=p.c*0; |
---|
| 569 | x(removethese)=p.lb(removethese); |
---|
| 570 | x(~removethese)=output.Primal; |
---|
| 571 | output.Primal = x; |
---|
| 572 | cost = output.Primal'*pp.Q*output.Primal + pp.c'*output.Primal + p.f; |
---|
| 573 | end |
---|
| 574 | end |
---|
| 575 | |
---|
| 576 | function output = solveupper(p,p_original,x,options,uppersolver) |
---|
| 577 | |
---|
| 578 | % The bounds and relaxed solutions have been computed w.r.t to the relaxed |
---|
| 579 | % bilinear model. We only need the original bounds and variables. |
---|
| 580 | p.lb = p.lb(1:length(p_original.c)); |
---|
| 581 | p.ub = p.ub(1:length(p_original.c)); |
---|
| 582 | x = x(1:length(p_original.c)); |
---|
| 583 | |
---|
| 584 | p_upper = p_original; |
---|
| 585 | |
---|
| 586 | % ...expand the current node just slightly |
---|
| 587 | p_upper.lb = p.lb; |
---|
| 588 | p_upper.ub = p.ub; |
---|
| 589 | fixed = find(abs([p.lb-p.ub]) < 1e-5); |
---|
| 590 | p_upper.lb(fixed) = (p.lb(fixed) +p.ub(fixed) )/2; |
---|
| 591 | p_upper.ub(fixed) = (p.lb(fixed) +p.ub(fixed) )/2; |
---|
| 592 | |
---|
| 593 | % Pick an initial point (this can be a bit tricky...) |
---|
| 594 | % Use relaxed point, shifted towards center of box |
---|
| 595 | if all(x<=p.ub) & all(x>=p.lb) |
---|
| 596 | p_upper.x0 = 0.1*x + 0.9*(p.lb+p.ub)/2; |
---|
| 597 | else |
---|
| 598 | p_upper.x0 = (p.lb+p.ub)/2; |
---|
| 599 | end |
---|
| 600 | % Shift towards interior for variables with unbounded lower or upper |
---|
| 601 | lbinfbounds = find(isinf(p.lb)); |
---|
| 602 | ubinfbounds = find(isinf(p.ub)); |
---|
| 603 | p_upper.x0(ubinfbounds) = x(ubinfbounds)+0.01; |
---|
| 604 | p_upper.x0(lbinfbounds) = x(lbinfbounds)-0.01; |
---|
| 605 | ublbinfbounds = find(isinf(p.lb) & isinf(p.ub)); |
---|
| 606 | p_upper.x0(ublbinfbounds) = x(ublbinfbounds); |
---|
| 607 | |
---|
| 608 | change_these_lb = setdiff(1:length(p.lb),fixed); |
---|
| 609 | change_these_lb = setdiff(change_these_lb,lbinfbounds); |
---|
| 610 | change_these_ub = setdiff(1:length(p.lb),fixed); |
---|
| 611 | change_these_ub = setdiff(change_these_ub,lbinfbounds); |
---|
| 612 | |
---|
| 613 | p_upper.lb(change_these_lb) = 0.99*p.lb(change_these_lb)+p_original.lb(change_these_lb)*0.01; |
---|
| 614 | p_upper.ub(change_these_ub) = 0.99*p.ub(change_these_ub)+p_original.ub(change_these_ub)*0.01; |
---|
| 615 | p_upper.lb(isinf(p_original.lb)) = p_upper.lb(isinf(p_original.lb)) - 0.001; |
---|
| 616 | p_upper.ub(isinf(p_original.ub)) = p_upper.ub(isinf(p_original.ub)) + 0.001; |
---|
| 617 | |
---|
| 618 | p_upper.options.saveduals = 0; |
---|
| 619 | |
---|
| 620 | ub = p_upper.ub ; |
---|
| 621 | lb = p_upper.lb ; |
---|
| 622 | |
---|
| 623 | % Remove redundant equality constraints (important for fmincon) |
---|
| 624 | if p_upper.K.f > 0 |
---|
| 625 | Aeq = -p_upper.F_struc(1:1:p_upper.K.f,2:end); |
---|
| 626 | beq = p_upper.F_struc(1:1:p_upper.K.f,1); |
---|
| 627 | redundant = find(((Aeq>0).*Aeq*(p_upper.ub-p_upper.lb) - (beq-Aeq*p_upper.lb) <1e-6)); |
---|
| 628 | p_upper.F_struc(redundant,:) = []; |
---|
| 629 | p_upper.K.f = p_upper.K.f - length(redundant); |
---|
| 630 | end |
---|
| 631 | |
---|
| 632 | % Solve upper bounding problem |
---|
| 633 | p_upper.options.usex0 = 1; |
---|
| 634 | |
---|
| 635 | output = feval(uppersolver,p_upper); |
---|
| 636 | % Project into the box (numerical issue) |
---|
| 637 | output.Primal(output.Primal<p_upper.lb) = p_upper.lb(output.Primal<p_upper.lb); |
---|
| 638 | output.Primal(output.Primal>p_upper.ub) = p_upper.ub(output.Primal>p_upper.ub); |
---|
| 639 | |
---|
| 640 | |
---|
| 641 | |
---|
| 642 | |
---|
| 643 | |
---|
| 644 | % ************************************************************************* |
---|
| 645 | % Heuristics from relaxed |
---|
| 646 | % Basically nothing coded yet. Just check feasibility... |
---|
| 647 | % ************************************************************************* |
---|
| 648 | function [upper,x_min,cost,info_text,numglobals] = heuristics_from_relaxed(p_upper,x,upper,x_min,cost,numglobals) |
---|
| 649 | x(p_upper.binary_variables) = round(x(p_upper.binary_variables)); |
---|
| 650 | x(p_upper.integer_variables) = round(x(p_upper.integer_variables)); |
---|
| 651 | |
---|
| 652 | z = evaluate_nonlinear(p_upper,x); |
---|
| 653 | |
---|
| 654 | relaxed_residual = constraint_residuals(p_upper,z); |
---|
| 655 | |
---|
| 656 | eq_ok = all(relaxed_residual(1:p_upper.K.f)>=-p_upper.options.bmibnb.eqtol); |
---|
| 657 | iq_ok = all(relaxed_residual(1+p_upper.K.f:end)>=p_upper.options.bmibnb.pdtol); |
---|
| 658 | |
---|
| 659 | relaxed_feasible = eq_ok & iq_ok; |
---|
| 660 | info_text = ''; |
---|
| 661 | if relaxed_feasible |
---|
| 662 | this_upper = p_upper.f+p_upper.c'*z+z'*p_upper.Q*z; |
---|
| 663 | if (this_upper < (1-1e-5)*upper) & (this_upper < upper - 1e-5) |
---|
| 664 | x_min = x; |
---|
| 665 | upper = this_upper; |
---|
| 666 | info_text = 'Improved solution'; |
---|
| 667 | cost = cost-1e-10; % Otherwise we'll fathome! |
---|
| 668 | numglobals = numglobals + 1; |
---|
| 669 | end |
---|
| 670 | end |
---|
| 671 | |
---|
| 672 | % ************************************************************************* |
---|
| 673 | % Solve local upper bound problem |
---|
| 674 | % ************************************************************************* |
---|
| 675 | function [upper,x_min,info_text,numglobals] = solve_upper_in_node(p,p_upper,x,upper,x_min,uppersolver,info_text,numglobals); |
---|
| 676 | |
---|
| 677 | output = solveupper(p,p_upper,x,p.options,uppersolver); |
---|
| 678 | output.Primal(p_upper.integer_variables) = round(output.Primal(p_upper.integer_variables)); |
---|
| 679 | output.Primal(p_upper.binary_variables) = round(output.Primal(p_upper.binary_variables)); |
---|
| 680 | |
---|
| 681 | xu = evaluate_nonlinear(p_upper,output.Primal); |
---|
| 682 | |
---|
| 683 | upper_residual = constraint_residuals(p_upper,xu); |
---|
| 684 | feasible = ~isempty(xu) & ~any(isnan(xu)) & all(upper_residual(1:p_upper.K.f)>=-p.options.bmibnb.eqtol) & all(upper_residual(1+p_upper.K.f:end)>=p.options.bmibnb.pdtol); |
---|
| 685 | |
---|
| 686 | if feasible |
---|
| 687 | this_upper = p_upper.f+p_upper.c'*xu+xu'*p_upper.Q*xu; |
---|
| 688 | if (this_upper < (1-1e-5)*upper) & (this_upper < upper - 1e-5) |
---|
| 689 | x_min = xu; |
---|
| 690 | upper = this_upper; |
---|
| 691 | info_text = 'Improved solution'; |
---|
| 692 | numglobals = numglobals + 1; |
---|
| 693 | end |
---|
| 694 | end |
---|
| 695 | |
---|
| 696 | % ************************************************************************* |
---|
| 697 | % Detect redundant constraints |
---|
| 698 | % ************************************************************************* |
---|
| 699 | function p = remove_redundant(p); |
---|
| 700 | |
---|
| 701 | b = p.F_struc(1+p.K.f:p.K.l+p.K.f,1); |
---|
| 702 | A = -p.F_struc(1+p.K.f:p.K.l+p.K.f,2:end); |
---|
| 703 | |
---|
| 704 | redundant = find(((A>0).*A*(p.ub-p.lb) - (b-A*p.lb) <-1e-2)); |
---|
| 705 | |
---|
| 706 | if length(redundant)>1 |
---|
| 707 | p.InequalityConstraintState(redundant) = inf; |
---|
| 708 | end |
---|
| 709 | |
---|
| 710 | if p.options.bmibnb.lpreduce |
---|
| 711 | b = p.lpcuts(:,1); |
---|
| 712 | A = -p.lpcuts(:,2:end); |
---|
| 713 | redundant = find(((A>0).*A*(p.ub-p.lb) - (b-A*p.lb) <-1e-2)); |
---|
| 714 | if length(redundant)>1 |
---|
| 715 | p.lpcuts(redundant,:) = []; |
---|
| 716 | p.cutState(redundant) = []; |
---|
| 717 | end |
---|
| 718 | end |
---|
| 719 | |
---|
| 720 | if p.K.f > 0 |
---|
| 721 | b = p.F_struc(1:p.K.f,1); |
---|
| 722 | A = -p.F_struc(1:p.K.f,2:end); |
---|
| 723 | s1 = ((A>0).*A*(p.ub-p.lb) - (b-A*p.lb) <1e-6); |
---|
| 724 | s2 = ((-A>0).*(-A)*(p.ub-p.lb) - ((-b)-(-A)*p.lb) <1e-6); |
---|
| 725 | redundant = find(s1 & s2); |
---|
| 726 | if length(redundant)>1 |
---|
| 727 | p.EqualityConstraintState(redundant) = inf; |
---|
| 728 | end |
---|
| 729 | end |
---|
| 730 | |
---|
| 731 | % ************************************************************************* |
---|
| 732 | % Clean the upper bound model |
---|
| 733 | % Remove cut constraints, and |
---|
| 734 | % possibly unused variables not used |
---|
| 735 | % ************************************************************************* |
---|
| 736 | function p = cleanuppermodel(p); |
---|
| 737 | |
---|
| 738 | % We might have created a bilinear model from an original polynomial model. |
---|
| 739 | % We should use the original model when we solve the upper bound problem. |
---|
| 740 | p_bilinear = p; |
---|
| 741 | p = p.originalModel; |
---|
| 742 | |
---|
| 743 | % Remove cuts |
---|
| 744 | p.F_struc(p.K.f+p.KCut.l,:)=[]; |
---|
| 745 | p.K.l = p.K.l - length(p.KCut.l); |
---|
| 746 | n_start = length(p.c); |
---|
| 747 | |
---|
| 748 | % Quadratic mode, and quadratic aware solver? |
---|
| 749 | bilinear_variables = find(p.variabletype == 1 | p.variabletype == 2); |
---|
| 750 | if ~isempty(bilinear_variables) |
---|
| 751 | used_in_c = find(p.c); |
---|
| 752 | quadraticterms = used_in_c(find(ismember(used_in_c,bilinear_variables))); |
---|
| 753 | if ~isempty(quadraticterms) & p.solver.uppersolver.objective.quadratic.nonconvex |
---|
| 754 | usedinquadratic = zeros(1,length(p.c)); |
---|
| 755 | for i = 1:length(quadraticterms) |
---|
| 756 | Qij = p.c(quadraticterms(i)); |
---|
| 757 | power_index = find(p.monomtable(quadraticterms(i),:)); |
---|
| 758 | if length(power_index) == 1 |
---|
| 759 | p.Q(power_index,power_index) = Qij; |
---|
| 760 | else |
---|
| 761 | p.Q(power_index(1),power_index(2)) = Qij/2; |
---|
| 762 | p.Q(power_index(2),power_index(1)) = Qij/2; |
---|
| 763 | end |
---|
| 764 | p.c(quadraticterms(i)) = 0; |
---|
| 765 | end |
---|
| 766 | end |
---|
| 767 | end |
---|
| 768 | |
---|
| 769 | % Remove SDP cuts |
---|
| 770 | if length(p.KCut.s)>0 |
---|
| 771 | starts = p.K.f+p.K.l + [1 1+cumsum((p.K.s).^2)]; |
---|
| 772 | remove_these = []; |
---|
| 773 | for i = 1:length(p.KCut.s) |
---|
| 774 | j = p.KCut.s(i); |
---|
| 775 | remove_these = [remove_these;(starts(j):starts(j+1)-1)']; |
---|
| 776 | end |
---|
| 777 | p.F_struc(remove_these,:)=[]; |
---|
| 778 | p.K.s(p.KCut.s) = []; |
---|
| 779 | end |
---|
| 780 | p.lb = p_bilinear.lb(1:length(p.c)); |
---|
| 781 | p.ub = p_bilinear.ub(1:length(p.c)); |
---|
| 782 | p.bilinears = []; |
---|