monomialreduction.m

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

function exponent_m = monomialreduction(exponent_m,exponent_p,options,csclasses)
%MONOMIALREDUCTION  Internal function for monomial reduction in SOS programs

% Author Johan L鰂berg
% $Id: monomialreduction.m,v 1.20 2005/05/27 16:19:13 joloef Exp $


% **********************************************
% TRIVIAL REDUCTIONS (stupid initial generation)
% **********************************************
mindegrees = min(exponent_p,[],1);
maxdegrees = max(exponent_p,[],1);
mindeg = min(sum(exponent_p,2));
maxdeg = max(sum(exponent_p,2));
if size(exponent_m{1},2)==0 % Stupid case : set(sos(parametric))
   if options.verbose>0;disp('Initially 1 monomials in R^0');end
else
    if options.verbose>0;disp(['Initially ' num2str(sum(cellfun('prodofsize',exponent_m)/size(exponent_m{1},2))) ' monomials in R^' num2str(size(exponent_p,2))]);end
end
    
for i = 1:length(csclasses)  
    t = cputime;
    % THE CODE BELOW IS MESSY TO HANDLE SEVERAL BUGS IN MATLAB
    %too_large_term = any(exponent_m-repmat(maxdegrees/2,size(exponent_m,1),1)>0,2);% DOES NOT HANDLE ODD
    % POLYNIMIALS CORRECTLY
    a1 = full(ceil((1+maxdegrees)/2)); % 6.5.1 in linux freaks on sparse stuff...
    if isempty(a1)
        a1 = zeros(size(maxdegrees));
    end
    a2 = full(size(exponent_m{i},1));
    too_large_term = any(exponent_m{i}-repmat(a1,a2,1)>0,2);
    %too_small_term = any(exponent_m-repmat(mindegrees/2,size(exponent_m,1),1)<0,2);
    a1 = full(floor(mindegrees/2));
     if isempty(a1)
        a1 = zeros(size(mindegrees));
    end
    a2 = full(size(exponent_m{i},1));
    too_small_term = any(exponent_m{i}-repmat(a1,a2,1)<0,2);%x^2+xz
    %too_large_sum = any(sum(exponent_m,2)-maxdeg/2>0,2); % DOES NOT HANDLE ODD
    % POLYNIMIALS CORRECTLY
    too_large_sum = any(sum(exponent_m{i},2)-ceil((1+maxdeg)/2)>0,2);
    too_small_sum = any(sum(exponent_m{i},2)-mindeg/2<0,2);
    keep = setdiff1D((1:size(exponent_m{i},1)),find(too_large_term | too_small_term | too_large_sum | too_small_sum));
    exponent_m{i} = exponent_m{i}(keep,:);
    t = cputime-t;    
end
if options.verbose>1;disp(['Removing large/small............Keeping ' num2str(sum(cellfun('prodofsize',exponent_m)/size(exponent_m{1},2))) ' monomials (' num2str(t) 'sec)']);end

% ************************************************
% Homogenuous?
% ************************************************
if all(sum(exponent_p,2)==sum(exponent_p(1,:)))
    for i = 1:length(csclasses)
        j = csclasses{i};
        t = cputime;
        exponent_m{i} = exponent_m{i}(sum(exponent_m{i},2)==sum(exponent_p(1,:))/2,:);
        t = cputime-t;        
    end
    if options.verbose>1;disp(['Homogenuous form!...............Keeping ' num2str(sum(cellfun('prodofsize',exponent_m)/size(exponent_m{1},2))) ' monomials (' num2str(t) 'sec)']);end
end

% ************************************************
% DIAGONAL CONSISTENCY : MONOMIAL ONLY IN
% DIAGONAL, CONSTRAINED TO BE ZER0, CAN BE REMOVED
% ************************************************
if (options.sos.inconsistent==1) &  ~options.sos.csp
    t = cputime;
    keep = consistent(exponent_m{1},exponent_p);
    exponent_m{1} = exponent_m{1}(keep,:);
    t = cputime-t;
    if options.verbose>0;disp(['Diagonal inconsistensies........Keeping ' num2str(size(exponent_m{1},1)) ' monomials (' num2str(t) 'sec)']);end
end

