solvemoment.m

国外专家做的求解LMI鲁棒控制的工具箱,可以相对高效的解决LMI问题

function [sol,x_extract,momentsstructure,sosout] = solvemoment(F,obj,options,k)
%SOLVEMOMENT Application of Lasserre's moment-method for polynomial programming
%  
%   min        h(x)
%   subject to
%              F(x) > 0, 
%  
%   [DIAGNOSTIC,X,MOMENT] = SOLVEMOMENT(F,h,options,k)
%  
%   diagnostic : Struct with diagnostics
%   x          : Extracted global solutions
%   moment     : Structure with various variables needed to recover solution
%   sos        : SOS decomposition {max a s.t h-a = p0+sum(pi*Fi), pi = vi'*Qi*vi}
%   
%   Input
%      F       : SET object with polynomial inequalities and equalities.
%      h       : SDPVAR object describing the polynomial h(x). 
%      options : solver options from SDPSETTINGS.
%      k       : Level of relaxation. If empty or not given, smallest possible applied.
%
%   The behaviour of the moment relaxation can be controlled
%   using the fields 'moment' in SDPSETTINGS
%
%      moment.refine      : Perform #refine Newton iterations in extracation of global solutions.
%                           This can improve numerical accuracy of extracted solutions in some cases.
%      moment.extractrank : Try (forcefully) to extract #extractrank global solutions.
%                           This feature should normally not be used and is then set to 0.
%      moment.rceftol     : Tolerance during Gaussian elimination used in extraction of global solutions.
%                           Default is -1 which means heuristic choice by YALMIP.
%
%   Some of the fields are only used when the moment relaxation is called
%   indirectly from SOLVESDP.
%
%      moment.solver : SDP solver used in moment relxation. Default ''
%      moment.order  : Order of relxation. Default [] meaning lowest possible.
%
%   See also SDPVAR, SET, SDPSETTINGS, SOLVESDP

% Author Johan L鰂berg
% $Id: solvemoment.m,v 1.28 2005/07/19 13:57:40 joloef Exp $

if nargin ==0
    help solvemoment
    return
end

if nargin<2
    obj=[];
end

if (nargin>=3) & (isa(options,'double') & ~isempty(options))
    help solvemoment
    error('Order of arguments have changed in solvemoment. Update code');   
end

if nargin<3 | (isempty(options))
    options = sdpsettings;
end

% Relaxation-order given?
if nargin<4
    k = [];
end

% Check for wrong syntax
if ~isempty(F) & ~isa(F,'lmi')
    error('First argument should be a SET object')
end

% Get all element-wise constraints, and put them in a vector
% Furthermore, gather the other constraints and place them 
% in a new LMI object.
% Additionally, we find out the variables on which we perform
% the relaxation.
vecConstraints = [];
sdpConstraints = [];
isinequality = [];
xvars = [];
Fnew = set([]);
for i = 1:length(F)
    if is(F(i),'elementwise') 
        X = sdpvar(F(i));
        vecConstraints = [vecConstraints;X(:)];
        isinequality = [isinequality ones(1,length(X))];
        xvars = [xvars depends(X(:))];
    elseif is(F(i),'equality') 
        X = sdpvar(F(i));
        vecConstraints = [vecConstraints;X(:)];
        isinequality = [isinequality zeros(1,length(X))];            
        xvars = [xvars depends(X(:))];
    elseif is(F(i),'sdp')        
        sdpConstraints{end+1} = sdpvar(F(i));
        xvars = [xvars depends(F(i))];
    else
        Fnew = Fnew+F(i); % Should only be SOCP constraints
    end
end

% Recover the involved variables
x = recover(unique([depends(obj) xvars]));
n = length(x);

% Check degrees of constraints
deg = [];
for i = 1:length(vecConstraints)
   deg(end+1) = degree(vecConstraints(i));
end
for i = 1:length(sdpConstraints)
   deg(end+1) = degree(sdpConstraints{i});
end
if isempty(deg)
    deg = 0;
end

% Perform Schur complements if possible to reduce the
% dimension of the SDP constraints
% for i = 1:length(sdpConstraints)
%     Fi = sdpConstraints{i};
%     j = 1;
%     while j<=length(Fi) & (length(Fi)>1)
%         if isa(Fi(j,j),'double')
%             ks = 1:length(Fi);ks(j)=[];
%             v = Fi(ks,j);
%             vv = v*v'/Fi(j,j);
%             if degree(vv)<=max(deg)
%                 Fi = Fi(ks,ks) - vv;
%             end
%         else
%             j = j+1;
%         end
%     end
% end

% Create lowest possible relaxation if k=[]
% k_min = floor((max(degree(obj)+1,max(deg))+1)/2); % why did I use this?
d = ceil((max(degree(obj),max(deg)))/2);
k_min = d;
if isempty(k)
    k = k_min;
else
    if k<k_min
        error('Higher order relaxation needed')
    end
end

% Generate monomials of order k
u{k} = monolist(x,k);

% Largest moment matrix. NOTE SHIFT M{k+1} = M_k.
M{k+1}=u{k}*u{k}';
% Moment matrices easily generated with this trick
% The matrices will NOT be rank-1 since the products
% generate the relaxed variables

% rt = [];
% MMM = M{k+1};
% setsdpvar(recover(getvariables(MMM)),0+0*getvariables(MMM)')
% Fb = set([]);
% if 1
%     involved = [];       
%     any_remove = 0;
%     removed = [];
%  
%     vars = getvariables(M{k+1});
%     for i = 1:length(vars)
%         B = getbasematrix(M{k+1},vars(i));        
%        v = eig(reshape(B,length(M{k+1}),length(M{k+1})));
%        if all(v>=0)
%            v
%        end
%     end
%     
%     while any_remove == 1
%         i = 1;
%         any_remove = 0;
%         bajs = [];
%         while i<=length(M{k+1})
%             Mij = M{k+1}(i,i);
%             if  isa(Mij,'sdpvar')
%                 B = getbasematrix(M{k+1},getvariables(Mij));
%                 if nnz(B) == 1
%                     if isempty(intersect(getvariables(Mij),getvariables(obj)))
%                         if isempty(intersect(getvariables(Mij),involved))
%                             
%                             Fb = Fb + set(M{k+1}(i,i) == 0);
%                             any_remove  = 1;
%                             bajs = [bajs i];
%                         end
%                     end
%                 end
%             end
%             i = i+1;
%         end
%     end
% 
%     if 1
%         for sd = 1:1:2
%             [exponent_p,p_base] = getexponentbase(obj+sum(sum(rt)),recover(depends(obj)));
%             Msparsity = zeros(size(M{k+1},1));
%             vars = recover(depends(obj));
%             for i = 1:size(M{k+1},1)
%                 for j = i:size(M{k+1},1)
%                     [exponent_Mij,Mij_base] = getexponentbase(M{k+1}(i,j),vars);
%                     if ~isempty(findrows(exponent_p,exponent_Mij))
%                         Msparsity(i,j) = 1;
%                         Msparsity(j,i) = 1;
%                     end
%                 end
%             end
% 
%             Mblocked = [];
%             [p,q,r,s] = dmperm(Msparsity+eye(length(Msparsity)));
%             Mrs = M{k+1}(p,p);
%             r(2:14)=[];
%             MM = ones(length(M{k+1}));
%             blocks = zeros(1,length(M{k+1}));
%             for i = 1:length(r)-1
%                 blocks(r(i+1)-r(i)) = blocks(r(i+1)-r(i)) + 1;
%                 MM = blkdiag(MM,ones(r(i+1)-r(i)));
%                 Mblocked = blkdiag(Mblocked,Mrs(r(i):(r(i+1)-1),r(i):(r(i+1)-1)));
%             end
%             string = 'Blocks : ';
%             for i = 1:length(blocks)
%                 if blocks(i)>0
%                     string = [string num2str(i) 'x' num2str(i) '(' num2str(blocks(i)) ') ' ];
%                 end
%             end
%             disp(string)
%             rt = Mblocked(:);
%         end
%         M{k+1} = Mblocked;
%     end
% end

% ... and lower degree localization matrices
M{1} = 1;
for i = 1:1:k-1;
    n_i = factorial(n+k-i)/(factorial(n)*factorial(k-i));
    M{k-i+1} = M{k+1}(1:n_i,1:n_i);
end

% Lasserres relaxation (Lasserre, SIAM J. OPTIM, 11(3) 796-817)
Fmoments = set(M{k+1}>0);
for i = 1:length(vecConstraints)
    v_k = floor((degree(vecConstraints(i))+1)/2);
    Localizer = vecConstraints(i)*M{k-v_k+1};
    if isinequality(i)
        Fmoments = Fmoments+set(Localizer>0);
    else
        indicies = find(triu(ones(length(Localizer))));
        Fmoments = Fmoments+set(Localizer(indicies)==0);
    end
end
for i = 1:length(sdpConstraints)
    v_k = floor((degree(sdpConstraints{i})+1)/2);       
    Fmoments = Fmoments+set(kron(M{k-v_k+1},sdpConstraints{i})>0);
end

% Add them all
Fnew = Fnew + unblkdiag(Fmoments);

% No objective, minimize trace on moment-matrix instead
if isempty(obj)
    obj = trace(M{k+1});
end

% Solve 
sol = solvesdp(Fnew,obj,sdpsettings(options,'relax',1));%,'shift',1e-8));

