robustify_sdp_conic.m

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

function F = robustify_sdp_conic(F_xw,Zmodel,x,w)

if length(F_xw) == 0
    F = [];
    return;
else
    Ftemp = [];
    Fvars = getvariables(F_xw);
    wvars = getvariables(w);
    [mt,vt] = yalmip('monomtable');
    if any(sum(mt(Fvars,wvars),2)>1)
        for i = 1:length(F_xw)
            Fi = sdpvar(F_xw(i));
        %    if sum(degree(Fi(:),w)) > 1
            if any(sum(mt(getvariables(Fi),wvars),2)>1)                
                [BilinearizeringConstraints,failure] = deriveBilinearizing(Fi,w);
                if failure
                    error('Cannot get rid of nonlinear uncertainty in uncertain SDP')
                else
                    Ftemp = Ftemp + BilinearizeringConstraints;
                end
            end
        end
        F_xw = F_xw + Ftemp;
    end

    if any(Zmodel.K.q) | any(Zmodel.K.s)
        error('Only polytope uncertainty supported for uncertain SDPs');
    else
        % FIX : Assumes all uncertainty in all constraints
        K = Zmodel.K;
        A = -Zmodel.F_struc((1+K.f):(K.f + K.l),2:end);
        b =  Zmodel.F_struc((1+K.f):(K.f + K.l),1);

        try
            % some preprocessing to extract bounds from equality
            % constraints in order to make the uncertainty polytope bounded
            [lo,up] = findulb(Zmodel.F_struc(1:K.f + K.l,:),K);
            Zmodel.lb = lo;Zmodel.ub = up;
            Zmodel = presolve_bounds_from_equalities(Zmodel);
            up = Zmodel.ub;
            lo = Zmodel.lb;
            upfi = find(~isinf(up));
            lofi = find(~isinf(lo));
            aux = Zmodel;
            aux.F_struc = [aux.F_struc;-lo(lofi) sparse(1:length(lofi),lofi,1,length(lofi),size(A,2))];
            aux.F_struc = [aux.F_struc;up(upfi) -sparse(1:length(upfi),upfi,1,length(upfi),size(A,2))] ;                       
            aux.K.l = aux.K.l + length(lofi) + length(upfi);                        
            K = aux.K;
            A = -aux.F_struc((1+K.f):(K.f + K.l),2:end);
            b =  aux.F_struc((1+K.f):(K.f + K.l),1);
        
            vertices = extreme(polytope(A,b))';
        catch
            disp('You probably need to install MPT (needed for vertex enumeration)')
            disp('http://control.ee.ethz.ch/~joloef/wiki/pmwiki.php?n=Solvers.MPT')
            disp('Alternatively, you need to add bounds on your uncertainty.')
            error('MPT missing');
        end
        if K.f > 0
            Aeq = -Zmodel.F_struc(1:K.f,2:end);
            beq =  Zmodel.F_struc(1:K.f,1);
            feasible = sum(abs(Aeq*vertices - repmat(beq,1,size(vertices,2))),1) < 1e-6;
            vertices = vertices(:,feasible);
            if isempty(feasible)
                error('The uncertainty space is infeasible.')                
            end
        end

        F_xw_lp = F_xw(find(is(F_xw,'elementwise')));
        F_xw_socp_sdp = F_xw -  F_xw_lp;
        F = set([]);
        
        % Do all LPs in one shot
        if length(F_xw_lp)>0
            rLP = [];
            %rLP = vec_replace(sdpvar(F_xw_lp),w,vertices);
            for i = 1:size(vertices,2)
                rLP = [rLP;replace(sdpvar(F_xw_lp),w,vertices(:,i),0)];
           %     F = F + set(replace(sdpvar(F_xw_lp),w,vertices(:,i),0));   
            end
            % Easily genrates redundnat constraints, so use simple check
            %rLP = sdpvar(F); 
            [aux,index] = uniquesafe(getbase(rLP),'rows');
            F = set(rLP(index) >= 0);
        end
        
        % Remaining conic stuff
        for j = 1:length(F_xw_socp_sdp)
            for i = 1:size(vertices,2)
                F = F + set(replace(F_xw_socp_sdp(j),w,vertices(:,i),0));
            end
        end
    end
end