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 |
---|