% Construct SOS decompositions if the user wants these
if nargout >= 4 
    sosout.Q0 = dual(Fnew(1));
    sosout.v0 = u{end};
    sosout.p0 = u{end}'*dual(Fnew(1))*u{end};
    for i=1:length(Fnew)-1
        sosout.Qi{i} = dual(Fnew(i+1));
        sosout.vi{i} = u{end}(1:length(sosout.Qi{i}));
        sosout.pi{i} = sosout.vi{i}'*sosout.Qi{i}*sosout.vi{i};
    end
end
    
    
% Get the moment matrices
for i = 1:k+1
    moments{i} = relaxdouble(M{i});    
end

% Extract solutions if possible (at-least fesible and unbounded)
momentsstructure.moment  = moments;
momentsstructure.x       = x;
momentsstructure.monomials = u{k};
momentsstructure.n       = n;
momentsstructure.d       = max(1,ceil(max(deg)/2));
x_extract = {};
if nargout>=2 & ~(sol.problem == 1 | sol.problem == 2)
    momentsstructure.moment  = moments;
    momentsstructure.x       = x;
    momentsstructure.monomials = u{k};
    momentsstructure.n       = n;
    momentsstructure.d       = max(1,ceil(max(deg)/2));
    x_extract = extractsolution(momentsstructure,options);
end