solvesos.m

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

            V(r,i)=1;
            ai = kron(V,ai);
            newAi = [newAi;ai(:)'];
        end
    end
end
        

function [Z,Q1,R] = sparsenull(A)

[Q,R] = qr(A');
n = max(find(sum(abs(R),2)));
Q1 = Q(:,1:n);
R = R(1:n,:);
Z = Q(:,n+1:end); % New basis


function [F,obj,BlockedQ,Primal_matrices,Free_variables] = create_kernelmodel(BlockedA,Blockedb,F_parametric,parobj,options,sparsityPattern);

% To get the primal kernel representation, we simply use
% the built-in dualization module!
% First, write the problem in primal kernel format
traceobj = 0;
dotraceobj = options.sos.traceobj;
F = F_parametric;
for i = 1:length(Blockedb)


    sizematrixSOS = sqrt(size(Blockedb{i},2));
    for k = 1:sizematrixSOS
        for r = k:sizematrixSOS
            res{(k-1)*sizematrixSOS+r} = 0;
        end
    end

    for j = 1:length(BlockedA{i})
        n = sqrt(size(BlockedA{i}{j},2));
        BlockedQ{i}{j} = sdpvar(n*sizematrixSOS,n*sizematrixSOS);
        F = F + set(BlockedQ{i}{j});
        if sizematrixSOS>0
            % Matrix valued sum of sqaures
            % Loop over all elements
            starttop = 1;
            for k = 1:sizematrixSOS
                startleft = 1;
                for r = 1:sizematrixSOS
                    if k<=r
                        Qkr = BlockedQ{i}{j}(starttop:starttop+n-1,startleft:startleft+n-1);
                        res{(k-1)*sizematrixSOS+r} = res{(k-1)*sizematrixSOS+r} + BlockedA{i}{j}*reshape(Qkr,n^2,1);
                    end
                    startleft = startleft + n;
                end
                starttop = starttop + n;
            end
        else
            % Standard case
            res{1} = res{1} + BlockedA{i}{j}*reshape(BlockedQ{i}{j},n^2,1);
        end
        if dotraceobj
            traceobj = traceobj + trace(BlockedQ{i}{j});
        end
    end
    for k = 1:sizematrixSOS
        for r = k:sizematrixSOS
            F = F + set(res{(k-1)*sizematrixSOS+r} == Blockedb{i}(:,(k-1)*sizematrixSOS+r));
        end
    end
end

% % Constrain elements according to a desired sparsity
if ~isempty(sparsityPattern)
    res = [];
    for i = 1:length(BlockedQ)
        for j = 1:length(BlockedQ{i})
            if ~isempty(sparsityPattern{i}{j})
                H = spalloc(length(BlockedQ{i}{j}),length(BlockedQ{i}{j}),length(sparsityPattern{i}{j}));
                H(sparsityPattern{i}{j}) = 1;
                k = find(triu(H));
                res = [res;BlockedQ{i}{j}(k)];
            end
        end
    end
    F = F + set(0 == res);
end

% And get the primal model of this
if isempty(parobj)
    if options.sos.traceobj
        [F,obj,Primal_matrices,Free_variables] = dualize(F,traceobj,1,options.sos.extlp);
    else
        [F,obj,Primal_matrices,Free_variables] = dualize(F,[],1,options.sos.extlp);
    end
else
    [F,obj,Primal_matrices,Free_variables] = dualize(F,parobj,1,options.sos.extlp);
end
% In dual mode, we maximize
obj = -obj;


function [F,obj,BlockedQ,sol] = create_imagemodel(BlockedA,Blockedb,F_parametric,parobj,options);


% *********************************
% IMAGE REPRESENTATION
% Needed for nonlinearly parameterized problems
% More efficient in some cases
% *********************************
g = [];
vecQ = [];
sol.problem = 0;
for i = 1:length(BlockedA)
    q = [];
    A = [];
    for j = 1:length(BlockedA{i})
        n = sqrt(size(BlockedA{i}{j},2));
        BlockedQ{i}{j} = sdpvar(n,n);
        q = [q;reshape(BlockedQ{i}{j},n^2,1)];
        A = [A BlockedA{i}{j}];
    end
    vecQ = [vecQ;q];
    g = [g;Blockedb{i}-A*q];
end

g_vars = getvariables(g);
q_vars = getvariables(vecQ);
x_vars = setdiff(g_vars,q_vars);


base = getbase(g);
if isempty(x_vars)
    A = base(:,1);base = base(:,2:end);
    B = (base(:,ismember(g_vars,q_vars)));
    Bnull = sparsenull(B);
    t = sdpvar(size(Bnull,2),1);
    imQ = -B\A+Bnull*t;
else
    A = base(:,1);base = base(:,2:end);
    C = base(:,ismember(g_vars,x_vars));
    B = (base(:,ismember(g_vars,q_vars)));
    if isa(B,'intval')
        if nnz(rad(B))>0
            error('Strange uncertainty, report this!');
        else
            B = mid(B);
        end
    end
    [Bnull,Q1,R1] = sparsenull(B);Bnull(abs(Bnull) < 1e-12) = 0;
    t = sdpvar(size(Bnull,2),1);
    if isa(A,'intval')
        % Limitation in INTLAB
        imQ = -Q1*(full(R1')\(A+C*recover(x_vars)))+Bnull*t;
    else
        imQ = -Q1*(R1'\(A+C*recover(x_vars)))+Bnull*t;
    end
end
notUsed = find(sum(abs(B),2)==0);
if ~isempty(notUsed)
    ff=g(notUsed);
    if isa(ff,'double')
        if nnz(ff)>0
            sol.yalmiptime = 0; % FIX
            sol.solvertime = 0;
            sol.problem = 2;
            sol.info = yalmiperror(1,'YALMIP');
            F = [];
            obj = [];
            if options.verbose > 0
                disp(' ');
                disp('-> During construction of data, I encountered a situation')
                disp('-> situation that tells me that the problem is trivially infeasible!')
                disp('-> Have you forgotten to define some parametric variables?,')
                disp('-> or perhaps you have a parametric problem where the highest')
                disp('-> power in some of the independent variables is odd for any choice')
                disp('-> of parametric variables, such as x^8+x^7+t*y^3')
                disp('-> Anyway, take a good look at your model again...');                
            end
            return
            %            error('You seem to have a strange model. Have you forgotten to define some parametric variable?');
        end
    else
       F_parametric = F_parametric + set(g(notUsed)==0);
    end
end
F_sos = set([]);
obj = 0;
for i = 1:length(BlockedQ)
    for j = 1:size(BlockedQ{i},2)
        Q_old = BlockedQ{i}{j};
        Q_old_vars = getvariables(Q_old);
        Q_old_base = getbase(Q_old);
        in_this = [];
        for k = 1:length(Q_old_vars)
            in_this = [in_this;find(Q_old_vars(k)==q_vars)];
        end
        Q_new = Q_old_base(:,1) + Q_old_base(:,2:end)*imQ(in_this);
        Q_new = reshape(Q_new,length(BlockedQ{i}{j}),length(BlockedQ{i}{j}));
        obj = obj+trace(Q_new);
        if ~isa(Q_new,'double')
            F_sos = F_sos + set(Q_new);
        elseif min(eig(Q_new))<-1e-8
            sol.yalmiptime = 0; % FIX
            sol.solvertime = 0;
            sol.problem = 2;
            sol.info = yalmiperror(1,'YALMIP');
            F = [];
            obj = [];
            error('Problem is trivially infeasible. After block-diagonalizing, I found constant negative definite blocks!');
            return
        end
        BlockedQ{i}{j} = Q_new;
    end
end

F = F_parametric + F_sos;

if isa(obj,'double') | (options.sos.traceobj == 0)
    obj = [];
end

if ~isempty(parobj)
    obj = parobj;
end


function   [F,obj,options] = create_lrmodel(BlockedA,Blockedb,F_parametric,parobj,options,ParametricVariables)
% Some special code for ther low-rank model in SDPLR
% Experimental code, not official yet
allb = [];
allA = [];
K.s = [];
for i = 1:length(Blockedb)
    allb = [allb;Blockedb{i}];
    Ai = [];
    for j = 1:size(BlockedA{i},2)
        Ai = [Ai BlockedA{i}{j}];
        K.s = [K.s sqrt(size(BlockedA{i}{j},2))];
    end
    %blkdiag bug in 7.0...
    [n1,m1] = size(allA);
    [n2,m2] = size(Ai);
    allA = [allA spalloc(n1,m2,0);spalloc(n2,m1,0) Ai];
end
options.solver = 'sdplr';
z = recover(ParametricVariables)
start = size(BlockedA,2)+1;
Mi = [];
for i = 1:length(allb)
    if isa(allb(i),'sdpvar')
        [Qi,ci,fi,xi,infoi] = quaddecomp(allb(i),z);
    else
        Qi = zeros(length(z));
        ci = zeros(length(z),1);
        fi = allb(i);
    end
    Z = -[fi ci'/2;ci/2 Qi];
    Mi = [Mi;Z(:)'];
end
K.s = [K.s length(z)+1];
zeroRow = zeros(1,size(allA,2));
allA = [allA Mi;zeroRow 1 zeros(1,K.s(end)^2-1)];
b = zeros(size(allA,1),1);b(end) = 1;
y = sdpvar(length(b),1);
CminusAy = -allA'*y;
start = 1;

% Get the cost, expressed in Z
[Qi,ci,fi,xi,infoi] = quaddecomp(parobj,z);
C = [fi ci'/2;ci/2 Qi];
F = set([]);
for i = 1:length(K.s)
    if i<length(K.s)
        F = F + set(reshape(CminusAy(start:start+K.s(i)^2-1),K.s(i),K.s(i)));
    else
        F = F + set(reshape(C(:) + CminusAy(start:start+K.s(i)^2-1),K.s(i),K.s(i)));
    end
    start = start + K.s(i)^2;
end
obj = -b'*y;

options.sdplr.forcerank = ceil(K.s/2);
options.sdplr.forcerank(end) = 1;
options.sdplr.feastol = 1e-7;



function [exponent_m2,N_unique] = expandmatrix(exponent_m2,N_unique,n);


function [sol,m,Q,residuals,everything] = solvesos_find_blocks(F,obj,options,params,candidateMonomials)

tol = options.sos.numblkdg;
if tol > 1e-2
    disp(' ');
    disp('-> Are you sure you meant to have a tolerance in numblk that big!')
    disp('-> The options numblkdiag controls the tolerance, it is not a 0/1 switch.')
    disp(' ');
end
options.sos.numblkdg = 0;
[sol,m,Q,residuals,everything] = solvesos(F,obj,options,params,candidateMonomials);

% Save old structure to find out when we have stalled
for i = 1:length(Q)
    oldlengths{i} = length(Q{i});
end

go_on = (sol.problem == 0 | sol.problem == 4);
while go_on

    for sosfun = 1:length(Q)
        Qtemp = Q{sosfun};
        keep = diag(Qtemp)>tol;
        Qtemp(:,find(~keep)) = [];
        Qtemp(find(~keep),:) = [];

        m{sosfun} = m{sosfun}(find(keep));

        Qtemp(abs(Qtemp) < tol) = 0;
        [v1,dummy1,r1,dummy3]=dmperm(Qtemp+eye(length(Qtemp)));
        lengths{sosfun} = [];
        n{sosfun} = {};
        for blocks = 1:length(r1)-1
            i1 = r1(blocks);
            i2 = r1(blocks+1)-1;
            if i2>i1
                n{sosfun}{blocks} = m{sosfun}(v1(i1:i2));
            else
                n{sosfun}{blocks} = m{sosfun}(v1(i1));
            end
            lengths{sosfun} =  [lengths{sosfun}  length(n{sosfun}{blocks})];
        end
        lengths{sosfun} = sort(lengths{sosfun});
    end
    go_on = ~isequal(lengths,oldlengths);
    oldlengths = lengths;
    if go_on
        [sol,m,Q,residuals,everything] = solvesos(F,obj,options,params,n);
        go_on = go_on & (sol.problem == 0 | sol.problem == 4);
        if sol.problem == 1
            disp('-> Feasibility was lost during the numerical block-diagonalization.')
            disp('-> The setting sos.numblkdiag is probably too big')
        end
    end
end


function sol =  solveranksos(F,obj,options,ranks,BlockedQ)

Frank = set([]);
for i = 1:length(ranks)
    if ~isinf(ranks(i))
        Frank = Frank + set(rank(BlockedQ{i}{1}) <= ranks(i));
    end
end
% rank adds the pos.def constraints again!!, so we remove them
check = ones(length(F),1);
keep  = ones(length(F),1);
for i = 1:length(BlockedQ)
    for j = 1:length(BlockedQ{i})
        Qij = BlockedQ{i}{j};
        for k = find(check)'
            if isequal(Qij,sdpvar(F(k)))
                keep(k)  = 0;
                check(k) = 0;
            end
        end
    end
end
% Let's hope LMIRANK is there
sol =  solvesdp(F(find(keep)) + Frank,[],sdpsettings(options,'solver',''));