solvesos.m

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

function [sol,m,Q,residuals,everything] = solvesos(F,obj,options,params,candidateMonomials)
%SOLVESOS Sum of squares decomposition
%
%   [sol,m,B] = solvesos(F,h,options,params,monomials) is used for finding
%   SOS-decompositions of polynomials.
%
%    The coefficients of the polynomials are aussumed linear w.r.t a set of decision
%    variables 'params' and polynomial with respect to a variable x.
%
%    An extension with a nonlinear parameterization in params is possible.
%    Note though that this gives BMIs or PMIs, solvable (locally) only if
%    PENBMI is installed, or by specifying 'moment' as solver to try to
%    solve the nonconvex semidefinite programming problem using a
%    semidefinite relaxation based on moments.
%
%    The SOS problem can be formulated as
%
%              min h(params)
%
%       subject to  F(i) >(=) 0 or F(i) is SOS w.r.t x
%
%   INPUT
%    F         : SET object with SOS-constrained polynomials and constraints on variables params
%    h         : scalar SDPVAR object (can be [])
%    options   : options structure obtained from SDPSETTINGS (can be [])
%    params    : SDPVAR object defining parametric variables (can be [])
%    monomials : SDPVAR object with user-specified monomials for decomposition (can be [])
%
%   OUTPUT
%    sol       : Solution diagnostic from SDP problem
%    v         : Cell with monomials used in decompositions
%    Q         : Cell with Gram matrices, p = v{i}'*Q{i}*v{i}, where p is the ith SOS polynomial in your model.
%    residuals : Mismatch between p and decompositions. Same values (modulo numerical issue) as checkset(find(is(F,'sos')))
%
%   EXAMPLE
%    x = sdpvar(1);solvesos(set(sos(x^4+x^3+1)));                    % Simple decompositions
%    x = sdpvar(1);t = sdpvar(1);solvesos(set(sos(x^4+x^3+1-t)),-t); % Lower bound by maximizing t
%
%   NOTES
%
%    Variables not part of params, but part of non-SOS constraints in F
%    or objective h will automatically be appended to the params list.
%
%    To extract SOS decomposition, use command SOSD
%
%    The field options.sos can be used to tune the SOS-calculations. See HTML help for details
%
%     sos.model          - Kernel or image representation of SOS problem [0|1|2 (0)]
%     sos.newton         - Use Newton polytope to reduce size [0|1 (1)]
%     sos.congruence     - Block-diagonalize using congruence classes [0|1|2 (2)]
%     sos.scale          - Scale polynomial [0|1 (1)]
%     sos.extlp          - Extract simple translated LP cones when performing dualization [0|1 (1)]
%     sos.inconsistent   - Remove diagonal-inconsistent monomials [0|1|2 (0)]
%     sos.clean          - Remove monomials with coefficients < clean [real > 0 (1e-4)]
%     sos.traceobj       - Minimize trace of Gram matrix in problems without objective function [0|1 (0)]
%
% See also SOSD, SDPSETTINGS, SOLVEMOMENT, SDPVAR, SDISPLAY

%% Time YALMIP
yalmip_time = clock;

% ************************************************
%% Check #inputs
% ************************************************
if nargin<5
    candidateMonomials = [];
    if nargin<4
        params = [];
        if nargin<3
            options = sdpsettings;
            if nargin<2
                obj = [];
                if nargin<1
                    help solvesos
                    return
                end
            end
        end
    end
end
% Lazy syntax (not official...)
if nargin==1 & isa(F,'sdpvar')
    F = set(sos(F));
end

% ****************************************************
%% Extract all SOS constraints and candidate monomials
% ****************************************************
if ~any(is(F,'sos'))
    error('At-least one constraint should be an SOS constraints!');
end
p = [];
ranks = [];
for i = 1:length(F)
    if is(F(i),'sos')
        pi = sdpvar(F(i));
        p{end+1} = pi;
        ranks(end+1) = getsosrank(pi); % Desired rank : Experimental code
    end    
end
if isempty(candidateMonomials)
    for i = 1:length(F)
        candidateMonomials{i}=[];
    end
elseif isa(candidateMonomials,'sdpvar')
    cM=candidateMonomials;
    candidateMonomials={};
    for i = 1:length(p)
        candidateMonomials{i}=cM;
    end
elseif isa(candidateMonomials,'cell')
    if length(p)~=length(candidateMonomials)
        error('Dimension mismatch between the candidate monomials and the number of SOS constraints');
    end
end

% ************************************************
%% Get the parametric constraints
% ************************************************
F_original = F;
F_parametric = F(find(~is(F,'sos')));
if isempty(F_parametric)
    F_parametric = set([]);
