bnb.m

国外专家做的求解LMI鲁棒控制的工具箱,可以相对高效的解决LMI问题

function output = bnb(p)
%BNB          General branch-and-bound scheme for conic programs
%
% BNB applies a branch-and-bound scheme to solve mixed integer
% conic programs (LP, QP, SOCP, SDP) and mixed integer geometric programs.
%
% BNB is never called by the user directly, but is called by
% YALMIP from SOLVESDP, by choosing the solver tag 'bnb' in sdpsettings.
%
% BNB is used if no other mixed integer solver is found, and
% is only useful for very small problems, due to its simple
% and naive implementation.
%
% The behaviour of BNB can be altered using the fields
% in the field 'bnb' in SDPSETTINGS
%
% bnb.branchrule   Deceides on what variable to branch
%                   'max'     : Variable furthest away from being integer
%                   'min'     : Variable closest to be being integer
%                   'first'   : First variable (lowest variable index in YALMIP)
%                   'last'    : Last variable (highest variable index in YALMIP)
%
% bnb.method       Branching strategy
%                   'depth'   : Depth first
%                   'breadth' : Breadth first
%                   'best'    : Expand branch with lowest lower bound
%                   'depthX'  : Depth until integer solution found, then X (e.g 'depthbest')
%
% solver           Solver for the relaxed problems (standard solver tag, see SDPSETTINGS)
%
% inttol           Tolerance for declaring a variable as integer
%
% feastol          Tolerance for declaring constraints as feasible
%
% gaptol           Exit when (upper bound-lower bound)/(1e-3+abs(lower bound)) < gaptol
%
% round            Round variables smaller than bnb.inttol
%
%
% See also SOLVESDP, BINVAR, INTVAR, BINARY, INTEGER, LMI

% Author Johan L鰂berg
% $Id: bnb.m,v 1.40 2005/06/24 16:04:35 joloef Exp $

% ********************************
%% INITIALIZE DIAGNOSTICS IN YALMIP
% ********************************
bnbsolvertime = clock;
showprogress('Branch and bound started',p.options.showprogress);

% ********************************
%% We might have a GP : pre-calc
% ********************************
p.nonlinear = find(~(sum(p.monomtable~=0,2)==1 & sum(p.monomtable,2)==1));

% ********************************
%% Define infinite bounds
% ********************************
if isempty(p.ub)
    p.ub = repmat(inf,length(p.c),1);
end
if isempty(p.lb)
    p.lb = repmat(-inf,length(p.c),1);
end

% ********************************
%% Extract bounds from model
% ********************************
if ~isempty(p.F_struc)
    [lb,ub,used_rows] = findulb(p.F_struc,p.K);
    if ~isempty(used_rows)
        used_rows = used_rows(~any(full(p.F_struc(used_rows,1+p.nonlinear)),2));
        if ~isempty(used_rows)
            lower_defined = find(~isinf(lb));
            if ~isempty(lower_defined)
                p.lb(lower_defined) = max(p.lb(lower_defined),lb(lower_defined));
            end
            upper_defined = find(~isinf(ub));
            if ~isempty(upper_defined)
                p.ub(upper_defined) = min(p.ub(upper_defined),ub(upper_defined));
            end
            p.F_struc(p.K.f+used_rows,:)=[];
            p.K.l = p.K.l - length(used_rows);
        end
    end
end

% ********************************
%% ADD CONSTRAINTS 0<x<1 FOR BINARY
% ********************************
if ~isempty(p.binary_variables)
    p.ub(p.binary_variables) =  min(p.ub(p.binary_variables),1);
    p.lb(p.binary_variables) =  max(p.lb(p.binary_variables),0);
    
    godown = find(p.ub(p.binary_variables) < 1-p.options.bnb.inttol);
    goup   = find(p.lb(p.binary_variables) > p.options.bnb.inttol);
    p.ub(p.binary_variables(godown)) = 0;
    p.lb(p.binary_variables(goup)) = 1;
end

