congruenceblocks.m

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

function N = congruenceblocks(exponent_m,exponent_p,options,csclasses)
%CONGRUENCEBLOCKS Internal function to reduce monomials in SOS problem

% Author Johan L鰂berg
% $Id: congruenceblocks.m,v 1.8 2005/02/28 20:59:14 johanl Exp $

if ~isempty(exponent_m{1}) & options.sos.congruence>0 & ((size(exponent_p,2)<=16)  | options.sos.congruence==1)

    % **********************************************
    % DEFINE CONGRUENCE CLASSES
    % **********************************************
    if options.verbose>0;fprintf('Finding symmetries..............');end;
    n = size(exponent_p,2);
    t = cputime;
    switch options.sos.congruence
        case 1
            Htemp = eye(n);                 % CHEAP VERSION; ONLY CHECK IF IT IS EVEN WRT x_i
        case 2
            Htemp = dec2decbin(1:2^n-1,n)'; % ALL POSSIBLE COMBINATIONS
        otherwise
            error('sos.congruence should be 0, 1 or 2')
    end
    H = Htemp(:,find(~any(rem(exponent_p*Htemp,2),1))); % Find "even" rows
    if isempty(H)
        N = exponent_m;
        if options.verbose>0;disp(['Found no symmetries (' num2str(cputime-t) 'sec)']);end
        return
    end

    t = cputime-t;
    if size(H,2)>=1
        if options.verbose>0
            if size(H,2)>1
                disp(['Found ' num2str(size(H,2)) ' symmetries  (' num2str(t) 'sec)']);
            else
                disp(['Found ' num2str(size(H,2)) ' symmetry  (' num2str(t) 'sec)']);
            end
        end
    else
        if options.verbose>0;disp(['Found no symmetries  (' num2str(t) 'sec)']);end;
    end

    % **********************************************
    % CLASSIFY MONMS ACCORDING TO CONGRUENCE CLASSES
    % **********************************************
    if size(H,2)>=1
        the_text = 'Partitioning using symmetry.....';
    end
    N = cell(0,1);
    for cs = 1:length(csclasses)
        [ur,j,k]=uniquesafe(mod(exponent_m{cs}*H,2),'rows');
        Ntemp = cell(size(ur,1),1);
        temp = [];
        for i = 1:length(k)
            Ntemp{k(i),1} = [Ntemp{k(i)};exponent_m{cs}(i,:)];
            temp = [temp;size(exponent_m{cs}(i,:),1)];
        end

        for i = 1:length(Ntemp)
            N{end+1,1} = Ntemp{i};
        end
    end
    % **********************************************
    % PRINT SOME RESULTS
    % **********************************************
    if size(H,2)>=1
        [uu,ii,oo] = uniquesafe(cellfun('prodofsize',N)/size(N{1},2));
        for i = 1:length(uu)
            n_this = length(find(oo==i));

            the_text = [the_text num2str(uu(i)) 'x' num2str(uu(i)) '(' num2str(n_this) ')' ' '];
        end
    end
    if options.verbose>0;;disp(the_text);end;
else
    N = exponent_m;
end