1 | function [sol,m,Q,residuals,everything] = solvesos(F,obj,options,params,candidateMonomials) |
---|
2 | %SOLVESOS Sum of squares decomposition |
---|
3 | % |
---|
4 | % [sol,m,B,residual] = solvesos(F,h,options,params,monomials) is used |
---|
5 | % for finding SOS decompositions of polynomials. |
---|
6 | % |
---|
7 | % The coefficients of the polynomials are assumed linear w.r.t a set of |
---|
8 | % decision variables params and polynomial with respect to a variable x. |
---|
9 | % |
---|
10 | % An extension with a nonlinear parameterization in params is possible. |
---|
11 | % Note though that this gives BMIs or PMIs, solvable (locally) only if |
---|
12 | % PENBMI is installed, or by specifying 'moment' as solver to try to |
---|
13 | % solve the nonconvex semidefinite programming problem using a |
---|
14 | % semidefinite relaxation based on moments. |
---|
15 | % |
---|
16 | % The SOS problem can be formulated as |
---|
17 | % |
---|
18 | % min h(params) |
---|
19 | % |
---|
20 | % subject to F(i) >(=) 0 or F(i) is SOS w.r.t x |
---|
21 | % |
---|
22 | % INPUT |
---|
23 | % F : SET object with SOS constrained polynomials and constraints on variables params |
---|
24 | % h : scalar SDPVAR object (can be []) |
---|
25 | % options : options structure obtained from SDPSETTINGS (can be []) |
---|
26 | % params : SDPVAR object defining parametric variables (can be []) |
---|
27 | % monomials : SDPVAR object with user-specified monomials for decomposition (can be []) |
---|
28 | % |
---|
29 | % OUTPUT |
---|
30 | % sol : Solution diagnostic from SDP problem |
---|
31 | % v : Cell with monomials used in decompositions |
---|
32 | % Q : Cell with Gram matrices, p = v{i}'*Q{i}*v{i}, where p is the ith SOS polynomial in your model. |
---|
33 | % residuals : Mismatch between p and decompositions. Same values (modulo numerical issue) as checkset(find(is(F,'sos'))) |
---|
34 | % Warning, these residuals are not computed on matrix sos-of-squares |
---|
35 | % |
---|
36 | % EXAMPLE |
---|
37 | % x = sdpvar(1);solvesos(set(sos(x^4+x^3+1))); % Simple decompositions |
---|
38 | % x = sdpvar(1);t = sdpvar(1);solvesos(set(sos(x^4+x^3+1-t)),-t); % Lower bound by maximizing t |
---|
39 | % |
---|
40 | % NOTES |
---|
41 | % |
---|
42 | % Variables not part of params, but part of non-SOS constraints in F |
---|
43 | % or objective h will automatically be appended to the params list. |
---|
44 | % |
---|
45 | % To extract SOS decomposition, use the command SOSD (or compute from use v and Q) |
---|
46 | % |
---|
47 | % If the 5th input argument is used, no additional monomial reduction is |
---|
48 | % performed (Newton, inconstency, congruence). It is thus assumed that |
---|
49 | % the supplied candidate monomials constitute a sufficient basis. |
---|
50 | % |
---|
51 | % The field options.sos can be used to tune the SOS-calculations. See HTML help for details |
---|
52 | % |
---|
53 | % sos.model - Kernel or image representation of SOS problem [0|1|2 (0)] |
---|
54 | % sos.newton - Use Newton polytope to reduce size [0|1 (1)] |
---|
55 | % sos.congruence - Block-diagonalize using congruence classes [0|1|2 (2)] |
---|
56 | % sos.scale - Scale polynomial [0|1 (1)] |
---|
57 | % sos.numblkdg - Try to perform a-posteriori block-diagonalization [real (0)] |
---|
58 | % sos.inconsistent - Remove diagonal-inconsistent monomials [0|1|2 (0)] |
---|
59 | % sos.clean - Remove monomials with coefficients < clean [real > 0 (1e-4)] |
---|
60 | % sos.traceobj - Minimize trace of Gram matrix in problems without objective function [0|1 (0)] |
---|
61 | % sos.extlp - Extract simple translated LP cones when performing dualization [0|1 (1)] |
---|
62 | % |
---|
63 | % See also SOS, SOSD, SDPSETTINGS, SOLVEMOMENT, SDPVAR, SDISPLAY |
---|
64 | |
---|
65 | %% Time YALMIP |
---|
66 | yalmip_time = clock; |
---|
67 | |
---|
68 | % ************************************************ |
---|
69 | %% Check #inputs |
---|
70 | % ************************************************ |
---|
71 | if nargin<5 |
---|
72 | candidateMonomials = []; |
---|
73 | if nargin<4 |
---|
74 | params = []; |
---|
75 | if nargin<3 |
---|
76 | options = sdpsettings; |
---|
77 | if nargin<2 |
---|
78 | obj = []; |
---|
79 | if nargin<1 |
---|
80 | help solvesos |
---|
81 | return |
---|
82 | end |
---|
83 | end |
---|
84 | end |
---|
85 | end |
---|
86 | end |
---|
87 | |
---|
88 | if isempty(options) |
---|
89 | options = sdpsettings; |
---|
90 | end |
---|
91 | |
---|
92 | % Lazy syntax (not official...) |
---|
93 | if nargin==1 & isa(F,'sdpvar') |
---|
94 | F = set(sos(F)); |
---|
95 | end |
---|
96 | |
---|
97 | if ~isempty(options) |
---|
98 | if options.sos.numblkdg |
---|
99 | [sol,m,Q,residuals,everything] = solvesos_find_blocks(F,obj,options,params,candidateMonomials); |
---|
100 | return |
---|
101 | end |
---|
102 | end |
---|
103 | |
---|
104 | % ************************************************************************* |
---|
105 | %% Extract all SOS constraints and candidate monomials |
---|
106 | % ************************************************************************* |
---|
107 | if ~any(is(F,'sos')) |
---|
108 | error('At-least one constraint should be an SOS constraints!'); |
---|
109 | end |
---|
110 | p = []; |
---|
111 | ranks = []; |
---|
112 | for i = 1:length(F) |
---|
113 | if is(F(i),'sos') |
---|
114 | pi = sdpvar(F(i)); |
---|
115 | p{end+1} = pi; |
---|
116 | ranks(end+1) = getsosrank(pi); % Desired rank : Experimental code |
---|
117 | end |
---|
118 | end |
---|
119 | if isempty(candidateMonomials) |
---|
120 | for i = 1:length(F) |
---|
121 | candidateMonomials{i}=[]; |
---|
122 | end |
---|
123 | elseif isa(candidateMonomials,'sdpvar') |
---|
124 | cM=candidateMonomials; |
---|
125 | candidateMonomials={}; |
---|
126 | for i = 1:length(p) |
---|
127 | candidateMonomials{i}=cM; |
---|
128 | end |
---|
129 | elseif isa(candidateMonomials,'cell') |
---|
130 | if length(p)~=length(candidateMonomials) |
---|
131 | error('Dimension mismatch between the candidate monomials and the number of SOS constraints'); |
---|
132 | end |
---|
133 | end |
---|
134 | |
---|
135 | % ************************************************************************* |
---|
136 | %% Get the parametric constraints |
---|
137 | % ************************************************************************* |
---|
138 | F_original = F; |
---|
139 | F_parametric = F(find(~is(F,'sos'))); |
---|
140 | if isempty(F_parametric) |
---|
141 | F_parametric = set([]); |
---|
142 | end |
---|
143 | |
---|
144 | % ************************************************************************* |
---|
145 | %% Expand the parametric constraints |
---|
146 | % ************************************************************************* |
---|
147 | if ~isempty(yalmip('extvariables')) |
---|
148 | [F_parametric,failure] = expandmodel(F_parametric,obj,options); |
---|
149 | F_parametric = expanded(F_parametric,1); |
---|
150 | obj = expanded(obj,1); |
---|
151 | if failure |
---|
152 | Q{1} = [];m{1} = [];residuals = [];everything = []; |
---|
153 | sol.yalmiptime = etime(clock,yalmip_time); |
---|
154 | sol.solvertime = 0; |
---|
155 | sol.info = yalmiperror(14,'YALMIP'); |
---|
156 | sol.problem = 14; |
---|
157 | end |
---|
158 | end |
---|
159 | |
---|
160 | if ~isempty(params) |
---|
161 | if ~isa(params,'sdpvar') |
---|
162 | error('Fourth argment should be a SDPVAR variable or empty') |
---|
163 | end |
---|
164 | end |
---|
165 | |
---|
166 | % ************************************************************************* |
---|
167 | % Collect all possible parametric variables |
---|
168 | % ************************************************************************* |
---|
169 | ParametricVariables = uniquestripped([depends(obj) depends(F_parametric) depends(params)]); |
---|
170 | |
---|
171 | if options.verbose>0; |
---|
172 | disp('-------------------------------------------------------------------------'); |
---|
173 | disp('YALMIP SOS module started...'); |
---|
174 | disp('-------------------------------------------------------------------------'); |
---|
175 | end |
---|
176 | |
---|
177 | % ************************************************************************* |
---|
178 | %% INITIALIZE SOS-DECOMPOSITIONS SDP CONSTRAINTS |
---|
179 | % ************************************************************************* |
---|
180 | F_sos = set([]); |
---|
181 | |
---|
182 | % ************************************************************************* |
---|
183 | %% FIGURE OUT ALL USED PARAMETRIC VARIABLES |
---|
184 | % ************************************************************************* |
---|
185 | AllVariables = uniquestripped([depends(obj) depends(F_original) depends(F_parametric)]); |
---|
186 | ParametricVariables = intersect(ParametricVariables,AllVariables); |
---|
187 | MonomVariables = setdiff(AllVariables,ParametricVariables); |
---|
188 | params = recover(ParametricVariables); |
---|
189 | if isempty(MonomVariables) |
---|
190 | error('No independent variables? Perhaps you added a constraint set(p(x)) when you meant set(sos(p(x)))'); |
---|
191 | end |
---|
192 | if options.verbose>0;disp(['Detected ' num2str(length(ParametricVariables)) ' parametric variables and ' num2str(length(MonomVariables)) ' independent variables.']);end |
---|
193 | |
---|
194 | % ************************************************ |
---|
195 | %% ANY BMI STUFF |
---|
196 | % ************************************************ |
---|
197 | NonLinearParameterization = 0; |
---|
198 | if ~isempty(ParametricVariables) |
---|
199 | monomtable = yalmip('monomtable'); |
---|
200 | ParametricMonomials = monomtable(uniquestripped([getvariables(obj) getvariables(F_original)]),ParametricVariables); |
---|
201 | if any(sum(abs(ParametricMonomials),2)>1) |
---|
202 | NonLinearParameterization = 1; |
---|
203 | end |
---|
204 | end |
---|
205 | |
---|
206 | % ************************************************ |
---|
207 | %% ANY INTEGER DATA |
---|
208 | % ************************************************ |
---|
209 | IntegerData = 0; |
---|
210 | if ~isempty(ParametricVariables) |
---|
211 | globalInteger = [yalmip('binvariables') yalmip('intvariables')]; |
---|
212 | integerVariables = getvariables(F_parametric(find(is(F_parametric,'binary') | is(F_parametric,'integer')))); |
---|
213 | integerVariables = [integerVariables intersect(ParametricVariables,globalInteger)]; |
---|
214 | integerVariables = intersect(integerVariables,ParametricVariables); |
---|
215 | IntegerData = ~isempty(integerVariables); |
---|
216 | end |
---|
217 | |
---|
218 | % ************************************************ |
---|
219 | %% ANY UNCERTAIN DATA |
---|
220 | % ************************************************ |
---|
221 | UncertainData = 0; |
---|
222 | if ~isempty(ParametricVariables) |
---|
223 | UncertainData = any(is(F_parametric,'uncertain')); |
---|
224 | end |
---|
225 | |
---|
226 | % ************************************************ |
---|
227 | %% DISPLAY WHAT WE FOUND |
---|
228 | % ************************************************ |
---|
229 | if options.verbose>0 & ~isempty(F_parametric) |
---|
230 | nLP = 0; |
---|
231 | nEQ = 0; |
---|
232 | nLMI = sum(full(is(F_parametric,'lmi')) & full(~is(F_parametric,'element-wise'))); %FULL due to bug in ML 7.0.1 |
---|
233 | for i = 1:length(F_parametric) |
---|
234 | if is(F_parametric,'element-wise') |
---|
235 | nLP = nLP + prod(size(F_parametric(i))); |
---|
236 | end |
---|
237 | if is(F_parametric,'equality') |
---|
238 | nEQ = nEQ + prod(size(F_parametric(i))); |
---|
239 | end |
---|
240 | end |
---|
241 | disp(['Detected ' num2str(full(nLP)) ' linear inequalities, ' num2str(full(nEQ)) ' equality constraints and ' num2str(full(nLMI)) ' LMIs.']); |
---|
242 | end |
---|
243 | |
---|
244 | % ************************************************ |
---|
245 | %% IMAGE OR KERNEL REPRESENTATION? |
---|
246 | % ************************************************ |
---|
247 | noRANK = all(isinf(ranks)); |
---|
248 | switch options.sos.model |
---|
249 | case 0 |
---|
250 | constraint_classes = constraintclass(F); |
---|
251 | noCOMPLICATING = ~any(ismember([7 8 9 10 12 13 14 15],constraint_classes)); |
---|
252 | if noCOMPLICATING & ~NonLinearParameterization & noRANK & ~IntegerData |
---|
253 | options.sos.model = 1; |
---|
254 | if options.verbose>0;disp('Using kernel representation (options.sos.model=1).');end |
---|
255 | else |
---|
256 | if NonLinearParameterization |
---|
257 | if options.verbose>0;disp('Using image representation (options.sos.model=2). Nonlinear parameterization found');end |
---|
258 | elseif ~noRANK |
---|
259 | if options.verbose>0;disp('Using image representation (options.sos.model=2). SOS-rank constraint was found.');end |
---|
260 | elseif IntegerData |
---|
261 | if options.verbose>0;disp('Using image representation (options.sos.model=2). Integrality constraint was found.');end |
---|
262 | elseif UncertainData |
---|
263 | if options.verbose>0;disp('Using image representation (options.sos.model=2). Uncertain data was found.');end |
---|
264 | else |
---|
265 | if options.verbose>0;disp('Using image representation (options.sos.model=2). Integer data, KYPs or similar was found.');end |
---|
266 | end |
---|
267 | options.sos.model = 2; |
---|
268 | end |
---|
269 | case 1 |
---|
270 | if NonLinearParameterization |
---|
271 | if options.verbose>0;disp('Switching to image model due to nonlinear parameterization (not supported in kernel model).');end |
---|
272 | options.sos.model = 2; |
---|
273 | end |
---|
274 | if ~noRANK |
---|
275 | if options.verbose>0;disp('Switching to image model due to SOS-rank constraints (not supported in kernel model).');end |
---|
276 | options.sos.model = 2; |
---|
277 | end |
---|
278 | if IntegerData |
---|
279 | if options.verbose>0;disp('Switching to image model due to integrality constraints (not supported in kernel model).');end |
---|
280 | options.sos.model = 2; |
---|
281 | end |
---|
282 | case 3 |
---|
283 | otherwise |
---|
284 | end |
---|
285 | |
---|
286 | if ~isempty(yalmip('extvariables')) & options.sos.model == 2 & nargin<4 |
---|
287 | disp(' ') |
---|
288 | disp('**Using nonlinear operators in SOS problems can cause problems.') |
---|
289 | disp('**Please specify all parametric variables using the fourth argument'); |
---|
290 | disp(' '); |
---|
291 | end |
---|
292 | |
---|
293 | % ************************************************ |
---|
294 | %% SKIP DIAGONAL INCONSISTENCY FOR PARAMETRIC MODEL |
---|
295 | % ************************************************ |
---|
296 | if ~isempty(params) & options.sos.inconsistent |
---|
297 | if options.verbose>0;disp('Turning off inconsistency based reduction (not supported in parametric models).');end |
---|
298 | options.sos.inconsistent = 0; |
---|
299 | end |
---|
300 | |
---|
301 | % ************************************************ |
---|
302 | %% INITIALIZE OBJECTIVE |
---|
303 | % ************************************************ |
---|
304 | if ~isempty(obj) |
---|
305 | options.sos.traceobj = 0; |
---|
306 | end |
---|
307 | parobj = obj; |
---|
308 | obj = []; |
---|
309 | |
---|
310 | % ************************************************ |
---|
311 | %% SCALE SOS CONSTRAINTS |
---|
312 | % ************************************************ |
---|
313 | if options.sos.scale |
---|
314 | for constraint = 1:length(p) |
---|
315 | normp(constraint) = sqrt(norm(full(getbase(p{constraint})))); |
---|
316 | p{constraint} = p{constraint}/normp(constraint); |
---|
317 | sizep(constraint) = size(p{constraint},1); |
---|
318 | end |
---|
319 | else |
---|
320 | normp = ones(length(p),1); |
---|
321 | end |
---|
322 | |
---|
323 | % ************************************************ |
---|
324 | %% Some stuff not supported for |
---|
325 | % matrix valued SOS yet, turn off for safety |
---|
326 | % ************************************************ |
---|
327 | for constraint = 1:length(p) |
---|
328 | sizep(constraint) = size(p{constraint},1); |
---|
329 | end |
---|
330 | if any(sizep>1) |
---|
331 | options.sos.postprocess = 0; |
---|
332 | options.sos.reuse = 0; |
---|
333 | end |
---|
334 | |
---|
335 | % ************************************************ |
---|
336 | %% SKIP CONGRUENCE REDUCTION WHEN SOS-RANK |
---|
337 | % ************************************************ |
---|
338 | if ~all(isinf(ranks)) |
---|
339 | options.sos.congruence = 0; |
---|
340 | end |
---|
341 | |
---|
342 | % ************************************************ |
---|
343 | %% Create an LP model to speed up things in Newton |
---|
344 | % polytope reduction |
---|
345 | % ************************************************ |
---|
346 | if options.sos.newton |
---|
347 | temp=sdpvar(1,1); |
---|
348 | tempops = options; |
---|
349 | tempops.solver = 'cdd,glpk,*'; % CDD is generally robust on these problems |
---|
350 | tempops.verbose = 0; |
---|
351 | tempops.saveduals = 0; |
---|
352 | [aux1,aux2,aux3,LPmodel] = export(set(temp>0),temp,tempops); |
---|
353 | else |
---|
354 | LPmodel = []; |
---|
355 | end |
---|
356 | |
---|
357 | |
---|
358 | % ************************************************ |
---|
359 | %% LOOP THROUGH ALL SOS CONSTRAINTS |
---|
360 | % ************************************************ |
---|
361 | for constraint = 1:length(p) |
---|
362 | % ********************************************* |
---|
363 | %% FIND THE VARIABLES IN p, SORT, UNIQUE ETC |
---|
364 | % ********************************************* |
---|
365 | if options.verbose>1;disp(['Creating SOS-description ' num2str(constraint) '/' num2str(length(p)) ]);end |
---|
366 | |
---|
367 | pVariables = depends(p{constraint}); |
---|
368 | AllVariables = uniquestripped([pVariables ParametricVariables]); |
---|
369 | MonomVariables = setdiff1D(pVariables,ParametricVariables); |
---|
370 | x = recover(MonomVariables); |
---|
371 | z = recover(AllVariables); |
---|
372 | MonomIndicies = find(ismember(AllVariables,MonomVariables)); |
---|
373 | ParametricIndicies = find(ismember(AllVariables,ParametricVariables)); |
---|
374 | |
---|
375 | if isempty(MonomIndicies) |
---|
376 | % This is the case set(sos(t)) where t is a parametric (matrix) variable |
---|
377 | % This used to create an error message befgore to avoid some silly |
---|
378 | % bug in the model generation. Creating this error message is |
---|
379 | % stupid, but at the same time I can not remember where the bug was |
---|
380 | % and I have no regression test for this case. To avoid |
---|
381 | % introducing same bug again by mistake, I create all data |
---|
382 | % specifically for this case |
---|
383 | previous_exponent_p_monoms = [];%exponent_p_monoms; |
---|
384 | n = length(p{constraint}); |
---|
385 | A_basis = getbase(sdpvar(n,n,'full'));d = find(triu(ones(n)));A_basis = A_basis(d,2:end); |
---|
386 | BlockedA{constraint} = {A_basis}; |
---|
387 | Blockedb{constraint} = p{constraint}(d); |
---|
388 | BlockedN{constraint} = {zeros(1,0)}; |
---|
389 | Blockedx{constraint} = x; |
---|
390 | Blockedvarchange{constraint}=zeros(1,0); |
---|
391 | continue |
---|
392 | % error('You have constraints of the type set(sos(f(parametric_variables))). Please use set(f(parametric_variables) > 0) instead') |
---|
393 | end |
---|
394 | |
---|
395 | % ********************************************* |
---|
396 | %% Express p in monimials and coefficients |
---|
397 | % ********************************************* |
---|
398 | [exponent_p,p_base] = getexponentbase(p{constraint},z); |
---|
399 | |
---|
400 | % ********************************************* |
---|
401 | %% Powers for user defined candidate monomials |
---|
402 | % (still experimental) |
---|
403 | % ********************************************* |
---|
404 | if ~all(cellfun('isempty',candidateMonomials)) |
---|
405 | exponent_c = []; |
---|
406 | if isa(candidateMonomials{constraint},'cell') |
---|
407 | for i = 1:length(candidateMonomials{constraint}) |
---|
408 | exponent_c{i} = getexponentbase(candidateMonomials{constraint}{i},z); |
---|
409 | exponent_c{i} = exponent_c{i}(:,MonomIndicies); |
---|
410 | end |
---|
411 | else |
---|
412 | exponent_c{1} = getexponentbase(candidateMonomials{constraint},z); |
---|
413 | exponent_c{1} = exponent_c{1}(:,MonomIndicies); |
---|
414 | end |
---|
415 | else |
---|
416 | exponent_c = []; |
---|
417 | end |
---|
418 | |
---|
419 | % ********************************************* |
---|
420 | %% STUPID PROBLEM WITH ODD HIGHEST POWER?... |
---|
421 | % ********************************************* |
---|
422 | if isempty(ParametricIndicies) |
---|
423 | max_degrees = max(exponent_p(:,MonomIndicies),[],1); |
---|
424 | bad_max = any(max_degrees-fix((max_degrees/2))*2); |
---|
425 | if bad_max |
---|
426 | for i = 1:length(p) |
---|
427 | Q{i}=[]; |
---|
428 | m{i}=[]; |
---|
429 | end |
---|
430 | residuals=[]; |
---|
431 | everything = []; |
---|
432 | sol.yalmiptime = etime(clock,yalmip_time); |
---|
433 | sol.solvertime = 0; |
---|
434 | sol.info = yalmiperror(1,'YALMIP'); |
---|
435 | sol.problem = 2; |
---|
436 | return |
---|
437 | end |
---|
438 | end |
---|
439 | |
---|
440 | % ********************************************* |
---|
441 | %% Can we make a smart variable change (no code) |
---|
442 | % ********************************************* |
---|
443 | exponent_p_monoms = exponent_p(:,MonomIndicies); |
---|
444 | varchange = ones(1,size(MonomIndicies,2)); |
---|
445 | |
---|
446 | % ********************************************* |
---|
447 | %% Unique monoms (copies due to parametric terms) |
---|
448 | % ********************************************* |
---|
449 | exponent_p_monoms = uniquesafe(exponent_p_monoms,'rows'); |
---|
450 | |
---|
451 | if options.sos.reuse & constraint > 1 & isequal(previous_exponent_p_monoms,exponent_p_monoms) |
---|
452 | % We don't have to do anything, candidate monomials can be-used |
---|
453 | if options.verbose>1;disp(['Re-using all candidate monomials (same problem structure)']);end |
---|
454 | else |
---|
455 | |
---|
456 | |
---|
457 | % ********************************************* |
---|
458 | %% PRUNE W.R.T USER DEFINED |
---|
459 | %********************************************* |
---|
460 | % if ~isempty(exponent_c) |
---|
461 | % hash = randn(size(exponent_m{1},2),1); |
---|
462 | % for i = 1:length(exponent_m) |
---|
463 | % hash = randn(size(exponent_m{i},2),1); |
---|
464 | % keep = find(ismember(exponent_m{i}*hash,exponent_c{1}*hash)); |
---|
465 | % if length(keep)>0 |
---|
466 | % exponent_m{i} = exponent_m{i}(keep,:); |
---|
467 | % else |
---|
468 | % exponent_m{i}=[]; |
---|
469 | % end |
---|
470 | % end |
---|
471 | % end |
---|
472 | |
---|
473 | % ********************************************* |
---|
474 | % User has supplied the whole candidate structure |
---|
475 | % Don't process this |
---|
476 | % ********************************************* |
---|
477 | if ~isempty(exponent_c) |
---|
478 | exponent_m{1} = []; |
---|
479 | N = {}; |
---|
480 | for i = 1:length(exponent_c) |
---|
481 | exponent_m{i} = [exponent_m{1};exponent_c{i}]; |
---|
482 | N{i,1} = exponent_c{i}; |
---|
483 | end |
---|
484 | else |
---|
485 | % ********************************************* |
---|
486 | %% CORRELATIVE SPARSITY PATTERN |
---|
487 | % ********************************************* |
---|
488 | [C,csclasses] = corrsparsity(exponent_p_monoms,options); |
---|
489 | |
---|
490 | % ********************************************* |
---|
491 | %% GENERATE MONOMIALS |
---|
492 | % ********************************************* |
---|
493 | exponent_m = monomialgeneration(exponent_p_monoms,csclasses); |
---|
494 | |
---|
495 | % ********************************************* |
---|
496 | %% REDUCE #of MONOMIALS |
---|
497 | % ********************************************* |
---|
498 | % Fix for matrix case, perform newton w.r.t |
---|
499 | % diagonal polynomials only. This can be |
---|
500 | % improved, but for now, keep it simple... |
---|
501 | n = length(p{constraint});diag_elements = 1:(n+1):n^2;used_diagonal = find(any(p_base(diag_elements,:),1)); |
---|
502 | exponent_p_monoms_diag = exponent_p(used_diagonal,MonomIndicies); |
---|
503 | exponent_m = monomialreduction(exponent_m,exponent_p_monoms_diag,options,csclasses,LPmodel); |
---|
504 | |
---|
505 | % ********************************************* |
---|
506 | %% BLOCK PARTITION THE MONOMIALS BY CONGRUENCE |
---|
507 | % ********************************************* |
---|
508 | N = congruenceblocks(exponent_m,exponent_p_monoms,options,csclasses); |
---|
509 | % ********************************************* |
---|
510 | %% REDUCE FURTHER BY EXPLOITING BLOCK-STRUCTURE |
---|
511 | % ********************************************* |
---|
512 | N = blockmonomialreduction(exponent_p_monoms_diag,N,options); |
---|
513 | |
---|
514 | end |
---|
515 | |
---|
516 | |
---|
517 | % ********************************************* |
---|
518 | %% PREPARE FOR SDP FORMULATION BY CALCULATING ALL |
---|
519 | % POSSIBLE MONOMIAL PRODUCS IN EACH BLOCK |
---|
520 | % ********************************************* |
---|
521 | [exponent_m2,N_unique] = monomialproducts(N); |
---|
522 | |
---|
523 | % ********************************************* |
---|
524 | %% CHECK FOR BUG/IDIOT PROBLEMS IN FIXED PROBLEM |
---|
525 | % ********************************************* |
---|
526 | if isempty(ParametricIndicies) |
---|
527 | if ~isempty(setdiff(exponent_p_monoms,N_unique(:,3:end),'rows')) |
---|
528 | for i = 1:length(p) |
---|
529 | Q{i} = []; |
---|
530 | m{i} = []; |
---|
531 | end |
---|
532 | residuals = [];everything = []; |
---|
533 | sol.problem = 2; |
---|
534 | sol.info = yalmiperror(1,'YALMIP'); |
---|
535 | warning('Problem is trivially infeasible (odd highest power?)'); |
---|
536 | return |
---|
537 | end |
---|
538 | end |
---|
539 | end |
---|
540 | |
---|
541 | previous_exponent_p_monoms = exponent_p_monoms; |
---|
542 | |
---|
543 | % ********************************************* |
---|
544 | %% GENERATE DATA FOR SDP FORMULATIONS |
---|
545 | % ********************************************* |
---|
546 | p_base_parametric = []; |
---|
547 | n = length(p{constraint}); |
---|
548 | for i=1:length(p{constraint}) |
---|
549 | for j = 1:length(p{constraint}) |
---|
550 | p_base_parametric = [p_base_parametric parameterizedbase(p{constraint}(i,j),z,params,ParametricIndicies,exponent_p,p_base((i-1)*n+j,:))]; |
---|
551 | end |
---|
552 | end |
---|
553 | [BlockedA{constraint},Blockedb{constraint}] = generate_kernel_representation_data(N,N_unique,exponent_m2,exponent_p,p{constraint},options,p_base_parametric,ParametricIndicies,MonomIndicies,constraint==1); |
---|
554 | |
---|
555 | % SAVE FOR LATER |
---|
556 | BlockedN{constraint} = N; |
---|
557 | Blockedx{constraint} = x; |
---|
558 | Blockedvarchange{constraint}=varchange; |
---|
559 | end |
---|
560 | |
---|
561 | % ********************************************* |
---|
562 | %% And now get the SDP formulations |
---|
563 | % |
---|
564 | % The code above has generated matrices A and b |
---|
565 | % in AQ == b(parametric) |
---|
566 | % |
---|
567 | % We use these to generate kernel or image models |
---|
568 | % ********************************************* |
---|
569 | sol.problem = 0; |
---|
570 | switch options.sos.model |
---|
571 | case 1 |
---|
572 | % Kernel model |
---|
573 | [F,obj,BlockedQ,Primal_matrices,Free_variables] = create_kernelmodel(BlockedA,Blockedb,F_parametric,parobj,options,[]); |
---|
574 | case 2 |
---|
575 | % Image model |
---|
576 | [F,obj,BlockedQ,sol] = create_imagemodel(BlockedA,Blockedb,F_parametric,parobj,options); |
---|
577 | |
---|
578 | case 3 |
---|
579 | % Un-official model to solve bilinearly parameterized SOS using SDPLR |
---|
580 | [F,obj,options] = create_lrmodel(BlockedA,Blockedb,F_parametric,parobj,options,ParametricVariables); |
---|
581 | |
---|
582 | otherwise |
---|
583 | end |
---|
584 | |
---|
585 | % Unofficial fifth output with pseudo-numerical data |
---|
586 | everything.BlockedA = BlockedA; |
---|
587 | everything.Blockedb = Blockedb; |
---|
588 | everything.F = F; |
---|
589 | everything.obj = obj; |
---|
590 | |
---|
591 | % % ********************************************** |
---|
592 | % %% SOLVE SDP |
---|
593 | % % ********************************************** |
---|
594 | if sol.problem == 0 |
---|
595 | if options.verbose > 0 |
---|
596 | disp(' '); |
---|
597 | end |
---|
598 | if all(isinf(ranks)) |
---|
599 | % Standard case |
---|
600 | %sol = solvesdp(F+set(-10<recover(depends(F)) < 10),obj,sdpsettings('bmibnb.maxit',200,'solver','bmibnb','bmibnb.lpred',1)) |
---|
601 | sol = solvesdp(F,obj,options); |
---|
602 | else |
---|
603 | % We have to alter the problem slightly if there are rank |
---|
604 | % constraints on the decompositions |
---|
605 | sol = solveranksos(F,obj,options,ranks,BlockedQ); |
---|
606 | end |
---|
607 | end |
---|
608 | |
---|
609 | % ********************************************** |
---|
610 | %% Recover solution (and rescale model+solution) |
---|
611 | % ********************************************** |
---|
612 | for constraint = 1:length(p) |
---|
613 | for i = 1:length(BlockedQ{constraint}) |
---|
614 | doubleQ{constraint}{i} = normp(constraint)*double(BlockedQ{constraint}{i}); |
---|
615 | end |
---|
616 | doubleb{constraint}=normp(constraint)*double(Blockedb{constraint}); |
---|
617 | end |
---|
618 | |
---|
619 | % ********************************************** |
---|
620 | %% Minor post-process |
---|
621 | % ********************************************** |
---|
622 | if all(sizep<=1) |
---|
623 | [doubleQ,residuals] = postprocesssos(BlockedA,doubleb,doubleQ,[],options); |
---|
624 | else |
---|
625 | residuals = 0; |
---|
626 | end |
---|
627 | |
---|
628 | % ********************************************** |
---|
629 | %% Safety check for bad solvers... |
---|
630 | % ********************************************** |
---|
631 | if any(residuals > 1e-3) & (sol.problem == 0) & options.verbose>0 |
---|
632 | disp(' '); |
---|
633 | disp('-> Although the solver indicates no problems,') |
---|
634 | disp('-> the residuals in the problem are really bad.') |
---|
635 | disp('-> My guess: the problem is probably infeasible.') |
---|
636 | disp('-> Make sure to check how well your decomposition') |
---|
637 | disp('-> matches your polynomial (see manual)') |
---|
638 | disp('-> You can also try to change the option sos.model') |
---|
639 | disp('-> or use another SDP solver.') |
---|
640 | disp(' '); |
---|
641 | end |
---|
642 | |
---|
643 | % ********************************************** |
---|
644 | %% Confused users. Primal dual kernel image?... |
---|
645 | % ********************************************** |
---|
646 | if options.verbose > 0 |
---|
647 | if sol.problem == 1 |
---|
648 | if options.sos.model == 1 |
---|
649 | disp(' ') |
---|
650 | disp('-> Solver reported infeasible dual problem.') |
---|
651 | disp('-> Your SOS problem is probably unbounded.') |
---|
652 | elseif options.sos.model == 2 |
---|
653 | disp(' ') |
---|
654 | disp('-> Solver reported infeasible primal problem.') |
---|
655 | disp('-> Your SOS problem is probably infeasible.') |
---|
656 | end |
---|
657 | elseif sol.problem == 2 |
---|
658 | if options.sos.model == 1 |
---|
659 | disp(' ') |
---|
660 | disp('-> Solver reported unboundness of the dual problem.') |
---|
661 | disp('-> Your SOS problem is probably infeasible.') |
---|
662 | elseif options.sos.model == 2 |
---|
663 | disp(' ') |
---|
664 | disp('-> Solver reported unboundness of the primal problem.') |
---|
665 | disp('-> Your SOS problem is probably unbounded.') |
---|
666 | end |
---|
667 | elseif sol.problem == 12 |
---|
668 | disp(' ') |
---|
669 | disp('-> Solver reported unboundness or infeasibility of the primal problem.') |
---|
670 | disp('-> Your SOS problem is probably unbounded.') |
---|
671 | end |
---|
672 | end |
---|
673 | |
---|
674 | % ********************************************** |
---|
675 | %% De-block |
---|
676 | % ********************************************** |
---|
677 | for constraint = 1:length(p) |
---|
678 | Qtemp = []; |
---|
679 | for i = 1:length(BlockedQ{constraint}) |
---|
680 | Qtemp = blkdiag(Qtemp,doubleQ{constraint}{i}); |
---|
681 | end |
---|
682 | Q{constraint} = full(Qtemp); |
---|
683 | end |
---|
684 | |
---|
685 | % ********************************************** |
---|
686 | %% Experimental code for declaring sparsity |
---|
687 | % ********************************************** |
---|
688 | if options.sos.impsparse == 1 |
---|
689 | somesmall = 0; |
---|
690 | for i = 1:length(BlockedQ) |
---|
691 | for j = 1:length(BlockedQ{i}) |
---|
692 | small = find(abs(double(BlockedQ{i}{j}))<options.sos.sparsetol); |
---|
693 | nullThese{i}{j} = small; |
---|
694 | if ~isempty(small) |
---|
695 | somesmall = 1; |
---|
696 | end |
---|
697 | end |
---|
698 | end |
---|
699 | if somesmall |
---|
700 | [F,obj,BlockedQ,Primal_matrices,Free_variables] = create_kernelmodel(BlockedA,Blockedb,F_parametric,parobj,options,nullThese); |
---|
701 | sol = solvesdp(F,obj,options); |
---|
702 | for constraint = 1:length(p) |
---|
703 | for i = 1:length(BlockedQ{constraint}) |
---|
704 | doubleQ{constraint}{i} = normp(constraint)*double(BlockedQ{constraint}{i}); |
---|
705 | end |
---|
706 | doubleb{constraint}=normp(constraint)*double(Blockedb{constraint}); |
---|
707 | end |
---|
708 | % for i = 1:length(doubleQ) |
---|
709 | % for j = 1:length(doubleQ{i}) |
---|
710 | % doubleQ{i}{j}(nullThese{i}{j}) = 0; |
---|
711 | % end |
---|
712 | % end |
---|
713 | % ********************************************** |
---|
714 | %% Post-process |
---|
715 | % ********************************************** |
---|
716 | [doubleQ,residuals] = postprocesssos(BlockedA,doubleb,doubleQ,nullThese,options); |
---|
717 | for constraint = 1:length(p) |
---|
718 | Qtemp = []; |
---|
719 | for i = 1:length(BlockedQ{constraint}) |
---|
720 | Qtemp = blkdiag(Qtemp,doubleQ{constraint}{i}); |
---|
721 | end |
---|
722 | Q{constraint} = Qtemp; |
---|
723 | end |
---|
724 | end |
---|
725 | end |
---|
726 | |
---|
727 | % ********************************************* |
---|
728 | %% EXTRACT DECOMPOSITION |
---|
729 | % ********************************************* |
---|
730 | switch sol.problem |
---|
731 | case {0,1,2,3,4,5} % Well, it didn't f**k up completely at-least |
---|
732 | |
---|
733 | % % If SDPLR was sucessful, the last dual should have rank 1 |
---|
734 | % Z = dual(F(end)); |
---|
735 | % i = 1e-16; |
---|
736 | % [R,fail] = chol(Z); |
---|
737 | % while fail & i<10 |
---|
738 | % i = i*100; |
---|
739 | % [R,fail] = chol(Z+i*eye(length(Z))); |
---|
740 | % end |
---|
741 | % Z = R(1,2:end)'; |
---|
742 | % assign(recover(ParametricVariables),Z); |
---|
743 | % start = 1; |
---|
744 | % for constraint = 1:length(p) |
---|
745 | % Qi = []; |
---|
746 | % for j = 1:length(BlockedA{constraint}) |
---|
747 | % Qi = blkdiag(Qi,dual(F(start))); |
---|
748 | % start = start + 1; |
---|
749 | % end |
---|
750 | % Q{constraint} = Qi; |
---|
751 | % end |
---|
752 | |
---|
753 | % ********************************************* |
---|
754 | %% GENERATE MONOMIALS IN SOS-DECOMPOSITION |
---|
755 | % ********************************************* |
---|
756 | for constraint = 1:length(p) |
---|
757 | |
---|
758 | if constraint > 1 & isequal(BlockedN{constraint},BlockedN{constraint-1}) & isequal(Blockedx{constraint},Blockedx{constraint-1}) & isequal(Blockedvarchange{constraint},Blockedvarchange{constraint-1}) & isequal(sizep(constraint),sizep(constraint-1)) |
---|
759 | monoms{constraint} = monoms{constraint-1}; |
---|
760 | else |
---|
761 | monoms{constraint} = []; |
---|
762 | totalN{constraint} = []; |
---|
763 | N = BlockedN{constraint}; |
---|
764 | x = Blockedx{constraint}; |
---|
765 | for i = 1:length(N) |
---|
766 | % Original variable |
---|
767 | for j = 1:size(N{i},1) |
---|
768 | N{i}(j,:)=N{i}(j,:).*Blockedvarchange{constraint}; |
---|
769 | end |
---|
770 | if isempty(N{i}) |
---|
771 | monoms{constraint} = [monoms{constraint};[]]; |
---|
772 | else |
---|
773 | mi = kron(eye(sizep(constraint)),recovermonoms(N{i},x)); |
---|
774 | monoms{constraint} = [monoms{constraint};mi]; |
---|
775 | end |
---|
776 | end |
---|
777 | if isempty(monoms{constraint}) |
---|
778 | monoms{constraint}=1; |
---|
779 | end |
---|
780 | end |
---|
781 | |
---|
782 | % For small negative eigenvalues |
---|
783 | % this is a good quick'n'dirty approximation |
---|
784 | % Improve...use shifted eigenvalues and chol or what ever... |
---|
785 | if ~any(any(isnan(Q{constraint}))) |
---|
786 | if isempty(Q{constraint}) |
---|
787 | Q{constraint}=0; |
---|
788 | h{constraint}=0; |
---|
789 | else |
---|
790 | usedVariables = find(any(Q{constraint},2)); |
---|
791 | if length(usedVariables)<length(Q{constraint}) |
---|
792 | Qpart = Q{constraint}(usedVariables,usedVariables); |
---|
793 | [U,S,V]=svd(Qpart); |
---|
794 | R = sqrt(S)*V'; |
---|
795 | h0 = R*monoms{constraint}(usedVariables); |
---|
796 | if isa(h0,'sdpvar') |
---|
797 | h{constraint} = clean(R*monoms{constraint}(usedVariables),options.sos.clean); |
---|
798 | h{constraint} = h{constraint}(find(h{constraint})); |
---|
799 | else |
---|
800 | h{constraint} = h0; |
---|
801 | end |
---|
802 | else |
---|
803 | [U,S,V]=svd(Q{constraint}); |
---|
804 | R = sqrt(S)*V'; |
---|
805 | h0 = R*monoms{constraint}; |
---|
806 | |
---|
807 | if isa(h0,'sdpvar') |
---|
808 | h{constraint} = clean(R*monoms{constraint},options.sos.clean); |
---|
809 | h{constraint} = h{constraint}(find(sum(h{constraint},2)),:); |
---|
810 | else |
---|
811 | h{constraint} = h0; |
---|
812 | end |
---|
813 | end |
---|
814 | end |
---|
815 | if isempty(ParametricVariables) |
---|
816 | ParametricVariables = []; |
---|
817 | end |
---|
818 | setsos(p{constraint},h{constraint},ParametricVariables,Q{constraint},monoms{constraint}); |
---|
819 | else |
---|
820 | if options.verbose>0; |
---|
821 | if UncertainData |
---|
822 | disp(' '); |
---|
823 | disp('-> Only partial decomposition is returned (since you have uncertain data).'); |
---|
824 | disp(''); |
---|
825 | else |
---|
826 | disp(' '); |
---|
827 | disp('-> FAILURE : SOS decomposition not available.'); |
---|
828 | disp('-> The reason is probably that you are using a solver that does not deliver a dual (LMILAB)'); |
---|
829 | disp('-> Use sdsettings(''sos.model'',2) to circumvent this, or use another solver (SDPT3, SEDUMI,...)'); |
---|
830 | disp(''); |
---|
831 | disp('-> An alternative reason is that YALMIP detected infeasibility during the compilation phase.'); |
---|
832 | end |
---|
833 | end |
---|
834 | end |
---|
835 | end |
---|
836 | |
---|
837 | m = monoms; |
---|
838 | |
---|
839 | otherwise |
---|
840 | Q = []; |
---|
841 | m = []; |
---|
842 | end |
---|
843 | |
---|
844 | % Don't need these outside |
---|
845 | yalmip('cleardual') |
---|
846 | |
---|
847 | % Done with YALMIP, this is the time it took, minus solver |
---|
848 | if ~isfield(sol,'solvertime') |
---|
849 | sol.solvertime = 0; |
---|
850 | end |
---|
851 | |
---|
852 | sol.yalmiptime = etime(clock,yalmip_time)-sol.solvertime; |
---|
853 | |
---|
854 | |
---|
855 | function p_base_parametric = parameterizedbase(p,z, params,ParametricIndicies,exponent_p,p_base) |
---|
856 | |
---|
857 | % Check for linear parameterization |
---|
858 | parametric_basis = exponent_p(:,ParametricIndicies); |
---|
859 | if all(sum(parametric_basis,2)==0) |
---|
860 | p_base_parametric = full(p_base(:)); |
---|
861 | return |
---|
862 | end |
---|
863 | if all(sum(parametric_basis,2)<=1) |
---|
864 | p_base_parametric = full(p_base(:)); |
---|
865 | n = length(p_base_parametric); |
---|
866 | |
---|
867 | |
---|
868 | if 1 |
---|
869 | [ii,vars] = find(parametric_basis); |
---|
870 | ii = ii(:)'; |
---|
871 | vars = vars(:)'; |
---|
872 | else |
---|
873 | ii = []; |
---|
874 | vars = []; |
---|
875 | js = sum(parametric_basis,1); |
---|
876 | indicies = find(js); |
---|
877 | for i = indicies |
---|
878 | if js(i) |
---|
879 | j = find(parametric_basis(:,i)); |
---|
880 | ii = [ii j(:)']; |
---|
881 | vars = [vars repmat(i,1,js(i))]; |
---|
882 | end |
---|
883 | end |
---|
884 | end |
---|
885 | |
---|
886 | k = setdiff1D(1:n,ii); |
---|
887 | if isempty(k) |
---|
888 | p_base_parametric = p_base_parametric.*sparse(ii,repmat(1,1,n),params(vars)); |
---|
889 | else |
---|
890 | pp = params(vars); % Must do this, bug in ML 6.1 (x=sparse(1);x([1 1]) gives different result in 6.1 and 7.0!) |
---|
891 | p_base_parametric = p_base_parametric.*sparse([ii k(:)'],repmat(1,1,n),[pp(:)' ones(1,1,length(k))]); |
---|
892 | end |
---|
893 | else |
---|
894 | % Bummer, nonlinear parameterization sucks... |
---|
895 | for i = 1:length(p_base) |
---|
896 | j = find(exponent_p(i,ParametricIndicies)); |
---|
897 | if ~isempty(j) |
---|
898 | temp = p_base(i); |
---|
899 | |
---|
900 | for k = 1:length(j) |
---|
901 | if exponent_p(i,ParametricIndicies(j(k)))==1 |
---|
902 | temp = temp*params(j(k)); |
---|
903 | else |
---|
904 | temp = temp*params(j(k))^exponent_p(i,ParametricIndicies(j(k))); |
---|
905 | end |
---|
906 | end |
---|
907 | xx{i} = temp; |
---|
908 | else |
---|
909 | xx{i} = p_base(i); |
---|
910 | end |
---|
911 | end |
---|
912 | p_base_parametric = stackcell(sdpvar(1,1),xx)'; |
---|
913 | end |
---|
914 | |
---|
915 | |
---|
916 | |
---|
917 | function [A,b] = generate_kernel_representation_data(N,N_unique,exponent_m2,exponent_p,p,options,p_base_parametric,ParametricIndicies,MonomIndicies,FirstRun) |
---|
918 | |
---|
919 | persistent saveData |
---|
920 | |
---|
921 | exponent_p_parametric = exponent_p(:,ParametricIndicies); |
---|
922 | exponent_p_monoms = exponent_p(:,MonomIndicies); |
---|
923 | pcoeffs = getbase(p); |
---|
924 | if any(exponent_p_monoms(1,:)) |
---|
925 | pcoeffs=pcoeffs(:,2:end); % No constant term in p |
---|
926 | end |
---|
927 | b = []; |
---|
928 | |
---|
929 | parametric = full((~isempty(ParametricIndicies) & any(any(exponent_p_parametric)))); |
---|
930 | |
---|
931 | % For problems with a lot of similar cones, this saves some time |
---|
932 | reuse = 0; |
---|
933 | if ~isempty(saveData) & isequal(saveData.N,N) & ~FirstRun |
---|
934 | n = saveData.n; |
---|
935 | ind = saveData.ind; |
---|
936 | if isequal(saveData.N_unique,N_unique) & isequal(saveData.exponent_m2,exponent_m2)% & isequal(saveData.epm,exponent_p_monoms) |
---|
937 | reuse = 1; |
---|
938 | end |
---|
939 | else |
---|
940 | % Congruence partition sizes |
---|
941 | for k = 1:size(N,1) |
---|
942 | n(k) = size(N{k},1); |
---|
943 | end |
---|
944 | % Save old SOS definition |
---|
945 | saveData.N = N; |
---|
946 | saveData.n = n; |
---|
947 | saveData.N_unique = N_unique; |
---|
948 | saveData.exponent_m2 = exponent_m2; |
---|
949 | saveData.N_unique = N_unique; |
---|
950 | end |
---|
951 | |
---|
952 | if reuse & options.sos.reuse |
---|
953 | % Get old stuff |
---|
954 | if size(exponent_m2{1},2)==2 % Stupid set(sos(parametric)) case |
---|
955 | ind = spalloc(1,1,0); |
---|
956 | ind(1)=1; |
---|
957 | allj = 1:size(exponent_p_monoms,1); |
---|
958 | used_in_p = ones(size(exponent_p_monoms,1),1); |
---|
959 | else |
---|
960 | allj = []; |
---|
961 | used_in_p = zeros(size(exponent_p_monoms,1),1); |
---|
962 | hash = randn(size(exponent_p_monoms,2),1); |
---|
963 | exponent_p_monoms_hash = exponent_p_monoms*hash; |
---|
964 | for i = 1:size(N_unique,1) |
---|
965 | monom = sparse(N_unique(i,3:end)); |
---|
966 | j = find(exponent_p_monoms_hash == (monom*hash)); |
---|
967 | |
---|
968 | if isempty(j) |
---|
969 | b = [b 0]; |
---|
970 | allj(end+1,1) = 0; |
---|
971 | else |
---|
972 | used_in_p(j) = 1; |
---|
973 | allj(end+1,1:length(j)) = j(:)'; |
---|
974 | end |
---|
975 | end |
---|
976 | ind = saveData.ind; |
---|
977 | end |
---|
978 | else |
---|
979 | allj = []; |
---|
980 | used_in_p = zeros(size(exponent_p_monoms,1),1); |
---|
981 | if size(exponent_m2{1},2)==2 % Stupid set(sos(parametric)) case |
---|
982 | ind = spalloc(1,1,0); |
---|
983 | ind(1)=1; |
---|
984 | allj = 1:size(exponent_p_monoms,1); |
---|
985 | used_in_p = ones(size(exponent_p_monoms,1),1); |
---|
986 | else |
---|
987 | % To speed up some searching, we random-hash data |
---|
988 | hash = randn(size(exponent_p_monoms,2),1); |
---|
989 | for k = 1:length(exponent_m2) |
---|
990 | if isempty(exponent_m2{k}) |
---|
991 | exp_hash{k}=[]; |
---|
992 | else |
---|
993 | exp_hash{k} = sparse((exponent_m2{k}(:,3:end)))*hash; % SPARSE NEEDED DUE TO STRANGE NUMERICS IN MATLAB ON 0s (the stuff will differ on last bit in hex format) |
---|
994 | end |
---|
995 | end |
---|
996 | |
---|
997 | p_hash = exponent_p_monoms*hash; |
---|
998 | ind = spalloc(size(N_unique,1),sum(n.^2),0); |
---|
999 | |
---|
1000 | for i = 1:size(N_unique,1) |
---|
1001 | monom = N_unique(i,3:end); |
---|
1002 | monom_hash = sparse(monom)*hash; |
---|
1003 | LHS = 0; |
---|
1004 | start = 0; |
---|
1005 | for k = 1:size(N,1) |
---|
1006 | j = find(exp_hash{k} == monom_hash); |
---|
1007 | if ~isempty(j) |
---|
1008 | pos=exponent_m2{k}(j,1:2); |
---|
1009 | nss = pos(:,1); |
---|
1010 | mss = pos(:,2); |
---|
1011 | indicies = nss+(mss-1)*n(k); |
---|
1012 | ind(i,indicies+start) = ind(i,indicies+start) + 1; |
---|
1013 | end |
---|
1014 | start = start + (n(k))^2; |
---|
1015 | % start = start + (matrixSOSsize*n(k))^2; |
---|
1016 | end |
---|
1017 | |
---|
1018 | j = find(p_hash == monom_hash); |
---|
1019 | |
---|
1020 | if isempty(j) |
---|
1021 | allj(end+1,1) = 0; |
---|
1022 | else |
---|
1023 | used_in_p(j) = 1; |
---|
1024 | allj(end+1,1:length(j)) = j(:)'; |
---|
1025 | end |
---|
1026 | end |
---|
1027 | end |
---|
1028 | end |
---|
1029 | saveData.ind = ind; |
---|
1030 | |
---|
1031 | % Some parametric terms in p(x,t) do not appear in v'Qv |
---|
1032 | % So these have to be added 0*Q = b |
---|
1033 | not_dealt_with = find(used_in_p==0); |
---|
1034 | while ~isempty(not_dealt_with) |
---|
1035 | j = findrows(exponent_p_monoms,exponent_p_monoms(not_dealt_with(1),:)); |
---|
1036 | allj(end+1,1:length(j)) = j(:)'; |
---|
1037 | used_in_p(j) = 1; |
---|
1038 | not_dealt_with = find(used_in_p==0); |
---|
1039 | ind(end+1,1)=0; |
---|
1040 | end |
---|
1041 | |
---|
1042 | matrixSOSsize = length(p); |
---|
1043 | if parametric |
---|
1044 | % Inconsistent behaviour in MATLAB |
---|
1045 | if size(allj,1)==1 |
---|
1046 | uu = [0;p_base_parametric]; |
---|
1047 | b = sum(uu(allj+1))'; |
---|
1048 | else |
---|
1049 | b = []; |
---|
1050 | for i = 1:matrixSOSsize |
---|
1051 | for j = i:matrixSOSsize |
---|
1052 | if i~=j |
---|
1053 | uu = [0;2*p_base_parametric(:,(i-1)*matrixSOSsize+j)]; |
---|
1054 | else |
---|
1055 | uu = [0;p_base_parametric(:,(i-1)*matrixSOSsize+j)]; |
---|
1056 | end |
---|
1057 | b = [b sum(uu(allj+1),2)']; |
---|
1058 | end |
---|
1059 | end |
---|
1060 | end |
---|
1061 | else |
---|
1062 | if matrixSOSsize == 1 |
---|
1063 | uu = [zeros(size(pcoeffs,1),1) pcoeffs]'; |
---|
1064 | b = sum(uu(allj+1,:),2)'; |
---|
1065 | else |
---|
1066 | b = []; |
---|
1067 | for i = 1:matrixSOSsize |
---|
1068 | for j = i:matrixSOSsize |
---|
1069 | if i~=j |
---|
1070 | uu = [0;2*pcoeffs((i-1)*matrixSOSsize+j,:)']; |
---|
1071 | else |
---|
1072 | uu = [0;pcoeffs((i-1)*matrixSOSsize+j,:)']; |
---|
1073 | end |
---|
1074 | b = [b;sum(uu(allj+1,:),2)']; |
---|
1075 | end |
---|
1076 | end |
---|
1077 | end |
---|
1078 | % uu = [0;pcoeffs(:)]; |
---|
1079 | % b = sum(uu(allj+1),2)'; |
---|
1080 | end |
---|
1081 | |
---|
1082 | b = b'; |
---|
1083 | dualbase = ind; |
---|
1084 | |
---|
1085 | j = 1; |
---|
1086 | A = cell(size(N,1),1); |
---|
1087 | for k = 1:size(N,1) |
---|
1088 | if matrixSOSsize==1 |
---|
1089 | A{k} = dualbase(:,j:j+n(k)^2-1); |
---|
1090 | else |
---|
1091 | % Quick fix for matrix SOS case, should be optimized |
---|
1092 | A{k} = inflate(dualbase(:,j:j+n(k)^2-1),matrixSOSsize,n(k)); |
---|
1093 | end |
---|
1094 | j = j + n(k)^2; |
---|
1095 | end |
---|
1096 | b = b(:); |
---|
1097 | |
---|
1098 | |
---|
1099 | |
---|
1100 | function newAi = inflate(Ai,matrixSOSsize,n); |
---|
1101 | % Quick fix for matrix SOS case, should be optimized |
---|
1102 | newAi = []; |
---|
1103 | for i = 1:matrixSOSsize |
---|
1104 | for r = i:matrixSOSsize |
---|
1105 | for m = 1:size(Ai,1) |
---|
1106 | ai = reshape(Ai(m,:),n,n); |
---|
1107 | V = zeros(matrixSOSsize,matrixSOSsize); |
---|
1108 | V(i,r)=1; |
---|
1109 | V(r,i)=1; |
---|
1110 | ai = kron(V,ai); |
---|
1111 | newAi = [newAi;ai(:)']; |
---|
1112 | end |
---|
1113 | end |
---|
1114 | end |
---|
1115 | |
---|
1116 | |
---|
1117 | function [Z,Q1,R] = sparsenull(A) |
---|
1118 | |
---|
1119 | [Q,R] = qr(A'); |
---|
1120 | n = max(find(sum(abs(R),2))); |
---|
1121 | Q1 = Q(:,1:n); |
---|
1122 | R = R(1:n,:); |
---|
1123 | Z = Q(:,n+1:end); % New basis |
---|
1124 | |
---|
1125 | |
---|
1126 | function [F,obj,BlockedQ,Primal_matrices,Free_variables] = create_kernelmodel(BlockedA,Blockedb,F_parametric,parobj,options,sparsityPattern); |
---|
1127 | |
---|
1128 | % To get the primal kernel representation, we simply use |
---|
1129 | % the built-in dualization module! |
---|
1130 | % First, write the problem in primal kernel format |
---|
1131 | traceobj = 0; |
---|
1132 | dotraceobj = options.sos.traceobj; |
---|
1133 | F = F_parametric; |
---|
1134 | for i = 1:length(Blockedb) |
---|
1135 | |
---|
1136 | |
---|
1137 | sizematrixSOS = sqrt(size(Blockedb{i},2)); |
---|
1138 | for k = 1:sizematrixSOS |
---|
1139 | for r = k:sizematrixSOS |
---|
1140 | res{(k-1)*sizematrixSOS+r} = 0; |
---|
1141 | end |
---|
1142 | end |
---|
1143 | |
---|
1144 | for j = 1:length(BlockedA{i}) |
---|
1145 | n = sqrt(size(BlockedA{i}{j},2)); |
---|
1146 | BlockedQ{i}{j} = sdpvar(n*sizematrixSOS,n*sizematrixSOS); |
---|
1147 | F = F + set(BlockedQ{i}{j}); |
---|
1148 | if sizematrixSOS>0 |
---|
1149 | % Matrix valued sum of sqaures |
---|
1150 | % Loop over all elements |
---|
1151 | starttop = 1; |
---|
1152 | for k = 1:sizematrixSOS |
---|
1153 | startleft = 1; |
---|
1154 | for r = 1:sizematrixSOS |
---|
1155 | if k<=r |
---|
1156 | Qkr = BlockedQ{i}{j}(starttop:starttop+n-1,startleft:startleft+n-1); |
---|
1157 | res{(k-1)*sizematrixSOS+r} = res{(k-1)*sizematrixSOS+r} + BlockedA{i}{j}*reshape(Qkr,n^2,1); |
---|
1158 | end |
---|
1159 | startleft = startleft + n; |
---|
1160 | end |
---|
1161 | starttop = starttop + n; |
---|
1162 | end |
---|
1163 | else |
---|
1164 | % Standard case |
---|
1165 | res{1} = res{1} + BlockedA{i}{j}*reshape(BlockedQ{i}{j},n^2,1); |
---|
1166 | end |
---|
1167 | if dotraceobj |
---|
1168 | traceobj = traceobj + trace(BlockedQ{i}{j}); |
---|
1169 | end |
---|
1170 | end |
---|
1171 | for k = 1:sizematrixSOS |
---|
1172 | for r = k:sizematrixSOS |
---|
1173 | F = F + set(res{(k-1)*sizematrixSOS+r} == Blockedb{i}(:,(k-1)*sizematrixSOS+r)); |
---|
1174 | end |
---|
1175 | end |
---|
1176 | end |
---|
1177 | |
---|
1178 | % % Constrain elements according to a desired sparsity |
---|
1179 | if ~isempty(sparsityPattern) |
---|
1180 | res = []; |
---|
1181 | for i = 1:length(BlockedQ) |
---|
1182 | for j = 1:length(BlockedQ{i}) |
---|
1183 | if ~isempty(sparsityPattern{i}{j}) |
---|
1184 | H = spalloc(length(BlockedQ{i}{j}),length(BlockedQ{i}{j}),length(sparsityPattern{i}{j})); |
---|
1185 | H(sparsityPattern{i}{j}) = 1; |
---|
1186 | k = find(triu(H)); |
---|
1187 | res = [res;BlockedQ{i}{j}(k)]; |
---|
1188 | end |
---|
1189 | end |
---|
1190 | end |
---|
1191 | F = F + set(0 == res); |
---|
1192 | end |
---|
1193 | |
---|
1194 | % And get the primal model of this |
---|
1195 | if isempty(parobj) |
---|
1196 | if options.sos.traceobj |
---|
1197 | [F,obj,Primal_matrices,Free_variables] = dualize(F,traceobj,1,options.sos.extlp); |
---|
1198 | else |
---|
1199 | [F,obj,Primal_matrices,Free_variables] = dualize(F,[],1,options.sos.extlp); |
---|
1200 | end |
---|
1201 | else |
---|
1202 | [F,obj,Primal_matrices,Free_variables] = dualize(F,parobj,1,options.sos.extlp); |
---|
1203 | end |
---|
1204 | % In dual mode, we maximize |
---|
1205 | obj = -obj; |
---|
1206 | |
---|
1207 | |
---|
1208 | function [F,obj,BlockedQ,sol] = create_imagemodel(BlockedA,Blockedb,F_parametric,parobj,options); |
---|
1209 | |
---|
1210 | |
---|
1211 | % ********************************* |
---|
1212 | % IMAGE REPRESENTATION |
---|
1213 | % Needed for nonlinearly parameterized problems |
---|
1214 | % More efficient in some cases |
---|
1215 | % ********************************* |
---|
1216 | g = []; |
---|
1217 | vecQ = []; |
---|
1218 | sol.problem = 0; |
---|
1219 | for i = 1:length(BlockedA) |
---|
1220 | q = []; |
---|
1221 | A = []; |
---|
1222 | for j = 1:length(BlockedA{i}) |
---|
1223 | n = sqrt(size(BlockedA{i}{j},2)); |
---|
1224 | BlockedQ{i}{j} = sdpvar(n,n); |
---|
1225 | q = [q;reshape(BlockedQ{i}{j},n^2,1)]; |
---|
1226 | A = [A BlockedA{i}{j}]; |
---|
1227 | end |
---|
1228 | vecQ = [vecQ;q]; |
---|
1229 | g = [g;Blockedb{i}-A*q]; |
---|
1230 | end |
---|
1231 | |
---|
1232 | g_vars = getvariables(g); |
---|
1233 | q_vars = getvariables(vecQ); |
---|
1234 | x_vars = setdiff(g_vars,q_vars); |
---|
1235 | |
---|
1236 | base = getbase(g); |
---|
1237 | if isempty(x_vars) |
---|
1238 | A = base(:,1);base = base(:,2:end); |
---|
1239 | B = (base(:,ismember(g_vars,q_vars))); |
---|
1240 | Bnull = sparsenull(B); |
---|
1241 | t = sdpvar(size(Bnull,2),1); |
---|
1242 | imQ = -B\A+Bnull*t; |
---|
1243 | else |
---|
1244 | A = base(:,1);base = base(:,2:end); |
---|
1245 | C = base(:,ismember(g_vars,x_vars)); |
---|
1246 | B = (base(:,ismember(g_vars,q_vars))); |
---|
1247 | [Bnull,Q1,R1] = sparsenull(B);Bnull(abs(Bnull) < 1e-12) = 0; |
---|
1248 | t = sdpvar(size(Bnull,2),1); |
---|
1249 | imQ = -Q1*(R1'\(A+C*recover(x_vars)))+Bnull*t; |
---|
1250 | end |
---|
1251 | notUsed = find(sum(abs(B),2)==0); |
---|
1252 | if ~isempty(notUsed) |
---|
1253 | ff=g(notUsed); |
---|
1254 | if isa(ff,'double') |
---|
1255 | if nnz(ff)>0 |
---|
1256 | sol.yalmiptime = 0; % FIX |
---|
1257 | sol.solvertime = 0; |
---|
1258 | sol.problem = 2; |
---|
1259 | sol.info = yalmiperror(1,'YALMIP'); |
---|
1260 | F = []; |
---|
1261 | obj = []; |
---|
1262 | if options.verbose > 0 |
---|
1263 | disp(' '); |
---|
1264 | disp('-> During construction of data, I encountered a situation') |
---|
1265 | disp('-> situation that tells me that the problem is trivially infeasible!') |
---|
1266 | disp('-> Have you forgotten to define some parametric variables?,') |
---|
1267 | disp('-> or perhaps you have a parametric problem where the highest') |
---|
1268 | disp('-> power in some of the independent variables is odd for any choice') |
---|
1269 | disp('-> of parametric variables, such as x^8+x^7+t*y^3') |
---|
1270 | disp('-> Anyway, take a good look at your model again...'); |
---|
1271 | end |
---|
1272 | return |
---|
1273 | % error('You seem to have a strange model. Have you forgotten to define some parametric variable?'); |
---|
1274 | end |
---|
1275 | else |
---|
1276 | F_parametric = F_parametric + set(g(notUsed)==0); |
---|
1277 | end |
---|
1278 | end |
---|
1279 | F_sos = set([]); |
---|
1280 | obj = 0; |
---|
1281 | for i = 1:length(BlockedQ) |
---|
1282 | for j = 1:size(BlockedQ{i},2) |
---|
1283 | Q_old = BlockedQ{i}{j}; |
---|
1284 | Q_old_vars = getvariables(Q_old); |
---|
1285 | Q_old_base = getbase(Q_old); |
---|
1286 | in_this = []; |
---|
1287 | for k = 1:length(Q_old_vars) |
---|
1288 | in_this = [in_this;find(Q_old_vars(k)==q_vars)]; |
---|
1289 | end |
---|
1290 | Q_new = Q_old_base(:,1) + Q_old_base(:,2:end)*imQ(in_this); |
---|
1291 | Q_new = reshape(Q_new,length(BlockedQ{i}{j}),length(BlockedQ{i}{j})); |
---|
1292 | obj = obj+trace(Q_new); |
---|
1293 | if ~isa(Q_new,'double') |
---|
1294 | F_sos = F_sos + set(Q_new); |
---|
1295 | elseif min(eig(Q_new))<-1e-8 |
---|
1296 | sol.yalmiptime = 0; % FIX |
---|
1297 | sol.solvertime = 0; |
---|
1298 | sol.problem = 2; |
---|
1299 | sol.info = yalmiperror(1,'YALMIP'); |
---|
1300 | F = []; |
---|
1301 | obj = []; |
---|
1302 | error('Problem is trivially infeasible. After block-diagonalizing, I found constant negative definite blocks!'); |
---|
1303 | return |
---|
1304 | end |
---|
1305 | BlockedQ{i}{j} = Q_new; |
---|
1306 | end |
---|
1307 | end |
---|
1308 | |
---|
1309 | F = F_parametric + F_sos; |
---|
1310 | |
---|
1311 | if isa(obj,'double') | (options.sos.traceobj == 0) |
---|
1312 | obj = []; |
---|
1313 | end |
---|
1314 | |
---|
1315 | if ~isempty(parobj) |
---|
1316 | obj = parobj; |
---|
1317 | end |
---|
1318 | |
---|
1319 | |
---|
1320 | function [F,obj,options] = create_lrmodel(BlockedA,Blockedb,F_parametric,parobj,options,ParametricVariables) |
---|
1321 | % Some special code for ther low-rank model in SDPLR |
---|
1322 | % Experimental code, not official yet |
---|
1323 | allb = []; |
---|
1324 | allA = []; |
---|
1325 | K.s = []; |
---|
1326 | for i = 1:length(Blockedb) |
---|
1327 | allb = [allb;Blockedb{i}]; |
---|
1328 | Ai = []; |
---|
1329 | for j = 1:size(BlockedA{i},2) |
---|
1330 | Ai = [Ai BlockedA{i}{j}]; |
---|
1331 | K.s = [K.s sqrt(size(BlockedA{i}{j},2))]; |
---|
1332 | end |
---|
1333 | %blkdiag bug in 7.0... |
---|
1334 | [n1,m1] = size(allA); |
---|
1335 | [n2,m2] = size(Ai); |
---|
1336 | allA = [allA spalloc(n1,m2,0);spalloc(n2,m1,0) Ai]; |
---|
1337 | end |
---|
1338 | options.solver = 'sdplr'; |
---|
1339 | z = recover(ParametricVariables) |
---|
1340 | start = size(BlockedA,2)+1; |
---|
1341 | Mi = []; |
---|
1342 | for i = 1:length(allb) |
---|
1343 | if isa(allb(i),'sdpvar') |
---|
1344 | [Qi,ci,fi,xi,infoi] = quaddecomp(allb(i),z); |
---|
1345 | else |
---|
1346 | Qi = zeros(length(z)); |
---|
1347 | ci = zeros(length(z),1); |
---|
1348 | fi = allb(i); |
---|
1349 | end |
---|
1350 | Z = -[fi ci'/2;ci/2 Qi]; |
---|
1351 | Mi = [Mi;Z(:)']; |
---|
1352 | end |
---|
1353 | K.s = [K.s length(z)+1]; |
---|
1354 | zeroRow = zeros(1,size(allA,2)); |
---|
1355 | allA = [allA Mi;zeroRow 1 zeros(1,K.s(end)^2-1)]; |
---|
1356 | b = zeros(size(allA,1),1);b(end) = 1; |
---|
1357 | y = sdpvar(length(b),1); |
---|
1358 | CminusAy = -allA'*y; |
---|
1359 | start = 1; |
---|
1360 | |
---|
1361 | % Get the cost, expressed in Z |
---|
1362 | [Qi,ci,fi,xi,infoi] = quaddecomp(parobj,z); |
---|
1363 | C = [fi ci'/2;ci/2 Qi]; |
---|
1364 | F = set([]); |
---|
1365 | for i = 1:length(K.s) |
---|
1366 | if i<length(K.s) |
---|
1367 | F = F + set(reshape(CminusAy(start:start+K.s(i)^2-1),K.s(i),K.s(i))); |
---|
1368 | else |
---|
1369 | F = F + set(reshape(C(:) + CminusAy(start:start+K.s(i)^2-1),K.s(i),K.s(i))); |
---|
1370 | end |
---|
1371 | start = start + K.s(i)^2; |
---|
1372 | end |
---|
1373 | obj = -b'*y; |
---|
1374 | |
---|
1375 | options.sdplr.forcerank = ceil(K.s/2); |
---|
1376 | options.sdplr.forcerank(end) = 1; |
---|
1377 | options.sdplr.feastol = 1e-7; |
---|
1378 | |
---|
1379 | |
---|
1380 | |
---|
1381 | function [exponent_m2,N_unique] = expandmatrix(exponent_m2,N_unique,n); |
---|
1382 | |
---|
1383 | |
---|
1384 | function [sol,m,Q,residuals,everything] = solvesos_find_blocks(F,obj,options,params,candidateMonomials) |
---|
1385 | |
---|
1386 | tol = options.sos.numblkdg; |
---|
1387 | if tol > 1e-2 |
---|
1388 | disp(' '); |
---|
1389 | disp('-> Are you sure you meant to have a tolerance in numblk that big!') |
---|
1390 | disp('-> The options numblkdiag controls the tolerance, it is not a 0/1 switch.') |
---|
1391 | disp(' '); |
---|
1392 | end |
---|
1393 | options.sos.numblkdg = 0; |
---|
1394 | [sol,m,Q,residuals,everything] = solvesos(F,obj,options,params,candidateMonomials); |
---|
1395 | |
---|
1396 | % Save old structure to find out when we have stalled |
---|
1397 | for i = 1:length(Q) |
---|
1398 | oldlengths{i} = length(Q{i}); |
---|
1399 | end |
---|
1400 | |
---|
1401 | go_on = (sol.problem == 0 | sol.problem == 4); |
---|
1402 | while go_on |
---|
1403 | |
---|
1404 | for sosfun = 1:length(Q) |
---|
1405 | Qtemp = Q{sosfun}; |
---|
1406 | keep = diag(Qtemp)>tol; |
---|
1407 | Qtemp(:,find(~keep)) = []; |
---|
1408 | Qtemp(find(~keep),:) = []; |
---|
1409 | |
---|
1410 | m{sosfun} = m{sosfun}(find(keep)); |
---|
1411 | |
---|
1412 | Qtemp(abs(Qtemp) < tol) = 0; |
---|
1413 | [v1,dummy1,r1,dummy3]=dmperm(Qtemp+eye(length(Qtemp))); |
---|
1414 | lengths{sosfun} = []; |
---|
1415 | n{sosfun} = {}; |
---|
1416 | for blocks = 1:length(r1)-1 |
---|
1417 | i1 = r1(blocks); |
---|
1418 | i2 = r1(blocks+1)-1; |
---|
1419 | if i2>i1 |
---|
1420 | n{sosfun}{blocks} = m{sosfun}(v1(i1:i2)); |
---|
1421 | else |
---|
1422 | n{sosfun}{blocks} = m{sosfun}(v1(i1)); |
---|
1423 | end |
---|
1424 | lengths{sosfun} = [lengths{sosfun} length(n{sosfun}{blocks})]; |
---|
1425 | end |
---|
1426 | lengths{sosfun} = sort(lengths{sosfun}); |
---|
1427 | end |
---|
1428 | go_on = ~isequal(lengths,oldlengths); |
---|
1429 | oldlengths = lengths; |
---|
1430 | if go_on |
---|
1431 | [sol,m,Q,residuals,everything] = solvesos(F,obj,options,params,n); |
---|
1432 | go_on = go_on & (sol.problem == 0 | sol.problem == 4); |
---|
1433 | if sol.problem == 1 |
---|
1434 | disp('-> Feasibility was lost during the numerical block-diagonalization.') |
---|
1435 | disp('-> The setting sos.numblkdiag is probably too big') |
---|
1436 | end |
---|
1437 | end |
---|
1438 | end |
---|
1439 | |
---|
1440 | |
---|
1441 | function sol = solveranksos(F,obj,options,ranks,BlockedQ) |
---|
1442 | |
---|
1443 | Frank = set([]); |
---|
1444 | for i = 1:length(ranks) |
---|
1445 | if ~isinf(ranks(i)) |
---|
1446 | Frank = Frank + set(rank(BlockedQ{i}{1}) <= ranks(i)); |
---|
1447 | end |
---|
1448 | end |
---|
1449 | % rank adds the pos.def constraints again!!, so we remove them |
---|
1450 | check = ones(length(F),1); |
---|
1451 | keep = ones(length(F),1); |
---|
1452 | for i = 1:length(BlockedQ) |
---|
1453 | for j = 1:length(BlockedQ{i}) |
---|
1454 | Qij = BlockedQ{i}{j}; |
---|
1455 | for k = find(check)' |
---|
1456 | if isequal(Qij,sdpvar(F(k))) |
---|
1457 | keep(k) = 0; |
---|
1458 | check(k) = 0; |
---|
1459 | end |
---|
1460 | end |
---|
1461 | end |
---|
1462 | end |
---|
1463 | % Let's hope LMIRANK is there |
---|
1464 | sol = solvesdp(F(find(keep)) + Frank,[],sdpsettings(options,'solver','')); |
---|
1465 | |
---|
1466 | |
---|