% *******************************
%% PRE-SOLVE (nothing fancy coded)
% *******************************
pss=[];
if isempty(find(isinf([p.ub;p.lb]))) & p.K.l>0 & isempty(p.nonlinear)
    [p.lb,p.ub,redundant,pss] = tightenbounds(-p.F_struc(1+p.K.f:p.K.f+p.K.l,2:end),p.F_struc(1+p.K.f:p.K.f+p.K.l,1),p.lb,p.ub,p.integer_variables,p.binary_variables,ones(length(p.lb),1));
    if length(redundant)>0
        pss.AL0A(redundant,:)=[];
        pss.AG0A(redundant,:)=[];
        p.F_struc(p.K.f+redundant,:)=[];
        p.K.l = p.K.l - length(redundant);  
    end
end

% *******************************
%% PERTURBATION OF LINEAR COST
% *******************************
p.corig = p.c;
if nnz(p.Q)~=0
    g = randn('seed');
    randn('state',1253); %For my testing, I keep this the same...
    % This perturbation has to be better. Crucial for many real LP problems
    p.c = (p.c).*(1+randn(length(p.c),1)*1e-4);
    randn('seed',g);
end

% *******************************
%% We don't need this
% *******************************
p.options.savesolverinput  = 0;
p.options.savesolveroutput = 0;
p.options.saveduals = 0;
p.options.dimacs = 0;

% *******************************
%% Display logics
% 0 : Silent
% 1 : Display branching
% 2 : Display node solver prints
% *******************************
switch max(min(p.options.verbose,3),0)
    case 0
        p.options.bnb.verbose = 0;
    case 1
        p.options.bnb.verbose = 1;
        p.options.verbose = 0;
    case 2
        p.options.bnb.verbose = 2;
        p.options.verbose = 0;
    case 3
        p.options.bnb.verbose = 2;
        p.options.verbose = 1;
    otherwise
        p.options.bnb.verbose = 0;
        p.options.verbose = 0;
end

% *******************************
%% START BRANCHING
% *******************************
[x_min,solved_nodes,lower,upper] = branch_and_bound(p,pss);
%% --

% **********************************
%% CREATE SOLUTION
% **********************************
output.problem = 0;
if isinf(upper)
    output.problem = 1;
end
if isinf(-lower)
    output.problem = 2;
end
if solved_nodes == p.options.bnb.maxiter
    output.problem = 3;
end
output.solved_nodes = solved_nodes;
output.Primal      = x_min;
output.Dual        = [];
output.Slack       = [];
output.infostr      = yalmiperror(output.problem,'BNB');
output.solverinput  = 0;
output.solveroutput =[];
output.solvertime   = etime(clock,bnbsolvertime);
%% --

function [x_min,solved_nodes,lower,upper] = branch_and_bound(p,pss)

% *******************************
%% SET-UP ROOT PROBLEM
% *******************************
p.depth = 0;
p.lower = NaN;
if p.options.usex0
    [x_min,upper] = initializesolution(p);
    p.x0    = x_min;
else
    upper   = inf;
    x_min   = zeros(length(p.c),1);    
    p.x0    = zeros(length(p.c),1);
end

% *******************************
%% Global stuff
% *******************************
lower   = NaN;
stack   = [];

% *******************************
%% Create function handle to solver
% *******************************
lowersolver = p.solver.lower.call;
uppersolver = p.options.bnb.uppersolver;

% *******************************
%% INVARIANT PROBLEM DATA
% *******************************
c = p.corig;
Q = p.Q;
f = p.f;
integer_variables = p.integer_variables;
solved_nodes = 0;

gap = inf;
node = 1;

