[37] | 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 |
---|