bmibnb.m

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

function output = bmibnb(p)
%BMIBNB          Branch-and-bound scheme for bilinear programs
%
% BMIBNB is never called by the user directly, but is called by
% YALMIP from SOLVESDP, by choosing the solver tag 'bmibnb' in sdpsettings
%
% The behaviour of BMIBNB can be altered using the fields
% in the field 'bmibnb' in SDPSETTINGS
%s
% WARNING: THIS IS EXPERIMENTAL CODE
%
% bmibnb.lowersolver    - Solver for lower bound [standard solver tag ('')]
% bmibnb.uppersolver    - Solver for upper bound [standard solver tag ('')]
% bmibnb.lpsolver       - Solver for LP bound tightening [standard solver tag ('')]
% bmibnb.branchmethod   - Branch strategy ['maxvol' | 'best' ('best')]
% bmibnb.branchrule     - Branch position ['omega' | 'bisect' ('omega')]
% bmibnb.lpreduce       - Improve variable bounds using LP [ real [0,1] (0 means no reduction, 1 means all variables)
% bmibnb.lowrank        - partition variables into two disjoint sets and branch on smallest [ 0|1 (0)]
% bmibnb.target         - Exit if upper found<target [double (-inf)]
% bmibnb.roottight      - Improve variable bounds in root using full problem [ 0|1 (1)]
% bmibnb.vartol         - Cut tree when x_U-x_L < vartol on all branching variables
% bmibnb.relgaptol      - Tolerance on relative objective error (UPPER-LOWER)/(1+|UPPER|) [real (0.01)]
% bmibnb.absgaptol      - Tolerance on objective error (UPPER-LOWER) [real (0.01)]
% bmibnb.pdtol          - A number is declared non-negative if larger than...[ double (-1e-6)]
% bmibnb.eqtol          - A number is declared zero if abs(x) smaller than...[ double (1e-6)]
% bmibnb.maxiter        - Maximum number of solved nodes [int (100)]
% bmibnb.maxtime        - Maximum CPU time (sec.) [int (3600)]

% Author Johan L鰂berg
% $Id: bmibnb.m,v 1.48 2005/06/20 15:26:48 joloef Exp $

% *************************************************************************
% INITIALIZE DIAGNOSTICS ETC
% *************************************************************************
bnbsolvertime = clock;
showprogress('Branch and bound started',p.options.showprogress);
switch max(min(p.options.verbose,3),0)
case 0
    p.options.bmibnb.verbose = 0;
case 1
    p.options.bmibnb.verbose = 1;
    p.options.verbose = 0;
case 2
    p.options.bmibnb.verbose = 2;
    p.options.verbose = 0;
case 3
    p.options.bmibnb.verbose = 2;
    p.options.verbose = 1;
otherwise
    p.options.bmibnb.verbose = 0;
    p.options.verbose = 0;
end

% *******************************
% Assume feasible
% *******************************
p.feasible = 1;

% *******************************
% Pre-calc linear/bilinear/nonlinear variables
% (assumes polynomial program)
% *******************************
[p.linears,p.bilinears,p.nonlinears] = compile_nonlinear_table(p);

% *******************************
% Select branch variables
% *******************************
p.branch_variables = decide_branch_variables(p);

% ********************************
% Extract bounds from model
% ********************************
p = preprocess_bounds(p);

% ********************************
% Is the initial value feasible
% ********************************
[p,x_min,upper] = initializesolution(p);

% ********************************
% Simple pre-solve (empty rows etc)
% Bound strenghtening
% The loop is needed for variables
% defined as w = x*y, x = t*u,y=..
% ********************************
start = [p.lb;p.ub]+1;i = 0;;goon = 1;
while goon
    start = [p.lb;p.ub];
    i = i+1;
    p = simple_presolve(p);
    p = updatenonlinearbounds(p);
    p = propagateequalities(p);
    p = updatenonlinearbounds(p);   
    p = updatenonlinearbounds(p);    
    i = i+1;
    goon = (norm(start-[p.lb;p.ub],'inf') > 1e-1) & i < 25;
    start = [p.lb;p.ub];
end

% ********************************
% Default root bound improvement
% ********************************
p = root_node_tighten(p,upper);
p = updatenonlinearbounds(p);
p = propagateequalities(p);
p = updatenonlinearbounds(p);

if p.feasible
    
    % *******************************
    % Bounded domain?
    % *******************************
    if ~isempty(p.bilinears)
        involved = unique(p.bilinears(:,2:3));
        if any(isinf(p.lb(involved))) | any(isinf(p.ub(involved)))
            output.Primal       = x_min;      
            output.Dual      = [];
            output.Slack     = [];            
            output.problem = -6;
            output.infostr = yalmiperror(-6);
            output.solved_nodes = 0;
            output.solverinput  = 0;
            output.solveroutput = [];
            output.solvertime   = 0;
            return
        end
    end
    
    % *******************************
    % Save time & memory
    % *******************************
    p.options.savesolverinput  = 0;
    p.options.savesolveroutput = 0;
    p.options.saveduals = 0;
    p.options.dimacs = 0;
    
    % *******************************
    % RUN BRANCH & BOUND
    % *******************************
    [x_min,solved_nodes,lower,upper] = branch_and_bound(p,x_min,upper);
    
    % **********************************
    % CREATE SOLUTION
    % **********************************
    output.problem = 0;
    if isinf(upper)
        output.problem = 1;
    end
    if isinf(lower) & (lower<0)
        output.problem = 2;
    end
    if isinf(lower) & (lower>0)
        output.problem = 1;
    end
    
    if solved_nodes == p.options.bnb.maxiter
        output.problem = 3;
    end
else
    output.problem = 1;
    x_min = repmat(nan,length(p.c),1);
    solved_nodes = 0;
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,x_min,upper)