if p.options.bnb.presolve
    savec = p.c;
    saveQ = p.Q;
    p.Q = p.Q*0;

    n = length(p.c);
    saveBinary = p.binary_variables;
    saveInteger = p.integer_variables;
    p.binary_variables = [];
    p.integer_variables = [];;

    for i = 1:length(c)
        p.c = eyev(n,i);
        output = feval(lowersolver,p);
        if output.problem == 0
            p.lb(i) = max(p.lb(i),output.Primal(i));
        end
        p.c = -eyev(n,i);
        output = feval(lowersolver,p);
        if output.problem == 0
            p.ub(i) = min(p.ub(i),output.Primal(i));
        end
        p.lb(saveBinary) = ceil(p.lb(saveBinary)-1e-3);
        p.ub(saveBinary) = floor(p.ub(saveBinary)+1e-3);
    end
    p.binary_variables = saveBinary;
    p.integer_variables = saveInteger;

    p.Q = saveQ;
    p.c = savec;
end

% ************************************************
% Some hacks to speed up solver calls
% ************************************************
p.getsolvertime = 0;

% *******************************
%% DISPLAY HEADER
% *******************************
originalDiscrete = [p.integer_variables(:);p.binary_variables(:)];
originalBinary   = p.binary_variables(:);

if nnz(Q)==0 & (nnz(p.c-fix(p.c))==0)    
    can_use_ceil_lower = all(ismember(find(p.c),originalDiscrete));
else
    can_use_ceil_lower = 0;
end

if p.options.bnb.verbose

    pc = p.problemclass;
    non_convex_obj = pc.objective.quadratic.nonconvex | pc.objective.polynomial;

    possiblynonconvex = non_convex_obj;
    if ~isequal(p.solver.lower.version,'')
        p.solver.lower.tag = [p.solver.lower.tag '-' p.solver.lower.version];
    end

    disp('* Starting YALMIP integer branch & bound.');
    disp(['* Lower solver   : ' p.solver.lower.tag]);
    disp(['* Upper solver   : ' p.options.bnb.uppersolver]);
    disp(['* Max iterations : ' num2str(p.options.bnb.maxiter)]);

    if possiblynonconvex
        disp(' ');
        disp('Warning : The relaxed problem may be nonconvex. This means ');
        disp('that the branching process not is guaranteed to find a');
        disp('globally optimal solution, since the lower bound can be');
        disp('invalid. Hence, do not trust the bound or the gap...')
    end
end
if p.options.bnb.verbose;            disp(' Node       Upper       Gap(%)      Lower    Open');end;

if nnz(Q)==0 & nnz(c)==1
    p.simplecost = 1;
else
    p.simplecost = 0;
end
poriginal = p;
p.cuts = [];

%% MAIN LOOP
p.options.rounding = [1 1 1 1];

