[37] | 1 | function [problem,integer_variables,binary_variables,parametric_variables,uncertain_variables,quad_info] = categorizeproblem(F,P,h,relax,parametric,evaluation) |
---|
| 2 | %categorizeproblem Internal function: tries to determine the type of optimization problem |
---|
| 3 | |
---|
| 4 | % Author Johan Löfberg |
---|
| 5 | % $Id: categorizeproblem.m,v 1.12 2006/09/22 11:57:38 joloef Exp $ |
---|
| 6 | |
---|
| 7 | Counter = size(F.clauses,2); |
---|
| 8 | Ftype = zeros(Counter,1); |
---|
| 9 | |
---|
| 10 | real_data = 1; |
---|
| 11 | int_data = 0; |
---|
| 12 | bin_data = 0; |
---|
| 13 | par_data = 0; |
---|
| 14 | |
---|
| 15 | poly_constraint = 0; |
---|
| 16 | bilin_constraint = 0; |
---|
| 17 | sigm_constraint = 0; |
---|
| 18 | rank_constraint = 0; |
---|
| 19 | rank_objective = 0; |
---|
| 20 | |
---|
| 21 | integer_variables = []; |
---|
| 22 | binary_variables = []; |
---|
| 23 | parametric_variables = []; |
---|
| 24 | kyp_prob = 0; |
---|
| 25 | |
---|
| 26 | % *********************************************** |
---|
| 27 | % Setup an empty problem definition |
---|
| 28 | % *********************************************** |
---|
| 29 | problem.objective.linear = 0; |
---|
| 30 | problem.objective.quadratic.convex = 0; |
---|
| 31 | problem.objective.quadratic.nonconvex = 0; |
---|
| 32 | problem.objective.polynomial = 0; |
---|
| 33 | problem.objective.maxdet = 0; |
---|
| 34 | problem.objective.sigmonial = 0; |
---|
| 35 | |
---|
| 36 | problem.constraint.equalities.linear = 0; |
---|
| 37 | problem.constraint.equalities.quadratic = 0; |
---|
| 38 | problem.constraint.equalities.polynomial = 0; |
---|
| 39 | problem.constraint.equalities.sigmonial = 0; |
---|
| 40 | |
---|
| 41 | problem.constraint.inequalities.elementwise.linear = 0; |
---|
| 42 | problem.constraint.inequalities.elementwise.quadratic.convex = 0; |
---|
| 43 | problem.constraint.inequalities.elementwise.quadratic.nonconvex = 0; |
---|
| 44 | problem.constraint.inequalities.elementwise.sigmonial = 0; |
---|
| 45 | problem.constraint.inequalities.elementwise.polynomial = 0; |
---|
| 46 | |
---|
| 47 | problem.constraint.inequalities.semidefinite.linear = 0; |
---|
| 48 | problem.constraint.inequalities.semidefinite.quadratic = 0; |
---|
| 49 | problem.constraint.inequalities.semidefinite.polynomial = 0; |
---|
| 50 | problem.constraint.inequalities.semidefinite.sigmonial = 0; |
---|
| 51 | |
---|
| 52 | problem.constraint.inequalities.rank = 0; |
---|
| 53 | |
---|
| 54 | problem.constraint.inequalities.secondordercone = []; |
---|
| 55 | problem.constraint.inequalities.rotatedsecondordercone = []; |
---|
| 56 | |
---|
| 57 | problem.constraint.integer = 0; |
---|
| 58 | problem.constraint.binary = 0; |
---|
| 59 | |
---|
| 60 | problem.complex = 0; |
---|
| 61 | problem.parametric = parametric; |
---|
| 62 | problem.evaluation = evaluation; |
---|
| 63 | |
---|
| 64 | % ******************************************************** |
---|
| 65 | % Make a list of all globally available discrete variables |
---|
| 66 | % ******************************************************** |
---|
| 67 | integer_variables = yalmip('intvariables'); |
---|
| 68 | binary_variables = yalmip('binvariables'); |
---|
| 69 | uncertain_variables = yalmip('uncvariables'); |
---|
| 70 | for i = 1:Counter |
---|
| 71 | switch F.clauses{i}.type |
---|
| 72 | case 7 |
---|
| 73 | integer_variables = union(integer_variables,getvariables(F.clauses{i}.data)); |
---|
| 74 | case 8 |
---|
| 75 | binary_variables = union(binary_variables,getvariables(F.clauses{i}.data)); |
---|
| 76 | case 13 |
---|
| 77 | parametric_variables = union(parametric_variables,getvariables(F.clauses{i}.data)); |
---|
| 78 | otherwise |
---|
| 79 | end |
---|
| 80 | end |
---|
| 81 | |
---|
| 82 | % ******************************************************** |
---|
| 83 | % Logarithmic objective? |
---|
| 84 | % ******************************************************** |
---|
| 85 | problem.objective.maxdet = ~isempty(P); |
---|
| 86 | |
---|
| 87 | % ******************************************************** |
---|
| 88 | % Rank variables |
---|
| 89 | % ******************************************************** |
---|
| 90 | rank_variables = yalmip('rankvariables'); |
---|
| 91 | any_rank_variables = ~isempty(rank_variables); |
---|
| 92 | |
---|
| 93 | % ******************************************************** |
---|
| 94 | % Make a list of all globally available nonlinear variables |
---|
| 95 | % ******************************************************** |
---|
| 96 | [monomtable,variabletype] = yalmip('monomtable'); |
---|
| 97 | |
---|
| 98 | linear_variables = find(variabletype==0); |
---|
| 99 | nonlinear_variables = find(variabletype~=0); |
---|
| 100 | sigmonial_variables = find(variabletype==4); |
---|
| 101 | |
---|
| 102 | allvars = getvariables(F); |
---|
| 103 | any_nonlinear_variables =~isempty(find(ismembc(nonlinear_variables,allvars))); |
---|
| 104 | any_discrete_variables = ~isempty(integer_variables) | ~isempty(binary_variables); |
---|
| 105 | |
---|
| 106 | for i = 1:Counter |
---|
| 107 | |
---|
| 108 | Fi = F.clauses{i}; |
---|
| 109 | % Only real-valued data? |
---|
| 110 | real_data = real_data & isreal(F.clauses{i}.data); |
---|
| 111 | |
---|
| 112 | % Any discrete variables used |
---|
| 113 | if any_discrete_variables |
---|
| 114 | Fvar = getvariables(Fi.data); |
---|
| 115 | int_data = int_data | any(ismembc(Fvar,integer_variables)); |
---|
| 116 | bin_data = bin_data | any(ismembc(Fvar,binary_variables)); |
---|
| 117 | par_data = par_data | any(ismembc(Fvar,parametric_variables)); |
---|
| 118 | end |
---|
| 119 | |
---|
| 120 | if any_rank_variables |
---|
| 121 | rank_constraint = rank_constraint | any(ismember(getvariables(Fi.data),rank_variables)); |
---|
| 122 | end |
---|
| 123 | |
---|
| 124 | if ~any_nonlinear_variables % No nonlinearly parameterized constraints |
---|
| 125 | |
---|
| 126 | switch Fi.type |
---|
| 127 | case {1,9} |
---|
| 128 | problem.constraint.inequalities.semidefinite.linear = 1; |
---|
| 129 | case 2 |
---|
| 130 | problem.constraint.inequalities.elementwise.linear = 1; |
---|
| 131 | case 3 |
---|
| 132 | problem.constraint.equalities.linear = 1; |
---|
| 133 | case 4 |
---|
| 134 | problem.constraint.inequalities.secondordercone = 1; |
---|
| 135 | case 5 |
---|
| 136 | problem.constraint.inequalities.rotatedsecondordercone = 1; |
---|
| 137 | otherwise |
---|
| 138 | end |
---|
| 139 | else |
---|
| 140 | % Can be nonlinear stuff |
---|
| 141 | vars = getvariables(Fi.data); |
---|
| 142 | usednonlins = find(ismembc(nonlinear_variables,vars)); |
---|
| 143 | if ~isempty(usednonlins) |
---|
| 144 | usedsigmonials = find(ismember(sigmonial_variables,vars)); |
---|
| 145 | if ~isempty(usedsigmonials) |
---|
| 146 | switch Fi.type |
---|
| 147 | case 1 |
---|
| 148 | problem.constraint.inequalities.semidefinite.sigmonial = 1; |
---|
| 149 | case 2 |
---|
| 150 | problem.constraint.inequalities.elementwise.sigmonial = 1; |
---|
| 151 | case 3 |
---|
| 152 | problem.constraint.equalities.sigmonial = 1; |
---|
| 153 | case 4 |
---|
| 154 | error('Sigmonial SOCP not supported'); |
---|
| 155 | case 5 |
---|
| 156 | error('Sigmonial RSOCP not supported'); |
---|
| 157 | otherwise |
---|
| 158 | error('Report bug in problem classification (sigmonial constraint)'); |
---|
| 159 | end |
---|
| 160 | else |
---|
| 161 | deg = degree(Fi.data); |
---|
| 162 | switch deg |
---|
| 163 | |
---|
| 164 | case 1 |
---|
| 165 | switch Fi.type |
---|
| 166 | case 1 |
---|
| 167 | problem.constraint.inequalities.semidefinite.linear = 1; |
---|
| 168 | case 2 |
---|
| 169 | problem.constraint.inequalities.elementwise.linear = 1; |
---|
| 170 | case 3 |
---|
| 171 | problem.constraint.equalities.linear = 1; |
---|
| 172 | case 4 |
---|
| 173 | problem.constraint.inequalities.secondordercone = 1; |
---|
| 174 | case 5 |
---|
| 175 | problem.constraint.inequalities.rotatedsecondordercone = 1; |
---|
| 176 | otherwise |
---|
| 177 | error('Report bug in problem classification (linear constraint)'); |
---|
| 178 | end |
---|
| 179 | case 2 |
---|
| 180 | switch Fi.type |
---|
| 181 | case 1 |
---|
| 182 | problem.constraint.inequalities.semidefinite.quadratic = 1; |
---|
| 183 | case 2 |
---|
| 184 | % FIX : This should be re-used from |
---|
| 185 | % convertconvexquad |
---|
| 186 | convex = 1; |
---|
| 187 | f = Fi.data;f = f(:); |
---|
| 188 | ii = 1; |
---|
| 189 | while convex & ii<=length(f) |
---|
| 190 | [Q,c,f,x,info] = quaddecomp(f(ii)); |
---|
| 191 | |
---|
| 192 | if info | any(eig(Q) > 0) |
---|
| 193 | convex = 0; |
---|
| 194 | end |
---|
| 195 | ii= ii + 1; |
---|
| 196 | end |
---|
| 197 | if convex |
---|
| 198 | problem.constraint.inequalities.elementwise.quadratic.convex = 1; |
---|
| 199 | else |
---|
| 200 | problem.constraint.inequalities.elementwise.quadratic.nonconvex = 1; |
---|
| 201 | end |
---|
| 202 | case 3 |
---|
| 203 | problem.constraint.equalities.quadratic = 1; |
---|
| 204 | case 4 |
---|
| 205 | error |
---|
| 206 | case 5 |
---|
| 207 | error |
---|
| 208 | otherwise |
---|
| 209 | error('Report bug in problem classification (quadratic constraint)'); |
---|
| 210 | end |
---|
| 211 | otherwise |
---|
| 212 | switch Fi.type |
---|
| 213 | case 1 |
---|
| 214 | problem.constraint.inequalities.semidefinite.polynomial = 1; |
---|
| 215 | case 2 |
---|
| 216 | problem.constraint.inequalities.elementwise.polynomial = 1; |
---|
| 217 | case 3 |
---|
| 218 | problem.constraint.equalities.polynomial = 1; |
---|
| 219 | case 4 |
---|
| 220 | % problem.constraint.inequalities.secondordercone = 1; |
---|
| 221 | case 5 |
---|
| 222 | % problem.constraint.inequalities.rotatedsecondordercone = 1; |
---|
| 223 | otherwise |
---|
| 224 | error('Report bug in problem classification (polynomial constraint)'); |
---|
| 225 | end |
---|
| 226 | |
---|
| 227 | end |
---|
| 228 | end |
---|
| 229 | else |
---|
| 230 | switch Fi.type |
---|
| 231 | case 1 |
---|
| 232 | problem.constraint.inequalities.semidefinite.linear = 1; |
---|
| 233 | case 2 |
---|
| 234 | problem.constraint.inequalities.elementwise.linear = 1; |
---|
| 235 | case 3 |
---|
| 236 | problem.constraint.equalities.linear = 1; |
---|
| 237 | case 4 |
---|
| 238 | problem.constraint.inequalities.secondordercone = 1; |
---|
| 239 | case 5 |
---|
| 240 | problem.constraint.inequalities.rotatedsecondordercone = 1; |
---|
| 241 | case 7 |
---|
| 242 | problem.constraint.integer = 1; |
---|
| 243 | case 8 |
---|
| 244 | problem.constraint.binary = 1; |
---|
| 245 | otherwise |
---|
| 246 | error('Report bug in problem classification (linear constraint)'); |
---|
| 247 | end |
---|
| 248 | end |
---|
| 249 | end |
---|
| 250 | end |
---|
| 251 | |
---|
| 252 | if int_data |
---|
| 253 | problem.constraint.integer = 1; |
---|
| 254 | end |
---|
| 255 | if bin_data |
---|
| 256 | problem.constraint.binary = 1; |
---|
| 257 | end |
---|
| 258 | if ~real_data |
---|
| 259 | problem.complex = 1; |
---|
| 260 | end |
---|
| 261 | if rank_constraint |
---|
| 262 | problem.constraint.inequalities.rank = 1; |
---|
| 263 | end |
---|
| 264 | if ~isempty(uncertain_variables) |
---|
| 265 | problem.uncertain = 1; |
---|
| 266 | end |
---|
| 267 | |
---|
| 268 | if (relax==1) | (relax==2) |
---|
| 269 | problem.constraint.integer = 0; |
---|
| 270 | problem.constraint.binary = 0; |
---|
| 271 | int_data = 0; |
---|
| 272 | bin_data = 0; |
---|
| 273 | % integer_variables = []; |
---|
| 274 | % binary_variables = []; |
---|
| 275 | end |
---|
| 276 | if (relax==1) | (relax==3) |
---|
| 277 | problem.constraint.equalities.linear = problem.constraint.equalities.linear | problem.constraint.equalities.quadratic | problem.constraint.equalities.polynomial | problem.constraint.equalities.sigmonial; |
---|
| 278 | problem.constraint.equalities.quadratic = 0; |
---|
| 279 | problem.constraint.equalities.polynomial = 0; |
---|
| 280 | problem.constraint.equalities.sigmonial = 0; |
---|
| 281 | |
---|
| 282 | problem.constraint.inequalities.elementwise.linear = problem.constraint.inequalities.elementwise.linear | problem.constraint.inequalities.elementwise.quadratic.convex | problem.constraint.inequalities.elementwise.quadratic.nonconvex | problem.constraint.inequalities.elementwise.sigmonial | problem.constraint.inequalities.elementwise.polynomial; |
---|
| 283 | problem.constraint.inequalities.elementwise.quadratic.convex = 0; |
---|
| 284 | problem.constraint.inequalities.elementwise.quadratic.nonconvex = 0; |
---|
| 285 | problem.constraint.inequalities.elementwise.sigmonial = 0; |
---|
| 286 | problem.constraint.inequalities.elementwise.polynomial = 0; |
---|
| 287 | |
---|
| 288 | problem.constraint.inequalities.semidefinite.linear = problem.constraint.inequalities.semidefinite.linear | problem.constraint.inequalities.semidefinite.quadratic | problem.constraint.inequalities.semidefinite.polynomial | problem.constraint.inequalities.semidefinite.sigmonial; |
---|
| 289 | problem.constraint.inequalities.semidefinite.quadratic = 0; |
---|
| 290 | problem.constraint.inequalities.semidefinite.polynomial = 0; |
---|
| 291 | problem.constraint.inequalities.semidefinite.sigmonial = 0; |
---|
| 292 | |
---|
| 293 | poly_constraint = 0; |
---|
| 294 | bilin_constraint = 0; |
---|
| 295 | sigm_constraint = 0; |
---|
| 296 | end |
---|
| 297 | |
---|
| 298 | |
---|
| 299 | % Analyse the objective function |
---|
| 300 | quad_info = []; |
---|
| 301 | if (~isempty(h)) & ~is(h,'linear') &~(relax==1) &~(relax==3) |
---|
| 302 | if any(ismember(getvariables(h),sigmonial_variables)) |
---|
| 303 | problem.objective.sigmonial = 1; |
---|
| 304 | else |
---|
| 305 | [Q,c,f,x,info] = quaddecomp(h); |
---|
| 306 | if ~isreal(Q) % Numerical noise common on imaginary parts |
---|
| 307 | Qr = real(Q); |
---|
| 308 | Qi = imag(Q); |
---|
| 309 | Qr(abs(Qr)<1e-10) = 0; |
---|
| 310 | Qi(abs(Qi)<1e-10) = 0; |
---|
| 311 | cr = real(c); |
---|
| 312 | ci = imag(c); |
---|
| 313 | cr(abs(cr)<1e-10) = 0; |
---|
| 314 | ci(abs(ci)<1e-10) = 0; |
---|
| 315 | Q = Qr + sqrt(-1)*Qi; |
---|
| 316 | c = cr + sqrt(-1)*ci; |
---|
| 317 | end |
---|
| 318 | if info==0 |
---|
| 319 | [R,p]=chol(Q); |
---|
| 320 | if p~=0 |
---|
| 321 | Q = full(Q); |
---|
| 322 | if min(eig(Q))>=-1e-10 |
---|
| 323 | p=0; |
---|
| 324 | try |
---|
| 325 | [U,S,V]=svd(Q); |
---|
| 326 | catch |
---|
| 327 | [U,S,V]=svd(full(Q)); |
---|
| 328 | end |
---|
| 329 | i = find(diag(S)>1e-10); |
---|
| 330 | R = sqrt(S(1:max(i),1:max(i)))*V(:,1:max(i))'; |
---|
| 331 | end |
---|
| 332 | end |
---|
| 333 | if p==0 |
---|
| 334 | problem.objective.quadratic.convex = 1; |
---|
| 335 | else |
---|
| 336 | problem.objective.quadratic.nonconvex = 1; |
---|
| 337 | end |
---|
| 338 | quad_info.Q = Q; |
---|
| 339 | quad_info.c = c; |
---|
| 340 | quad_info.f = f; |
---|
| 341 | quad_info.x = x; |
---|
| 342 | quad_info.R = R; |
---|
| 343 | quad_info.p = p; |
---|
| 344 | else |
---|
| 345 | problem.objective.polynomial = 1; |
---|
| 346 | end |
---|
| 347 | end |
---|
| 348 | else |
---|
| 349 | problem.objective.linear = ~isempty(h); |
---|
| 350 | end |
---|