[37] | 1 | function [Fdual,objdual,X,t,err] = dualize(F,obj,auto,extlp,extend) |
---|
| 2 | % DUALIZE Create the dual of an SDP given in primal form |
---|
| 3 | % |
---|
| 4 | % [Fd,objd,X,t,err] = dualize(F,obj,auto) |
---|
| 5 | % |
---|
| 6 | % Input |
---|
| 7 | % F : Primal constraint in form AX=b+dt, X>0, t free. |
---|
| 8 | % obj : Primal cost CX+ct |
---|
| 9 | % auto : If set to 0, YALMIP will not automatically handle variables |
---|
| 10 | % and update variable values when the dual problem is solved. |
---|
| 11 | % extlp: If set to 0, YALMIP will not try to perform variables changes in |
---|
| 12 | % order to convert simple translated LP cones (as in x>1) to |
---|
| 13 | % standard unit cone constraints (x>0) |
---|
| 14 | % |
---|
| 15 | % Output |
---|
| 16 | % Fd : Dual constraints in form C-A'y>0, c-dy==0 |
---|
| 17 | % obj : Dual cost b'y (to be MAXIMIZED!) |
---|
| 18 | % X : The detected primal cone variables |
---|
| 19 | % t : The detected primal free variables |
---|
| 20 | % err : Error status (returns 0 if no problems) |
---|
| 21 | % |
---|
| 22 | % Example |
---|
| 23 | % See the HTML help. |
---|
| 24 | % |
---|
| 25 | % See also DUAL, SOLVESDP, PRIMALIZE |
---|
| 26 | |
---|
| 27 | % Author Johan Löfberg |
---|
| 28 | % $Id: dualize.m,v 1.51 2006/10/18 07:59:09 joloef Exp $ |
---|
| 29 | |
---|
| 30 | % Check for unsupported problems |
---|
| 31 | |
---|
| 32 | if isempty(F) |
---|
| 33 | F = set([]); |
---|
| 34 | end |
---|
| 35 | |
---|
| 36 | if nargin < 2 |
---|
| 37 | obj = []; |
---|
| 38 | end |
---|
| 39 | |
---|
| 40 | err = 0; |
---|
| 41 | p1 = ~(isreal(F) & isreal(obj)); |
---|
| 42 | p2 = ~(islinear(F) & islinear(obj)); |
---|
| 43 | p3 = any(is(F,'integer')) | any(is(F,'binary')); |
---|
| 44 | if p1 | p2 | p3 |
---|
| 45 | if nargout == 5 |
---|
| 46 | Fdual = set([]);objdual = [];y = []; X = []; t = []; err = 1; |
---|
| 47 | else |
---|
| 48 | problems = {'Cannot dualize complex-valued problems','Cannot dualize nonlinear problems','Cannot dualize discrete problems'}; |
---|
| 49 | error(problems{min(find([p1 p2 p3]))}); |
---|
| 50 | end |
---|
| 51 | end |
---|
| 52 | |
---|
| 53 | if nargin<5 |
---|
| 54 | extend = 1; |
---|
| 55 | end |
---|
| 56 | |
---|
| 57 | if extend |
---|
| 58 | options.dualize = 1; |
---|
| 59 | options.allowmilp = 0; |
---|
| 60 | options.solver = ''; |
---|
| 61 | [F,failure,cause] = expandmodel(F,obj,options); |
---|
| 62 | if failure |
---|
| 63 | error('Failed during convexity propagation. Avoid nonlinear operators when applying dualization.'); |
---|
| 64 | end |
---|
| 65 | end |
---|
| 66 | |
---|
| 67 | if nargin<3 |
---|
| 68 | auto = 1; |
---|
| 69 | end |
---|
| 70 | |
---|
| 71 | if nargin<4 |
---|
| 72 | extlp = 1; |
---|
| 73 | end |
---|
| 74 | |
---|
| 75 | % Cones and equalities |
---|
| 76 | F_AXb = set([]); |
---|
| 77 | |
---|
| 78 | |
---|
| 79 | % Shiftmatrix is a bit messy at the moment. |
---|
| 80 | % We want to be able to allow cones X>shift |
---|
| 81 | % by using a new variable X-shift = Z |
---|
| 82 | shiftMatrix = {}; |
---|
| 83 | X={}; |
---|
| 84 | |
---|
| 85 | % First, get variables in initial SDP cones |
---|
| 86 | % We need this to avoid adding the same variable twice |
---|
| 87 | % when we add simple LP constraints (as in P>0, P(1,3)>0) |
---|
| 88 | varSDP = []; |
---|
| 89 | SDPset = zeros(length(F),1); |
---|
| 90 | isSDP = is(F,'sdp'); |
---|
| 91 | for i = 1:length(F) |
---|
| 92 | if isSDP(i); |
---|
| 93 | Fi = sdpvar(F(i)); |
---|
| 94 | if is(Fi,'shiftsdpcone') |
---|
| 95 | vars = getvariables(Fi); |
---|
| 96 | if isempty(findrows(varSDP,[vars(1) vars(end)])) |
---|
| 97 | SDPset(i) = 1; |
---|
| 98 | varSDP = [varSDP;vars(1) vars(end)]; |
---|
| 99 | shiftMatrix{end+1} = getbasematrix(Fi,0); |
---|
| 100 | X{end+1}=Fi; |
---|
| 101 | end |
---|
| 102 | end |
---|
| 103 | end |
---|
| 104 | end |
---|
| 105 | F_SDP = F(find(SDPset)); |
---|
| 106 | F = F(find(~SDPset)); |
---|
| 107 | |
---|
| 108 | % Same thing for second order cones |
---|
| 109 | % However, we must not add any SOC cones |
---|
| 110 | % that we already defined as SDP cones |
---|
| 111 | varSOC = []; |
---|
| 112 | SOCset = zeros(length(F),1); |
---|
| 113 | isSOCP = is(F,'socp'); |
---|
| 114 | for i = 1:length(F) |
---|
| 115 | if isSOCP(i);%is(F(i),'socp') |
---|
| 116 | Fi = sdpvar(F(i)); |
---|
| 117 | if is(Fi,'socone') |
---|
| 118 | vars = getvariables(Fi); |
---|
| 119 | % Make sure these variables are not SDP cone variables |
---|
| 120 | % This can actually only happen for set(X>0) + set(Xcone((2:end,1),X(1))) |
---|
| 121 | if ~isempty(varSDP) |
---|
| 122 | inSDP = any(varSDP(:,1)<=vars(1)& vars(1) <=varSDP(:,2)) | any(varSDP(:,1)<=vars(end)& vars(end) <=varSDP(:,2)); |
---|
| 123 | else |
---|
| 124 | inSDP = 0; |
---|
| 125 | end |
---|
| 126 | if ~inSDP |
---|
| 127 | SOCset(i) = 1; |
---|
| 128 | vars = getvariables(Fi); |
---|
| 129 | varSOC = [varSOC;vars(1) vars(end)]; |
---|
| 130 | shiftMatrix{end+1} = getbasematrix(Fi,0); |
---|
| 131 | X{end+1}=Fi; |
---|
| 132 | end |
---|
| 133 | end |
---|
| 134 | end |
---|
| 135 | end |
---|
| 136 | F_SOC = F(find(SOCset)); |
---|
| 137 | F = F(find(~SOCset)); |
---|
| 138 | |
---|
| 139 | % Merge SDP and SOC data |
---|
| 140 | varCONE = [varSDP;varSOC]; |
---|
| 141 | F_CONE = F_SDP + F_SOC; |
---|
| 142 | |
---|
| 143 | % Find all LP constraints, add slacks and extract simple cones |
---|
| 144 | % to speed up things, we treat LP cone somewhat different |
---|
| 145 | % compared to the conceptually similiar SOCP/SDP cones |
---|
| 146 | % This code is pretty messy, since there are a lot off odd |
---|
| 147 | % cases to take care of (x>0 and x>1 etc etc) |
---|
| 148 | elementwise = is(F,'element-wise'); |
---|
| 149 | elementwise_index = find(elementwise); |
---|
| 150 | if ~isempty(elementwise_index) |
---|
| 151 | |
---|
| 152 | % Find element-wise inequalities |
---|
| 153 | Flp = F(elementwise_index); |
---|
| 154 | F = F(find(~elementwise)); % remove these LPs |
---|
| 155 | % Find LP cones |
---|
| 156 | lpconstraint = []; |
---|
| 157 | for i = 1:length(Flp) |
---|
| 158 | temp = sdpvar(Flp(i)); |
---|
| 159 | if min(size(temp))>1 |
---|
| 160 | temp = temp(:); |
---|
| 161 | end |
---|
| 162 | lpconstraint = [lpconstraint reshape(temp,1,length(temp))]; |
---|
| 163 | end |
---|
| 164 | |
---|
| 165 | % Find all constraints of type a_i+x_i >0 and extract the unique and |
---|
| 166 | % most constraining inequalities (i.e. remove redundant lower bounds) |
---|
| 167 | base = getbase(lpconstraint); |
---|
| 168 | candidates = find(sum(base(:,2:end)~=0,2)==1); |
---|
| 169 | if length(candidates)>0 |
---|
| 170 | % The other ones... |
---|
| 171 | alwayskeep = find(sum(base(:,2:end)~=0,2)~=1); |
---|
| 172 | w1 = lpconstraint(alwayskeep); |
---|
| 173 | w2 = lpconstraint(candidates); |
---|
| 174 | |
---|
| 175 | % Find unique rows |
---|
| 176 | base = getbase(w2); |
---|
| 177 | [i,uniquerows,k] = unique(base(:,2:end)*randn(size(base,2)-1,1)); |
---|
| 178 | aUniqueRow=k(:)'; |
---|
| 179 | keep = []; |
---|
| 180 | rhsLP = base(:,1); |
---|
| 181 | rr = histc(k,[1:max(k)]); |
---|
| 182 | if all(rr==1) |
---|
| 183 | lpconstraint = [w1 w2]; |
---|
| 184 | else |
---|
| 185 | for i=1:length(k) |
---|
| 186 | sameRow=find(k==k(i)); |
---|
| 187 | if length(sameRow)==1 |
---|
| 188 | keep = [keep sameRow]; |
---|
| 189 | else |
---|
| 190 | rhs=base(sameRow,1); |
---|
| 191 | [val,pos]=min(rhsLP(sameRow)); |
---|
| 192 | keep = [keep sameRow(pos)]; |
---|
| 193 | end |
---|
| 194 | end |
---|
| 195 | lpconstraint = [w1 w2(unique(keep))]; |
---|
| 196 | end |
---|
| 197 | |
---|
| 198 | end |
---|
| 199 | |
---|
| 200 | % LP cone will be saved in a vector for speed |
---|
| 201 | x = []; |
---|
| 202 | |
---|
| 203 | % Pure cones of the type x>0 |
---|
| 204 | base = getbase(lpconstraint); |
---|
| 205 | purelpcones = (base(:,1)==0) & (sum(base(:,2:end),2)==1) & (sum(base(:,2:end)==1,2)==1); |
---|
| 206 | if ~isempty(purelpcones) |
---|
| 207 | x = [x lpconstraint(purelpcones)]; |
---|
| 208 | lpconstraint = lpconstraint(find(~purelpcones)); |
---|
| 209 | end |
---|
| 210 | |
---|
| 211 | |
---|
| 212 | % Translated cones x>k, k positive |
---|
| 213 | % User does not want to make variable changes based on k |
---|
| 214 | % But if k>=0, we can at-least mark x as a simple LP cone variable and |
---|
| 215 | % thus avoid a free variable. |
---|
| 216 | if ~extlp & ~isempty(lpconstraint) |
---|
| 217 | base = getbase(lpconstraint); |
---|
| 218 | lpcones = (base(:,1)<0) & (sum(base(:,2:end),2)==1) & (sum(base(:,2:end)==1,2)==1); |
---|
| 219 | if ~isempty(find(lpcones)) |
---|
| 220 | s = recover(getvariables(lpconstraint(find(lpcones)))); |
---|
| 221 | x = [x reshape(s,1,length(s))]; |
---|
| 222 | end |
---|
| 223 | end |
---|
| 224 | |
---|
| 225 | % Translated cones x>k |
---|
| 226 | % Extract these and perform the associated variable change y=x-k |
---|
| 227 | if ~isempty(lpconstraint)%Avoid warning in 5.3.1 |
---|
| 228 | base = getbase(lpconstraint); |
---|
| 229 | lpcones = (sum(base(:,2:end),2)==1) & (sum(base(:,2:end)==1,2)==1); |
---|
| 230 | if ~isempty(lpcones) & extlp |
---|
| 231 | x = [x lpconstraint(find(lpcones))]; |
---|
| 232 | nlp = lpconstraint(find(~lpcones)); |
---|
| 233 | if ~isempty(nlp) |
---|
| 234 | s = sdpvar(1,length(nlp)); |
---|
| 235 | F_AXb = F_AXb + set(nlp-s==0); |
---|
| 236 | x = [x s]; |
---|
| 237 | end |
---|
| 238 | elseif length(lpconstraint) > 0 |
---|
| 239 | s = sdpvar(1,length(lpconstraint)); |
---|
| 240 | x = [x s]; % New LP cones |
---|
| 241 | F_AXb = F_AXb + set(lpconstraint-s==0); |
---|
| 242 | end |
---|
| 243 | end |
---|
| 244 | |
---|
| 245 | |
---|
| 246 | % Sort asccording to variable index |
---|
| 247 | % (Code below assumes x is sorted in increasing variables indicies) |
---|
| 248 | base = getbase(x);base = base(:,2:end);[i,j,k] = find(base); |
---|
| 249 | x = x(i); |
---|
| 250 | xv = getvariables(x); |
---|
| 251 | |
---|
| 252 | % For mixed LP/SDP problems, we must ensure that LP cone variables are |
---|
| 253 | % not actually an element in a SDP cone variable |
---|
| 254 | if ~isempty(varCONE) |
---|
| 255 | keep = zeros(length(x),1); |
---|
| 256 | for i = 1:length(xv) |
---|
| 257 | if any(varCONE(:,1)<=xv(i) & xv(i) <=varCONE(:,2)) |
---|
| 258 | else |
---|
| 259 | keep(i) = 1; |
---|
| 260 | end |
---|
| 261 | end |
---|
| 262 | if ~all(keep) |
---|
| 263 | % We need to add some explicit constraints on some elements and |
---|
| 264 | % remove the x variables since they are already in a cone |
---|
| 265 | % variable |
---|
| 266 | xcone = x(find(~keep)); |
---|
| 267 | s = sdpvar(1,length(xcone)); |
---|
| 268 | F_AXb = F_AXb + set(xcone-s==0); |
---|
| 269 | x = x(find(keep)); |
---|
| 270 | x = [x s]; |
---|
| 271 | end |
---|
| 272 | end |
---|
| 273 | else |
---|
| 274 | x = []; |
---|
| 275 | end |
---|
| 276 | |
---|
| 277 | % Check for mixed cones, ie terms C-A'y > 0. |
---|
| 278 | keep = ones(length(F),1); |
---|
| 279 | isSDP = is(F,'sdp'); |
---|
| 280 | isSOCP = is(F,'socp'); |
---|
| 281 | |
---|
| 282 | % Pre-allocate all SDP slacks in one call |
---|
| 283 | % This is a lot faster |
---|
| 284 | if nnz(isSDP) > 0 |
---|
| 285 | SDPindicies = find(isSDP)'; |
---|
| 286 | for i = 1:length(SDPindicies)%find(isSDP)' |
---|
| 287 | Fi = sdpvar(F(SDPindicies(i))); |
---|
| 288 | ns(i) = size(Fi,1); |
---|
| 289 | ms(i) = ns(i); |
---|
| 290 | end |
---|
| 291 | Slacks = sdpvar(ns,ms); |
---|
| 292 | if ~isa(Slacks,'cell') |
---|
| 293 | Slacks = {Slacks}; |
---|
| 294 | end |
---|
| 295 | end |
---|
| 296 | prei = 1; |
---|
| 297 | for i = 1:length(F) |
---|
| 298 | if isSDP(i) |
---|
| 299 | % Semidefinite dual-form cone |
---|
| 300 | Fi = sdpvar(F(i)); |
---|
| 301 | n = size(Fi,1); |
---|
| 302 | % S = sdpvar(n,n); |
---|
| 303 | S = Slacks{prei};prei = prei + 1; |
---|
| 304 | slack = Fi-S; |
---|
| 305 | ind = find(triu(reshape(1:n^2,n,n))); |
---|
| 306 | F_AXb = F_AXb + set(slack(ind)==0); |
---|
| 307 | F_CONE = F_CONE + lmi(S,[],[],[],1); |
---|
| 308 | shiftMatrix{end+1} = spalloc(n,n,0); |
---|
| 309 | X{end+1}=S; |
---|
| 310 | keep(i)=0; |
---|
| 311 | elseif isSOCP(i) |
---|
| 312 | % SOC dual-form cone |
---|
| 313 | Fi = sdpvar(F(i)); |
---|
| 314 | n = size(Fi,1); |
---|
| 315 | S = sdpvar(n,1); |
---|
| 316 | % S = Slacks{i}; |
---|
| 317 | slack = Fi-S; |
---|
| 318 | F_AXb = F_AXb + set(slack==0); |
---|
| 319 | F_CONE = F_CONE + set(cone(S(2:end),S(1))); |
---|
| 320 | shiftMatrix{end+1} = spalloc(n,1,0); |
---|
| 321 | X{end+1}=S; |
---|
| 322 | keep(i)=0; |
---|
| 323 | end |
---|
| 324 | end |
---|
| 325 | |
---|
| 326 | % Now, remove all mixed cones... |
---|
| 327 | F = F(find(keep)); |
---|
| 328 | |
---|
| 329 | % Get the equalities |
---|
| 330 | AXbset = is(F,'equality'); |
---|
| 331 | if any(AXbset) |
---|
| 332 | % Get the constraints |
---|
| 333 | F_AXb = F_AXb + F(find(AXbset)); |
---|
| 334 | F = F(find(~AXbset)); |
---|
| 335 | end |
---|
| 336 | |
---|
| 337 | % Is thee something we missed in our tests? |
---|
| 338 | if length(F)>0 |
---|
| 339 | error('DUALIZE can only treat standard SDPs (and LPs) at the moment.') |
---|
| 340 | end |
---|
| 341 | |
---|
| 342 | |
---|
| 343 | % Sort the SDP cone variables X according to YALMIP |
---|
| 344 | % This is just to simplify some indexing later |
---|
| 345 | ns = []; |
---|
| 346 | first_var = []; |
---|
| 347 | for i = 1:length(F_CONE) |
---|
| 348 | ns = [ns length(X{i})]; |
---|
| 349 | first_var = [first_var min(getvariables(X{i}))]; |
---|
| 350 | end |
---|
| 351 | [sorted,index] = sort(first_var); |
---|
| 352 | X={X{index}}; |
---|
| 353 | shiftMatrix={shiftMatrix{index}}; |
---|
| 354 | |
---|
| 355 | shift = []; |
---|
| 356 | for i = 1:length(F_CONE) |
---|
| 357 | ns(i) = length(X{i}); |
---|
| 358 | if size(X{i},2)==1 |
---|
| 359 | % SOCP |
---|
| 360 | shift = [shift;shiftMatrix{i}(:)]; |
---|
| 361 | else |
---|
| 362 | % SDP |
---|
| 363 | ind = find(tril(reshape(1:ns(i)^2,ns(i),ns(i)))); |
---|
| 364 | shift = [shift;shiftMatrix{i}(ind)]; |
---|
| 365 | end |
---|
| 366 | end |
---|
| 367 | |
---|
| 368 | % Free variables (here called t) is everything except the cone variables |
---|
| 369 | X_variables = getvariables(F_CONE); |
---|
| 370 | x_variables = getvariables(x); |
---|
| 371 | Xx_variables = [X_variables x_variables]; |
---|
| 372 | |
---|
| 373 | other_variables = [getvariables(obj) getvariables(F_AXb)]; |
---|
| 374 | % For quadratic case |
---|
| 375 | %other_variables = [depends(obj) getvariables(F_AXb)]; |
---|
| 376 | |
---|
| 377 | all_variables = uniquestripped([other_variables Xx_variables]); |
---|
| 378 | |
---|
| 379 | % Avoid set-diff |
---|
| 380 | if isequal(all_variables,Xx_variables) |
---|
| 381 | t_variables = []; |
---|
| 382 | t_in_all = []; |
---|
| 383 | t = []; |
---|
| 384 | else |
---|
| 385 | t_variables = setdiff(all_variables,Xx_variables); |
---|
| 386 | ind = ismembc(all_variables,t_variables); |
---|
| 387 | t_in_all = find(ind); |
---|
| 388 | t = recover(t_variables); |
---|
| 389 | end |
---|
| 390 | |
---|
| 391 | ind = ismembc(all_variables,x_variables); |
---|
| 392 | x_in_all = find(ind); |
---|
| 393 | ind = ismembc(all_variables,X_variables); |
---|
| 394 | X_in_all = find(ind); |
---|
| 395 | |
---|
| 396 | vecF1 = []; |
---|
| 397 | nvars = length(all_variables); |
---|
| 398 | for i = 1:length(F_AXb) |
---|
| 399 | AXb = sdpvar(F_AXb(i)); |
---|
| 400 | mapper = find(ismembc(all_variables,getvariables(AXb))); |
---|
| 401 | |
---|
| 402 | [n,m] = size(AXb); |
---|
| 403 | data = getbase(AXb); |
---|
| 404 | [iF,jF,sF] = find(data); |
---|
| 405 | if 1 % New |
---|
| 406 | smapper = [1 1+mapper(:)']; |
---|
| 407 | F_structemp = sparse(iF,smapper(jF),sF,n*m,1+nvars); |
---|
| 408 | else |
---|
| 409 | F_structemp = spalloc(n*m,1+nvars,nnz(data)); |
---|
| 410 | F_structemp(:,[1 1+mapper(:)'])= data; |
---|
| 411 | end |
---|
| 412 | vecF1 = [vecF1;F_structemp]; |
---|
| 413 | end |
---|
| 414 | |
---|
| 415 | %Remove trivially redundant constraints |
---|
| 416 | h = 1+rand(size(vecF1,2),1); |
---|
| 417 | h = vecF1*h; |
---|
| 418 | [dummy,uniquerows] = uniquesafe(h); |
---|
| 419 | if length(uniquerows) < length(h) |
---|
| 420 | % Sort to ensure run-to-run consistency |
---|
| 421 | vecF1 = vecF1((sort(uniquerows)),:); |
---|
| 422 | end |
---|
| 423 | |
---|
| 424 | if isempty(obj) |
---|
| 425 | vecF1(end+1,1) = 0; |
---|
| 426 | else |
---|
| 427 | if is(obj,'linear') |
---|
| 428 | mapper = find(ismembc(all_variables,getvariables(obj))); |
---|
| 429 | [n,m] = size(obj); |
---|
| 430 | data = getbase(obj); |
---|
| 431 | [iF,jF,sF] = find(data); |
---|
| 432 | if 1 |
---|
| 433 | smapper = [1 1+mapper(:)']; |
---|
| 434 | F_structemp = sparse(iF,smapper(jF),sF,n*m,1+nvars); |
---|
| 435 | else |
---|
| 436 | F_structemp = spalloc(n*m,1+nvars,nnz(data)); |
---|
| 437 | F_structemp(:,[1 1+mapper(:)'])= data; |
---|
| 438 | end |
---|
| 439 | vecF1 = [vecF1;F_structemp]; |
---|
| 440 | else |
---|
| 441 | % FIX: Generalize to QP duality |
---|
| 442 | % min c'x+0.5x'Qx, Ax==b, x>=0 |
---|
| 443 | % max b'y-0.5x'Qx, c-A'y+Qx >=0 |
---|
| 444 | [Q,c,xreally,info] = quaddecomp(obj,recover(all_variables)) |
---|
| 445 | mapper = find(ismembc(all_variables,getvariables(c'*xreally))); |
---|
| 446 | [n,m] = size(c'*xreally); |
---|
| 447 | data = getbase(c'*xreally); |
---|
| 448 | F_structemp = spalloc(n*m,1+nvars,nnz(data)); |
---|
| 449 | F_structemp(:,[1 1+mapper(:)'])= data; |
---|
| 450 | vecF1 = [vecF1;F_structemp] |
---|
| 451 | end |
---|
| 452 | |
---|
| 453 | end |
---|
| 454 | |
---|
| 455 | vecF1(end+1,1) = 0; |
---|
| 456 | Fbase = vecF1; |
---|
| 457 | |
---|
| 458 | %Fbase = unique(Fbase','rows')'; |
---|
| 459 | |
---|
| 460 | b = Fbase(1:end-2,1); |
---|
| 461 | |
---|
| 462 | Fbase = -Fbase(1:end-1,2:end); |
---|
| 463 | vecA = []; |
---|
| 464 | Fbase_t = Fbase(:,t_in_all); |
---|
| 465 | Fbase_x = Fbase(:,x_in_all); |
---|
| 466 | Fbase_X = Fbase; |
---|
| 467 | %Fbase_X(:,unionstripped(t_in_all,x_in_all)) = []; |
---|
| 468 | if 1 |
---|
| 469 | removethese = unique([t_in_all x_in_all]); |
---|
| 470 | if length(removethese) > 0.5*size(Fbase_X,2) |
---|
| 471 | Fbase_X = Fbase_X(:,setdiff(1:size(Fbase_X,2),removethese)); |
---|
| 472 | else |
---|
| 473 | Fbase_X(:,[t_in_all x_in_all]) = []; |
---|
| 474 | end |
---|
| 475 | else |
---|
| 476 | removecols = uniquestripped([t_in_all x_in_all]); |
---|
| 477 | if ~isempty(removecols) |
---|
| 478 | [i,j,k] = find(Fbase_X); |
---|
| 479 | keep = find(~ismember(j,removecols)); |
---|
| 480 | i = i(keep); |
---|
| 481 | k = k(keep); |
---|
| 482 | j = j(keep); |
---|
| 483 | map = find(1:length(unique(j)),j); |
---|
| 484 | |
---|
| 485 | end |
---|
| 486 | end |
---|
| 487 | |
---|
| 488 | % Shift due to translated dual cones X = Z+shift |
---|
| 489 | if length(shift > 0) |
---|
| 490 | b = b + Fbase_X(1:end-1,:)*shift; |
---|
| 491 | end |
---|
| 492 | if length(x)>0 |
---|
| 493 | % Arrgh |
---|
| 494 | base = getbase(x); |
---|
| 495 | constant = base(:,1); |
---|
| 496 | base = base(:,2:end);[i,j,k] = find(base); |
---|
| 497 | b = b + Fbase_x(1:end-1,:)*constant(i); |
---|
| 498 | end |
---|
| 499 | |
---|
| 500 | start = 0; |
---|
| 501 | n_cones = length(ns); |
---|
| 502 | % All LPs in one cone |
---|
| 503 | if length(x)>0 |
---|
| 504 | n_lp = 1; |
---|
| 505 | else |
---|
| 506 | n_lp = 0; |
---|
| 507 | end |
---|
| 508 | n_free = length(t_variables); |
---|
| 509 | |
---|
| 510 | % SDP cones |
---|
| 511 | for j = 1:1:n_cones |
---|
| 512 | |
---|
| 513 | if size(X{j},1)==size(X{j},2) |
---|
| 514 | % SDP cone |
---|
| 515 | ind = reshape(1:ns(j)^2,ns(j),ns(j)); |
---|
| 516 | ind = find(tril(ind)); |
---|
| 517 | |
---|
| 518 | % Get non-symmetric constraint AiX=b |
---|
| 519 | Fi = Fbase_X(1:length(b),start+(1:length(ind)))'/2; |
---|
| 520 | |
---|
| 521 | if 1 |
---|
| 522 | [iF,jF,sF] = find(Fi); |
---|
| 523 | iA = ind(iF); |
---|
| 524 | jA = jF; |
---|
| 525 | sA = sF; |
---|
| 526 | the_col = 1+floor((iA-1)/ns(j)); |
---|
| 527 | the_row = iA-(the_col-1)*ns(j); |
---|
| 528 | these_must_be_transposed = find(the_row > the_col); |
---|
| 529 | if ~isempty(these_must_be_transposed) |
---|
| 530 | new_rowcols = the_col(these_must_be_transposed) + (the_row(these_must_be_transposed)-1)*ns(j); |
---|
| 531 | iA = [iA;new_rowcols]; |
---|
| 532 | jA = [jA;jA(these_must_be_transposed)]; |
---|
| 533 | sA = [sA;sA(these_must_be_transposed)]; |
---|
| 534 | end |
---|
| 535 | % Fix diagonal term |
---|
| 536 | diags = find(diag(1:ns(j))); |
---|
| 537 | id = find(ismembc(iA,diags)); |
---|
| 538 | sA(id) = 2*sA(id); |
---|
| 539 | Ai = sparse(iA,jA,sA,ns(j)^2,length(b)); |
---|
| 540 | |
---|
| 541 | else % Old slooooooow version |
---|
| 542 | Ai = spalloc(ns(j)^2,length(b),nnz(Fi)); |
---|
| 543 | Ai(ind,:) = Fi; |
---|
| 544 | % Symmetrize |
---|
| 545 | [rowcols,varindicies,vals]=find(Ai); |
---|
| 546 | the_col = 1+floor((rowcols-1)/ns(j)); |
---|
| 547 | the_row = rowcols-(the_col-1)*ns(j); |
---|
| 548 | these_must_be_transposed = find(the_row > the_col); |
---|
| 549 | if ~isempty(these_must_be_transposed) |
---|
| 550 | new_rowcols = the_col(these_must_be_transposed) + (the_row(these_must_be_transposed)-1)*ns(j); |
---|
| 551 | Ai(sub2ind(size(Ai),new_rowcols,varindicies(these_must_be_transposed))) = vals(these_must_be_transposed); |
---|
| 552 | end |
---|
| 553 | |
---|
| 554 | % Fix diagonal term |
---|
| 555 | diags = find(diag(1:ns(j))); |
---|
| 556 | Ai(diags,:) = 2*Ai(diags,:); |
---|
| 557 | end |
---|
| 558 | |
---|
| 559 | % if norm(Ai-Ai2,inf)>1e-12 |
---|
| 560 | % error |
---|
| 561 | % end |
---|
| 562 | |
---|
| 563 | |
---|
| 564 | vecA{j} = Ai; |
---|
| 565 | |
---|
| 566 | start = start + length(ind); |
---|
| 567 | else |
---|
| 568 | % Second order cone |
---|
| 569 | ind = 1:ns(j); |
---|
| 570 | |
---|
| 571 | % Get constraint AiX=b |
---|
| 572 | Fi = Fbase_X(1:length(b),start+(1:length(ind)))'; |
---|
| 573 | Ai = spalloc(ns(j),length(b),nnz(Fi)); |
---|
| 574 | Ai(ind,:) = Fi; |
---|
| 575 | vecA{j} = Ai; |
---|
| 576 | start = start + length(ind); |
---|
| 577 | end |
---|
| 578 | |
---|
| 579 | end |
---|
| 580 | % LP Cone |
---|
| 581 | if n_lp>0 |
---|
| 582 | Alp=Fbase_x(1:length(b),:)'; |
---|
| 583 | end |
---|
| 584 | % FREE VARIABLES |
---|
| 585 | start = 0; |
---|
| 586 | if n_free>0 |
---|
| 587 | Afree=Fbase_t(1:length(b),:)'; |
---|
| 588 | end |
---|
| 589 | |
---|
| 590 | % COST MATRIX |
---|
| 591 | % SDP CONE |
---|
| 592 | start = 0; |
---|
| 593 | for j = 1:1:n_cones |
---|
| 594 | if size(X{j},1)==size(X{j},2) |
---|
| 595 | %ind = reshape(1:ns(j)^2,ns(j),ns(j)); |
---|
| 596 | %ind = find(tril(ind)); |
---|
| 597 | %C{j} = spalloc(ns(j),ns(j),0); |
---|
| 598 | %C{j}(ind) = -Fbase_X(end,start+(1:length(ind))); |
---|
| 599 | %C{j} = (C{j}+ C{j}')/2; |
---|
| 600 | %start = start + length(ind); |
---|
| 601 | |
---|
| 602 | ind = reshape(1:ns(j)^2,ns(j),ns(j)); |
---|
| 603 | [indi,indj] = find(tril(ind)); |
---|
| 604 | C{j} = sparse(indi,indj,-Fbase_X(end,start+(1:length(indi))),ns(j),ns(j)); |
---|
| 605 | C{j} = (C{j}+ C{j}')/2; |
---|
| 606 | start = start + length(indi); |
---|
| 607 | else |
---|
| 608 | ind = 1:ns(j); |
---|
| 609 | C{j} = spalloc(ns(j),1,0); |
---|
| 610 | C{j}(ind) = -Fbase_X(end,start+(1:length(ind))); |
---|
| 611 | start = start + length(ind); |
---|
| 612 | end |
---|
| 613 | end |
---|
| 614 | % LP CONE |
---|
| 615 | for j = 1:1:n_lp |
---|
| 616 | Clp = -Fbase_x(end,:)'; |
---|
| 617 | end |
---|
| 618 | % FREE CONE |
---|
| 619 | if n_free>0 |
---|
| 620 | Cfree = -Fbase_t(end,1:end)'; |
---|
| 621 | end |
---|
| 622 | |
---|
| 623 | % Create dual |
---|
| 624 | y = sdpvar(length(b),1); |
---|
| 625 | yvars = getvariables(y); |
---|
| 626 | cf = []; |
---|
| 627 | Af = []; |
---|
| 628 | Fdual = set([]); |
---|
| 629 | for j = 1:n_cones |
---|
| 630 | if size(X{j},1)==size(X{j},2) |
---|
| 631 | % Yep, this is optimized... |
---|
| 632 | S = sdpvar(ns(j),ns(j),[],yvars,[C{j}(:) -vecA{j}]); |
---|
| 633 | % Fast call avoids symmetry check |
---|
| 634 | Fdual = Fdual + lmi(S,[],[],[],1); |
---|
| 635 | else |
---|
| 636 | Ay = reshape(vecA{j}*y,ns(j),1); |
---|
| 637 | S = C{j}-Ay; |
---|
| 638 | Fdual = Fdual + lmi(cone(S(2:end),S(1))); |
---|
| 639 | end |
---|
| 640 | end |
---|
| 641 | if n_lp > 0 |
---|
| 642 | keep = any(Alp,2); |
---|
| 643 | if ~all(keep) |
---|
| 644 | % Fix for unused primal cones |
---|
| 645 | preset=find(~keep); |
---|
| 646 | xfix = x(preset); |
---|
| 647 | assign(xfix(:),Clp(preset(:))); |
---|
| 648 | end |
---|
| 649 | keep = find(keep); |
---|
| 650 | if ~isempty(keep) |
---|
| 651 | z = Clp(keep)-Alp(keep,:)*y; |
---|
| 652 | if isa(z,'double') |
---|
| 653 | assign(x(:),z(:)); |
---|
| 654 | else |
---|
| 655 | Fdual = Fdual + set(z); |
---|
| 656 | x =x(keep); |
---|
| 657 | X{end+1} = x(:); % We have worked with a row for performance reasons |
---|
| 658 | end |
---|
| 659 | end |
---|
| 660 | end |
---|
| 661 | if n_free > 0 |
---|
| 662 | CfreeAfreey = Cfree-Afree*y; |
---|
| 663 | if isa(CfreeAfreey,'double') |
---|
| 664 | if nnz(CfreeAfreey)>0 |
---|
| 665 | error('Trivially unbounded!'); |
---|
| 666 | end |
---|
| 667 | else |
---|
| 668 | Fdual = Fdual + set(0 == CfreeAfreey); |
---|
| 669 | end |
---|
| 670 | end |
---|
| 671 | |
---|
| 672 | objdual = b'*y; |
---|
| 673 | if auto |
---|
| 674 | for i = 1:length(X) |
---|
| 675 | yalmip('associatedual',getlmiid(Fdual(i)),X{i}); |
---|
| 676 | end |
---|
| 677 | if n_free>0 |
---|
| 678 | yalmip('associatedual',getlmiid(Fdual(end)),t); |
---|
| 679 | end |
---|
| 680 | end |
---|