end

if isempty(options)
    options = sdpsettings;
end

% ************************************************
%% Get model for nonlinear operators if we are
% going for an kernel model
% ************************************************
expanded_nonlinear_operators = 0;
if ~isempty(yalmip('extvariables')) & options.sos.model == 1
    [F_parametric,failure] = convexitypropagation(F_parametric,obj);  
    expanded_nonlinear_operators = 1;
    if failure
        Q{1} = [];
        m{1} = [];
        residuals = [];
        everything = [];
        sol.yalmiptime = etime(clock,yalmip_time);
        sol.solvertime = 0;
        sol.info = yalmiperror(14,'YALMIP');
        sol.problem = 14;
    end
end

if ~isempty(params)
    if ~isa(params,'sdpvar')
        error('Fourth argment should be a SDPVAR variable or empty')
    end
end

if options.verbose>0;
    disp('-------------------------------------------------------------------------');
    disp('YALMIP SOS module started...');
    disp('-------------------------------------------------------------------------');
end

% ************************************************
%% INITIALIZE SOS-DECOMPOSITIONS SDP CONSTRAINTS
% ************************************************
F_sos = set([]);

% ************************************************
%% FIGURE OUT ALL PARAMETRIC VARIABLES
% ************************************************
AllVariables =  uniquestripped([depends(obj) depends(F_original)]);
ParametricVariables = uniquestripped([depends(obj) depends(F_parametric) depends(params)]);
ParametricVariables = intersect(ParametricVariables,AllVariables);
MonomVariables = setdiff(AllVariables,ParametricVariables);
params = recover(ParametricVariables);
if isempty(MonomVariables)
    error('No independent variables? Perhaps you added a constraint set(p(x)) when you meant set(sos(p(x)))');
end
if options.verbose>0;disp(['Detected ' num2str(length(ParametricVariables)) ' parametric variables and ' num2str(length(MonomVariables)) ' independent variables.']);end

% ************************************************
%% ANY BMI STUFF
% ************************************************
NonLinearParameterization = 0;
if ~isempty(ParametricVariables)
    monomtable = yalmip('monomtable');
    ParametricMonomials = monomtable(uniquestripped([getvariables(obj) getvariables(F_original)]),ParametricVariables);
    if any(sum(abs(ParametricMonomials),2)>1)
        NonLinearParameterization = 1;
    end
end

% ************************************************
%% DISPLAY WHAT WE FOUND
% ************************************************
if options.verbose>0 & ~isempty(F_parametric)
    nLP = 0;
    nEQ = 0;
    nLMI = sum(full(is(F_parametric,'lmi')) &  full(~is(F_parametric,'element-wise'))); %FULL due to bug in ML 7.0.1
    for i = 1:length(F_parametric)
        if is(F_parametric,'element-wise')
            nLP = nLP + prod(size(F_parametric(i)));
        end
        if is(F_parametric,'equality')
            nEQ = nEQ + prod(size(F_parametric(i)));
        end
    end
    disp(['Detected ' num2str(full(nLP)) ' linear inequalities, ' num2str(full(nEQ)) ' equality constraints and ' num2str(full(nLMI)) ' LMIs.']);
end

% ************************************************
%% IMAGE OR KERNEL REPRESENTATION?
% ************************************************
noRANK = all(isinf(ranks));
switch options.sos.model
case 0
    constraint_classes = constraintclass(F);
    noCOMPLICATING = ~any(ismember([7 8 9 10 12 13 14],constraint_classes));    
    if noCOMPLICATING & ~NonLinearParameterization & noRANK
        options.sos.model = 1;
        if options.verbose>0;disp('Using kernel representation (model 1).');end
    else
        if NonLinearParameterization
            if options.verbose>0;disp('Using image representation (model 2). Nonlinear parameterization found');end
        elseif ~noRANK
            if options.verbose>0;disp('Using image representation (model 2). SOS-rank constraint was found.');end
        else
            if options.verbose>0;disp('Using image representation (model 2). Integer data, KYPs or similar was found.');end
        end                
        options.sos.model = 2;            
    end
case 1
    if NonLinearParameterization
        if options.verbose>0;disp('Switching to image model due to nonlinear parameterization (not supported in kernel model).');end
        options.sos.model = 2;
    end
    if ~noRANK
        if options.verbose>0;disp('Switching to image model due to SOS-rank constraints (not supported in kernel model).');end
        options.sos.model = 2;
    end    
case 3
    
otherwise
end

