[37] | 1 | function output = cutsdp(p) |
---|
| 2 | % CUTSDP |
---|
| 3 | % |
---|
| 4 | % See also SOLVESDP, BNB, BINVAR, INTVAR, BINARY, INTEGER, LMI |
---|
| 5 | |
---|
| 6 | % Author Johan Löfberg |
---|
| 7 | % $Id: cutsdp.m,v 1.5 2006/05/23 21:55:59 joloef Exp $ |
---|
| 8 | |
---|
| 9 | % ************************************************************************* |
---|
| 10 | %% INITIALIZE DIAGNOSTICS IN YALMIP |
---|
| 11 | % ************************************************************************* |
---|
| 12 | bnbsolvertime = clock; |
---|
| 13 | showprogress('Cutting plane solver started',p.options.showprogress); |
---|
| 14 | |
---|
| 15 | % ************************************************************************* |
---|
| 16 | %% If we want duals, we may not extract bounds. However, bounds must be |
---|
| 17 | % extracted in discrete problems. |
---|
| 18 | % ************************************************************************* |
---|
| 19 | if p.options.cutsdp.recoverdual |
---|
| 20 | warning('Dual recovery not implemented yet in CUTSDP') |
---|
| 21 | end |
---|
| 22 | |
---|
| 23 | % ************************************************************************* |
---|
| 24 | %% Define infinite bounds |
---|
| 25 | % ************************************************************************* |
---|
| 26 | if isempty(p.ub) |
---|
| 27 | p.ub = repmat(inf,length(p.c),1); |
---|
| 28 | end |
---|
| 29 | if isempty(p.lb) |
---|
| 30 | p.lb = repmat(-inf,length(p.c),1); |
---|
| 31 | end |
---|
| 32 | |
---|
| 33 | % ************************************************************************* |
---|
| 34 | %% ADD CONSTRAINTS 0<x<1 FOR BINARY |
---|
| 35 | % ************************************************************************* |
---|
| 36 | if ~isempty(p.binary_variables) |
---|
| 37 | p.ub(p.binary_variables) = min(p.ub(p.binary_variables),1); |
---|
| 38 | p.lb(p.binary_variables) = max(p.lb(p.binary_variables),0); |
---|
| 39 | end |
---|
| 40 | |
---|
| 41 | % ************************************************************************* |
---|
| 42 | %% Extract better bounds from model |
---|
| 43 | % ************************************************************************* |
---|
| 44 | if ~isempty(p.F_struc) |
---|
| 45 | [lb,ub,used_rows] = findulb(p.F_struc,p.K); |
---|
| 46 | if ~isempty(used_rows) |
---|
| 47 | lower_defined = find(~isinf(lb)); |
---|
| 48 | if ~isempty(lower_defined) |
---|
| 49 | p.lb(lower_defined) = max(p.lb(lower_defined),lb(lower_defined)); |
---|
| 50 | end |
---|
| 51 | upper_defined = find(~isinf(ub)); |
---|
| 52 | if ~isempty(upper_defined) |
---|
| 53 | p.ub(upper_defined) = min(p.ub(upper_defined),ub(upper_defined)); |
---|
| 54 | end |
---|
| 55 | p.F_struc(p.K.f+used_rows,:)=[]; |
---|
| 56 | p.K.l = p.K.l - length(used_rows); |
---|
| 57 | end |
---|
| 58 | end |
---|
| 59 | |
---|
| 60 | % ************************************************************************* |
---|
| 61 | %% ADD CONSTRAINTS 0<x<1 FOR BINARY |
---|
| 62 | % ************************************************************************* |
---|
| 63 | if ~isempty(p.binary_variables) |
---|
| 64 | p.ub(p.binary_variables) = min(p.ub(p.binary_variables),1); |
---|
| 65 | p.lb(p.binary_variables) = max(p.lb(p.binary_variables),0); |
---|
| 66 | end |
---|
| 67 | |
---|
| 68 | p.ub = min(p.ub,p.options.cutsdp.variablebound'); |
---|
| 69 | p.lb = max(p.lb,-p.options.cutsdp.variablebound'); |
---|
| 70 | |
---|
| 71 | % ************************************************************************* |
---|
| 72 | %% PRE-SOLVE (nothing fancy coded) |
---|
| 73 | % ************************************************************************* |
---|
| 74 | if isempty(find(isinf([p.ub;p.lb]))) & p.K.l>0 |
---|
| 75 | [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,p.integer_variables); |
---|
| 76 | end |
---|
| 77 | |
---|
| 78 | % ************************************************************************* |
---|
| 79 | %% PERTURBATION OF LINEAR COST |
---|
| 80 | % ************************************************************************* |
---|
| 81 | p.corig = p.c; |
---|
| 82 | if nnz(p.Q)~=0 |
---|
| 83 | g = randn('seed'); |
---|
| 84 | randn('state',1253); %For my testing, I keep this the same... |
---|
| 85 | % This perturbation has to be better. Crucial for many real LP problems |
---|
| 86 | p.c = (p.c).*(1+randn(length(p.c),1)*1e-4); |
---|
| 87 | randn('seed',g); |
---|
| 88 | end |
---|
| 89 | |
---|
| 90 | % ************************************************************************* |
---|
| 91 | %% We don't need this |
---|
| 92 | % ************************************************************************* |
---|
| 93 | p.options.savesolverinput = 0; |
---|
| 94 | p.options.savesolveroutput = 0; |
---|
| 95 | |
---|
| 96 | % ************************************************************************* |
---|
| 97 | %% Display logics |
---|
| 98 | % 0 : Silent |
---|
| 99 | % 1 : Display cut progress |
---|
| 100 | % 2 : Display node solver prints |
---|
| 101 | % ************************************************************************* |
---|
| 102 | switch max(min(p.options.verbose,3),0) |
---|
| 103 | case 0 |
---|
| 104 | p.options.cutsdp.verbose = 0; |
---|
| 105 | case 1 |
---|
| 106 | p.options.cutsdp.verbose = 1; |
---|
| 107 | p.options.verbose = 0; |
---|
| 108 | case 2 |
---|
| 109 | p.options.cutsdp.verbose = 2; |
---|
| 110 | p.options.verbose = 0; |
---|
| 111 | case 3 |
---|
| 112 | p.options.cutsdp.verbose = 2; |
---|
| 113 | p.options.verbose = 1; |
---|
| 114 | otherwise |
---|
| 115 | p.options.cutsdp.verbose = 0; |
---|
| 116 | p.options.verbose = 0; |
---|
| 117 | end |
---|
| 118 | |
---|
| 119 | % ************************************************************************* |
---|
| 120 | %% START CUTTING |
---|
| 121 | % ************************************************************************* |
---|
| 122 | [x_min,solved_nodes,lower,feasible,D_struc] = cutting(p); |
---|
| 123 | %% -- |
---|
| 124 | |
---|
| 125 | % ************************************************************************* |
---|
| 126 | %% CREATE SOLUTION |
---|
| 127 | % ************************************************************************* |
---|
| 128 | output.problem = 0; |
---|
| 129 | if ~feasible |
---|
| 130 | output.problem = 1; |
---|
| 131 | end |
---|
| 132 | if solved_nodes == p.options.cutsdp.maxiter |
---|
| 133 | output.problem = 3; |
---|
| 134 | end |
---|
| 135 | output.solved_nodes = solved_nodes; |
---|
| 136 | output.Primal = x_min; |
---|
| 137 | output.Dual = D_struc; |
---|
| 138 | output.Slack = []; |
---|
| 139 | output.solverinput = 0; |
---|
| 140 | output.solveroutput =[]; |
---|
| 141 | output.solvertime = etime(clock,bnbsolvertime); |
---|
| 142 | %% -- |
---|
| 143 | |
---|
| 144 | function [x,solved_nodes,lower,feasible,D_struc] = cutting(p) |
---|
| 145 | |
---|
| 146 | % ************************************************************************* |
---|
| 147 | %% Sanity check |
---|
| 148 | % ************************************************************************* |
---|
| 149 | if any(p.lb>p.ub) |
---|
| 150 | x = zeros(length(p.c),1); |
---|
| 151 | solved_nodes = 0; |
---|
| 152 | lower = inf; |
---|
| 153 | feasible = 0; |
---|
| 154 | D_struc = []; |
---|
| 155 | return |
---|
| 156 | end |
---|
| 157 | |
---|
| 158 | % ************************************************************************* |
---|
| 159 | %% Create function handle to solver |
---|
| 160 | % ************************************************************************* |
---|
| 161 | cutsolver = p.solver.cutsolver.call; |
---|
| 162 | |
---|
| 163 | % ************************************************************************* |
---|
| 164 | %% Create copy of model without |
---|
| 165 | % the SDP part |
---|
| 166 | % ************************************************************************* |
---|
| 167 | p_lp = p; |
---|
| 168 | p_lp.F_struc = p_lp.F_struc(1:p.K.l+p.K.f,:); |
---|
| 169 | p_lp.K.s = 0; |
---|
| 170 | |
---|
| 171 | % ************************************************************************* |
---|
| 172 | %% DISPLAY HEADER |
---|
| 173 | % ************************************************************************* |
---|
| 174 | if p.options.cutsdp.verbose |
---|
| 175 | disp('* Starting YALMIP cutting plane for MISDP based on MILP'); |
---|
| 176 | disp(['* Lower solver : ' p.solver.cutsolver.tag]); |
---|
| 177 | disp(['* Max iterations : ' num2str(p.options.cutsdp.maxiter)]); |
---|
| 178 | end |
---|
| 179 | |
---|
| 180 | if p.options.bnb.verbose; disp(' Node Infeasibility. Lower LP cuts');end; |
---|
| 181 | |
---|
| 182 | %% Initialize diagnostic |
---|
| 183 | infeasibility = -inf; |
---|
| 184 | solved_nodes = 0; |
---|
| 185 | feasible = 1; |
---|
| 186 | lower = -inf; |
---|
| 187 | saveduals = 1; |
---|
| 188 | |
---|
| 189 | % ************************************************************************* |
---|
| 190 | %% Add diagonal cuts to begin with |
---|
| 191 | % ************************************************************************* |
---|
| 192 | savedCuts = []; |
---|
| 193 | savedIndicies = []; |
---|
| 194 | if p.K.s(1)>0 |
---|
| 195 | top = p.K.f+p.K.l+1; |
---|
| 196 | for i = 1:length(p.K.s) |
---|
| 197 | n = p.K.s(i); |
---|
| 198 | newF=[]; |
---|
| 199 | for m = 1:p.K.s(i) |
---|
| 200 | d = eyev(p.K.s(i),m); |
---|
| 201 | index = (1+(m-1)*(p.K.s(i)+1)); |
---|
| 202 | newF = [newF;p.F_struc(top+index-1,:)]; |
---|
| 203 | end |
---|
| 204 | % Clean |
---|
| 205 | newF(abs(newF)<1e-12) = 0; |
---|
| 206 | keep=find(any(newF(:,2:end),2)); |
---|
| 207 | newF = newF(keep,:); |
---|
| 208 | |
---|
| 209 | p_lp.F_struc = [p_lp.F_struc;newF]; |
---|
| 210 | p_lp.K.l = p_lp.K.l + size(newF,1); |
---|
| 211 | top = top+n^2; |
---|
| 212 | end |
---|
| 213 | end |
---|
| 214 | |
---|
| 215 | goon = 1; |
---|
| 216 | rr = p_lp.integer_variables; |
---|
| 217 | rb = p_lp.binary_variables; |
---|
| 218 | only_solvelp = 0; |
---|
| 219 | pplot = 0; |
---|
| 220 | |
---|
| 221 | % ************************************************************************* |
---|
| 222 | % Crude lower bound |
---|
| 223 | % FIX for quadratic case |
---|
| 224 | % ************************************************************************* |
---|
| 225 | lower = 0; |
---|
| 226 | if nnz(p.Q) == 0 |
---|
| 227 | for i = 1:length(p.c) |
---|
| 228 | if p.c(i)>0 |
---|
| 229 | if isinf(p.lb(i)) |
---|
| 230 | lower = -inf; |
---|
| 231 | break |
---|
| 232 | else |
---|
| 233 | lower = lower + p.c(i)*p.lb(i); |
---|
| 234 | end |
---|
| 235 | elseif p.c(i)<0 |
---|
| 236 | if isinf(p.ub(i)) |
---|
| 237 | lower = -inf; |
---|
| 238 | break |
---|
| 239 | else |
---|
| 240 | lower = lower + p.c(i)*p.ub(i); |
---|
| 241 | end |
---|
| 242 | end |
---|
| 243 | end |
---|
| 244 | end |
---|
| 245 | %lower = sum(sign(p.c).*(p.lb)); |
---|
| 246 | if isinf(lower) | nnz(p.Q)~=0 |
---|
| 247 | lower = -1e6; |
---|
| 248 | end |
---|
| 249 | |
---|
| 250 | % ************************************************************************* |
---|
| 251 | % Experimental stuff for variable fixing |
---|
| 252 | % ************************************************************************* |
---|
| 253 | if p.options.cutsdp.nodefix & (p.K.s(1)>0) |
---|
| 254 | top=1+p.K.f+p.K.l+sum(p.K.q); |
---|
| 255 | for i=1:length(p.K.s) |
---|
| 256 | n=p.K.s(i); |
---|
| 257 | for j=1:size(p.F_struc,2)-1; |
---|
| 258 | X=full(reshape(p.F_struc(top:top+n^2-1,j+1),p.K.s(i),p.K.s(i))); |
---|
| 259 | X=(X+X')/2; |
---|
| 260 | v=real(eig(X+sqrt(eps)*eye(length(X)))); |
---|
| 261 | if all(v>=0) |
---|
| 262 | sdpmonotinicity(i,j)=-1; |
---|
| 263 | elseif all(v<=0) |
---|
| 264 | sdpmonotinicity(i,j)=1; |
---|
| 265 | else |
---|
| 266 | sdpmonotinicity(i,j)=nan; |
---|
| 267 | end |
---|
| 268 | end |
---|
| 269 | top=top+n^2; |
---|
| 270 | end |
---|
| 271 | else |
---|
| 272 | sdpmonotinicity=[]; |
---|
| 273 | end |
---|
| 274 | |
---|
| 275 | while goon |
---|
| 276 | |
---|
| 277 | % Add lower bound |
---|
| 278 | if ~isinf(lower) |
---|
| 279 | p_lp.F_struc = [p_lp.F_struc;-lower p_lp.c']; |
---|
| 280 | p_lp.K.l = p_lp.K.l + 1; |
---|
| 281 | end |
---|
| 282 | |
---|
| 283 | if p.options.cutsdp.nodetight |
---|
| 284 | % Extract LP part Ax<=b |
---|
| 285 | A = -p_lp.F_struc(p_lp.K.f + (1:p_lp.K.l),2:end); |
---|
| 286 | b = p_lp.F_struc(p_lp.K.f + (1:p_lp.K.l),1); |
---|
| 287 | c = p_lp.c; |
---|
| 288 | % Tighten bounds and find redundant constraints |
---|
| 289 | [p_lp.lb,p_lp.ub,redundant,pss] = milppresolve(A,b,p_lp.lb,p_lp.ub,p.integer_variables,p.binary_variables,ones(length(p.lb),1)); |
---|
| 290 | A(redundant,:) = []; |
---|
| 291 | b(redundant) = []; |
---|
| 292 | p_lp.F_struc(p_lp.K.f+redundant,:) = []; |
---|
| 293 | p_lp.K.l = p_lp.K.l-length(redundant); |
---|
| 294 | end |
---|
| 295 | |
---|
| 296 | if p.options.cutsdp.nodefix |
---|
| 297 | % Try to find variables to fix w.l.o.g |
---|
| 298 | [fix_up,fix_down] = presolve_fixvariables(A,b,c,p_lp.lb,p_lp.ub,sdpmonotinicity); |
---|
| 299 | p_lp.lb(fix_up) = p_lp.ub(fix_up); |
---|
| 300 | p_lp.ub(fix_down) = p_lp.lb(fix_down); |
---|
| 301 | while ~(isempty(fix_up) & isempty(fix_down)) |
---|
| 302 | [p_lp.lb,p_lp.ub,redundant,pss] = milppresolve(A,b,p_lp.lb,p_lp.ub,p.integer_variables,p.binary_variables,ones(length(p.lb),1)); |
---|
| 303 | A(redundant,:) = []; |
---|
| 304 | b(redundant) = []; |
---|
| 305 | p_lp.F_struc(p_lp.K.f+redundant,:) = []; |
---|
| 306 | p_lp.K.l = p_lp.K.l-length(redundant); |
---|
| 307 | fix_up = []; |
---|
| 308 | fix_down = []; |
---|
| 309 | % Try to find variables to fix w.l.o.g |
---|
| 310 | [fix_up,fix_down] = presolve_fixvariables(A,b,c,p_lp.lb,p_lp.ub,sdpmonotinicity); |
---|
| 311 | p_lp.lb(fix_up) = p_lp.ub(fix_up); |
---|
| 312 | p_lp.ub(fix_down) = p_lp.lb(fix_down); |
---|
| 313 | end |
---|
| 314 | end |
---|
| 315 | |
---|
| 316 | |
---|
| 317 | output = feval(cutsolver,p_lp); |
---|
| 318 | |
---|
| 319 | % Remove lower bound (avoid accumulating them) |
---|
| 320 | if ~isinf(lower) |
---|
| 321 | p_lp.K.l = p_lp.K.l - 1; |
---|
| 322 | p_lp.F_struc = p_lp.F_struc(1:end-1,:); |
---|
| 323 | end |
---|
| 324 | |
---|
| 325 | if output.problem == 1 | output.problem == 12 |
---|
| 326 | % LP relaxation was infeasible, hence problem is infeasible |
---|
| 327 | feasible = 0; |
---|
| 328 | lower = inf; |
---|
| 329 | goon = 0; |
---|
| 330 | x = zeros(length(p.c),1); |
---|
| 331 | lower = inf; |
---|
| 332 | else |
---|
| 333 | % Relaxed solution |
---|
| 334 | x = output.Primal; |
---|
| 335 | lower = p.f+p.c'*x+x'*p.Q*x; |
---|
| 336 | |
---|
| 337 | infeasibility = 0; |
---|
| 338 | if p.K.s(1)>0 |
---|
| 339 | % Add cuts |
---|
| 340 | top = p.K.f+p.K.l+1; |
---|
| 341 | for i = 1:1:length(p.K.s) |
---|
| 342 | n = p.K.s(i); |
---|
| 343 | X = p.F_struc(top:top+n^2-1,:)*[1;x]; |
---|
| 344 | X = reshape(X,n,n); |
---|
| 345 | Y = randn(n,n); |
---|
| 346 | newcuts = 1; |
---|
| 347 | newF = zeros(n,size(p.F_struc,2)); |
---|
| 348 | [d,v] = eig(X); |
---|
| 349 | infeasibility = min(infeasibility,min(diag(v))); |
---|
| 350 | dummy=[]; |
---|
| 351 | newF = []; |
---|
| 352 | if infeasibility<0 |
---|
| 353 | [ii,jj] = sort(diag(v)); |
---|
| 354 | for m = jj(1:min(length(jj),p.options.cutsdp.cutlimit))'%find(diag(v<0))%1:1%length(v) |
---|
| 355 | if v(m,m)<0 |
---|
| 356 | bA = d(:,m)'*(kron(d(:,m),speye(n))'*p.F_struc(top:top+n^2-1,:)); |
---|
| 357 | b = bA(:,1); |
---|
| 358 | A = -bA(:,2:end); |
---|
| 359 | newF = [newF;b -A]; |
---|
| 360 | %f = b-floor(b); |
---|
| 361 | %fj = A - floor(A); |
---|
| 362 | %A = floor(A) + max((fj-f),0)./(1-f); |
---|
| 363 | %b = floor(b); |
---|
| 364 | %newF = [newF;b -A]; |
---|
| 365 | newcuts = newcuts + 1; |
---|
| 366 | end |
---|
| 367 | end |
---|
| 368 | end |
---|
| 369 | newF(abs(newF)<1e-12) = 0; |
---|
| 370 | keep=find(any(newF(:,2:end),2)); |
---|
| 371 | newF = newF(keep,:); |
---|
| 372 | if size(newF,1)>0 |
---|
| 373 | p_lp.F_struc = [p_lp.F_struc;newF]; |
---|
| 374 | p_lp.K.l = p_lp.K.l + size(newF,1); |
---|
| 375 | [i,j] = sort(p_lp.F_struc*[1;x]); |
---|
| 376 | end |
---|
| 377 | top = top+n^2; |
---|
| 378 | end |
---|
| 379 | end |
---|
| 380 | % res = find(p_lp.F_struc*[1;x] > mean(abs(p_lp.F_struc*[1;x]))); |
---|
| 381 | % p_lp.F_struc(res,:) = []; |
---|
| 382 | % p_lp.K.l = p_lp.K.l-length(res); |
---|
| 383 | goon = infeasibility <= p.options.cutsdp.feastol; |
---|
| 384 | goon = goon & feasible; |
---|
| 385 | goon = goon & (solved_nodes < p.options.cutsdp.maxiter-1); |
---|
| 386 | end |
---|
| 387 | |
---|
| 388 | solved_nodes = solved_nodes + 1; |
---|
| 389 | if p.options.cutsdp.verbose;fprintf(' %4.0f : %12.3E %12.3E %2.0f\n',solved_nodes,infeasibility,lower,p_lp.K.l-p.K.l);end |
---|
| 390 | end |
---|
| 391 | |
---|
| 392 | D_struc = []; |
---|