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