cutsdp.m

matlab波形优化算法经常要用到的matlab toolbox工具箱:yalmip

goon = 1;
rr = p_lp.integer_variables;
rb = p_lp.binary_variables;
only_solvelp = 0;
pplot = 0;

% *************************************************************************
% Crude lower bound
% FIX for quadratic case
% *************************************************************************
lower = 0;
if nnz(p.Q) == 0
    for i = 1:length(p.c)
        if p.c(i)>0
            if isinf(p.lb(i))
                lower = -inf;
                break
            else
                lower = lower + p.c(i)*p.lb(i);
            end
        elseif p.c(i)<0
            if isinf(p.ub(i))
                lower = -inf;
                break
            else
                lower = lower + p.c(i)*p.ub(i);
            end
        end
    end
end
%lower = sum(sign(p.c).*(p.lb));
if isinf(lower) | nnz(p.Q)~=0
    lower = -1e6;
end

% *************************************************************************
% Experimental stuff for variable fixing
% *************************************************************************
if p.options.cutsdp.nodefix & (p.K.s(1)>0)
    top=1+p.K.f+p.K.l+sum(p.K.q);
    for i=1:length(p.K.s)
        n=p.K.s(i);
        for j=1:size(p.F_struc,2)-1;
            X=full(reshape(p.F_struc(top:top+n^2-1,j+1),p.K.s(i),p.K.s(i)));
            X=(X+X')/2;
            v=real(eig(X+sqrt(eps)*eye(length(X))));
            if all(v>=0)
                sdpmonotinicity(i,j)=-1;
            elseif all(v<=0)
                sdpmonotinicity(i,j)=1;
            else
                sdpmonotinicity(i,j)=nan;
            end
        end
        top=top+n^2;
    end
else
    sdpmonotinicity=[];
end

hist_infeasibility = [];
while goon

    % Add lower bound
    if ~isinf(lower)
        p_lp.F_struc = [p_lp.F_struc;-lower p_lp.c'];
        p_lp.K.l = p_lp.K.l + 1;
    end

    if p.options.cutsdp.nodetight
        % Extract LP part Ax<=b
        A = -p_lp.F_struc(p_lp.K.f + (1:p_lp.K.l),2:end);
        b = p_lp.F_struc(p_lp.K.f + (1:p_lp.K.l),1);
        c = p_lp.c;
        % Tighten bounds and find redundant constraints
        [p_lp.lb,p_lp.ub,redundant,pss] = milppresolve(A,b,p_lp.lb,p_lp.ub,p.integer_variables,p.binary_variables,ones(length(p.lb),1));
        A(redundant,:) = [];
        b(redundant) = [];
        p_lp.F_struc(p_lp.K.f+redundant,:) = [];
        p_lp.K.l = p_lp.K.l-length(redundant);
    end

    if p.options.cutsdp.nodefix
        % Try to find variables to fix w.l.o.g
        [fix_up,fix_down] = presolve_fixvariables(A,b,c,p_lp.lb,p_lp.ub,sdpmonotinicity);
        p_lp.lb(fix_up)   = p_lp.ub(fix_up);
        p_lp.ub(fix_down) = p_lp.lb(fix_down);
        while ~(isempty(fix_up) & isempty(fix_down))
            [p_lp.lb,p_lp.ub,redundant,pss] = milppresolve(A,b,p_lp.lb,p_lp.ub,p.integer_variables,p.binary_variables,ones(length(p.lb),1));
            A(redundant,:) = [];
            b(redundant) = [];
            p_lp.F_struc(p_lp.K.f+redundant,:) = [];
            p_lp.K.l = p_lp.K.l-length(redundant);
            fix_up = [];
            fix_down = [];
            % Try to find variables to fix w.l.o.g
            [fix_up,fix_down] = presolve_fixvariables(A,b,c,p_lp.lb,p_lp.ub,sdpmonotinicity);
            p_lp.lb(fix_up)   = p_lp.ub(fix_up);
            p_lp.ub(fix_down) = p_lp.lb(fix_down);
        end
    end


    output = feval(cutsolver,p_lp);

    % Remove lower bound (avoid accumulating them)
    if ~isinf(lower)
        p_lp.K.l = p_lp.K.l - 1;
        p_lp.F_struc = p_lp.F_struc(1:end-1,:);
    end

    if output.problem == 1 | output.problem == 12
        % LP relaxation was infeasible, hence problem is infeasible
        feasible = 0;
        lower = inf;
        goon = 0;
        x = zeros(length(p.c),1);
        lower = inf;
    else
        % Relaxed solution
        x = output.Primal;
        lower = p.f+p.c'*x+x'*p.Q*x;

        infeasibility = 0;

        [p_lp,infeasibility] = add_socp_cut(p,p_lp,x,infeasibility);
        [p_lp,infeasibility] =  add_sdp_cut(p,p_lp,x,infeasibility);

        %    res = find(p_lp.F_struc*[1;x] > mean(abs(p_lp.F_struc*[1;x])));
        %    p_lp.F_struc(res,:) = [];
        %    p_lp.K.l = p_lp.K.l-length(res);
        goon = infeasibility <= p.options.cutsdp.feastol;
        goon = goon & feasible;
        goon = goon & (solved_nodes < p.options.cutsdp.maxiter-1);
    end

    solved_nodes = solved_nodes + 1;
    if p.options.cutsdp.verbose;fprintf(' %4.0f :      %12.3E      %12.3E      %2.0f\n',solved_nodes,infeasibility,lower,p_lp.K.l-p.K.l);end
end

D_struc = [];



function [p_lp,infeasibility] = add_sdp_cut(p,p_lp,x,infeasibility);

if p.K.s(1)>0
    % Add cuts
    top = p.K.f+p.K.l+p.K.q+1;
    for i = 1:1:length(p.K.s)
        n = p.K.s(i);
        X = p.F_struc(top:top+n^2-1,:)*[1;x];
        X = reshape(X,n,n);
        Y = randn(n,n);
        newcuts = 1;
        newF = zeros(n,size(p.F_struc,2));
        [d,v] = eig(X);
        infeasibility = min(infeasibility,min(diag(v)));
        dummy=[];
        newF = [];
        if infeasibility<0
            [ii,jj] = sort(diag(v));
            for m = jj(1:min(length(jj),p.options.cutsdp.cutlimit))'%find(diag(v<0))%1:1%length(v)
                if v(m,m)<0
                    bA =  d(:,m)'*(kron(d(:,m),speye(n)).'*p.F_struc(top:top+n^2-1,:));
                    b = bA(:,1);
                    A = -bA(:,2:end);
                    newF = real([newF;[b -A]]);
                    %f = b-floor(b);
                    %fj = A - floor(A);
                    %A = floor(A) + max((fj-f),0)./(1-f);
                    %b = floor(b);
                    %newF = [newF;b -A];
                    newcuts = newcuts + 1;
                end
            end
        end
        newF(abs(newF)<1e-12) = 0;
        keep=find(any(newF(:,2:end),2));
        newF = newF(keep,:);
        if size(newF,1)>0
            p_lp.F_struc = [p_lp.F_struc;newF];
            p_lp.K.l = p_lp.K.l + size(newF,1);
            [i,j] = sort(p_lp.F_struc*[1;x]);
        end
        top = top+n^2;
    end
end





function [p_lp,infeasibility] = add_socp_cut(p,p_lp,x,infeasibility);

if p.K.q(1)>0
    % Add cuts
    top = p.K.f+p.K.l+1;
    for i = 1:1:length(p.K.s)
        n = p.K.q(i);
        X = p.F_struc(top:top+n-1,:)*[1;x];
        X = [X(1) X(2:end)';X(2:end) eye(n-1)*X(1)];
        Y = randn(n,n);
        newcuts = 1;
        newF = zeros(n,size(p.F_struc,2));
        [d,v] = eig(X);
        infeasibility = min(infeasibility,min(diag(v)));
        dummy=[];
        newF = [];
        if infeasibility<0
            [ii,jj] = sort(diag(v));
            for m = jj(1:min(length(jj),p.options.cutsdp.cutlimit))'%find(diag(v<0))%1:1%length(v)
                if v(m,m)<0
                    v1 = d(1,m);v2 = d(2:end,m);
                    newF = [newF;p.F_struc(top,:) + 2*v1*v2'*p.F_struc(top+1:top+n-1,:)];
                    newcuts = newcuts + 1;
                end
            end
        end
        newF(abs(newF)<1e-12) = 0;
        keep=find(any(newF(:,2:end),2));
        newF = newF(keep,:);
        if size(newF,1)>0
            p_lp.F_struc = [p_lp.F_struc;newF];
            p_lp.K.l = p_lp.K.l + size(newF,1);
            [i,j] = sort(p_lp.F_struc*[1;x]);
        end
        top = top+n;
    end
end