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 = []; |
---|