1 | function [sol,x_extract,momentsstructure,sosout] = solvemoment(F,obj,options,k) |
---|
2 | %SOLVEMOMENT Application of Lasserre's moment-method for polynomial programming |
---|
3 | % |
---|
4 | % min h(x) |
---|
5 | % subject to |
---|
6 | % F(x) > 0, |
---|
7 | % |
---|
8 | % [DIAGNOSTIC,X,MOMENT,SOS] = SOLVEMOMENT(F,h,options,k) |
---|
9 | % |
---|
10 | % diagnostic : Struct with diagnostics |
---|
11 | % x : Extracted global solutions |
---|
12 | % moment : Structure with various variables needed to recover solution |
---|
13 | % sos : SOS decomposition {max t s.t h-t = p0+sum(pi*Fi), pi = vi'*Qi*vi} |
---|
14 | % |
---|
15 | % Input |
---|
16 | % F : SET object with polynomial inequalities and equalities. |
---|
17 | % h : SDPVAR object describing the polynomial h(x). |
---|
18 | % options : solver options from SDPSETTINGS. |
---|
19 | % k : Level of relaxation. If empty or not given, smallest possible applied. |
---|
20 | % |
---|
21 | % The behaviour of the moment relaxation can be controlled |
---|
22 | % using the fields 'moment' in SDPSETTINGS |
---|
23 | % |
---|
24 | % moment.refine : Perform #refine Newton iterations in extracation of global solutions. |
---|
25 | % This can improve numerical accuracy of extracted solutions in some cases. |
---|
26 | % moment.extractrank : Try (forcefully) to extract #extractrank global solutions. |
---|
27 | % This feature should normally not be used and is default 0. |
---|
28 | % moment.rceftol : Tolerance during Gaussian elimination used in extraction of global solutions. |
---|
29 | % Default is -1 which means heuristic choice by YALMIP. |
---|
30 | % |
---|
31 | % Some of the fields are only used when the moment relaxation is called |
---|
32 | % indirectly from SOLVESDP. |
---|
33 | % |
---|
34 | % moment.solver : SDP solver used in moment relxation. Default '' |
---|
35 | % moment.order : Order of relxation. Default [] meaning lowest possible. |
---|
36 | % |
---|
37 | % See also SDPVAR, SET, SDPSETTINGS, SOLVESDP |
---|
38 | |
---|
39 | % Author Johan Löfberg |
---|
40 | % $Id: solvemoment.m,v 1.4 2006/10/17 14:02:02 joloef Exp $ |
---|
41 | |
---|
42 | if nargin ==0 |
---|
43 | help solvemoment |
---|
44 | return |
---|
45 | end |
---|
46 | |
---|
47 | if nargin<2 |
---|
48 | obj=[]; |
---|
49 | end |
---|
50 | |
---|
51 | if (nargin>=3) & (isa(options,'double') & ~isempty(options)) |
---|
52 | help solvemoment |
---|
53 | error('Order of arguments have changed in solvemoment. Update code'); |
---|
54 | end |
---|
55 | |
---|
56 | if nargin<3 | (isempty(options)) |
---|
57 | options = sdpsettings; |
---|
58 | end |
---|
59 | |
---|
60 | % Relaxation-order given? |
---|
61 | if nargin<4 |
---|
62 | k = []; |
---|
63 | end |
---|
64 | |
---|
65 | % Check for wrong syntax |
---|
66 | if ~isempty(F) & ~isa(F,'lmi') & ~isa(F,'constraint') |
---|
67 | error('First argument should be a SET object') |
---|
68 | end |
---|
69 | |
---|
70 | if isa(F,'constraint') |
---|
71 | F = set(F); |
---|
72 | end |
---|
73 | |
---|
74 | % Take care of rational expressions |
---|
75 | [F,failure] = expandmodel(F,obj); |
---|
76 | if failure |
---|
77 | error('Could not expand model (rational functions, min/max etc). Avoid nonlinear operators in moment problems.'); |
---|
78 | end |
---|
79 | |
---|
80 | % Get all element-wise constraints, and put them in a vector |
---|
81 | % Furthermore, gather the other constraints and place them |
---|
82 | % in a new LMI object. |
---|
83 | % Additionally, we find out the variables on which we perform |
---|
84 | % the relaxation. |
---|
85 | vecConstraints = []; |
---|
86 | sdpConstraints = []; |
---|
87 | isinequality = []; |
---|
88 | binaries = []; |
---|
89 | xvars = []; |
---|
90 | Fnew = set([]); |
---|
91 | for i = 1:length(F) |
---|
92 | if is(F(i),'elementwise') |
---|
93 | X = sdpvar(F(i)); |
---|
94 | vecConstraints = [vecConstraints;X(:)]; |
---|
95 | isinequality = [isinequality ones(1,prod(size(X)))]; |
---|
96 | xvars = [xvars depends(X(:))]; |
---|
97 | elseif is(F(i),'equality') |
---|
98 | X = sdpvar(F(i)); |
---|
99 | vecConstraints = [vecConstraints;X(:)]; |
---|
100 | isinequality = [isinequality zeros(1,prod(size(X)))]; |
---|
101 | xvars = [xvars depends(X(:))]; |
---|
102 | elseif is(F(i),'sdp') |
---|
103 | sdpConstraints{end+1} = sdpvar(F(i)); |
---|
104 | xvars = [xvars depends(F(i))]; |
---|
105 | elseif is(F(i),'binary') |
---|
106 | binaries = [binaries getvariables(F(i))]; |
---|
107 | else |
---|
108 | Fnew = Fnew+F(i); % Should only be SOCP constraints |
---|
109 | end |
---|
110 | end |
---|
111 | |
---|
112 | % Recover the involved variables |
---|
113 | x = recover(unique([depends(obj) xvars])); |
---|
114 | n = length(x); |
---|
115 | |
---|
116 | % Check degrees of constraints |
---|
117 | deg = []; |
---|
118 | for i = 1:length(vecConstraints) |
---|
119 | deg(end+1) = degree(vecConstraints(i)); |
---|
120 | end |
---|
121 | for i = 1:length(sdpConstraints) |
---|
122 | deg(end+1) = degree(sdpConstraints{i}); |
---|
123 | end |
---|
124 | if isempty(deg) |
---|
125 | deg = 0; |
---|
126 | end |
---|
127 | |
---|
128 | % Perform Schur complements if possible to reduce the |
---|
129 | % dimension of the SDP constraints |
---|
130 | % for i = 1:length(sdpConstraints) |
---|
131 | % Fi = sdpConstraints{i}; |
---|
132 | % j = 1; |
---|
133 | % while j<=length(Fi) & (length(Fi)>1) |
---|
134 | % if isa(Fi(j,j),'double') |
---|
135 | % ks = 1:length(Fi);ks(j)=[]; |
---|
136 | % v = Fi(ks,j); |
---|
137 | % vv = v*v'/Fi(j,j); |
---|
138 | % if degree(vv)<=max(deg) |
---|
139 | % Fi = Fi(ks,ks) - vv; |
---|
140 | % end |
---|
141 | % else |
---|
142 | % j = j+1; |
---|
143 | % end |
---|
144 | % end |
---|
145 | % end |
---|
146 | |
---|
147 | % Create lowest possible relaxation if k=[] |
---|
148 | % k_min = floor((max(degree(obj)+1,max(deg))+1)/2); % why did I use this? |
---|
149 | d = ceil((max(degree(obj),max(deg)))/2); |
---|
150 | k_min = d; |
---|
151 | if isempty(k) |
---|
152 | k = k_min; |
---|
153 | else |
---|
154 | if k<k_min |
---|
155 | error('Higher order relaxation needed') |
---|
156 | end |
---|
157 | end |
---|
158 | |
---|
159 | % Generate monomials of order k |
---|
160 | u{k} = monolist(x,k); |
---|
161 | |
---|
162 | % Largest moment matrix. NOTE SHIFT M{k+1} = M_k. |
---|
163 | M{k+1}=u{k}*u{k}'; |
---|
164 | % Moment matrices easily generated with this trick |
---|
165 | % The matrices will NOT be rank-1 since the products |
---|
166 | % generate the relaxed variables |
---|
167 | |
---|
168 | % ... and lower degree localization matrices |
---|
169 | M{1} = 1; |
---|
170 | for i = 1:1:k-1; |
---|
171 | n_i = round(factorial(n+k-i)/(factorial(n)*factorial(k-i))); |
---|
172 | M{k-i+1} = M{k+1}(1:n_i,1:n_i); |
---|
173 | end |
---|
174 | |
---|
175 | % %Moment structure exploition |
---|
176 | % M_original = M{end}; |
---|
177 | % setsdpvar(recover(getvariables(M_original)),0*getvariables(M_original)'); |
---|
178 | % if 1%options.moment.blockdiag & isempty(F) |
---|
179 | % Ms = blockdiagmoment(obj,M); |
---|
180 | % end |
---|
181 | |
---|
182 | % Lasserres relaxation (Lasserre, SIAM J. OPTIM, 11(3) 796-817) |
---|
183 | Fmoments = set(M{k+1}>0); |
---|
184 | for i = 1:length(vecConstraints) |
---|
185 | v_k = floor((degree(vecConstraints(i))+1)/2); |
---|
186 | Localizer = vecConstraints(i)*M{k-v_k+1}; |
---|
187 | if isinequality(i) |
---|
188 | Fmoments = Fmoments+set(Localizer>0); |
---|
189 | else |
---|
190 | indicies = find(triu(ones(length(Localizer)))); |
---|
191 | Fmoments = Fmoments+set(Localizer(indicies)==0); |
---|
192 | end |
---|
193 | end |
---|
194 | for i = 1:length(sdpConstraints) |
---|
195 | v_k = floor((degree(sdpConstraints{i})+1)/2); |
---|
196 | Fmoments = Fmoments+set(kron(M{k-v_k+1},sdpConstraints{i})>0); |
---|
197 | end |
---|
198 | |
---|
199 | % Add them all |
---|
200 | Fnew = Fnew + Fmoments;%unblkdiag(Fmoments); |
---|
201 | |
---|
202 | % No objective, minimize trace on moment-matrix instead |
---|
203 | if isempty(obj) |
---|
204 | obj = trace(M{k+1}); |
---|
205 | end |
---|
206 | |
---|
207 | % Get all binary and reduce problem |
---|
208 | binaries = union(binaries,yalmip('binvariables')); |
---|
209 | if ~isempty(binaries) |
---|
210 | obj = eliminateBinary(obj,binaries); |
---|
211 | for i = 1:length(Fmoments) |
---|
212 | Fnew(i) = eliminateBinary(Fnew(i),binaries); |
---|
213 | end |
---|
214 | for i = 2:1:k+1; |
---|
215 | M{i} = eliminateBinary(M{i},binaries); |
---|
216 | end |
---|
217 | end |
---|
218 | |
---|
219 | % Solve |
---|
220 | sol = solvesdp(Fnew,obj,sdpsettings(options,'relax',1)); |
---|
221 | |
---|
222 | % Construct SOS decompositions if the user wants these |
---|
223 | if nargout >= 4 |
---|
224 | sosout.t = relaxdouble(obj); |
---|
225 | sosout.Q0 = dual(Fnew(1)); |
---|
226 | sosout.v0 = u{end}; |
---|
227 | sosout.p0 = u{end}'*dual(Fnew(1))*u{end}; |
---|
228 | for i = 1:length(vecConstraints) |
---|
229 | if isinequality(i) |
---|
230 | sosout.Qi{i} = dual(Fnew(i+1)); |
---|
231 | sosout.vi{i} = u{end}(1:length(sosout.Qi{i})); |
---|
232 | sosout.pi{i} = sosout.vi{i}'*sosout.Qi{i}*sosout.vi{i}; |
---|
233 | else |
---|
234 | % For equalities, we need to reconstruct a matrix from the |
---|
235 | % upper triangle |
---|
236 | vecDuals = dual(Fnew(i+1)); |
---|
237 | n = -(1/2)+sqrt(1/4+length(vecDuals)*2); |
---|
238 | [ix,jx] = find(triu(ones(n))); |
---|
239 | temp = sparse(ix,jx,vecDuals); |
---|
240 | temp = temp + temp';%- diag(diag(temp)); |
---|
241 | sosout.Qi{i} = temp/2; |
---|
242 | sosout.vi{i} = u{end}(1:length(sosout.Qi{i})); |
---|
243 | sosout.pi{i} = sosout.vi{i}'*sosout.Qi{i}*sosout.vi{i}; |
---|
244 | end |
---|
245 | end |
---|
246 | end |
---|
247 | |
---|
248 | % Get the moment matrices |
---|
249 | % M{end} = M_original; |
---|
250 | for i = 1:k+1 |
---|
251 | moments{i} = relaxdouble(M{i}); |
---|
252 | end |
---|
253 | |
---|
254 | % Extract solutions if possible (at-least fesible and unbounded) |
---|
255 | momentsstructure.moment = moments; |
---|
256 | momentsstructure.x = x; |
---|
257 | momentsstructure.monomials = u{k}; |
---|
258 | momentsstructure.n = n; |
---|
259 | momentsstructure.d = max(1,ceil(max(deg)/2)); |
---|
260 | x_extract = {}; |
---|
261 | if nargout>=2 & ~(sol.problem == 1 | sol.problem == 2) |
---|
262 | momentsstructure.moment = moments; |
---|
263 | momentsstructure.x = x; |
---|
264 | momentsstructure.monomials = u{k}; |
---|
265 | momentsstructure.n = n; |
---|
266 | momentsstructure.d = max(1,ceil(max(deg)/2)); |
---|
267 | x_extract = extractsolution(momentsstructure,options); |
---|
268 | end |
---|