[37] | 1 | function [interfacedata,recoverdata,solver,diagnostic,F,Fremoved] = compileinterfacedata(F,aux_obsolete,P,h,options,findallsolvers,parametric) |
---|
| 2 | |
---|
| 3 | persistent CACHED_SOLVERS |
---|
| 4 | persistent EXISTTIME |
---|
| 5 | persistent NCHECKS |
---|
| 6 | |
---|
| 7 | diagnostic = []; |
---|
| 8 | interfacedata = []; |
---|
| 9 | recoverdata = []; |
---|
| 10 | solver = []; |
---|
| 11 | Fremoved = []; |
---|
| 12 | |
---|
| 13 | if nargin<7 |
---|
| 14 | parametric = 0; |
---|
| 15 | end |
---|
| 16 | |
---|
| 17 | if isequal(h,0) |
---|
| 18 | h = []; |
---|
| 19 | end |
---|
| 20 | |
---|
| 21 | % ************************************************************************* |
---|
| 22 | % EXTRACT LOW-RANK DESCRIPTION |
---|
| 23 | % ************************************************************************* |
---|
| 24 | lowrankdetails = getlrdata(F); |
---|
| 25 | if ~isempty(lowrankdetails) |
---|
| 26 | F = F(~is(F,'lowrank')); |
---|
| 27 | end |
---|
| 28 | |
---|
| 29 | % ************************************************************************* |
---|
| 30 | % PERTURB STRICT INEQULAITIES |
---|
| 31 | % ************************************************************************* |
---|
| 32 | if isa(options.shift,'sdpvar') | (options.shift~=0) |
---|
| 33 | F = shift(F,options.shift); |
---|
| 34 | end |
---|
| 35 | |
---|
| 36 | % ************************************************************************* |
---|
| 37 | % ADD RADIUS CONSTRAINT |
---|
| 38 | % ************************************************************************* |
---|
| 39 | if isa(options.radius,'sdpvar') | ~isinf(options.radius) |
---|
| 40 | x = recover(unique(union(depends(h),depends(F)))); |
---|
| 41 | if length(x)>1 |
---|
| 42 | F = F + set(cone(x,options.radius)); |
---|
| 43 | else |
---|
| 44 | F = F + set(-options.radius < x < options.radius); |
---|
| 45 | end |
---|
| 46 | end |
---|
| 47 | |
---|
| 48 | % ************************************************************************* |
---|
| 49 | % CONVERT LOGIC CONSTRAINTS |
---|
| 50 | % ************************************************************************* |
---|
| 51 | [F,changed] = convertlogics(F); |
---|
| 52 | if changed |
---|
| 53 | options.saveduals = 0; % Don't calculate duals since we changed the problem |
---|
| 54 | end |
---|
| 55 | |
---|
| 56 | % ************************************************************************* |
---|
| 57 | % Take care of the nonlinear operators by converting expressions such as |
---|
| 58 | % t = max(x,y) to standard conic models and mixed integer models |
---|
| 59 | % This part also adds parts from logical expressions and mpower terms |
---|
| 60 | % ************************************************************************* |
---|
| 61 | if options.expand |
---|
| 62 | % Experimental hack due to support for the PWQ function used for |
---|
| 63 | % quadratic dynamic programming with MPT. |
---|
| 64 | % FIX: Clean up and generalize |
---|
| 65 | try |
---|
| 66 | variables = uniquestripped([depends(h) getvariables(h)]); |
---|
| 67 | extendedvariables = yalmip('extvariables'); |
---|
| 68 | index_in_extended = find(ismembc(variables,extendedvariables)); |
---|
| 69 | if ~isempty(index_in_extended) |
---|
| 70 | extstruct = yalmip('extstruct',variables(index_in_extended)); |
---|
| 71 | if isequal(extstruct.fcn ,'pwq_yalmip') |
---|
| 72 | [Fz,properties,arguments]=model(extstruct.var,'milp',options,extstruct); |
---|
| 73 | gain = getbasematrix(h,getvariables(extstruct.var)); |
---|
| 74 | h = replace(h,extstruct.var,0); |
---|
| 75 | h = h + gain*properties.replacer; |
---|
| 76 | F = F + Fz; |
---|
| 77 | end |
---|
| 78 | end |
---|
| 79 | catch |
---|
| 80 | end |
---|
| 81 | [F,failure,cause] = expandmodel(F,h,options); |
---|
| 82 | if failure % Convexity propgation failed |
---|
| 83 | interfacedata = []; |
---|
| 84 | recoverdata = []; |
---|
| 85 | solver = ''; |
---|
| 86 | diagnostic.solvertime = 0; |
---|
| 87 | diagnostic.problem = 14; |
---|
| 88 | diagnostic.info = yalmiperror(14,cause); |
---|
| 89 | return |
---|
| 90 | end |
---|
| 91 | evalVariables = yalmip('evalVariables'); |
---|
| 92 | if isempty(evalVariables) |
---|
| 93 | evaluation_based = 0; |
---|
| 94 | else |
---|
| 95 | evaluation_based = ~isempty(intersect([getvariables(h) getvariables(F)],evalVariables)); |
---|
| 96 | end |
---|
| 97 | else |
---|
| 98 | evaluation_based = 0; |
---|
| 99 | end |
---|
| 100 | |
---|
| 101 | % ************************************************************************* |
---|
| 102 | % CONVERT CONVEX QUADRATIC CONSTRAINTS |
---|
| 103 | % We do not convert quadratic constraints to SOCPs if we have have |
---|
| 104 | % sigmonial terms (thus indicating a GP problem), if we have relaxed |
---|
| 105 | % nonlinear expressions, or if we have specified a nonlinear solver. |
---|
| 106 | % Why do we convert them already here? Don't remember, should be cleaned up |
---|
| 107 | % ************************************************************************* |
---|
| 108 | [monomtable,variabletype] = yalmip('monomtable'); |
---|
| 109 | % if any(variabletype(getvariables(F))==4) |
---|
| 110 | % do_not_convert = 1; |
---|
| 111 | % else |
---|
| 112 | % do_not_convert = 0; |
---|
| 113 | % end |
---|
| 114 | do_not_convert = any(variabletype(getvariables(F))==4); |
---|
| 115 | do_not_convert = do_not_convert | strcmpi(options.solver,'pennlp') | strcmpi(options.solver,'penbmi'); |
---|
| 116 | do_not_convert = do_not_convert | strcmpi(options.solver,'fmincon') | strcmpi(options.solver,'lindo'); |
---|
| 117 | do_not_convert = do_not_convert | strcmpi(options.solver,'fmincon-geometric') | strcmpi(options.solver,'fmincon-standard'); |
---|
| 118 | do_not_convert = do_not_convert | strcmpi(options.solver,'bmibnb'); |
---|
| 119 | do_not_convert = do_not_convert | (options.convertconvexquad == 0); |
---|
| 120 | do_not_convert = do_not_convert | (options.relax == 1); |
---|
| 121 | if ~do_not_convert |
---|
| 122 | [F,socp_changed] = convertquadratics(F); |
---|
| 123 | if socp_changed % changed holds the number of QC -> SOCC conversions |
---|
| 124 | options.saveduals = 0; % We cannot calculate duals since we changed the problem |
---|
| 125 | end |
---|
| 126 | else |
---|
| 127 | socp_changed = 0; |
---|
| 128 | end |
---|
| 129 | |
---|
| 130 | % CHEAT FOR QC |
---|
| 131 | if socp_changed>0 & length(find(is(F,'socc')))==socp_changed |
---|
| 132 | socp_are_really_qc = 1; |
---|
| 133 | else |
---|
| 134 | socp_are_really_qc = 0; |
---|
| 135 | end |
---|
| 136 | |
---|
| 137 | % ************************************************************************* |
---|
| 138 | % LOOK FOR AVAILABLE SOLVERS |
---|
| 139 | % Finding solvers can be very slow on some systems. To alleviate this |
---|
| 140 | % problem, YALMIP can cache the list of avaialble solvers. |
---|
| 141 | % ************************************************************************* |
---|
| 142 | if (options.cachesolvers==0) | isempty(CACHED_SOLVERS) |
---|
| 143 | getsolvertime = clock; |
---|
| 144 | solvers = getavailablesolvers(findallsolvers,options); |
---|
| 145 | getsolvertime = etime(clock,getsolvertime); |
---|
| 146 | % CODE TO INFORM USERS ABOUT SLOW NETWORKS! |
---|
| 147 | if isempty(EXISTTIME) |
---|
| 148 | EXISTTIME = getsolvertime; |
---|
| 149 | NCHECKS = 1; |
---|
| 150 | else |
---|
| 151 | EXISTTIME = [EXISTTIME getsolvertime]; |
---|
| 152 | NCHECKS = NCHECKS + 1; |
---|
| 153 | end |
---|
| 154 | if (options.cachesolvers==0) |
---|
| 155 | if ((NCHECKS >= 3 & (sum(EXISTTIME)/NCHECKS > 1)) | EXISTTIME(end)>2) |
---|
| 156 | if warningon |
---|
| 157 | info = 'Warning: YALMIP has detected that your drive or network is unusually slow.\nThis causes a severe delay in SOLVESDP when I try to find available solvers.\nTo avoid this, use the options CACHESOLVERS in SDPSETTINGS.\nSee the FAQ for more information.\n'; |
---|
| 158 | fprintf(info); |
---|
| 159 | end |
---|
| 160 | end |
---|
| 161 | end |
---|
| 162 | if length(EXISTTIME) > 5 |
---|
| 163 | EXISTTIME = EXISTTIME(end-4:end); |
---|
| 164 | NCHECKS = 5; |
---|
| 165 | end |
---|
| 166 | CACHED_SOLVERS = solvers; |
---|
| 167 | else |
---|
| 168 | solvers = CACHED_SOLVERS; |
---|
| 169 | end |
---|
| 170 | |
---|
| 171 | % ************************************************************************* |
---|
| 172 | % NO SOLVER AVAILABLE |
---|
| 173 | % ************************************************************************* |
---|
| 174 | if isempty(solvers) |
---|
| 175 | diagnostic.solvertime = 0; |
---|
| 176 | if isempty(options.solver) |
---|
| 177 | diagnostic.info = yalmiperror(-2,'YALMIP'); |
---|
| 178 | diagnostic.problem = -2; |
---|
| 179 | else |
---|
| 180 | diagnostic.info = yalmiperror(-3,'YALMIP'); |
---|
| 181 | diagnostic.problem = -3; |
---|
| 182 | end |
---|
| 183 | if warningon & options.warning & isempty(findstr(diagnostic.info,'No problems detected')) |
---|
| 184 | disp(['Warning: ' diagnostic.info]); |
---|
| 185 | end |
---|
| 186 | return |
---|
| 187 | end |
---|
| 188 | |
---|
| 189 | % ************************************************************************* |
---|
| 190 | % WHAT KIND OF PROBLEM DO WE HAVE NOW? |
---|
| 191 | % ************************************************************************* |
---|
| 192 | [ProblemClass,integer_variables,binary_variables,parametric_variables,uncertain_variables,quad_info] = categorizeproblem(F,P,h,options.relax,parametric,evaluation_based); |
---|
| 193 | |
---|
| 194 | % ************************************************************************* |
---|
| 195 | % SELECT SUITABLE SOLVER |
---|
| 196 | % ************************************************************************* |
---|
| 197 | [solver,problem] = selectsolver(options,ProblemClass,solvers,socp_are_really_qc); |
---|
| 198 | |
---|
| 199 | if isempty(solver) |
---|
| 200 | diagnostic.solvertime = 0; |
---|
| 201 | if problem == -4 |
---|
| 202 | diagnostic.info = yalmiperror(problem,options.solver); |
---|
| 203 | else |
---|
| 204 | diagnostic.info = yalmiperror(problem,'YALMIP'); |
---|
| 205 | end |
---|
| 206 | diagnostic.problem = problem; |
---|
| 207 | |
---|
| 208 | if warningon & options.warning |
---|
| 209 | disp(['Warning: ' diagnostic.info]); |
---|
| 210 | end |
---|
| 211 | return |
---|
| 212 | end |
---|
| 213 | if length(solver.version)>0 |
---|
| 214 | solver.tag = [solver.tag '-' solver.version]; |
---|
| 215 | end |
---|
| 216 | |
---|
| 217 | % ************************************************************************* |
---|
| 218 | % DID WE SELECT THE INTERNAL BNB SOLVER |
---|
| 219 | % IN THAT CASE, SELECT LOCAL SOLVER |
---|
| 220 | % (UNLESS ALREADY SPECIFIED IN OPTIONS.BNB) |
---|
| 221 | % ************************************************************************* |
---|
| 222 | localsolver.qc = 0; |
---|
| 223 | localsolver = solver; |
---|
| 224 | if strcmpi(solver.tag,'bnb') |
---|
| 225 | temp_options = options; |
---|
| 226 | temp_options.solver = options.bnb.solver; |
---|
| 227 | tempProblemClass = ProblemClass; |
---|
| 228 | tempProblemClass.constraint.binary = 0; |
---|
| 229 | tempProblemClass.constraint.integer = 0; |
---|
| 230 | localsolver = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc); |
---|
| 231 | if isempty(localsolver) | strcmpi(localsolver.tag,'bnb') |
---|
| 232 | diagnostic.solvertime = 0; |
---|
| 233 | diagnostic.info = yalmiperror(-2,'YALMIP'); |
---|
| 234 | diagnostic.problem = -2; |
---|
| 235 | return |
---|
| 236 | end |
---|
| 237 | solver.lower = localsolver; |
---|
| 238 | end |
---|
| 239 | |
---|
| 240 | % ************************************************************************* |
---|
| 241 | % DID WE SELECT THE INTERNAL EXPERIMENTAL KKT SOLVER |
---|
| 242 | % IN THAT CASE, SELECT SOLVER TO SOLVE THE MILP PROBLEM |
---|
| 243 | % ************************************************************************* |
---|
| 244 | localsolver.qc = 0; |
---|
| 245 | localsolver = solver; |
---|
| 246 | if strcmpi(solver.tag,'kktqp') |
---|
| 247 | temp_options = options; |
---|
| 248 | temp_options.solver = ''; |
---|
| 249 | tempProblemClass = ProblemClass; |
---|
| 250 | tempProblemClass.constraint.binary = 1; |
---|
| 251 | tempProblemClass.objective.quadratic.convex = 0; |
---|
| 252 | tempProblemClass.objective.quadratic.nonconvex = 0; |
---|
| 253 | localsolver = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc); |
---|
| 254 | if isempty(localsolver) | strcmpi(localsolver.tag,'bnb') | strcmpi(localsolver.tag,'kktqp') |
---|
| 255 | diagnostic.solvertime = 0; |
---|
| 256 | diagnostic.info = yalmiperror(-2,'YALMIP'); |
---|
| 257 | diagnostic.problem = -2; |
---|
| 258 | return |
---|
| 259 | end |
---|
| 260 | solver.lower = localsolver; |
---|
| 261 | end |
---|
| 262 | |
---|
| 263 | % ************************************************************************* |
---|
| 264 | % DID WE SELECT THE LMIRANK? |
---|
| 265 | % FIND SDP SOLVER FOR INITIAL SOLUTION |
---|
| 266 | % ************************************************************************* |
---|
| 267 | if strcmpi(solver.tag,'lmirank') |
---|
| 268 | temp_options = options; |
---|
| 269 | temp_options.solver = options.lmirank.solver; |
---|
| 270 | tempProblemClass = ProblemClass; |
---|
| 271 | tempProblemClass.constraint.inequalities.rank = 0; |
---|
| 272 | tempProblemClass.constraint.inequalities.semidefinite.linear = 1; |
---|
| 273 | tempProblemClass.objective.linear = 1; |
---|
| 274 | initialsolver = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc); |
---|
| 275 | if isempty(initialsolver) | strcmpi(initialsolver.tag,'lmirank') |
---|
| 276 | diagnostic.solvertime = 0; |
---|
| 277 | diagnostic.info = yalmiperror(-2,'YALMIP'); |
---|
| 278 | diagnostic.problem = -2; |
---|
| 279 | return |
---|
| 280 | end |
---|
| 281 | solver.initialsolver = initialsolver; |
---|
| 282 | end |
---|
| 283 | |
---|
| 284 | % ************************************************************************* |
---|
| 285 | % DID WE SELECT THE INTERNAL BMIBNB SOLVER? SELECT UPPER/LOWER SOLVERs |
---|
| 286 | % (UNLESS ALREADY SPECIFIED IN OPTIONS) |
---|
| 287 | % ************************************************************************* |
---|
| 288 | if strcmpi(solver.tag,'bmibnb') |
---|
| 289 | |
---|
| 290 | % Relax problem for lower solver |
---|
| 291 | tempProblemClass = ProblemClass; |
---|
| 292 | |
---|
| 293 | sdp = tempProblemClass.constraint.inequalities.semidefinite; |
---|
| 294 | tempProblemClass.constraint.inequalities.semidefinite.linear = sdp.linear | sdp.quadratic | sdp.polynomial; |
---|
| 295 | tempProblemClass.constraint.inequalities.semidefinite.quadratic = 0; |
---|
| 296 | tempProblemClass.constraint.inequalities.semidefinite.polynomial = 0; |
---|
| 297 | |
---|
| 298 | lp = tempProblemClass.constraint.inequalities.elementwise; |
---|
| 299 | tempProblemClass.constraint.inequalities.elementwise.linear = lp.linear | lp.quadratic.convex | lp.quadratic.nonconvex | sdp.polynomial; |
---|
| 300 | tempProblemClass.constraint.inequalities.elementwise.quadratic.convex = 0; |
---|
| 301 | tempProblemClass.constraint.inequalities.elementwise.quadratic.nonconvex = 0; |
---|
| 302 | tempProblemClass.constraint.inequalities.elementwise.polynomial = 0; |
---|
| 303 | |
---|
| 304 | equ = tempProblemClass.constraint.equalities; |
---|
| 305 | tempProblemClass.constraint.equalities.linear = equ.linear | equ.quadratic | equ.polynomial; |
---|
| 306 | tempProblemClass.constraint.equalities.quadratic = 0; |
---|
| 307 | tempProblemClass.constraint.equalities.polynomial = 0; |
---|
| 308 | |
---|
| 309 | % if nnz(quad_info.Q) > 0 |
---|
| 310 | % if all(eig(quad_info.Q) > -1e-12) |
---|
| 311 | % tempProblemClass.objective.quadratic.convex = 0; |
---|
| 312 | % tempProblemClass.objective.quadratic.nonconvex = 0; |
---|
| 313 | % else |
---|
| 314 | % tempProblemClass.objective.quadratic.convex = 0; |
---|
| 315 | tempProblemClass.objective.quadratic.nonconvex = 0; |
---|
| 316 | tempProblemClass.objective.polynomial = 0; |
---|
| 317 | % end |
---|
| 318 | % end |
---|
| 319 | |
---|
| 320 | tempProblemClass.constraint.inequalities.rank = 0; |
---|
| 321 | tempProblemClass.evaluation = 0; |
---|
| 322 | |
---|
| 323 | temp_options = options; |
---|
| 324 | temp_options.solver = options.bmibnb.lowersolver; |
---|
| 325 | |
---|
| 326 | % If the problem actually is quadratic, try to get a convex problem |
---|
| 327 | % this will typically allow us to solver better lower bounding problems |
---|
| 328 | % (we don't have to linearize the cost) |
---|
| 329 | [lowersolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc); |
---|
| 330 | if isempty(lowersolver)| strcmpi(lowersolver.tag,'bmibnb') | strcmpi(lowersolver.tag,'bnb') |
---|
| 331 | % No, probably non-convex cost. Pick a linear solver instead and go |
---|
| 332 | % for lower bound based on a complete "linearization" |
---|
| 333 | tempProblemClass.objective.quadratic.convex = 0; |
---|
| 334 | [lowersolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc); |
---|
| 335 | end |
---|
| 336 | |
---|
| 337 | if isempty(lowersolver)| strcmpi(lowersolver.tag,'bmibnb') | strcmpi(lowersolver.tag,'bnb') |
---|
| 338 | tempbinary = tempProblemClass.constraint.binary; |
---|
| 339 | tempinteger = tempProblemClass.constraint.integer; |
---|
| 340 | tempProblemClass.constraint.binary = 0; |
---|
| 341 | tempProblemClass.constraint.integer = 0; |
---|
| 342 | [lowersolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc); |
---|
| 343 | tempProblemClass.constraint.binary = tempbinary; |
---|
| 344 | tempProblemClass.constraint.integer = tempinteger; |
---|
| 345 | end |
---|
| 346 | |
---|
| 347 | if isempty(lowersolver) | strcmpi(lowersolver.tag,'bmibnb') | strcmpi(lowersolver.tag,'bnb') |
---|
| 348 | diagnostic.solvertime = 0; |
---|
| 349 | diagnostic.info = yalmiperror(-2,'YALMIP'); |
---|
| 350 | diagnostic.problem = -2; |
---|
| 351 | return |
---|
| 352 | end |
---|
| 353 | solver.lowercall = lowersolver.call; |
---|
| 354 | solver.lowersolver = lowersolver; |
---|
| 355 | |
---|
| 356 | temp_options = options; |
---|
| 357 | temp_options.solver = options.bmibnb.uppersolver; |
---|
| 358 | temp_ProblemClass = ProblemClass; |
---|
| 359 | temp_ProblemClass.constraint.binary = 0; |
---|
| 360 | temp_ProblemClass.constraint.integer = 0; |
---|
| 361 | temp_ProblemClass.objective.quadratic.convex = 0; |
---|
| 362 | temp_ProblemClass.objective.quadratic.nonconvex = 0; |
---|
| 363 | [uppersolver,problem] = selectsolver(temp_options,temp_ProblemClass,solvers,socp_are_really_qc); |
---|
| 364 | if strcmpi(uppersolver.tag,'bnb') |
---|
| 365 | temp_options.solver = 'none'; |
---|
| 366 | [uppersolver,problem] = selectsolver(temp_options,temp_ProblemClass,solvers,socp_are_really_qc); |
---|
| 367 | end |
---|
| 368 | if isempty(uppersolver) | strcmpi(uppersolver.tag,'bmibnb') |
---|
| 369 | diagnostic.solvertime = 0; |
---|
| 370 | diagnostic.info = yalmiperror(-2,'YALMIP'); |
---|
| 371 | diagnostic.problem = -2; |
---|
| 372 | return |
---|
| 373 | end |
---|
| 374 | solver.uppercall = uppersolver.call; |
---|
| 375 | solver.uppersolver = uppersolver; |
---|
| 376 | |
---|
| 377 | |
---|
| 378 | temp_options = options; |
---|
| 379 | temp_options.solver = options.bmibnb.lpsolver; |
---|
| 380 | tempProblemClass.constraint.inequalities.semidefinite.linear = 0; |
---|
| 381 | tempProblemClass.constraint.inequalities.semidefinite.quadratic = 0; |
---|
| 382 | tempProblemClass.constraint.inequalities.semidefinite.polynomial = 0; |
---|
| 383 | tempProblemClass.constraint.inequalities.secondordercone = 0; |
---|
| 384 | tempProblemClass.objective.quadratic.convex = 0; |
---|
| 385 | tempProblemClass.objective.quadratic.nonconvex = 0; |
---|
| 386 | tempProblemClass.objective.quadratic.nonconvex = 0; |
---|
| 387 | tempProblemClass.objective.polynomial = 0; |
---|
| 388 | |
---|
| 389 | [lpsolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc); |
---|
| 390 | |
---|
| 391 | if isempty(lowersolver)| strcmpi(lowersolver.tag,'bmibnb') |
---|
| 392 | tempbinary = tempProblemClass.constraint.binary; |
---|
| 393 | tempProblemClass.constraint.binary = 0; |
---|
| 394 | [lpsolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc); |
---|
| 395 | tempProblemClass.constraint.binary = tempbinary; |
---|
| 396 | end |
---|
| 397 | |
---|
| 398 | if isempty(lpsolver) | strcmpi(lpsolver.tag,'bmibnb') |
---|
| 399 | diagnostic.solvertime = 0; |
---|
| 400 | diagnostic.info = yalmiperror(-2,'YALMIP'); |
---|
| 401 | diagnostic.problem = -2; |
---|
| 402 | return |
---|
| 403 | end |
---|
| 404 | solver.lpsolver = lpsolver; |
---|
| 405 | solver.lpcall = lpsolver.call; |
---|
| 406 | end |
---|
| 407 | |
---|
| 408 | % ************************************************************************* |
---|
| 409 | % DID WE SELECT THE INTERNAL SDPMILP SOLVER |
---|
| 410 | % This solver solves MISDP problems by solving MILP problems and adding SDP |
---|
| 411 | % cuts based on the infasible MILP solution. |
---|
| 412 | % ************************************************************************* |
---|
| 413 | if strcmpi(solver.tag,'cutsdp') |
---|
| 414 | |
---|
| 415 | % Relax problem for lower solver |
---|
| 416 | tempProblemClass = ProblemClass; |
---|
| 417 | tempProblemClass.constraint.inequalities.elementwise.linear = tempProblemClass.constraint.inequalities.elementwise.linear | tempProblemClass.constraint.inequalities.semidefinite.linear | tempProblemClass.constraint.inequalities.secondordercone; |
---|
| 418 | tempProblemClass.constraint.inequalities.semidefinite.linear = 0; |
---|
| 419 | tempProblemClass.constraint.inequalities.secondordercone = 0; |
---|
| 420 | |
---|
| 421 | temp_options = options; |
---|
| 422 | temp_options.solver = options.cutsdp.solver; |
---|
| 423 | |
---|
| 424 | [cutsolver,problem] = selectsolver(temp_options,tempProblemClass,solvers,socp_are_really_qc); |
---|
| 425 | |
---|
| 426 | if isempty(cutsolver) | strcmpi(cutsolver.tag,'cutsdp') |strcmpi(cutsolver.tag,'bmibnb') | strcmpi(cutsolver.tag,'bnb') |
---|
| 427 | diagnostic.solvertime = 0; |
---|
| 428 | diagnostic.info = yalmiperror(-2,'YALMIP'); |
---|
| 429 | diagnostic.problem = -2; |
---|
| 430 | return |
---|
| 431 | end |
---|
| 432 | solver.cutsolver = cutsolver; |
---|
| 433 | end |
---|
| 434 | |
---|
| 435 | showprogress(['Solver chosen : ' solver.tag],options.showprogress); |
---|
| 436 | |
---|
| 437 | % ************************************************************************* |
---|
| 438 | % CONVERT MAXDET TO SDP USING GEOMEAN? |
---|
| 439 | % ************************************************************************* |
---|
| 440 | % MAXDET using geometric mean construction |
---|
| 441 | if ~isempty(P) & solver.objective.maxdet==0 |
---|
| 442 | t = sdpvar(1,1); |
---|
| 443 | if isequal(P,sdpvar(F(end))) |
---|
| 444 | F = F(1:end-1); |
---|
| 445 | end |
---|
| 446 | F = F + detset(t,P); |
---|
| 447 | if isempty(h) |
---|
| 448 | h = -t; |
---|
| 449 | else |
---|
| 450 | h = h-t; |
---|
| 451 | % Warn about logdet -> det^1/m |
---|
| 452 | if options.verbose>0 & options.warning>0 |
---|
| 453 | disp(' ') |
---|
| 454 | disp('Objective c''x-logdet(P) has been changed to c''x-det(P)^(1/(2^ceil(log2(length(X))))).') |
---|
| 455 | disp('See the MAXDET section in the manual for details.') |
---|
| 456 | disp(' ') |
---|
| 457 | end |
---|
| 458 | end |
---|
| 459 | end |
---|
| 460 | |
---|
| 461 | % ************************************************************************* |
---|
| 462 | % Change binary variables to integer? |
---|
| 463 | % ************************************************************************* |
---|
| 464 | old_binary_variables = binary_variables; |
---|
| 465 | if ~isempty(binary_variables) & (solver.constraint.binary==0) |
---|
| 466 | x_bin = recover(binary_variables(ismember(binary_variables,unique([getvariables(h) getvariables(F)])))); |
---|
| 467 | F = F + set(x_bin<1)+set(x_bin>0); |
---|
| 468 | integer_variables = union(binary_variables,integer_variables); |
---|
| 469 | binary_variables = []; |
---|
| 470 | end |
---|
| 471 | |
---|
| 472 | % ************************************************************************* |
---|
| 473 | % Model quadratics using SOCP? |
---|
| 474 | % Should not be done when using PENNLP or BMIBNB, or if we have relaxed the |
---|
| 475 | % monmial terms or...Ouch, need to clean up all special cases, this sucks. |
---|
| 476 | % ************************************************************************* |
---|
| 477 | convertQuadraticObjective = ~strcmpi(solver.tag,'pennlp-standard'); |
---|
| 478 | convertQuadraticObjective = convertQuadraticObjective & ~strcmpi(solver.tag,'bmibnb'); |
---|
| 479 | convertQuadraticObjective = convertQuadraticObjective & ( ~(options.relax==1 | options.relax==3) &(~isempty(quad_info) & solver.objective.quadratic.convex==0 | (~isempty(quad_info) & strcmp(solver.tag,'bnb') & localsolver.objective.quadratic.convex==0))); |
---|
| 480 | if strcmpi(solver.tag,'bnb') & ~isempty(quad_info) |
---|
| 481 | if solver.lower.objective.quadratic.convex==0 |
---|
| 482 | convertQuadraticObjective = 1; |
---|
| 483 | end |
---|
| 484 | end |
---|
| 485 | if convertQuadraticObjective |
---|
| 486 | t = sdpvar(1,1); |
---|
| 487 | x = quad_info.x; |
---|
| 488 | R = quad_info.R; |
---|
| 489 | if ~isempty(R) |
---|
| 490 | c = quad_info.c; |
---|
| 491 | f = quad_info.f; |
---|
| 492 | F = F + lmi(cone([2*R*x;1-(t-c'*x-f)],1+t-c'*x-f)); |
---|
| 493 | h = t; |
---|
| 494 | end |
---|
| 495 | quad_info = []; |
---|
| 496 | end |
---|
| 497 | if solver.constraint.inequalities.rotatedsecondordercone == 0 |
---|
| 498 | [F,changed] = convertlorentz(F); |
---|
| 499 | if changed |
---|
| 500 | options.saveduals = 0; % We cannot calculate duals since we change the problem |
---|
| 501 | end |
---|
| 502 | end |
---|
| 503 | % Whoa, horrible tests to find out when to convert SOCP to SDP |
---|
| 504 | % This should not be done if : |
---|
| 505 | % 1. Problem is actually a QCQP and solver supports this |
---|
| 506 | % 2. Problem is integer, local solver supports SOCC |
---|
| 507 | % 3. Solver supports SOCC |
---|
| 508 | if ~((solver.constraint.inequalities.elementwise.quadratic.convex == 1) & socp_are_really_qc) |
---|
| 509 | if ~(strcmp(solver.tag,'bnb') & socp_are_really_qc & localsolver.constraint.inequalities.elementwise.quadratic.convex==1 ) |
---|
| 510 | if ((solver.constraint.inequalities.secondordercone == 0) | (strcmpi(solver.tag,'bnb') & localsolver.constraint.inequalities.secondordercone==0)) |
---|
| 511 | [F,changed] = convertsocp(F); |
---|
| 512 | if changed |
---|
| 513 | options.saveduals = 0; % We cannot calculate duals since we change the problem |
---|
| 514 | end |
---|
| 515 | end |
---|
| 516 | end |
---|
| 517 | end |
---|
| 518 | |
---|
| 519 | % ************************************************************************* |
---|
| 520 | % Add logaritmic barrier cost/constraint for MAXDET and SDPT3-4. Note we |
---|
| 521 | % have to add it her ein order for a complex valued matrix to be converted. |
---|
| 522 | % ************************************************************************* |
---|
| 523 | if ~isempty(P) & solver.objective.maxdet==1 |
---|
| 524 | F = F + set(P > 0); |
---|
| 525 | end |
---|
| 526 | |
---|
| 527 | if ((solver.complex==0) & ProblemClass.complex) | ((strcmp(solver.tag,'bnb') & localsolver.complex==0) & ProblemClass.complex) |
---|
| 528 | showprogress('Converting to real constraints',options.showprogress) |
---|
| 529 | F = imag2reallmi(F); |
---|
| 530 | options.saveduals = 0; % We cannot calculate duals since we change the problem |
---|
| 531 | end |
---|
| 532 | |
---|
| 533 | % ************************************************************************* |
---|
| 534 | % CREATE OBJECTIVE FUNCTION c'*x+x'*Q*x |
---|
| 535 | % ************************************************************************* |
---|
| 536 | showprogress('Processing objective h(x)',options.showprogress); |
---|
| 537 | try |
---|
| 538 | if strcmpi(solver.tag,'bmibnb') | strcmpi(solver.tag,'pennlp-standard') |
---|
| 539 | tempoptions = options; |
---|
| 540 | tempoptions.relax = 1; |
---|
| 541 | [c,Q,f]=createobjective(h,P,tempoptions,quad_info); |
---|
| 542 | else |
---|
| 543 | [c,Q,f]=createobjective(h,P,options,quad_info); |
---|
| 544 | end |
---|
| 545 | catch |
---|
| 546 | error(lasterr) |
---|
| 547 | end |
---|
| 548 | |
---|
| 549 | % ************************************************************************* |
---|
| 550 | % Convert {F(x),G(x)} to a numerical SeDuMi-like format |
---|
| 551 | % ************************************************************************* |
---|
| 552 | showprogress('Processing F(x)',options.showprogress); |
---|
| 553 | [F_struc,K,KCut] = lmi2sedumistruct(F); |
---|
| 554 | % We add a field to remember the dimmesion of the logarithmic cost. |
---|
| 555 | % Actually, the actually value is not interesting, we know that the |
---|
| 556 | % logarithmic cost corresponds to the last SDP constraint anyway |
---|
| 557 | K.m = length(P); |
---|
| 558 | |
---|
| 559 | % ************************************************************************* |
---|
| 560 | % SOME HORRIBLE CODE TO DETERMINE USED VARIABLES |
---|
| 561 | % ************************************************************************* |
---|
| 562 | % Which sdpvar variables are actually in the problem |
---|
| 563 | used_variables_LMI = find(any(F_struc(:,2:end),1)); |
---|
| 564 | used_variables_obj = find(any(c',1) | any(Q)); |
---|
| 565 | used_variables = uniquestripped([used_variables_LMI used_variables_obj]); |
---|
| 566 | % The problem is that linear terms might be missing in problems with only |
---|
| 567 | % nonlinear expressions |
---|
| 568 | [monomtable,variabletype] = yalmip('monomtable'); |
---|
| 569 | if (options.relax==1)|(options.relax==3) |
---|
| 570 | monomtable = []; |
---|
| 571 | nonlinearvariables = []; |
---|
| 572 | linearvariables = used_variables; |
---|
| 573 | else |
---|
| 574 | nonlinearvariables = find(variabletype); |
---|
| 575 | linearvariables = used_variables(find(variabletype(used_variables)==0)); |
---|
| 576 | end |
---|
| 577 | needednonlinear = nonlinearvariables(ismembc(nonlinearvariables,used_variables)); |
---|
| 578 | linearinnonlinear = find(sum(abs(monomtable(needednonlinear,:)),1)); |
---|
| 579 | missinglinear = setdiff(linearinnonlinear(:),linearvariables); |
---|
| 580 | used_variables = uniquestripped([used_variables(:);missinglinear(:)]); |
---|
| 581 | |
---|
| 582 | % ************************************************************************* |
---|
| 583 | % So are we done now? No... What about variables hiding inside so called |
---|
| 584 | % evaluation variables. We detect these, and at the same time set up the |
---|
| 585 | % structures needed to support general functions such as exp, log, etc |
---|
| 586 | % NOTE : This is experimental code |
---|
| 587 | % FIX : Clean up... |
---|
| 588 | % ************************************************************************* |
---|
| 589 | [evalMap,evalVariables,used_variables,nonlinearvariables,linearvariables] = detectHiddenNonlinear(used_variables,options,nonlinearvariables,linearvariables); |
---|
| 590 | |
---|
| 591 | % ************************************************************************* |
---|
| 592 | % REMOVE UNNECESSARY VARIABLES FROM PROBLEM |
---|
| 593 | % ************************************************************************* |
---|
| 594 | if length(used_variables)<yalmip('nvars') |
---|
| 595 | c = c(used_variables); |
---|
| 596 | Q = Q(:,used_variables);Q = Q(used_variables,:); |
---|
| 597 | if ~isempty(F_struc) |
---|
| 598 | F_struc = sparse(F_struc(:,[1 1+used_variables])); |
---|
| 599 | end |
---|
| 600 | end |
---|
| 601 | |
---|
| 602 | % ************************************************************************* |
---|
| 603 | % Map variables and constraints in low-rank definition to local stuff |
---|
| 604 | % ************************************************************************* |
---|
| 605 | if ~isempty(lowrankdetails) |
---|
| 606 | % Identifiers of the SDP constraints |
---|
| 607 | lmiid = getlmiid(F); |
---|
| 608 | for i = 1:length(lowrankdetails) |
---|
| 609 | lowrankdetails{i}.id = find(ismember(lmiid,lowrankdetails{i}.id)); |
---|
| 610 | if ~isempty(lowrankdetails{i}.variables) |
---|
| 611 | index = ismember(used_variables,lowrankdetails{i}.variables); |
---|
| 612 | lowrankdetails{i}.variables = find(index); |
---|
| 613 | end |
---|
| 614 | end |
---|
| 615 | end |
---|
| 616 | |
---|
| 617 | % ************************************************************************* |
---|
| 618 | % SPECIAL VARIABLES |
---|
| 619 | % Relax = 1 : relax both integers and nonlinear stuff |
---|
| 620 | % Relax = 2 : relax integers |
---|
| 621 | % Relax = 3 : relax nonlinear stuff |
---|
| 622 | % ************************************************************************* |
---|
| 623 | if (options.relax==1) | (options.relax==3) |
---|
| 624 | nonlins = []; |
---|
| 625 | end |
---|
| 626 | if (options.relax) | (options.relax==2) |
---|
| 627 | integer_variables = []; |
---|
| 628 | binary_variables = []; |
---|
| 629 | old_binary_variables = find(ismember(used_variables,old_binary_variables)); |
---|
| 630 | else |
---|
| 631 | integer_variables = find(ismember(used_variables,integer_variables)); |
---|
| 632 | binary_variables = find(ismember(used_variables,binary_variables)); |
---|
| 633 | old_binary_variables = find(ismember(used_variables,old_binary_variables)); |
---|
| 634 | end |
---|
| 635 | parametric_variables = find(ismember(used_variables,parametric_variables)); |
---|
| 636 | extended_variables = find(ismember(used_variables,yalmip('extvariables'))); |
---|
| 637 | |
---|
| 638 | |
---|
| 639 | % ************************************************************************* |
---|
| 640 | % Equality constraints not supported or supposed to be removed |
---|
| 641 | % ************************************************************************* |
---|
| 642 | % We may save some data in order to reconstruct |
---|
| 643 | % dual variables related to equality constraints that |
---|
| 644 | % have been removed. |
---|
| 645 | oldF_struc = []; |
---|
| 646 | oldQ = []; |
---|
| 647 | oldc = []; |
---|
| 648 | oldK = K; |
---|
| 649 | Fremoved = []; |
---|
| 650 | if (K.f>0) |
---|
| 651 | % reduce if user explicitely says remove, or user says nothing but |
---|
| 652 | % solverdefinitions does, and there are no nonlinear variables |
---|
| 653 | if ((options.removeequalities==1 | options.removeequalities==2) & isempty(intersect(used_variables,nonlinearvariables))) | ((options.removeequalities==0) & (solver.constraint.equalities.linear==-1)) |
---|
| 654 | showprogress('Solving equalities',options.showprogress); |
---|
| 655 | [x_equ,H,A_equ,b_equ,factors] = solveequalities(F_struc,K,options.removeequalities==1); |
---|
| 656 | % Exit if no consistent solution exist |
---|
| 657 | if (norm(A_equ*x_equ-b_equ,'inf')>1e-5)%sqrt(eps)*size(A_equ,2)) |
---|
| 658 | diagnostic.solvertime = 0; |
---|
| 659 | diagnostic.info = yalmiperror(1,'YALMIP'); |
---|
| 660 | diagnostic.problem = 1; |
---|
| 661 | solution = diagnostic; |
---|
| 662 | solution.variables = used_variables(:); |
---|
| 663 | solution.optvar = x_equ; |
---|
| 664 | % And we are done! Save the result |
---|
| 665 | % sdpvar('setSolution',solution); |
---|
| 666 | return |
---|
| 667 | end |
---|
| 668 | % We dont need the rows for equalities anymore |
---|
| 669 | oldF_struc = F_struc; |
---|
| 670 | oldc = c; |
---|
| 671 | oldQ = Q; |
---|
| 672 | oldK = K; |
---|
| 673 | F_struc = F_struc(K.f+1:end,:); |
---|
| 674 | K.f = 0; |
---|
| 675 | Fold = F; |
---|
| 676 | [nlmi neq]=size(F); |
---|
| 677 | F = lmi; |
---|
| 678 | Fremoved = set([]); |
---|
| 679 | for i = 1:nlmi+neq |
---|
| 680 | if ~is(Fold(i),'equality') |
---|
| 681 | F = F + Fold(i); |
---|
| 682 | else |
---|
| 683 | Fremoved = Fremoved + Fold(i); |
---|
| 684 | end |
---|
| 685 | end |
---|
| 686 | |
---|
| 687 | % No variables left. Problem solved! |
---|
| 688 | if size(H,2)==0 |
---|
| 689 | diagnostic.solvertime = 0; |
---|
| 690 | diagnostic.info = yalmiperror(0,'YALMIP'); |
---|
| 691 | diagnostic.problem = 0; |
---|
| 692 | solution = diagnostic; |
---|
| 693 | solution.variables = used_variables(:); |
---|
| 694 | solution.optvar = x_equ; |
---|
| 695 | % And we are done! Save the result |
---|
| 696 | % Note, no dual is saved |
---|
| 697 | sdpvar('setSolution',solution); |
---|
| 698 | p = checkset(F); |
---|
| 699 | if any(p<1e-5) |
---|
| 700 | diagnostic.info = yalmiperror(1,'YALMIP'); |
---|
| 701 | diagnostic.problem = 1; |
---|
| 702 | end |
---|
| 703 | return |
---|
| 704 | end |
---|
| 705 | showprogress('Converting problem to new basis',options.showprogress) |
---|
| 706 | |
---|
| 707 | % objective in new basis |
---|
| 708 | f = f + x_equ'*Q*x_equ; |
---|
| 709 | c = H'*c + 2*H'*Q*x_equ; |
---|
| 710 | Q = H'*Q*H;Q=((Q+Q')/2); |
---|
| 711 | % LMI in new basis |
---|
| 712 | F_struc = [F_struc*[1;x_equ] F_struc(:,2:end)*H]; |
---|
| 713 | else |
---|
| 714 | % Solver does not support equality constraints and user specifies |
---|
| 715 | % double-sided inequalitis to remove them |
---|
| 716 | if (solver.constraint.equalities.linear==0 | options.removeequalities==-1) |
---|
| 717 | % Add equalities |
---|
| 718 | F_struc = [-F_struc(1:1:K.f,:);F_struc]; |
---|
| 719 | K.l = K.l+K.f*2; |
---|
| 720 | % Keep this in mind... |
---|
| 721 | K.fold = K.f; |
---|
| 722 | K.f = 0; |
---|
| 723 | end |
---|
| 724 | % For simpliciy we introduce a dummy coordinate change |
---|
| 725 | x_equ = 0; |
---|
| 726 | H = 1; |
---|
| 727 | factors = []; |
---|
| 728 | end |
---|
| 729 | else |
---|
| 730 | x_equ = 0; |
---|
| 731 | H = 1; |
---|
| 732 | factors = []; |
---|
| 733 | end |
---|
| 734 | |
---|
| 735 | |
---|
| 736 | % ************************************************************************* |
---|
| 737 | % Setup the initial solution |
---|
| 738 | % ************************************************************************* |
---|
| 739 | x0 = []; |
---|
| 740 | if options.usex0 |
---|
| 741 | if options.relax |
---|
| 742 | x0_used = relaxdouble(recover(used_variables)); |
---|
| 743 | else |
---|
| 744 | %FIX : Do directly using yalmip('solution') |
---|
| 745 | x0_used = double(recover(used_variables)); |
---|
| 746 | end |
---|
| 747 | x0 = zeros(sdpvar('nvars'),1); |
---|
| 748 | x0(used_variables) = x0_used(:); |
---|
| 749 | x0(isnan(x0))=0; |
---|
| 750 | end |
---|
| 751 | if ~isempty(x0) |
---|
| 752 | % Get a coordinate in the reduced space |
---|
| 753 | x0 = H\(x0(used_variables)-x_equ); |
---|
| 754 | end |
---|
| 755 | |
---|
| 756 | % Monomial table for nonlinear variables |
---|
| 757 | % FIX : Why here!!! mt handled above also |
---|
| 758 | [mt,variabletype] = yalmip('monomtable'); |
---|
| 759 | if size(mt,1)>size(mt,2) |
---|
| 760 | mt(size(mt,1),size(mt,1)) = 0; |
---|
| 761 | end |
---|
| 762 | % In local variables |
---|
| 763 | mt = mt(used_variables,used_variables); |
---|
| 764 | variabletype = variabletype(used_variables); |
---|
| 765 | if (options.relax)|(options.relax==3) |
---|
| 766 | mt = eye(length(used_variables)); |
---|
| 767 | end |
---|
| 768 | |
---|
| 769 | % FIX : Make sure these things work... |
---|
| 770 | lub = yalmip('getbounds',used_variables); |
---|
| 771 | lb = lub(:,1)-inf; |
---|
| 772 | ub = lub(:,2)+inf; |
---|
| 773 | lb(old_binary_variables) = max(lb(old_binary_variables),0); |
---|
| 774 | ub(old_binary_variables) = min(ub(old_binary_variables),1); |
---|
| 775 | |
---|
| 776 | % This does not work if we have used removeequalities, so we clear them for |
---|
| 777 | % safety. note that bounds are not guaranteed to be used according to the |
---|
| 778 | % manual, so this is allowed, although it might be a bit inconsistent to |
---|
| 779 | % some users. |
---|
| 780 | if ~isempty(oldc) |
---|
| 781 | lb = []; |
---|
| 782 | ub = []; |
---|
| 783 | end |
---|
| 784 | |
---|
| 785 | % ************************************************************************* |
---|
| 786 | % GENERAL DATA EXCHANGE WITH SOLVER |
---|
| 787 | % ************************************************************************* |
---|
| 788 | interfacedata.F_struc = F_struc; |
---|
| 789 | interfacedata.c = c; |
---|
| 790 | interfacedata.Q = Q; |
---|
| 791 | interfacedata.f = f; |
---|
| 792 | interfacedata.K = K; |
---|
| 793 | interfacedata.lb = lb; |
---|
| 794 | interfacedata.ub = ub; |
---|
| 795 | interfacedata.x0 = x0; |
---|
| 796 | interfacedata.options = options; |
---|
| 797 | interfacedata.solver = solver; |
---|
| 798 | interfacedata.monomtable = mt; |
---|
| 799 | interfacedata.variabletype = variabletype; |
---|
| 800 | interfacedata.integer_variables = integer_variables; |
---|
| 801 | interfacedata.binary_variables = binary_variables; |
---|
| 802 | interfacedata.uncertain_variables = []; |
---|
| 803 | interfacedata.parametric_variables= parametric_variables; |
---|
| 804 | interfacedata.extended_variables = extended_variables; |
---|
| 805 | interfacedata.used_variables = used_variables; |
---|
| 806 | interfacedata.lowrankdetails = lowrankdetails; |
---|
| 807 | interfacedata.problemclass = ProblemClass; |
---|
| 808 | interfacedata.KCut = KCut; |
---|
| 809 | interfacedata.getsolvertime = 1; |
---|
| 810 | % Data to be able to recover duals when model is reduced |
---|
| 811 | interfacedata.oldF_struc = oldF_struc; |
---|
| 812 | interfacedata.oldc = oldc; |
---|
| 813 | interfacedata.oldK = oldK; |
---|
| 814 | interfacedata.factors = factors; |
---|
| 815 | interfacedata.Fremoved = Fremoved; |
---|
| 816 | interfacedata.evalMap = evalMap; |
---|
| 817 | interfacedata.evalVariables = evalVariables; |
---|
| 818 | |
---|
| 819 | % ************************************************************************* |
---|
| 820 | % GENERAL DATA EXCANGE TO RECOVER SOLUTION AND UPDATE YALMIP VARIABLES |
---|
| 821 | % ************************************************************************* |
---|
| 822 | recoverdata.H = H; |
---|
| 823 | recoverdata.x_equ = x_equ; |
---|
| 824 | recoverdata.used_variables = used_variables; |
---|
| 825 | |
---|
| 826 | function yesno = warningon |
---|
| 827 | |
---|
| 828 | s = warning; |
---|
| 829 | if isa(s,'char') |
---|
| 830 | yesno = isequal(s,'on'); |
---|
| 831 | else |
---|
| 832 | yesno = isequal(s(1).state,'on'); |
---|
| 833 | end |
---|
| 834 | |
---|
| 835 | |
---|
| 836 | |
---|
| 837 | |
---|
| 838 | |
---|
| 839 | |
---|
| 840 | |
---|
| 841 | |
---|
| 842 | |
---|
| 843 | |
---|
| 844 | function [evalMap,evalVariables,used_variables,nonlinearvariables,linearvariables] = detectHiddenNonlinear(used_variables,options,nonlinearvariables,linearvariables) |
---|
| 845 | |
---|
| 846 | evalVariables = yalmip('evalVariables'); |
---|
| 847 | old_used_variables = used_variables; |
---|
| 848 | goon = 1; |
---|
| 849 | if ~isempty(evalVariables) |
---|
| 850 | while goon |
---|
| 851 | % Which used_variables are representing general functions |
---|
| 852 | evalVariables = yalmip('evalVariables'); |
---|
| 853 | usedEvalVariables = find(ismember(used_variables,evalVariables)); |
---|
| 854 | evalMap = yalmip('extstruct',used_variables(usedEvalVariables)); |
---|
| 855 | if ~isa(evalMap,'cell') |
---|
| 856 | evalMap = {evalMap}; |
---|
| 857 | end |
---|
| 858 | % Find all variables used in the arguments of these functions |
---|
| 859 | hidden = []; |
---|
| 860 | for i = 1:length(evalMap) |
---|
| 861 | if isequal(getbase(evalMap{i}.arg{1}),[0 1])% & is(evalMap{i}.arg{1},'linear') |
---|
| 862 | for j = 1:length(evalMap{i}.arg)-1 |
---|
| 863 | % The last argument is the help variable z in the |
---|
| 864 | % transformation from f(ax+b) to f(z),z==ax+b. We should not |
---|
| 865 | % use this transformation if the argument already is unitary |
---|
| 866 | hidden = [hidden getvariables(evalMap{i}.arg{j})]; |
---|
| 867 | end |
---|
| 868 | else |
---|
| 869 | for j = 1:length(evalMap{i}.arg) |
---|
| 870 | % The last argument is the help variable z in the |
---|
| 871 | % transformation from f(ax+b) to f(z),z==ax+b. We should not |
---|
| 872 | % use this transformation if the argument already is unitary |
---|
| 873 | hidden = [hidden getvariables(evalMap{i}.arg{j})]; |
---|
| 874 | end |
---|
| 875 | end |
---|
| 876 | end |
---|
| 877 | used_variables = union(used_variables,hidden); |
---|
| 878 | |
---|
| 879 | % The problem is that linear terms might be missing in problems with only |
---|
| 880 | % nonlinear expressions |
---|
| 881 | [monomtable,variabletype] = yalmip('monomtable'); |
---|
| 882 | if (options.relax==1)|(options.relax==3) |
---|
| 883 | monomtable = []; |
---|
| 884 | nonlinearvariables = []; |
---|
| 885 | linearvariables = used_variables; |
---|
| 886 | else |
---|
| 887 | nonlinearvariables = find(variabletype); |
---|
| 888 | linearvariables = used_variables(find(variabletype(used_variables)==0)); |
---|
| 889 | end |
---|
| 890 | needednonlinear = nonlinearvariables(ismembc(nonlinearvariables,used_variables)); |
---|
| 891 | linearinnonlinear = find(sum(abs(monomtable(needednonlinear,:)),1)); |
---|
| 892 | missinglinear = setdiff(linearinnonlinear(:),linearvariables); |
---|
| 893 | used_variables = uniquestripped([used_variables(:);missinglinear(:)]); |
---|
| 894 | |
---|
| 895 | |
---|
| 896 | usedEvalVariables = find(ismember(used_variables,evalVariables)); |
---|
| 897 | evalMap = yalmip('extstruct',used_variables(usedEvalVariables)); |
---|
| 898 | if ~isa(evalMap,'cell') |
---|
| 899 | evalMap = {evalMap}; |
---|
| 900 | end |
---|
| 901 | evalVariables = usedEvalVariables; |
---|
| 902 | |
---|
| 903 | for i = 1:length(evalMap) |
---|
| 904 | if isequal(getbase(evalMap{i}.arg{1}),[0 1]) |
---|
| 905 | index = ismember(used_variables,getvariables(evalMap{i}.arg{1})); |
---|
| 906 | evalMap{i}.variableIndex = find(index); |
---|
| 907 | else |
---|
| 908 | index = ismember(used_variables,getvariables(evalMap{i}.arg{end})); |
---|
| 909 | evalMap{i}.variableIndex = find(index); |
---|
| 910 | end |
---|
| 911 | end |
---|
| 912 | goon = ~isequal(used_variables,old_used_variables); |
---|
| 913 | old_used_variables = used_variables; |
---|
| 914 | end |
---|
| 915 | else |
---|
| 916 | evalMap = []; |
---|
| 917 | end |
---|
| 918 | |
---|