% ***********************************************************
% NEWTON POLYTOPE CHECK
% ***********************************************************
if options.sos.newton
    t = cputime;
    for j = 1:length(csclasses)               
        try
            % Basic idea : Try to find separating hyperplane between Newton polytope and candidate,
            % for all candidates.
            % Pro : Numerical stability, polynomial (does NOT calculate H-representation of Newton polytope)
            % Con : Slightly slower than explicit convex hull calculation in most cases.
            %
            % For speed, assuming GLPKLMEX or CDDMEX available (full-fledged YALMIP too slow)
            no_lp_solved = 0;
            if 1
                keep = ismember(exponent_m{j}*2,exponent_p,'rows');
                dontkeep = 0*keep;
                ii = 0;
                for i = 1:length(keep)
                    if ~keep(i) & ~dontkeep(i)
                        q = exponent_m{j}(i,:)'*2;
                        IN=struct('obj',full([-q' 1]),'A',full([exponent_p -ones(size(exponent_p,1),1);q' -1]),'B',[zeros(size(exponent_p,1),1);1]);
                        OUT=cddmex('solve_lp',IN);
                        no_lp_solved = no_lp_solved  + 1;
                        ii = ii + 1;
                        if OUT.how==5
                            % BUG in CDDMEX
                            warning('Resorting to YALMIP to solve LP (cddmex gave strange result). This can be slow');
                            xx = sdpvar(length(IN.obj),1);
                            sol = solvesdp(set(IN.A*xx < IN.B),IN.obj*xx,sdpsettings('verbose',0));
                            no_lp_solved = no_lp_solved  + 1;
                            OUT.how = sol.problem==0;
                            OUT.objlp = double(IN.obj*xx);
                            OUT.xopt = double(xx);
                        end
                        if (OUT.how == 1 & OUT.objlp<0)
                            a = OUT.xopt(1:end-1);
                            b = OUT.xopt(end);
                            u = find(a'*2*exponent_m{j}' - b > sqrt(eps));
                            dontkeep(u) = 1;
                        end
                    end
                    if keep(i)
                        q = exponent_m{j}(i,:)'*2;
                    end
                end
                keep = find(~dontkeep);
            else
                V1=struct('V',full(exponent_p));
                uu = cddmex('hull',V1);
                A = uu.A;
                b = uu.B;
                keep = zeros(1,length(exponent_m{j}));
                for i = 1:length(exponent_m{j})

                    if all(A*(2*exponent_m{j}(i,:))'<=b)
                        keep(i)=1;
                    else
                        keep(i)=0;
                    end
                end
                keep = find(keep);
            end
        catch
            try
                keep = ismember(exponent_m{j}*2,exponent_p,'rows');
                dontkeep = 0*keep;
                ops.msglev = 0;
                CTYPE = repmat('U',full(size(exponent_p,1)+1),1);
                VTYPE = repmat('C',full(size(exponent_p,2)+1),1);
                for i = 1:length(keep)
                    if ~keep(i) & ~dontkeep(i)
                        q = exponent_m{j}(i,:)'*2;
                        A = [exponent_p -ones(size(exponent_p,1),1);q' -1];
                        b = [zeros(size(exponent_p,1),1);1];
                        c = [-q(:);1];
                        [xopt,FMIN,STATUS]=glpkmex(1,c,A,b,CTYPE,[],[],VTYPE,ops);
                        no_lp_solved = no_lp_solved  + 1;
                        if (STATUS == 180) & (FMIN<0)
                            a = xopt(1:end-1);
                            b = xopt(end);
                            u = find(a'*2*exponent_m{j}' - b > sqrt(eps));
                            dontkeep(u) = 1;
                        end
                    end
                    if keep(i)
                        q = exponent_m{j}(i,:)'*2;
                    end
                end
                keep = find(~dontkeep);

            catch
                if options.verbose>0;warning('Resorting to LINPROG to solve LP (tried glpk and cdd but failed). This can be slow...');end

                try
                    keep = ismember(exponent_m{j}*2,exponent_p,'rows');
                    dontkeep = 0*keep;
                    ops = optimset('display','off');
                    for i = 1:length(keep)
                        if ~keep(i) & ~dontkeep(i)
                            q = exponent_m{j}(i,:)'*2;
                            A = [exponent_p -ones(size(exponent_p,1),1);q' -1];
                            b = [zeros(size(exponent_p,1),1);1];
                            c = [-q(:);1];
                            [xopt,FMIN,STATUS] = linprog(c,A,b,[],[],[],[],[],ops);ii = ii + 1;
                            no_lp_solved = no_lp_solved  + 1;
                            if (STATUS > 0) & (FMIN<0)
                                a = xopt(1:end-1);
                                b = xopt(end);
                                u = find(a'*2*exponent_m{j}' - b > sqrt(eps));
                                dontkeep(u) = 1;
                            end
                        end
                        if keep(i)
                            q = exponent_m{j}(i,:)'*2;
                        end
                    end
                    keep = find(~dontkeep);

                catch
                    disp(' ');
                    disp('Could not find cddmex, glpkmex or linprog (or it failed). Switching to crappy convhull for Newton polytope...');
                    disp(' ');
                    keep = newtonpolytope(exponent_m{j}*2,exponent_p);
                end
            end
        end

        changed = length(keep)<size(exponent_m{j},1);
        exponent_m{j} = exponent_m{j}(keep,:);       
    end
    t = cputime-t;
    %if options.verbose>1 | (options.verbose>0 & changed);disp(['Newton polytope.................Keeping ' num2str(size(exponent_m{j},1)) ' monomials (' num2str(t) 'sec)']);end
    
    if options.verbose>1 | (options.verbose>0 & 1+changed);
        info = ['Newton polytope (' num2str(no_lp_solved) ' LPs)'];
        info = [info repmat('.',1,32-length(info))];
        disp([info 'Keeping ' num2str(size(exponent_m{j},1)) ' monomials (' num2str(t) 'sec)']);
    end
end

if (options.sos.inconsistent==2) &  ~options.sos.csp
    t = cputime;
    keep = consistent(exponent_m{1},exponent_p);
    exponent_m{1} = exponent_m{1}(keep,:);
    t = cputime-t;
    if options.verbose>1;disp(['Diagonal inconsistensies........Keeping ' num2str(size(exponent_m,1)) ' monomials (' num2str(t) 'sec)']);end
end