bmibnb.m

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

        p.complementary = [p.complementary ;x y];
    end
end

% Some code to extract bounds in complementary conditions, when the
% complementary variables hasn't been defined
for i = 1:size(p.complementary,1)
    %    if 1%(isinf(p.lb(p.complementary(i,1))) |  isinf(p.ub(p.complementary(i,1)))) &
    if ismember(p.complementary(i,1),p.branch_variables)
        %x*y == 0 where x is unbounded. Hence x is either 0, or bounded by
        %the value obtained when y is 0
        pp=p;
        pp.ub(p.complementary(i,2))=0;
        pp = presolveloop(pp,upper);
        if pp.feasible
        p.ub(p.complementary(i,1)) = min(pp.ub(p.complementary(i,1)),p.ub(p.complementary(i,1)));
        p.lb(p.complementary(i,1)) = max(pp.lb(p.complementary(i,1)),p.lb(p.complementary(i,1)));
        else
            p.ub(p.complementary(i,1))=0;
            p.lb(p.complementary(i,1))=0;   
       
        end
    end
%     %    if 1%(isinf(p.lb(p.complementary(i,2))) |  isinf(p.ub(p.complementary(i,2)))) &
    if ismember(p.complementary(i,2),p.branch_variables)
        %x*y == 0 where y is unbounded. Hence y is either 0, or bounded by
        %the bounds obtained when x is 0
        pp=p;
        pp.ub(p.complementary(i,1))=0;
        pp = presolveloop(pp,upper);
        if pp.feasible
        p.ub(p.complementary(i,2)) = min(pp.ub(p.complementary(i,2)),p.ub(p.complementary(i,2)));
        p.lb(p.complementary(i,2)) = max(pp.lb(p.complementary(i,2)),p.lb(p.complementary(i,2)));
        else
            % This case is infeasible, hence the other term has to be zero
            p.ub(p.complementary(i,2))=0;
            p.lb(p.complementary(i,2))=0;       
        end
    end
end

if p.feasible
    
    % *******************************
    % Bounded domain?
    % *******************************
    if ~isempty(p.branch_variables)
        if any(isinf(p.lb(p.branch_variables))) | any(isinf(p.ub(p.branch_variables)))
            output = yalmip_default_output;
            output.Primal  = x_min;                             
            output.problem = -6;
            output.infostr = yalmiperror(-6);
            output.solved_nodes = 0;           
            return
        end
    end
    
    % *******************************
    % Save time & memory
    % *******************************
    p.options.savesolverinput  = 0;
    p.options.savesolveroutput = 0;
    p.options.saveduals = 0;
    p.options.dimacs = 0;    
    
    % *******************************
    % RUN BILINEAR BRANCH & BOUND
    % *******************************
    [x_min,solved_nodes,lower,upper] = branch_and_bound(p,x_min,upper);
    
    % ********************************
    % CREATE SOLUTION AND DIAGNOSTICS
    % ********************************
    problem = 0;
    if isinf(upper)
        problem = 1;
    end
    if isinf(lower) & (lower<0)
        problem = 2;
    end
    if isinf(lower) & (lower>0)
        problem = 1;
    end    
    if solved_nodes == p.options.bmibnb.maxiter
        problem = 3;
    end
else
    problem = 1;
    x_min = repmat(nan,length(p.c),1);
    solved_nodes = 0;
    lower = inf;
end

x_min = dediagonalize(p,x_min);

output = yalmip_default_output;
output.problem = problem;
output.solved_nodes = solved_nodes;
output.Primal        = zeros(length(p.kept),1);
output.Primal(p.kept(1:n_in))= x_min(1:n_in);
output.infostr      = yalmiperror(output.problem,'BNB');
output.solvertime   = etime(clock,bnbsolvertime);
output.lower = lower;

function pnew = diagonalize_quadratic_program(p);
pnew = p;
pnew.V = [];
pnew.diagonalized = 0;

