newcatsdp.m

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

function [problem,integer_variables,binary_variables,parametric_variables,quad_info] = catsdp(F,G,h,relax)
%catsdp          Internal function: tries to determine the type of
%optimization problem

% Author Johan L鰂berg 
% $Id: newcatsdp.m,v 1.10 2005/04/27 15:47:18 joloef Exp $

Counter = size(F.clauses,2);
Ftype = zeros(Counter,1);

real_data = 1;
int_data  = 0;
bin_data  = 0;
par_data  = 0;

poly_constraint = 0;
bilin_constraint = 0;
sigm_constraint = 0;
rank_constraint = 0;
rank_objective = 0;

integer_variables = [];
binary_variables = [];
parametric_variables = [];
kyp_prob  = 0;

% ***********************************************
% Setup an empty problem definition
% ***********************************************
problem.objective.linear = 0;
problem.objective.quadratic.convex = 0;
problem.objective.quadratic.nonconvex = 0;
problem.objective.polynomial = 0;
problem.objective.maxdet = 0;
problem.objective.sigmonial = 0;

problem.constraint.equalities.linear     = 0;
problem.constraint.equalities.quadratic  = 0;
problem.constraint.equalities.polynomial = 0;
problem.constraint.equalities.sigmonial  = 0;

problem.constraint.inequalities.elementwise.linear = 0;
problem.constraint.inequalities.elementwise.quadratic.convex = 0;
problem.constraint.inequalities.elementwise.quadratic.nonconvex = 0;
problem.constraint.inequalities.elementwise.sigmonial = 0;
problem.constraint.inequalities.elementwise.polynomial = 0;

problem.constraint.inequalities.semidefinite.linear = 0;
problem.constraint.inequalities.semidefinite.quadratic = 0;
problem.constraint.inequalities.semidefinite.polynomial = 0;
problem.constraint.inequalities.semidefinite.sigmonial = 0;

problem.constraint.inequalities.rank = 0;

problem.constraint.inequalities.secondordercone = [];
problem.constraint.inequalities.rotatedsecondordercone = [];

problem.constraint.integer = 0;
problem.constraint.binary = 0;

problem.complex = 0;
problem.parametric = 0;

% ********************************************************
% Make a list of all globally available discrete variables
% ********************************************************
integer_variables = yalmip('intvariables');
binary_variables  = yalmip('binvariables');
for i = 1:Counter
    switch F.clauses{i}.type    
    case 7
        integer_variables = union(integer_variables,getvariables(F.clauses{i}.data));
    case 8
        binary_variables = union(binary_variables,getvariables(F.clauses{i}.data));        
    case 13
        parametric_variables = union(parametric_variables,getvariables(F.clauses{i}.data));        
    otherwise
    end   
end

% ********************************************************
% Rank variables
% ********************************************************
rank_variables = yalmip('rankvariables');
any_rank_variables = ~isempty(rank_variables);

% ********************************************************
% Make a list of all globally available nonlinear variables
% ********************************************************
[monomtable,variabletype] = yalmip('monomtable');

%temp = (sum(abs(monomtable),2)==1) & (any(monomtable==1,2));
%linear_variables = find(temp);
%nonlinear_variables = find(~temp);
%sigmonial_variables = find(any(0>monomtable,2) | any(monomtable-fix(monomtable),2));
linear_variables = find(variabletype==0);
nonlinear_variables = find(variabletype~=0);
sigmonial_variables = find(variabletype==4);

allvars = getvariables(F);
any_nonlinear_variables =~isempty(find(ismembc(nonlinear_variables,allvars)));

%any_nonlinear_variables = ~isempty(find(ismembc(nonlinear_variables,allvars)));

any_discrete_variables = ~isempty(integer_variables) | ~isempty(binary_variables);