% ***************************************
% Create function handles to solvers
% ***************************************
try
    lowersolver = eval(['@' p.solver.lowersolver.call]); % Solver for relaxed problem
    uppersolver = eval(['@' p.solver.uppersolver.call]); % Local nonlinear solver
    lpsolver    = eval(['@' p.solver.lpsolver.call]);    % LP solver
catch
    disp(' ');
    disp('The internal branch & bound solver requires MATLAB 6')
    disp('I am too lazy too do the changes to make it compatible')
    disp('with MATLAB 5. If you really need it, contact me...');
    disp(' ');
    error(lasterr);
end

% ************************************************
% GLOBAL PROBLEM DATA
% ************************************************
c       = p.c;
Q       = p.Q;
f       = p.f;
K       = p.K;
options = p.options;

% ************************************************
% DEFINE ORIGINAL PROBLEM (used in local solver)
% ************************************************
p_upper = cleanuppermodel(p);

% ************************************************
% Active constraints in main model
% 0  : Inactive constraint
% 1  : Active constraint
% inf: Removed constraint
% ************************************************
p.constraintState = ones(p.K.l,1);
p.constraintState(p.KCut.l,1) = 0;

% ************************************************
% LPs ARE USED IN  BOX-REDUCTION
% ************************************************
p.lpcuts = p.F_struc(1+p.K.f:1:p.K.l+p.K.f,:);
p.cutState = ones(p.K.l,1);
p.cutState(p.KCut.l,1) = 0; % Don't use to begin with

% ************************************************
% Remove sdp cutting planes from local problem
% ************************************************
if length(p_upper.KCut.s)>0
    starts = p_upper.K.f+p_upper.K.l + [1 1+cumsum((p_upper.K.s).^2)];
    remove_these = [];
    for i = 1:length(p_upper.KCut.s)
        j = p_upper.KCut.s(i);
        remove_these = [remove_these;(starts(j):starts(j+1)-1)'];
    end
    p_upper.F_struc(remove_these,:)=[];
    p_upper.K.s(p_upper.KCut.s) = [];
end

% ************************************************
% Some messy code to reduce the size of the
% convex envelope model, based on the objective function
% ************************************************
if 0
    s = find(p.c);
    if all(ismember(s,p.bilinears(:,1)))
        p.r = p.bilinears(s,1);
        p.s = p.c(s);               
    end
else
    p.r = [];
    p.s = [];
end

% ************************************************
% INITIALITAZION
% ************************************************
p.depth = 0;        % depth in search tree
p.dpos  = 0;         % used for debugging
p.lower = NaN;
lower   = NaN;
gap     = inf;
stack   = [];
solved_nodes = 0;
numGlobalSolutions = 0;


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

% ************************************************
% Branch & bound loop
% ************************************************
if options.bmibnb.verbose>0
    disp('* Starting YALMIP bilinear branch & bound.');
    disp(['* Upper solver   : ' p.solver.uppersolver.tag]);    
    disp(['* Lower solver   : ' p.solver.lowersolver.tag]);
    if p.options.bmibnb.lpreduce
        disp(['* LP solver      : ' p.solver.lpsolver.tag]);
    end
    disp(' Node       Upper      Gap(%)       Lower    Open');
end

t_start = cputime;
go_on  = 1;

while go_on
    
    % ********************************************
    % ASSUME THAT WE WON'T FATHOME
    % ********************************************
    keep_digging = 1;
    
    % ********************************************
    % Reduce size of current box (bound tightening)
    % ********************************************
    p = propagatequadratics(p,upper,lower);  
    [p,feasible,vol_reduction] =  domain_reduction(p,upper,lower,lpsolver);
    
    % ********************************************
    % Detect redundant constraints
    % ********************************************
    p = remove_redundant(p);
        
    % ********************************************
    % SOLVE LOWER AND UPPER
    % ********************************************    
    if feasible      
        output = solvelower(p,options,lowersolver);
        
        info_text = '';
        switch output.problem
        case 1 % Infeasible              
            info_text = 'Infeasible';
            keep_digging = 0;
            cost = inf;
            feasible = 0;
            
        case 2 % Unbounded (should not happen!)
            cost = -inf;
            x = output.Primal;
            
        case {0,3,4} % No problems (disregard numerical problems)
            
            x = output.Primal;
            
            cost = f+c'*x+x'*Q*x;
            
           
            % UPDATE THE LOWER BOUND
            if isnan(lower)
                lower = cost;
            end
            if ~isempty(stack)
                lower = min(cost,min([stack.lower]));
            else
                lower = min(lower,cost);
            end
            
            if cost<upper-1e-5
                
                z = evaluate_nonlinear(p,x);
                
                % Manage cuts etc
                p = addsdpcut(p,z);
                p = addlpcuts(p,x);
                
                if numGlobalSolutions < p.options.bmibnb.numglobal
                    [upper,x_min,cost,info_text,numGlobalSolutions] = heuristics_from_relaxed(p_upper,x,upper,x_min,cost,numGlobalSolutions);
                    [upper,x_min,info_text,numGlobalSolutions] = solve_upper_in_node(p,p_upper,x,upper,x_min,uppersolver,info_text,numGlobalSolutions);
                end
            else
                keep_digging = 0;
                info_text = 'Poor bound';
            end
        otherwise
            cost = -inf;
            x = (p.lb+p.ub)/2;
        end