[37] | 1 | function output = callpenbmi(interfacedata); |
---|
| 2 | |
---|
| 3 | % Author Johan Löfberg |
---|
| 4 | % $Id: callpenbmim.m,v 1.21 2006/08/24 17:12:27 joloef Exp $ |
---|
| 5 | |
---|
| 6 | if any(interfacedata.variabletype > 2) |
---|
| 7 | % Polynomial problem, YALMIP has to bilienarize |
---|
| 8 | output = callpenbmi_with_bilinearization(interfacedata); |
---|
| 9 | else |
---|
| 10 | % Old standard code |
---|
| 11 | output = callpenbmi_without_bilinearization(interfacedata); |
---|
| 12 | end |
---|
| 13 | |
---|
| 14 | function output = callpenbmi_without_bilinearization(interfacedata); |
---|
| 15 | |
---|
| 16 | % Retrieve needed data |
---|
| 17 | clear penbmim |
---|
| 18 | |
---|
| 19 | options = interfacedata.options; |
---|
| 20 | F_struc = interfacedata.F_struc; |
---|
| 21 | c = interfacedata.c; |
---|
| 22 | Q = interfacedata.Q; |
---|
| 23 | K = interfacedata.K; |
---|
| 24 | x0 = interfacedata.x0; |
---|
| 25 | monomtable = interfacedata.monomtable; |
---|
| 26 | ub = interfacedata.ub; |
---|
| 27 | lb = interfacedata.lb; |
---|
| 28 | |
---|
| 29 | % Linear before bilinearize |
---|
| 30 | temp = sum(monomtable,2)>1; |
---|
| 31 | tempnonlinearindicies = find(temp); |
---|
| 32 | templinearindicies = find(~temp); |
---|
| 33 | |
---|
| 34 | % Any stupid constant>0 constraints |
---|
| 35 | % FIX : Recover duals afterwards |
---|
| 36 | % Better fix : Do this outside |
---|
| 37 | zrow = []; |
---|
| 38 | if K.l >0 |
---|
| 39 | zrow = find(any(F_struc(1:K.l+K.f,:),2)==0); |
---|
| 40 | if ~isempty(zrow) |
---|
| 41 | K.l = K.l - nnz(zrow>K.f); |
---|
| 42 | K.f = K.f - nnz(zrow<=K.f); |
---|
| 43 | F_struc(zrow,:) = []; |
---|
| 44 | end |
---|
| 45 | end |
---|
| 46 | |
---|
| 47 | % Bounded variables converted to constraints |
---|
| 48 | if ~isempty(ub) |
---|
| 49 | [F_struc,K] = addbounds(F_struc,K,ub,lb); |
---|
| 50 | end |
---|
| 51 | |
---|
| 52 | % This one only occurs if called from bmibnb |
---|
| 53 | if K.f>0 |
---|
| 54 | F_struc = [-F_struc(1:K.f,:);F_struc]; |
---|
| 55 | F_struc(1:K.f,1) = F_struc(1:K.f,1)+sqrt(eps); |
---|
| 56 | K.l = K.l + 2*K.f; |
---|
| 57 | K.f = 0; |
---|
| 58 | end |
---|
| 59 | |
---|
| 60 | if isempty(monomtable) |
---|
| 61 | monomtable = eye(length(c)); |
---|
| 62 | end |
---|
| 63 | |
---|
| 64 | temp = sum(monomtable,2)>1; |
---|
| 65 | nonlinearindicies = find(temp); |
---|
| 66 | linearindicies = find(~temp); |
---|
| 67 | c0 = c; |
---|
| 68 | c = c(linearindicies); |
---|
| 69 | Q = Q(linearindicies,linearindicies); |
---|
| 70 | nonlinear_scalars = []; |
---|
| 71 | |
---|
| 72 | % Any non-linear scalar inequalities? |
---|
| 73 | % Move these to the BMI part |
---|
| 74 | if K.l>0 |
---|
| 75 | nonlinear_scalars = find(any(full(F_struc(1:K.l,[nonlinearindicies(:)'+1])),2)); |
---|
| 76 | if ~isempty(nonlinear_scalars) |
---|
| 77 | Kold = K; |
---|
| 78 | % SETDIFF DONE FASTER |
---|
| 79 | aa = 1:K.l; |
---|
| 80 | bb = nonlinear_scalars; |
---|
| 81 | tf = ~(ismembc(aa,bb)); |
---|
| 82 | cc = aa(tf); |
---|
| 83 | cc = unique(cc); |
---|
| 84 | linear_scalars = cc; |
---|
| 85 | % linear_scalars = setdiff1(1:K.l,nonlinear_scalars); |
---|
| 86 | F_struc = [F_struc(linear_scalars,:);F_struc(nonlinear_scalars,:);F_struc(K.l+1:end,:)]; |
---|
| 87 | K.l = K.l-length(nonlinear_scalars); |
---|
| 88 | if (length(K.s)==1) & (K.s==0) |
---|
| 89 | K.s = [repmat(1,1,length(nonlinear_scalars))]; |
---|
| 90 | K.rank = repmat(1,1,length(nonlinear_scalars)); |
---|
| 91 | else |
---|
| 92 | K.s = [repmat(1,1,length(nonlinear_scalars)) K.s]; |
---|
| 93 | K.rank = [repmat(1,1,length(nonlinear_scalars)) K.rank]; |
---|
| 94 | end |
---|
| 95 | end |
---|
| 96 | end |
---|
| 97 | |
---|
| 98 | if ~isempty(F_struc) |
---|
| 99 | pen = sedumi2pen(F_struc(:,[1 linearindicies(:)'+1]),K,c,x0); |
---|
| 100 | else |
---|
| 101 | pen = sedumi2pen([],K,c,x0); |
---|
| 102 | end |
---|
| 103 | |
---|
| 104 | if ~isempty(nonlinearindicies) |
---|
| 105 | bmi = sedumi2pen(F_struc(:,[nonlinearindicies(:)'+1]),K,[],[]); |
---|
| 106 | pen.ki_dim = bmi.ai_dim; |
---|
| 107 | % Nonlinear index |
---|
| 108 | pen.ki_dim = bmi.ai_dim; |
---|
| 109 | pen.ki_row = bmi.ai_row; |
---|
| 110 | pen.ki_col = bmi.ai_col; |
---|
| 111 | pen.ki_nzs = bmi.ai_nzs; |
---|
| 112 | pen.ki_val = bmi.ai_val; |
---|
| 113 | if 0 |
---|
| 114 | for i = 1:length(bmi.ai_idx) |
---|
| 115 | nl = nonlinearindicies(1+bmi.ai_idx(i)); |
---|
| 116 | v = find(monomtable(nl,:)); |
---|
| 117 | if length(v)==1 |
---|
| 118 | v(2)=v(1); |
---|
| 119 | end |
---|
| 120 | pen.ki_idx(i)=find(linearindicies == v(1)); |
---|
| 121 | pen.kj_idx(i)=find(linearindicies == v(2)); |
---|
| 122 | end |
---|
| 123 | else |
---|
| 124 | top = 1; |
---|
| 125 | [ii,jj,kk] = find(monomtable(nonlinearindicies(1+bmi.ai_idx),:)'); |
---|
| 126 | pen.ki_idx = zeros(1,length(bmi.ai_idx)); |
---|
| 127 | pen.kj_idx = zeros(1,length(bmi.ai_idx)); |
---|
| 128 | for i = 1:length(bmi.ai_idx) |
---|
| 129 | if kk(top)==2 |
---|
| 130 | v(1) = ii(top); |
---|
| 131 | v(2) = ii(top); |
---|
| 132 | top = top + 1; |
---|
| 133 | else |
---|
| 134 | v(1) = ii(top); |
---|
| 135 | v(2) = ii(top+1); |
---|
| 136 | top = top + 2; |
---|
| 137 | end |
---|
| 138 | % FIX : precompute this map |
---|
| 139 | pen.ki_idx(i)=find(linearindicies == v(1)); |
---|
| 140 | pen.kj_idx(i)=find(linearindicies == v(2)); |
---|
| 141 | end |
---|
| 142 | end |
---|
| 143 | |
---|
| 144 | else |
---|
| 145 | pen.ki_dim = 0*pen.ai_dim; |
---|
| 146 | pen.ki_row = 0; |
---|
| 147 | pen.ki_col = 0; |
---|
| 148 | pen.ki_nzs = 0; |
---|
| 149 | pen.ki_idx = 0; |
---|
| 150 | pen.kj_idx = 0; |
---|
| 151 | pen.kj_val = 0; |
---|
| 152 | end |
---|
| 153 | |
---|
| 154 | if nnz(Q)>0 |
---|
| 155 | [row,col,vals] = find(triu(Q)); |
---|
| 156 | pen.q_nzs = length(row); |
---|
| 157 | pen.q_val = vals'; |
---|
| 158 | pen.q_col = col'-1; |
---|
| 159 | pen.q_row = row'-1; |
---|
| 160 | else |
---|
| 161 | pen.q_nzs = 0; |
---|
| 162 | pen.q_val = 0; |
---|
| 163 | pen.q_col = 0; |
---|
| 164 | pen.q_row = 0; |
---|
| 165 | end |
---|
| 166 | |
---|
| 167 | ops = struct2cell(options.penbmi);ops = [ops{1:end}]; |
---|
| 168 | pen.ioptions = ops(1:12); |
---|
| 169 | pen.foptions = ops(13:end); |
---|
| 170 | pen.ioptions(4) = max(0,min(3,options.verbose+1)); |
---|
| 171 | if pen.ioptions(4)==1 |
---|
| 172 | pen.ioptions(4)=0; |
---|
| 173 | end |
---|
| 174 | |
---|
| 175 | % **************************************** |
---|
| 176 | % UNCOMMENT THIS FOR PENBMI version 1 |
---|
| 177 | % **************************************** |
---|
| 178 | %pen.ioptions = pen.ioptions(1:8); |
---|
| 179 | %pen.foptions = pen.foptions(1:8); |
---|
| 180 | |
---|
| 181 | if ~isempty(x0) |
---|
| 182 | pen.x0(isnan(pen.x0)) = 0; |
---|
| 183 | pen.x0 = x0(linearindicies); |
---|
| 184 | pen.x0 = pen.x0(:)'; |
---|
| 185 | end |
---|
| 186 | |
---|
| 187 | if options.savedebug |
---|
| 188 | save penbmimdebug pen |
---|
| 189 | end |
---|
| 190 | |
---|
| 191 | if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end |
---|
| 192 | |
---|
| 193 | solvertime = clock; |
---|
| 194 | [f,xout,u,iflag,niter,feas] = penbmim(pen); |
---|
| 195 | solvertime = etime(clock,solvertime); |
---|
| 196 | |
---|
| 197 | if options.saveduals & isempty(zrow) |
---|
| 198 | |
---|
| 199 | % Get dual variable |
---|
| 200 | % First, get the nonlinear scalars treated as BMIs |
---|
| 201 | if ~isempty(nonlinear_scalars) |
---|
| 202 | n_orig_scalars = length(nonlinear_scalars)+K.l; |
---|
| 203 | linear_scalars = setdiff(1:n_orig_scalars,nonlinear_scalars); |
---|
| 204 | u_nonlinear=u(K.l+1:K.l+length(nonlinear_scalars)); |
---|
| 205 | u(K.l+1:K.l+length(nonlinear_scalars))=[]; |
---|
| 206 | u_linear = u(1:K.l); |
---|
| 207 | u_scalar = zeros(1,n_orig_scalars); |
---|
| 208 | u_scalar(linear_scalars)=u_linear; |
---|
| 209 | u_scalar(nonlinear_scalars)=u_nonlinear; |
---|
| 210 | u = [u_scalar u(1+K.l:end)]; |
---|
| 211 | K = Kold; |
---|
| 212 | end |
---|
| 213 | |
---|
| 214 | u = u(:); |
---|
| 215 | D_struc = u(1:1:K.l); |
---|
| 216 | if length(K.s)>0 |
---|
| 217 | if K.s(1)>0 |
---|
| 218 | pos = K.l+1; |
---|
| 219 | for i = 1:length(K.s) |
---|
| 220 | temp = zeros(K.s(i),K.s(i)); |
---|
| 221 | vecZ = u(pos:pos+0.5*K.s(i)*(K.s(i)+1)-1); |
---|
| 222 | top = 1; |
---|
| 223 | for j = 1:K.s(i) |
---|
| 224 | len = K.s(i)-j+1; |
---|
| 225 | temp(j:end,j)=vecZ(top:top+len-1);top=top+len; |
---|
| 226 | end |
---|
| 227 | temp = (temp+temp');j = find(speye(K.s(i)));temp(j)=temp(j)/2; |
---|
| 228 | D_struc = [D_struc;temp(:)]; |
---|
| 229 | pos = pos + (K.s(i)+1)*K.s(i)/2; |
---|
| 230 | end |
---|
| 231 | end |
---|
| 232 | end |
---|
| 233 | else |
---|
| 234 | D_struc = []; |
---|
| 235 | end |
---|
| 236 | |
---|
| 237 | %Recover solution |
---|
| 238 | if isempty(nonlinearindicies) |
---|
| 239 | x = xout(:); |
---|
| 240 | else |
---|
| 241 | x = zeros(length(interfacedata.c),1); |
---|
| 242 | for i = 1:length(templinearindicies) |
---|
| 243 | x(templinearindicies(i)) = xout(i); |
---|
| 244 | end |
---|
| 245 | end |
---|
| 246 | |
---|
| 247 | problem = 0; |
---|
| 248 | switch iflag |
---|
| 249 | case 0 |
---|
| 250 | problem = 0; % OK |
---|
| 251 | case {1,3} |
---|
| 252 | problem = 4; |
---|
| 253 | case 2 |
---|
| 254 | problem = 1; % INFEASIBLE |
---|
| 255 | case 4 |
---|
| 256 | problem = 3; % Numerics |
---|
| 257 | case 5 |
---|
| 258 | problem = 7; |
---|
| 259 | case {6,7} |
---|
| 260 | problem = 11; |
---|
| 261 | otherwise |
---|
| 262 | problem = -1; |
---|
| 263 | end |
---|
| 264 | infostr = yalmiperror(problem,interfacedata.solver.tag); |
---|
| 265 | |
---|
| 266 | if options.savesolveroutput |
---|
| 267 | solveroutput.f = f; |
---|
| 268 | solveroutput.xout = xout; |
---|
| 269 | solveroutput.u = u; |
---|
| 270 | solveroutput.iflag = iflag; |
---|
| 271 | solveroutput.niter = niter; |
---|
| 272 | solveroutput.feas = feas; |
---|
| 273 | else |
---|
| 274 | solveroutput = []; |
---|
| 275 | end |
---|
| 276 | |
---|
| 277 | if options.savesolverinput |
---|
| 278 | solverinput.pen = pen; |
---|
| 279 | else |
---|
| 280 | solverinput = []; |
---|
| 281 | end |
---|
| 282 | |
---|
| 283 | % Standard interface |
---|
| 284 | output.Primal = x(:); |
---|
| 285 | output.Dual = D_struc; |
---|
| 286 | output.Slack = []; |
---|
| 287 | output.problem = problem; |
---|
| 288 | output.infostr = infostr; |
---|
| 289 | output.solverinput = solverinput; |
---|
| 290 | output.solveroutput= solveroutput; |
---|
| 291 | output.solvertime = solvertime; |
---|
| 292 | |
---|
| 293 | function output = callpenbmi_with_bilinearization(interfacedata); |
---|
| 294 | |
---|
| 295 | % Bilinearize |
---|
| 296 | [p,changed] = bilinearize(interfacedata); |
---|
| 297 | |
---|
| 298 | % Convert bilinearizing equalities to inequalities |
---|
| 299 | if p.K.f>0 |
---|
| 300 | p.F_struc = [-p.F_struc(1:p.K.f,:);p.F_struc]; |
---|
| 301 | p.F_struc(1:p.K.f,1) = p.F_struc(1:p.K.f,1)+sqrt(eps); |
---|
| 302 | p.K.l = p.K.l + 2*p.K.f; |
---|
| 303 | p.K.f = 0; |
---|
| 304 | end |
---|
| 305 | |
---|
| 306 | % Solve bilinearized problem |
---|
| 307 | output = callpenbmi_without_bilinearization(p); |
---|
| 308 | |
---|
| 309 | % Get our original variables & duals |
---|
| 310 | output.Primal = output.Primal(1:length(interfacedata.c)); |
---|
| 311 | if ~isempty(output.Dual) |
---|
| 312 | n_equ = p.K.f - interfacedata.K.f; |
---|
| 313 | % First 2*n_eq are the duals for the new inequalities |
---|
| 314 | output.Dual = output.Dual(1+2*n_equ:end); |
---|
| 315 | end |
---|
| 316 | |
---|
| 317 | |
---|
| 318 | |
---|
| 319 | |
---|
| 320 | |
---|