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