solvesos.m

求解线性矩阵不等式简单方便--与LMI工具箱相比

        %         i = 1e-16;
        %         [R,fail] = chol(Z);
        %         while fail & i<10
        %             i = i*100;
        %             [R,fail] = chol(Z+i*eye(length(Z)));
        %         end
        %         Z = R(1,2:end)';
        %         assign(recover(ParametricVariables),Z);
        %         start = 1;
        %         for constraint = 1:length(p)
        %             Qi = [];
        %             for j = 1:length(BlockedA{constraint})
        %                 Qi = blkdiag(Qi,dual(F(start)));
        %                 start = start + 1;
        %             end
        %             Q{constraint} = Qi;
        %         end

        % *********************************************
        %% GENERATE MONOMIALS IN SOS-DECOMPOSITION
        % *********************************************
        for constraint = 1:length(p)

            if constraint > 1 & isequal(BlockedN{constraint},BlockedN{constraint-1}) & isequal(Blockedx{constraint},Blockedx{constraint-1}) & isequal(Blockedvarchange{constraint},Blockedvarchange{constraint-1}) & isequal(sizep(constraint),sizep(constraint-1))
                monoms{constraint} = monoms{constraint-1};
            else
                monoms{constraint} = [];
                totalN{constraint} = [];
                N = BlockedN{constraint};
                x = Blockedx{constraint};
                for i = 1:length(N)
                    % Original variable
                    for j = 1:size(N{i},1)
                        N{i}(j,:)=N{i}(j,:).*Blockedvarchange{constraint};
                    end
                    if isempty(N{i})
                        monoms{constraint} = [monoms{constraint};[]];
                    else
                        mi = kron(eye(sizep(constraint)),recovermonoms(N{i},x));
                        monoms{constraint} = [monoms{constraint};mi];
                    end
                end
                if isempty(monoms{constraint})
                    monoms{constraint}=1;
                end
            end

            % For small negative eigenvalues
            % this is a good quick'n'dirty approximation
            % Improve...use shifted eigenvalues and chol or what ever...
            if ~any(any(isnan(Q{constraint})))
                if isempty(Q{constraint})
                    Q{constraint}=0;
                    h{constraint}=0;
                else
                    usedVariables = find(any(Q{constraint},2));
                    if length(usedVariables)<length(Q{constraint})
                        Qpart = Q{constraint}(usedVariables,usedVariables);
                        [U,S,V]=svd(Qpart);
                        R = sqrt(S)*V';
                        h0 = R*monoms{constraint}(usedVariables);
                        if isa(h0,'sdpvar')
                            h{constraint} = clean(R*monoms{constraint}(usedVariables),options.sos.clean);
                            h{constraint} = h{constraint}(find(h{constraint}));
                        else
                            h{constraint} = h0;
                        end
                    else
                        [U,S,V]=svd(mid(Q{constraint}));
                        R = sqrt(S)*V';
                        h0 = R*monoms{constraint};

                        if isa(h0,'sdpvar')
                            h{constraint} = clean(R*monoms{constraint},options.sos.clean);
                            h{constraint} = h{constraint}(find(sum(h{constraint},2)),:);
                        else
                            h{constraint} = h0;
                        end
                    end
                end
                if isempty(ParametricVariables)
                    ParametricVariables = [];
                end
                setsos(p{constraint},h{constraint},ParametricVariables,Q{constraint},monoms{constraint});
            else
                if options.verbose>0;
                    if UncertainData
                        disp(' ');
                        disp('-> Only partial decomposition is returned (since you have uncertain data).');
                        disp('');
                    else
                        disp(' ');
                        disp('-> FAILURE : SOS decomposition not available.');
                        disp('-> The reason is probably that you are using a solver that does not deliver a dual (LMILAB)');
                        disp('-> Use sdsettings(''sos.model'',2) to circumvent this, or use another solver (SDPT3, SEDUMI,...)');
                        disp('');
                        disp('-> An alternative reason is that YALMIP detected infeasibility during the compilation phase.');
                    end
                end
            end
        end

        m = monoms;

    otherwise
        Q = [];
        m = [];
end

% Don't need these outside
yalmip('cleardual')

% Done with YALMIP, this is the time it took, minus solver
if ~isfield(sol,'solvertime')
    sol.solvertime = 0;
end

sol.yalmiptime = etime(clock,yalmip_time)-sol.solvertime;


