1 | function output = bmibnb(p) |
---|
2 | %BMIBNB Branch-and-bound scheme for bilinear programs |
---|
3 | % |
---|
4 | % BMIBNB is never called by the user directly, but is called by YALMIP from |
---|
5 | % SOLVESDP, by choosing the solver tag 'bmibnb' in sdpsettings |
---|
6 | % |
---|
7 | % The behaviour of BMIBNB can be altered using the fields in the field |
---|
8 | % 'bmibnb' in SDPSETTINGS |
---|
9 | % |
---|
10 | % bmibnb.lowersolver - Solver for lower bound [solver tag ('')] |
---|
11 | % bmibnb.uppersolver - Solver for upper bound [solver tag ('')] |
---|
12 | % bmibnb.lpsolver - Solver for LP bound tightening [solver tag ('')] |
---|
13 | % bmibnb.branchmethod - Branch strategy ['maxvol' | 'best' ('best')] |
---|
14 | % bmibnb.branchrule - Branch position ['omega' | 'bisect' ('omega')] |
---|
15 | % bmibnb.lpreduce - Improve variable bounds using LP [ real [0,1] (0 means no reduction, 1 means all variables) |
---|
16 | % bmibnb.lowrank - Partition variables into two disjoint sets and branch on smallest [ 0|1 (0)] |
---|
17 | % bmibnb.target - Exit if upper bound<target [double (-inf)] |
---|
18 | % bmibnb.roottight - Improve variable bounds in root using full problem [ 0|1 (1)] |
---|
19 | % bmibnb.vartol - Cut tree when x_U-x_L < vartol on all branching variables |
---|
20 | % bmibnb.relgaptol - Tolerance on relative objective error (UPPER-LOWER)/(1+|UPPER|) [real (0.01)] |
---|
21 | % bmibnb.absgaptol - Tolerance on objective error (UPPER-LOWER) [real (0.01)] |
---|
22 | % bmibnb.pdtol - A number is declared non-negative if larger than...[ double (-1e-6)] |
---|
23 | % bmibnb.eqtol - A number is declared zero if abs(x) smaller than...[ double (1e-6)] |
---|
24 | % bmibnb.maxiter - Maximum number nodes [int (100)] |
---|
25 | % bmibnb.maxtime - Maximum CPU time (sec.) [int (3600)] |
---|
26 | |
---|
27 | % Author Johan Löfberg |
---|
28 | % $Id: bmibnb.m,v 1.13 2006/09/20 12:43:45 joloef Exp $ |
---|
29 | |
---|
30 | bnbsolvertime = clock; |
---|
31 | showprogress('Branch and bound started',p.options.showprogress); |
---|
32 | |
---|
33 | % ************************************************************************* |
---|
34 | % INITIALIZE DIAGNOSTICS ETC |
---|
35 | % ************************************************************************* |
---|
36 | switch max(min(p.options.verbose,3),0) |
---|
37 | case 0 |
---|
38 | p.options.bmibnb.verbose = 0; |
---|
39 | case 1 |
---|
40 | p.options.bmibnb.verbose = 1; |
---|
41 | p.options.verbose = 0; |
---|
42 | case 2 |
---|
43 | p.options.bmibnb.verbose = 2; |
---|
44 | p.options.verbose = 0; |
---|
45 | case 3 |
---|
46 | p.options.bmibnb.verbose = 2; |
---|
47 | p.options.verbose = 1; |
---|
48 | otherwise |
---|
49 | p.options.bmibnb.verbose = 0; |
---|
50 | p.options.verbose = 0; |
---|
51 | end |
---|
52 | |
---|
53 | % ************************************************************************* |
---|
54 | % Assume feasible (this property can be changed in the presolve codes |
---|
55 | % ************************************************************************* |
---|
56 | p.feasible = 1; |
---|
57 | |
---|
58 | % ************************************************************************* |
---|
59 | % Extract bounds from model |
---|
60 | % ************************************************************************* |
---|
61 | [p.lb,p.ub,used_rows] = findulb(p.F_struc,p.K); |
---|
62 | |
---|
63 | % These two are needed for ex14_1_7 |
---|
64 | p = fixbounds(p); |
---|
65 | p = presolve_bounds_from_equalities(p); |
---|
66 | p = preprocess_eval_bounds(p); |
---|
67 | p = fixbounds(p); |
---|
68 | |
---|
69 | % ************************************************************************* |
---|
70 | % The bilinear solver does not support higher order polynomial (i.e. >2). |
---|
71 | % Hence, we need to bilinearize the problem. however, the upper bound |
---|
72 | % solver might support polynomial problems, so we should keep a copy of the |
---|
73 | % original problem also |
---|
74 | % ************************************************************************* |
---|
75 | n_in = length(p.c); |
---|
76 | [p_bilin,changed] = bilinearize(p); |
---|
77 | if (p.solver.uppersolver.constraint.equalities.polynomial & p.solver.uppersolver.objective.polynomial) |
---|
78 | p_bilin.originalModel = p; |
---|
79 | p = p_bilin; |
---|
80 | else |
---|
81 | p = p_bilin; |
---|
82 | p.originalModel = p_bilin; |
---|
83 | end |
---|
84 | |
---|
85 | % ************************************************************************* |
---|
86 | % Pre-calc lists of linear/bilinear/nonlinear variables |
---|
87 | % ************************************************************************* |
---|
88 | [p.linears,p.bilinears,p.nonlinears] = compile_nonlinear_table(p); |
---|
89 | |
---|
90 | % ************************************************************************* |
---|
91 | % Select branch variables. We should branch on all variables that are |
---|
92 | % involved in bilinear terms. |
---|
93 | % ************************************************************************* |
---|
94 | p.branch_variables = decide_branch_variables(p); |
---|
95 | |
---|
96 | % ************************************************************************* |
---|
97 | % Extract bounds from model |
---|
98 | % ************************************************************************* |
---|
99 | p = preprocess_bilinear_bounds(p); |
---|
100 | p = preprocess_eval_bounds(p); |
---|
101 | |
---|
102 | % ************************************************************************* |
---|
103 | % Now reduce the branch variables by removing bilinear terms that only have |
---|
104 | % been introduced due to a convex quadratic objective |
---|
105 | % ************************************************************************* |
---|
106 | p = reduce_bilinear_branching_variables(p); |
---|
107 | |
---|
108 | % ************************************************************************* |
---|
109 | % Is the given initial solution feasible, or is the zero-solution feasible? |
---|
110 | % ************************************************************************* |
---|
111 | [p,x_min,upper] = initializesolution(p); |
---|
112 | |
---|
113 | % ************************************************************************* |
---|
114 | % Simple pre-solve loop. The loop is needed for variables defined as w = |
---|
115 | % x*y, x = t*u,y=.. |
---|
116 | % ************************************************************************* |
---|
117 | i = 0; |
---|
118 | goon = 1; |
---|
119 | while goon |
---|
120 | start = [p.lb;p.ub]; |
---|
121 | i = i+1; |
---|
122 | p = updatenonlinearbounds(p); |
---|
123 | p = presolve_bounds_from_equalities(p); |
---|
124 | p = updatenonlinearbounds(p); |
---|
125 | p = updatenonlinearbounds(p); |
---|
126 | p = preprocess_eval_bounds(p); |
---|
127 | i = i+1; |
---|
128 | goon = (norm(start-[p.lb;p.ub],'inf') > 1e-1) & i < 50; |
---|
129 | start = [p.lb;p.ub]; |
---|
130 | end |
---|
131 | |
---|
132 | % ************************************************************************* |
---|
133 | % Default root bound improvement |
---|
134 | % ************************************************************************* |
---|
135 | %p.lb = (p.lb+1e5)-1e5; |
---|
136 | %p.ub = (p.ub+1e5)-1e5; |
---|
137 | close = find(abs(p.lb - p.ub) < 1e-12); |
---|
138 | p.lb(close) = (p.lb(close)+p.ub(close))/2; |
---|
139 | p.ub(close) = p.lb(close); |
---|
140 | p = root_node_tighten(p,upper); |
---|
141 | p = updatenonlinearbounds(p); |
---|
142 | p = preprocess_eval_bounds(p); |
---|
143 | p = presolve_bounds_from_equalities(p); |
---|
144 | p = updatenonlinearbounds(p); |
---|
145 | |
---|
146 | output = yalmip_default_output; |
---|
147 | if p.feasible |
---|
148 | |
---|
149 | % ******************************* |
---|
150 | % Bounded domain? |
---|
151 | % ******************************* |
---|
152 | if ~isempty(p.bilinears) |
---|
153 | involved = unique(p.bilinears(:,2:3)); |
---|
154 | if any(isinf(p.lb(involved))) | any(isinf(p.ub(involved))) |
---|
155 | output = yalmip_default_output; |
---|
156 | output.Primal = x_min; |
---|
157 | output.problem = -6; |
---|
158 | output.infostr = yalmiperror(-6); |
---|
159 | output.solved_nodes = 0; |
---|
160 | return |
---|
161 | end |
---|
162 | end |
---|
163 | |
---|
164 | % ******************************* |
---|
165 | % Save time & memory |
---|
166 | % ******************************* |
---|
167 | p.options.savesolverinput = 0; |
---|
168 | p.options.savesolveroutput = 0; |
---|
169 | p.options.saveduals = 0; |
---|
170 | p.options.dimacs = 0; |
---|
171 | |
---|
172 | % ******************************* |
---|
173 | % RUN BILINEAR BRANCH & BOUND |
---|
174 | % ******************************* |
---|
175 | [x_min,solved_nodes,lower,upper] = branch_and_bound(p,x_min,upper); |
---|
176 | |
---|
177 | % ******************************** |
---|
178 | % CREATE SOLUTION AND DIAGNOSTICS |
---|
179 | % ******************************** |
---|
180 | problem = 0; |
---|
181 | if isinf(upper) |
---|
182 | problem = 1; |
---|
183 | end |
---|
184 | if isinf(lower) & (lower<0) |
---|
185 | problem = 2; |
---|
186 | end |
---|
187 | if isinf(lower) & (lower>0) |
---|
188 | problem = 1; |
---|
189 | end |
---|
190 | if solved_nodes == p.options.bnb.maxiter |
---|
191 | problem = 3; |
---|
192 | end |
---|
193 | else |
---|
194 | problem = 1; |
---|
195 | x_min = repmat(nan,length(p.c),1); |
---|
196 | solved_nodes = 0; |
---|
197 | end |
---|
198 | output = yalmip_default_output; |
---|
199 | output.problem = problem; |
---|
200 | output.solved_nodes = solved_nodes; |
---|
201 | output.Primal = x_min(1:n_in); |
---|
202 | output.infostr = yalmiperror(output.problem,'BNB'); |
---|
203 | output.solvertime = etime(clock,bnbsolvertime); |
---|
204 | |
---|
205 | function model = fixbounds(model); |
---|
206 | polynomials = find(model.variabletype > 0); |
---|
207 | LU = max([abs(model.lb) abs(model.ub)],[],2); |
---|
208 | for i = 1:length(polynomials) |
---|
209 | j = polynomials(i); |
---|
210 | if j<=length(model.lb) |
---|
211 | vars = find(model.monomtable(j,:)); |
---|
212 | bound = 1; |
---|
213 | for k = 1:length(vars) |
---|
214 | bound = bound * LU(vars(k))^model.monomtable(j,vars(k)); |
---|
215 | end |
---|
216 | model.lb(j) = max(model.lb(j),-bound); |
---|
217 | model.ub(j) = min(model.ub(j),bound); |
---|
218 | end |
---|
219 | end |
---|
220 | |
---|
221 | |
---|