% Any polynomial terms or simple linear by some reason
if any(p.variabletype > 2) & ~all(p.variabletype == 0) | p.options.bmibnb.diagonalize==0
    return
end

if ~isempty(p.evalVariables)
    return
end

if ~isempty(p.binary_variables) | ~isempty(p.integer_variables)
    return
end


nonlinear = find(p.variabletype > 0);
linear = find(p.variabletype == 0);
if ~isempty(p.F_struc)
    % Nonlinear terms in constraints
    if nnz(p.F_struc(:,1 + nonlinear))> 0
        return
    end
end

if ~isempty(p.lb)
    if ~all(isinf(p.lb(nonlinear)))
        return
    end
end
if ~isempty(p.ub)
    if ~all(isinf(p.ub(nonlinear)))
        return
    end
end

% Find quadratic and linear terms
used_in_c = find(p.c);
quadraticterms = used_in_c(find(ismember(used_in_c,nonlinear)));
Q = zeros(length(p.c),length(p.c));
if ~isempty(quadraticterms)
    usedinquadratic = zeros(1,length(p.c));
    for i = 1:length(quadraticterms)
        Qij = p.c(quadraticterms(i));
        power_index = find(p.monomtable(quadraticterms(i),:));
        if length(power_index) == 1
            Q(power_index,power_index) = Qij;
        else
            Q(power_index(1),power_index(2)) = Qij/2;
            Q(power_index(2),power_index(1)) = Qij/2;
        end       
    end
end
Qlin = Q(linear,linear);
clin = p.c(linear);

% Decompose Q
[V,D] = eig(full(Qlin));
V = real(V);
D = real(D);
V(abs(V)<1e-11) = 0;
D(abs(D)<1e-11) = 0;
lb = p.lb(linear);
ub = p.ub(linear);
Z = V';
for i = 1:length(lb)
    z = Z(i,:);
    newub(i,1) = z(z>0)*ub(z>0) + z(z<0)*lb(z<0);
    newlb(i,1) = z(z>0)*lb(z>0) + z(z<0)*ub(z<0);
end


% Create new problem
clin = V'*clin;
%A = A*V;

n = length(linear);
pnew.original_linear = linear;
pnew.original_n = length(p.c);
pnew.V = V;
pnew.c = [clin;diag(D)];
pnew.Q = spalloc(2*n,2*n,0);

% find constraint polytope
if size(p.F_struc,1)>0
    A = -p.F_struc(:,1 + linear);
    b = p.F_struc(:,1);
    pnew.F_struc = [b -A*V zeros(length(b),n)];
end


pnew.variabletype = [zeros(1,n) ones(1,n)*2];
pnew.monomtable = [eye(n);2*eye(n)];
pnew.monomtable(2*n,2*n) = 0;
pnew.lb =-inf(2*n,1);
pnew.ub = inf(2*n,1);
pnew.lb(1:n) = newlb;
pnew.ub(1:n) = newub;
if length(pnew.x0)>0
    pnew.x0 = V*p.x0(1:2);
end
pnew.diagonalized = 1;

function x = dediagonalize(p,x);
if isempty(p.V)
    return
end
y = p.V*x(1:length(x)/2);
x = zeros(p.original_n,1);
x(p.original_linear) = y;


function p = presolveloop(p,upper)
i = 0;
goon = 1;
while goon
    start = [p.lb;p.ub];
    i = i+1;
    p = updateboundsfromupper(p,upper);
    p = propagatecomplementary(p);
    p = updatemonomialbounds(p);
    p = presolve_bounds_from_equalities(p);
    p = updatemonomialbounds(p);
    p = updatemonomialbounds(p);
    p = update_eval_bounds(p);
    i = i+1;
    goon = (norm(start-[p.lb;p.ub],'inf') > 1e-2) & i < 150;
    start = [p.lb;p.ub];
end