function p_base_parametric = parameterizedbase(p,z, params,ParametricIndicies,exponent_p,p_base)

% Check for linear parameterization
parametric_basis = exponent_p(:,ParametricIndicies);
if all(sum(parametric_basis,2)==0)
    p_base_parametric = full(p_base(:));
    return
end
if all(sum(parametric_basis,2)<=1)
    p_base_parametric = full(p_base(:));
    n = length(p_base_parametric);


    if 1
        [ii,vars] = find(parametric_basis);
        ii = ii(:)';
        vars = vars(:)';
    else
        ii = [];
        vars = [];
        js = sum(parametric_basis,1);
        indicies = find(js);
        for i = indicies
            if js(i)
                j = find(parametric_basis(:,i));
                ii = [ii j(:)'];
                vars = [vars repmat(i,1,js(i))];
            end
        end
    end

    k = setdiff1D(1:n,ii);
    if isempty(k)
        p_base_parametric = p_base_parametric.*sparse(ii,repmat(1,1,n),params(vars));
    else
        pp = params(vars); % Must do this, bug in ML 6.1 (x=sparse(1);x([1 1]) gives different result in 6.1 and 7.0!)
        p_base_parametric = p_base_parametric.*sparse([ii k(:)'],repmat(1,1,n),[pp(:)' ones(1,1,length(k))]);
    end
else
    % Bummer, nonlinear parameterization sucks...
    for i = 1:length(p_base)
        j = find(exponent_p(i,ParametricIndicies));
        if ~isempty(j)
            temp = p_base(i);

            for k = 1:length(j)
                if exponent_p(i,ParametricIndicies(j(k)))==1
                    temp = temp*params(j(k));
                else
                    temp = temp*params(j(k))^exponent_p(i,ParametricIndicies(j(k)));
                end
            end
            xx{i} = temp;
        else
            xx{i} = p_base(i);
        end
    end
    p_base_parametric = stackcell(sdpvar(1,1),xx)';
end



function [A,b] = generate_kernel_representation_data(N,N_unique,exponent_m2,exponent_p,p,options,p_base_parametric,ParametricIndicies,MonomIndicies,FirstRun)

persistent saveData

exponent_p_parametric = exponent_p(:,ParametricIndicies);
exponent_p_monoms = exponent_p(:,MonomIndicies);
pcoeffs = getbase(p);
if any(exponent_p_monoms(1,:))
    pcoeffs=pcoeffs(:,2:end); % No constant term in p
end
b = [];

parametric = full((~isempty(ParametricIndicies) & any(any(exponent_p_parametric))));

% For problems with a lot of similar cones, this saves some time
reuse = 0;
if ~isempty(saveData) & isequal(saveData.N,N) & ~FirstRun
    n = saveData.n;
    ind = saveData.ind;
    if  isequal(saveData.N_unique,N_unique) & isequal(saveData.exponent_m2,exponent_m2)% & isequal(saveData.epm,exponent_p_monoms)
        reuse = 1;
    end
else
    % Congruence partition sizes
    for k = 1:size(N,1)
        n(k) = size(N{k},1);
    end
    % Save old SOS definition
    saveData.N = N;
    saveData.n = n;
    saveData.N_unique = N_unique;
    saveData.exponent_m2 = exponent_m2;
    saveData.N_unique = N_unique;
end

if reuse & options.sos.reuse
    % Get old stuff
    if size(exponent_m2{1},2)==2 % Stupid set(sos(parametric)) case
        ind = spalloc(1,1,0);
        ind(1)=1;
        allj = 1:size(exponent_p_monoms,1);
        used_in_p = ones(size(exponent_p_monoms,1),1);
    else
        allj = [];
        used_in_p = zeros(size(exponent_p_monoms,1),1);
        hash = randn(size(exponent_p_monoms,2),1);
        exponent_p_monoms_hash = exponent_p_monoms*hash;
        for i = 1:size(N_unique,1)
            monom = sparse(N_unique(i,3:end));
            j = find(exponent_p_monoms_hash == (monom*hash));
          
            if isempty(j)
                b = [b 0];
                allj(end+1,1) = 0;
            else
                used_in_p(j) = 1;
                allj(end+1,1:length(j)) = j(:)';
            end
        end
        ind = saveData.ind;
    end
else
    allj = [];
    used_in_p = zeros(size(exponent_p_monoms,1),1);
    if size(exponent_m2{1},2)==2 % Stupid set(sos(parametric)) case
        ind = spalloc(1,1,0);
        ind(1)=1;
        allj = 1:size(exponent_p_monoms,1);
        used_in_p = ones(size(exponent_p_monoms,1),1);
    else
        % To speed up some searching, we random-hash data
        hash = randn(size(exponent_p_monoms,2),1);
        for k = 1:length(exponent_m2)
            if isempty(exponent_m2{k})
                exp_hash{k}=[];
            else
                exp_hash{k} = sparse((exponent_m2{k}(:,3:end)))*hash; % SPARSE NEEDED DUE TO STRANGE NUMERICS IN MATLAB ON 0s (the stuff will differ on last bit in hex format)
            end
        end

        p_hash = exponent_p_monoms*hash;
        ind = spalloc(size(N_unique,1),sum(n.^2),0);       

        for i = 1:size(N_unique,1)
            monom = N_unique(i,3:end);
            monom_hash = sparse(monom)*hash;
            LHS = 0;
            start = 0;
            for k = 1:size(N,1)
                j = find(exp_hash{k} == monom_hash);
                if ~isempty(j)
                    pos=exponent_m2{k}(j,1:2);
                    nss = pos(:,1);
                    mss = pos(:,2);
                    indicies = nss+(mss-1)*n(k);
                    ind(i,indicies+start) = ind(i,indicies+start) + 1;                         
                end
                start = start + (n(k))^2;
                %                start = start + (matrixSOSsize*n(k))^2;
            end

            j = find(p_hash == monom_hash);

            if isempty(j)
                allj(end+1,1) = 0;
            else
                used_in_p(j) = 1;
                allj(end+1,1:length(j)) = j(:)';
            end
        end
    end
end
saveData.ind = ind;

% Some parametric terms in p(x,t) do not appear in v'Qv
% So these have to be added 0*Q = b
not_dealt_with  = find(used_in_p==0);
while ~isempty(not_dealt_with)
    j = findrows(exponent_p_monoms,exponent_p_monoms(not_dealt_with(1),:));
    allj(end+1,1:length(j)) = j(:)';
    used_in_p(j) = 1;
    not_dealt_with  = find(used_in_p==0);
    ind(end+1,1)=0;
end

matrixSOSsize = length(p);
if parametric
    % Inconsistent behaviour in MATLAB
    if size(allj,1)==1
        uu = [0;p_base_parametric];
        b = sum(uu(allj+1))';
    else
        b = [];
        for i = 1:matrixSOSsize
            for j = i:matrixSOSsize
                if i~=j
                    uu = [0;2*p_base_parametric(:,(i-1)*matrixSOSsize+j)];
                else
                    uu = [0;p_base_parametric(:,(i-1)*matrixSOSsize+j)];
                end
                b = [b sum(uu(allj+1),2)'];
            end
        end
    end
else
    if matrixSOSsize == 1
        uu = [zeros(size(pcoeffs,1),1) pcoeffs]';
        b = sum(uu(allj+1,:),2)';
    else
        b = [];
        for i = 1:matrixSOSsize
            for j = i:matrixSOSsize
                if i~=j
                    uu = [0;2*pcoeffs((i-1)*matrixSOSsize+j,:)'];
                else
                    uu = [0;pcoeffs((i-1)*matrixSOSsize+j,:)'];
                end
                b = [b;sum(uu(allj+1,:),2)'];
            end
        end
    end
    % uu = [0;pcoeffs(:)];
    % b = sum(uu(allj+1),2)';
end

b = b';
dualbase = ind;

j = 1;
A = cell(size(N,1),1);
for k = 1:size(N,1)
    if matrixSOSsize==1
        A{k} = dualbase(:,j:j+n(k)^2-1);
    else
        % Quick fix for matrix SOS case, should be optimized
        A{k} = inflate(dualbase(:,j:j+n(k)^2-1),matrixSOSsize,n(k));
    end
    j = j + n(k)^2;
end
b = b(:);



function newAi = inflate(Ai,matrixSOSsize,n);
% Quick fix for matrix SOS case, should be optimized
newAi = [];
for i = 1:matrixSOSsize
    for r = i:matrixSOSsize
        for m = 1:size(Ai,1)
            ai = reshape(Ai(m,:),n,n);
            V = zeros(matrixSOSsize,matrixSOSsize);
            V(i,r)=1;