[37] | 1 | function [F,failure,cause] = expandmodel(F,h,options) |
---|
| 2 | % Author Johan Löfberg |
---|
| 3 | % $Id: expandmodel.m,v 1.66 2006/09/13 09:28:51 joloef Exp $ |
---|
| 4 | |
---|
| 5 | % FIX : Current code experimental, complex, conservative, has issues with |
---|
| 6 | % nonlinearities and is slow... |
---|
| 7 | |
---|
| 8 | % All extended variables in the problem. It is expensive to extract this |
---|
| 9 | % one so we will keep it and pass it along in the recursion |
---|
| 10 | extendedvariables = yalmip('extvariables'); |
---|
| 11 | |
---|
| 12 | % Assume success |
---|
| 13 | failure = 0; |
---|
| 14 | cause = ''; |
---|
| 15 | |
---|
| 16 | % Early bail out |
---|
| 17 | if isempty(extendedvariables) |
---|
| 18 | return |
---|
| 19 | end |
---|
| 20 | |
---|
| 21 | % Check if it already has ben expanded |
---|
| 22 | already_expanded = expanded(F); |
---|
| 23 | if all(already_expanded) |
---|
| 24 | if isempty(setdiff(getvariables(h),expanded(h))) |
---|
| 25 | return |
---|
| 26 | end |
---|
| 27 | end |
---|
| 28 | |
---|
| 29 | % Extract all simple bounds from the model, and update the internal bounds |
---|
| 30 | % in YALMIP. This is done in order to get as tighter big-M models |
---|
| 31 | if ~isempty(F) |
---|
| 32 | nv = yalmip('nvars'); |
---|
| 33 | yalmip('setbounds',1:nv,repmat(-inf,nv,1),repmat(inf,nv,1)); |
---|
| 34 | LU = getbounds(F); |
---|
| 35 | yalmip('setbounds',1:nv,LU(:,1),LU(:,2)); |
---|
| 36 | end |
---|
| 37 | |
---|
| 38 | % Temporary hack to deal with a bug in CPLEX. For the implies operator (and |
---|
| 39 | % some more) YALMIP creates a dummy variable x with set(x==1). Cplex fails |
---|
| 40 | % to solve problem with these stupid variables kept, hence we need to |
---|
| 41 | % remove these variables and constraints... |
---|
| 42 | global MARKER_VARIABLES |
---|
| 43 | MARKER_VARIABLES = []; |
---|
| 44 | |
---|
| 45 | % Temporary hack to deal with geometric programs. GPs are messy here, |
---|
| 46 | % becasue we can by mistake claim nonconvexity, since we may have no |
---|
| 47 | % sigmonial terms but indefinite quadratic term, but the whole problem is |
---|
| 48 | % meant to be solved using a GP solver. YES, globals suck, but this is |
---|
| 49 | % only temporary...hrm. |
---|
| 50 | global DUDE_ITS_A_GP |
---|
| 51 | DUDE_ITS_A_GP = 0; |
---|
| 52 | |
---|
| 53 | % Keep track of expressions that already have been modelled. Note that if a |
---|
| 54 | % graph-model already has been constructed but we now require a milp, for |
---|
| 55 | % numerical reasons, we should remove the old graph descriptions (important |
---|
| 56 | % for MPT models in particular) |
---|
| 57 | % FIX: Pre-parse the whole problem etc (solves the issues with GP also) |
---|
| 58 | global ALREADY_MODELLED |
---|
| 59 | global REMOVE_THESE_IN_THE_END |
---|
| 60 | ALREADY_MODELLED = {}; |
---|
| 61 | REMOVE_THESE_IN_THE_END = []; |
---|
| 62 | |
---|
| 63 | % Nonlinear operator variables are not allowed to be used in polynomial |
---|
| 64 | % expressions, except if they are exactly modelled, i.e. modelled using |
---|
| 65 | % MILP models. We will expand the model and collect variables that are in |
---|
| 66 | % polynomials, and check in the end if they are exaclty modelled |
---|
| 67 | global OPERATOR_IN_POLYNOM |
---|
| 68 | OPERATOR_IN_POLYNOM = []; |
---|
| 69 | |
---|
| 70 | % All variable indicies used in the problem |
---|
| 71 | v1 = getvariables(F); |
---|
| 72 | v2 = depends(F); |
---|
| 73 | v3 = getvariables(h); |
---|
| 74 | v4 = depends(h); |
---|
| 75 | |
---|
| 76 | % Speed-hack for LARGE-scale dualizations |
---|
| 77 | if isequal(v3,v4) & isequal(v1,v2) |
---|
| 78 | variables = uniquestripped([v1 v3]); |
---|
| 79 | else |
---|
| 80 | variables = uniquestripped([v1 v2 v3 v4]); |
---|
| 81 | end |
---|
| 82 | |
---|
| 83 | % Index to variables modeling operators |
---|
| 84 | extended = find(ismembc(variables,extendedvariables)); |
---|
| 85 | |
---|
| 86 | if nargin < 3 |
---|
| 87 | options = sdpsettings; |
---|
| 88 | end |
---|
| 89 | |
---|
| 90 | % This is a tweak to allow epxansion of bilinear terms in robust problems, |
---|
| 91 | % is expression such as abs(x*w) < 1 for all -1 < w < 1 |
---|
| 92 | % This field is set to 1 in robustify and tells YALMIP to skip checking for |
---|
| 93 | % polynomial nonconvexity in the propagation |
---|
| 94 | if ~isfield(options,'expandbilinear') |
---|
| 95 | options.expandbilinear = 0; |
---|
| 96 | end |
---|
| 97 | |
---|
| 98 | % Monomial information. Expensive to retrieve, so we pass this along |
---|
| 99 | [monomtable,variabletype] = yalmip('monomtable'); |
---|
| 100 | |
---|
| 101 | % Is this trivially a GP, or meant to be at least? |
---|
| 102 | if strcmpi(options.solver,'gpposy') | strcmpi(options.solver,'fmincon-geometric') | strcmpi(options.solver,'mosek-geometric') |
---|
| 103 | DUDE_ITS_A_GP = 1; |
---|
| 104 | else |
---|
| 105 | if ~isequal(options.solver,'fmincon') & ~isequal(options.solver,'') & ~isequal(options.solver,'mosek') |
---|
| 106 | % User has specified some other solver, which does not |
---|
| 107 | % support GPs, hence it cannot be intended to be a GP |
---|
| 108 | DUDE_ITS_A_GP = 0; |
---|
| 109 | else |
---|
| 110 | % Check to see if there are any sigmonial terms on top-level |
---|
| 111 | DUDE_ITS_A_GP = ~isempty(find(variabletype(variables) == 4)); |
---|
| 112 | end |
---|
| 113 | end |
---|
| 114 | |
---|
| 115 | % Constraints generated during recursive process to model operators |
---|
| 116 | F_expand = set([]); |
---|
| 117 | |
---|
| 118 | if isempty(F) |
---|
| 119 | F = set([]); |
---|
| 120 | end |
---|
| 121 | |
---|
| 122 | % First, check the objective |
---|
| 123 | variables = uniquestripped([depends(h) getvariables(h)]); |
---|
| 124 | monomtable = monomtable(:,extendedvariables); |
---|
| 125 | |
---|
| 126 | % However, some of the variables are already expanded (expand can be called |
---|
| 127 | % sequentially from solvemp and solverobust) |
---|
| 128 | variables = setdiff(variables,expanded(h)); |
---|
| 129 | |
---|
| 130 | % Determine if we should aim for MILP model directly |
---|
| 131 | if options.allowmilp == 2 |
---|
| 132 | method = 'milp'; |
---|
| 133 | else |
---|
| 134 | method = 'graph'; |
---|
| 135 | end |
---|
| 136 | |
---|
| 137 | if DUDE_ITS_A_GP == 1 |
---|
| 138 | options.allowmilp = 0; |
---|
| 139 | method = 'graph'; |
---|
| 140 | end |
---|
| 141 | |
---|
| 142 | % ************************************************************************* |
---|
| 143 | % OK, looks good. Apply recursive expansion on the objective |
---|
| 144 | % ************************************************************************* |
---|
| 145 | index_in_extended = find(ismembc(variables,extendedvariables)); |
---|
| 146 | if ~isempty(index_in_extended) |
---|
| 147 | extstruct = yalmip('extstruct',variables(index_in_extended)); |
---|
| 148 | if ~isa(extstruct,'cell') |
---|
| 149 | extstruct = {extstruct}; |
---|
| 150 | end |
---|
| 151 | [F_expand,failure,cause] = expand(index_in_extended,variables,h,F_expand,extendedvariables,monomtable,'objective',0,options,method,extstruct); |
---|
| 152 | end |
---|
| 153 | |
---|
| 154 | % ************************************************************************* |
---|
| 155 | % Continue with constraints |
---|
| 156 | % ************************************************************************* |
---|
| 157 | constraint = 1; |
---|
| 158 | all_extstruct = yalmip('extstruct'); |
---|
| 159 | while constraint <=length(F) & ~failure |
---|
| 160 | if ~already_expanded(constraint) |
---|
| 161 | variables = uniquestripped([depends(F(constraint)) getvariables(F(constraint))]); |
---|
| 162 | [ix,jx,kx] = find(monomtable(variables,:)); |
---|
| 163 | if ~isempty(jx) % Bug in 6.1 |
---|
| 164 | if any(kx>1) |
---|
| 165 | OPERATOR_IN_POLYNOM = [OPERATOR_IN_POLYNOM extendedvariables(jx(find(kx>1)))]; |
---|
| 166 | end |
---|
| 167 | end |
---|
| 168 | |
---|
| 169 | index_in_extended = find(ismembc(variables,extendedvariables)); |
---|
| 170 | if ~isempty(index_in_extended) |
---|
| 171 | global_index = variables(index_in_extended); |
---|
| 172 | local_index = []; |
---|
| 173 | for i = 1:length(global_index) |
---|
| 174 | local_index = [local_index find(global_index(i) == extendedvariables)]; |
---|
| 175 | end |
---|
| 176 | extstruct = num2cell(all_extstruct(local_index)); |
---|
| 177 | if is(F(constraint),'equality') |
---|
| 178 | if options.allowmilp |
---|
| 179 | [F_expand,failure,cause] = expand(index_in_extended,variables,-sdpvar(F(constraint)),F_expand,extendedvariables,monomtable,['constraint #' num2str(constraint)],0,options,'milp',extstruct); |
---|
| 180 | else |
---|
| 181 | failure = 1; |
---|
| 182 | cause = ['MILP model required for equality in constraint #' num2str(constraint)]; |
---|
| 183 | end |
---|
| 184 | else |
---|
| 185 | [F_expand,failure,cause] = expand(index_in_extended,variables,-sdpvar(F(constraint)),F_expand,extendedvariables,monomtable,['constraint #' num2str(constraint)],0,options,method,extstruct); |
---|
| 186 | end |
---|
| 187 | end |
---|
| 188 | end |
---|
| 189 | constraint = constraint+1; |
---|
| 190 | end |
---|
| 191 | |
---|
| 192 | % ************************************************************************* |
---|
| 193 | % Temporary hack to fix the implies operator (cplex has some problem on |
---|
| 194 | % these trivial models where a variable only is used in x==1 |
---|
| 195 | % FIX: Automatically support this type of nonlinear operators |
---|
| 196 | % ************************************************************************* |
---|
| 197 | if ~isempty(MARKER_VARIABLES) |
---|
| 198 | MARKER_VARIABLES = sort(MARKER_VARIABLES); |
---|
| 199 | equalities = find(is(F,'equality')); |
---|
| 200 | equalities = equalities(:)'; |
---|
| 201 | remove = []; |
---|
| 202 | for j = equalities |
---|
| 203 | v = getvariables(F(j)); |
---|
| 204 | if length(v)==1 |
---|
| 205 | if ismembc(v,MARKER_VARIABLES) |
---|
| 206 | remove = [remove j]; |
---|
| 207 | end |
---|
| 208 | end |
---|
| 209 | end |
---|
| 210 | if ~isempty(remove) |
---|
| 211 | F(remove) = []; |
---|
| 212 | end |
---|
| 213 | end |
---|
| 214 | |
---|
| 215 | F_expand = lifted(F_expand,1); |
---|
| 216 | % ************************************************************************* |
---|
| 217 | % We are done. We might have generated some stuff more than once, but |
---|
| 218 | % luckily we keep track of these mistakes and remove them in the end (this |
---|
| 219 | % happens if we have constraints like set(max(x)<1) + set(max(x)>0) where |
---|
| 220 | % the first constraint would genrate a graph-model but the second set |
---|
| 221 | % creates a milp model. |
---|
| 222 | % ************************************************************************* |
---|
| 223 | if ~failure |
---|
| 224 | F = F + F_expand; |
---|
| 225 | if length(REMOVE_THESE_IN_THE_END) > 0 |
---|
| 226 | F = F(find(~ismember(getlmiid(F),REMOVE_THESE_IN_THE_END))); |
---|
| 227 | end |
---|
| 228 | end |
---|
| 229 | |
---|
| 230 | % ************************************************************************* |
---|
| 231 | % Normally, operators are not allowed in polynomial expressions. We do |
---|
| 232 | % however allow this if the variable has been modelled with an exact MILP |
---|
| 233 | % model. |
---|
| 234 | % ************************************************************************* |
---|
| 235 | Final_model = {ALREADY_MODELLED{unique(OPERATOR_IN_POLYNOM)}}; |
---|
| 236 | for i = 1:length(Final_model) |
---|
| 237 | if ~(strcmp(Final_model{i}.method,'milp') | strcmp(Final_model{i}.method,'none') | options.allownonconvex) |
---|
| 238 | failure = 1; |
---|
| 239 | cause = 'Nonlinear operator in polynomial expression.'; |
---|
| 240 | return |
---|
| 241 | end |
---|
| 242 | end |
---|
| 243 | |
---|
| 244 | % declare this model as expanded |
---|
| 245 | F = expanded(F,1); |
---|
| 246 | |
---|
| 247 | function [F_expand,failure,cause] = expand(index_in_extended,variables,expression,F_expand,extendedvariables,monomtable,where,level,options,method,extstruct) |
---|
| 248 | global DUDE_ITS_A_GP ALREADY_MODELLED REMOVE_THESE_IN_THE_END OPERATOR_IN_POLYNOM |
---|
| 249 | |
---|
| 250 | % ************************************************************************* |
---|
| 251 | % Go through all parts of expression to check for convexity/concavity |
---|
| 252 | % First, a small gateway function before calling the recursive stuff |
---|
| 253 | % ************************************************************************* |
---|
| 254 | if ~DUDE_ITS_A_GP |
---|
| 255 | [ix,jx,kx] = find(monomtable(variables,:)); |
---|
| 256 | if ~isempty(jx) % Bug in 6.1 |
---|
| 257 | if any(kx>1) |
---|
| 258 | OPERATOR_IN_POLYNOM = [OPERATOR_IN_POLYNOM extendedvariables(jx(find(kx>1)))]; |
---|
| 259 | end |
---|
| 260 | end |
---|
| 261 | end |
---|
| 262 | |
---|
| 263 | failure = 0; |
---|
| 264 | j = 1; |
---|
| 265 | |
---|
| 266 | while j<=length(index_in_extended) & ~failure |
---|
| 267 | i = index_in_extended(j); |
---|
| 268 | basis = getbasematrix(expression,variables(i)); |
---|
| 269 | if all(basis >= 0) |
---|
| 270 | [F_expand,failure,cause] = expandrecursive(recover(variables(i)),F_expand,extendedvariables,monomtable,where,level+1,options,method,extstruct{j},'convex'); |
---|
| 271 | else |
---|
| 272 | [F_expand,failure,cause] = expandrecursive(recover(variables(i)),F_expand,extendedvariables,monomtable,where,level+1,options,method,extstruct{j},'concave'); |
---|
| 273 | end |
---|
| 274 | j=j+1; |
---|
| 275 | end |
---|