for i = 1:Counter
    
    Fi = F.clauses{i};
    % Only real-valued data?
    real_data = real_data & isreal(F.clauses{i}.data);
    
    % Any discrete variables used
    if any_discrete_variables
        int_data = int_data | any(ismember(getvariables(Fi.data),integer_variables));
        bin_data = bin_data | any(ismember(getvariables(Fi.data),binary_variables));
        par_data = par_data | any(ismember(getvariables(Fi.data),parametric_variables));
    end
    
    if any_rank_variables       
        rank_constraint = rank_constraint | any(ismember(getvariables(Fi.data),rank_variables));
    end  
    
    if ~any_nonlinear_variables % No nonlinearly parameterized constraints
        
        switch Fi.type
        case {1,9}
            problem.constraint.inequalities.semidefinite.linear = 1;
        case 2
            problem.constraint.inequalities.elementwise.linear = 1;
        case 3
            problem.constraint.equalities.linear = 1;
        case 4
            problem.constraint.inequalities.secondordercone = 1;
        case 5
            problem.constraint.inequalities.rotatedsecondordercone = 1;            
        otherwise
        end
    else
        % Can be nonlinear stuff
        vars = getvariables(Fi.data);
        usednonlins = find(ismembc(nonlinear_variables,vars));
        if ~isempty(usednonlins)
            usedsigmonials = find(ismember(sigmonial_variables,vars));
            if ~isempty(usedsigmonials)
                switch Fi.type
                case 1
                    problem.constraint.inequalities.semidefinite.sigmonial = 1;
                case 2
                    problem.constraint.inequalities.elementwise.sigmonial = 1;
                case 3
                    problem.constraint.equalities.sigmonial = 1;
                case 4
                    error('Sigmonial SOCP not supported');                     
                case 5
                    error('Sigmonial RSOCP not supported');                    
                otherwise
                    error('Report bug in problem classification (sigmonial constraint)');
                end
            else
                deg = degree(Fi.data);
                switch deg
                    
                case 1
                    switch Fi.type
                    case 1
                        problem.constraint.inequalities.semidefinite.linear = 1;
                    case 2
                        problem.constraint.inequalities.elementwise.linear = 1;
                    case 3
                        problem.constraint.equalities.linear = 1;
                    case 4
                        problem.constraint.inequalities.secondordercone = 1;
                    case 5
                        problem.constraint.inequalities.rotatedsecondordercone = 1;                                
                    otherwise
                        error('Report bug in problem classification (linear constraint)');
                    end
                case 2
                    switch Fi.type
                    case 1
                        problem.constraint.inequalities.semidefinite.quadratic = 1;
                    case 2
                        problem.constraint.inequalities.elementwise.quadratic.convex = 1;
                    case 3
                        problem.constraint.equalities.quadratic = 1;
                    case 4
                        error                                
                    case 5
                        error                                
                    otherwise
                        error('Report bug in problem classification (quadratic constraint)');
                    end
                otherwise
                    switch Fi.type
                    case 1
                        problem.constraint.inequalities.semidefinite.polynomial = 1;
                    case 2
                        problem.constraint.inequalities.elementwise.polynomial = 1;
                    case 3
                        problem.constraint.equalities.polynomial = 1;
                    case 4
                        %   problem.constraint.inequalities.secondordercone = 1;
                    case 5
                        %   problem.constraint.inequalities.rotatedsecondordercone = 1;
                    otherwise
                        error('Report bug in problem classification (polynomial constraint)');
                    end
                    
                end
            end
        else
            switch Fi.type
            case 1
                problem.constraint.inequalities.semidefinite.linear = 1;
            case 2
                problem.constraint.inequalities.elementwise.linear = 1;
            case 3
                problem.constraint.equalities.linear = 1;
            case 4
                problem.constraint.inequalities.secondordercone = 1;                    
            case 5
                problem.constraint.inequalities.rotatedsecondordercone = 1;                                                   
            case 7
                problem.constraint.integer = 1;
            case 8
                problem.constraint.binary = 1;                     
            otherwise
                error('Report bug in problem classification (linear constraint)');
            end
        end
    end
end

if int_data    
    problem.constraint.integer = 1;
end
if bin_data
    problem.constraint.binary = 1;
end
if ~real_data
    problem.complex = 1;
end
if rank_constraint
    problem.constraint.inequalities.rank = 1;
end

if (relax==1) | (relax==2)
    problem.constraint.integer = 0;
    problem.constraint.binary = 0;        
    int_data  = 0;
    bin_data  = 0;
    integer_variables = [];
    binary_variables = [];
end
if (relax==1) | (relax==3)    
    problem.constraint.equalities.linear = problem.constraint.equalities.linear | problem.constraint.equalities.quadratic | problem.constraint.equalities.polynomial | problem.constraint.equalities.sigmonial; 
    problem.constraint.equalities.quadratic = 0;
    problem.constraint.equalities.polynomial = 0;
    problem.constraint.equalities.sigmonial = 0;
    
    problem.constraint.inequalities.elementwise.linear = problem.constraint.inequalities.elementwise.linear | problem.constraint.inequalities.elementwise.quadratic.convex | problem.constraint.inequalities.elementwise.quadratic.nonconvex | problem.constraint.inequalities.elementwise.sigmonial | problem.constraint.inequalities.elementwise.polynomial;
    problem.constraint.inequalities.elementwise.quadratic.convex = 0;
    problem.constraint.inequalities.elementwise.quadratic.nonconvex = 0;
    problem.constraint.inequalities.elementwise.sigmonial   = 0;
    problem.constraint.inequalities.elementwise.polynomial  = 0;    
    
    problem.constraint.inequalities.semidefinite.linear =  problem.constraint.inequalities.semidefinite.linear | problem.constraint.inequalities.semidefinite.quadratic | problem.constraint.inequalities.semidefinite.polynomial | problem.constraint.inequalities.semidefinite.sigmonial;    
    problem.constraint.inequalities.semidefinite.quadratic  = 0;
    problem.constraint.inequalities.semidefinite.polynomial = 0;
    problem.constraint.inequalities.semidefinite.sigmonial  = 0;
    
    poly_constraint = 0;
    bilin_constraint = 0;
    sigm_constraint = 0;    
end


% Analyse the objective function
quad_info = [];
if (~isempty(h)) & ~is(h,'linear') &~(relax==1) &~(relax==3)
    if any(ismember(getvariables(h),sigmonial_variables))
        problem.objective.sigmonial = 1;
    else
        [Q,c,f,x,info] = quaddecomp(h);
        if info==0
            [R,p]=chol(Q);
            if p~=0
                if min(eig(Q))>=-1e-10
                    p=0;
                    [U,S,V]=svd(Q);
                    i = find(diag(S)>1e-10);
                    R = sqrt(S(1:max(i),1:max(i)))*V(:,1:max(i))';
                end
            end
            if p==0
                problem.objective.quadratic.convex = 1;
            else
                problem.objective.quadratic.nonconvex = 1;
            end
            quad_info.Q = Q;
            quad_info.c = c;
            quad_info.f = f;
            quad_info.x = x;
            quad_info.R = R;
            quad_info.p = p;
        else
            problem.objective.polynomial = 1;
        end
    end
else   
    problem.objective.linear = ~isempty(h);
end