% ************************************************
%% Damn, we switched model
% ************************************************
if ~isempty(yalmip('extvariables')) & options.sos.model == 1 & ~expanded_nonlinear_operators
    [F_parametric,failure] = convexitypropagation(F_parametric,obj);  
    expanded_nonlinear_operators = 1;
    if failure
        Q{1} = [];
        m{1} = [];
        residuals = [];
        everything = [];
        sol.yalmiptime = etime(clock,yalmip_time);
        sol.solvertime = 0;
        sol.info = yalmiperror(14,'YALMIP');
        sol.problem = 14;
    end
end

% ************************************************
%% SKIP DIAGONAL INCONSITENCY FOR PARAMETRIC MODEL
% ************************************************
if ~isempty(params) & options.sos.inconsistent
    if options.verbose>0;disp('Turning off inconsistency based reduction (not supported in parametric models).');end
    options.sos.inconsistent = 0;
end

% ************************************************
%% INITIALIZE OBJECTIVE
% ************************************************
if ~isempty(obj)
    options.sos.traceobj = 0;
end
parobj = obj;
obj    = [];

% ************************************************
%% SCALE SOS CONSTRAINTS
% ************************************************
if options.sos.scale
    for constraint = 1:length(p)
        normp(constraint) = sqrt(norm(getbase(p{constraint})));
        p{constraint} = p{constraint}/normp(constraint);
    end
else
    normp = ones(length(p),1);
end

% ************************************************
%% SKIP CONGRUENCE REDUCTION WHEN SOS-RANK 
% ************************************************
options.sos.congruence = options.sos.congruence & all(isinf(ranks));

% ************************************************
%% LOOP THROUGH ALL SOS CONSTRAINTS
% ************************************************
for constraint = 1:length(p)
    % *********************************************
    %% FIND THE VARIABLES IN p, SORT, GET UNIQUE ETC
    % *********************************************
    if options.verbose>1;disp(['Creating SOS-description ' num2str(constraint) '/' num2str(length(p)) ]);end
    
    pVariables = depends(p{constraint});
    AllVariables = uniquestripped([pVariables ParametricVariables]);
    MonomVariables = setdiff1D(pVariables,ParametricVariables);
    x = recover(MonomVariables);
    z = recover(AllVariables);
    MonomIndicies = find(ismember(AllVariables,MonomVariables));
    ParametricIndicies = find(ismember(AllVariables,ParametricVariables));
    
    if isempty(MonomIndicies)
        error('You have constraints of the type set(sos(f(parametric_variables))). Please use set(f(parametric_variables) > 0) instead')
    end
    
    % *********************************************
    %% Express p in monimials and coefficients
    % *********************************************
    [exponent_p,p_base] = getexponentbase(p{constraint},z);
    
    % *********************************************
    %% Powers for user defined candidate monomials
    % (still experimental)
    % *********************************************
    [exponent_c] = getexponentbase(candidateMonomials{constraint},z);    
    if ~isempty(exponent_c)
        exponent_c = exponent_c(:,MonomIndicies);
    end
    
    % *********************************************
    %% STUPID PROBLEM WITH ODD HIGHEST POWER?...
    % *********************************************
    if isempty(ParametricIndicies)
        max_degrees = max(exponent_p(:,MonomIndicies),[],1);
        bad_max = any(max_degrees-fix((max_degrees/2))*2);
        if bad_max
            for i = 1:length(p)
                Q{i}=[];
            end
            
            sol.yalmiptime = etime(clock,yalmip_time);
            sol.solvertime = 0;
            sol.info = yalmiperror(1,'YALMIP');
            sol.problem = 2;
            return
        end
    end
    
    % *********************************************
    %% Can we make a variable smart change (no code)
    % *********************************************
    exponent_p_monoms = exponent_p(:,MonomIndicies);
    varchange = ones(1,size(MonomIndicies,2));
    
    % *********************************************
    %% Unique monoms (copies due to parametric terms)
    % *********************************************
    exponent_p_monoms = uniquesafe(exponent_p_monoms,'rows');
    
    if options.sos.reuse & constraint > 1 & isequal(previous_exponent_p_monoms,exponent_p_monoms)
        % We don't have to do anything, candidate monomials can be-used
        if options.verbose>1;disp(['Re-using all candidate monomials (same problem structure)']);end
    else
        
        % *********************************************
        %% CORRELATIVE SPARSITY PATTERN
        % *********************************************
        [C,csclasses] = corrsparsity(exponent_p_monoms,options);
        
        % *********************************************
        %% GENERATE MONOMIALS
        % *********************************************
        exponent_m = monomialgeneration(exponent_p_monoms,csclasses);
        
        % *********************************************
        %% PRUNE W.R.T USER DEFINED
        % *********************************************