while ~isempty(node) & (solved_nodes < p.options.bnb.maxiter) & (isinf(lower) | gap>p.options.bnb.gaptol)

    % ********************************************
    % BINARY VARIABLES ARE FIXED ALONG THE PROCESS
    % ********************************************
    binary_variables  = p.binary_variables;

    % ********************************************
    % ASSUME THAT WE WON'T FATHOME
    % ********************************************
    keep_digging = 1;
    message = '';
  
    p = Updatecostbound(p,upper);     

    if 0
        if p.K.l>0
            [p.lb,p.ub] = tightenbounds(-p.F_struc(1+p.K.f:p.K.f+p.K.l,2:end),p.F_struc(1+p.K.f:p.K.f+p.K.l,1),p.lb,p.ub,originalDiscrete,originalBinary);
        end
    end

   % [p.lb p.ub]
    % *************************************
    % SOLVE NODE PROBLEM
    % *************************************
    if any(p.ub<p.lb)
        x = zeros(length(p.c),1);
        output.Primal = x;
        output.problem=1;
    else
        relaxed_p = p;
        relaxed_p.integer_variables = [];
        relaxed_p.binary_variables = [];
        output = solvelower(lowersolver,relaxed_p,upper,lower);

        x = setnonlinearvariables(p,output.Primal);
        
        % **************************************
        % Hmm, don't remember why this fix...
        % **************************************
        if (p.K.l>0) & any(p.F_struc(p.K.f+1:p.K.f+p.K.l,:)*[1;x]<-1e-5)
            output.problem = 1;
        end                    
    end

    solved_nodes = solved_nodes+1;

    % **************************************
    % THIS WILL BE INTIAL GUESS FOR CHILDREN
    % **************************************
    p.x0 = x;

    % *************************************
    % ROOT NODE HEURISTICS (NOTHING CODED)
    % *************************************
    if output.problem==0 | output.problem==3 | output.problem==4

        cost = f+c'*x+x'*Q*x;        
        
        if isnan(lower)
            lower = cost;   
        end

        [upper1,x_min1] = feval(uppersolver,poriginal,output);
        if upper1 < upper            
            x_min = x_min1;
            upper = upper1;
            [stack,stacklower] = prune(stack,upper,p.options,solved_nodes,p);            
            lower = min(lower,stacklower);
        end
    end
    p = adaptivestrategy(p,upper,solved_nodes);

    % *************************************
    % ANY INTEGERS? ROUND?
    % *************************************
    non_integer_binary = abs(x(binary_variables)-round(x(binary_variables)))>p.options.bnb.inttol;
    non_integer_integer = abs(x(integer_variables)-round(x(integer_variables)))>p.options.bnb.inttol;
    if p.options.bnb.round
        x(binary_variables(~non_integer_binary))   = round(x(binary_variables(~non_integer_binary)));
        x(integer_variables(~non_integer_integer)) = round(x(integer_variables(~non_integer_integer)));
    end
    non_integer_binary = find(non_integer_binary);
    non_integer_integer = find(non_integer_integer);

    % *************************************
    % CHECK FATHOMING POSSIBILITIES
    % *************************************
    feasible = 1;
    switch output.problem
        case 0          
             if can_use_ceil_lower 
                lower = ceil(lower);
             end
        case {1,12}            
            keep_digging = 0;
            cost = inf;
            feasible = 0;
        case 2            
            cost = -inf;
        otherwise            
            % This part has to be much more robust
            cost = f+c'*x+x'*Q*x;
    end

    % **************************************
    % YAHOO! INTEGER SOLUTION FOUND
    % **************************************
    if isempty(non_integer_binary) & isempty(non_integer_integer)
        if (cost<upper) & feasible           
            x_min = x;
            upper = cost;
            [stack,lower] = prune(stack,upper,p.options,solved_nodes,p);
        end
        p = adaptivestrategy(p,upper,solved_nodes);
        keep_digging = 0;
    end
    if cost>upper*(1-1e-6)        
        keep_digging = 0;
    end
    % **********************************
    % CONTINUE SPLITTING?
    % **********************************
    if keep_digging & (cost<upper)

        % **********************************
        % BRANCH VARIABLE
        % **********************************
        [index,whatsplit] = branchvariable(x,integer_variables,binary_variables,p.options,x_min);

        % **********************************
        % CREATE NEW PROBLEMS
        % **********************************
        switch whatsplit
            case 'binary'
                [p0,p1,index] = binarysplit(p,x,index,cost);
            case 'integer'
                [p0,p1] = integersplit(p,x,index,cost,x_min);
            otherwise
        end

        
        % **********************************
        % Incrementally improve bounds
        % **********************************
      %  A=[p.c';-p.F_struc(1+p.K.f:p.K.f+p.K.l,2:end)];
      %  b=[upper;p.F_struc(1+p.K.f:p.K.f+p.K.l,1)];
      %  [p1.lb,p1.ub]=tightenbounds(A,b,p1.lb,p1.ub,originalDiscrete,originalBinary,eyev(length(p.lb),index));
      %  [p0.lb,p0.ub]=tightenbounds(A,b,p0.lb,p0.ub,originalDiscrete,originalBinary,eyev(length(p.lb),index));
        % **********************************
        % Only save varying data in the tree
        % **********************************

        % **********************************
        % Only save varying data in the tree
        % **********************************
        node1.lb = p1.lb;
        node1.ub = p1.ub;
        node1.depth = p1.depth;
        node1.lower = p1.lower;
        node1.x0 = p1.x0;
        node1.binary_variables = p1.binary_variables;
        node0.lb = p0.lb;