callmpcvx.m

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

function output = mpcvx(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
%
% 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: callmpcvx.m,v 1.2 2005/05/07 13:53:20 joloef Exp $

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

% *******************************
% Display-logics
% *******************************
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

% *******************************
% Actual linear variables
% *******************************
p.linears = find(sum(p.monomtable,2)==1);

% *******************************
% PRE-CALCULATE INDICIES
% FOR BILNEAR VARIABLES
% *******************************
nonlinear = find(~(sum(p.monomtable~=0,2)==1 & sum(p.monomtable,2)==1));
nonlins   = [];
for i = 1:length(nonlinear)
    z = find(p.monomtable(nonlinear(i),:));
    if length(z)==1
        nonlins = [nonlins;nonlinear(i) z z];
    else
        nonlins = [nonlins;nonlinear(i) z(1) z(2)];
    end
end
p.nonlins = nonlins;

p.options.saveduals = 0;

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

% ********************************
% Extract bounds from model
% ********************************
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
if ~isempty(p.F_struc)
    [lb,ub,used_rows] = findulb(p.F_struc,p.K);
    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   
        % Remove linear bounds
        used_rows = used_rows(find(~any(p.F_struc(p.K.f+used_rows,1+nonlinear),2)));
        not_used_rows = setdiff(1:p.K.l,used_rows);
        for i = 1:length(p.KCut.l)
            p.KCut.l(i) = find(not_used_rows==p.KCut.l(i) );
        end
        if ~isempty(used_rows)
            p.F_struc(p.K.f+used_rows,:)=[];
            p.K.l = p.K.l - length(used_rows);
        end       
    end
end
p = clean_bounds(p);
p = updatenonlinearbounds(p);
feasible = all(p.lb<=p.ub);

% ********************************
% Remove empty linear rows
% ********************************
if p.K.l > 0
    empty_rows = find(~any(p.F_struc(p.K.f+1:p.K.f+p.K.l,2:end),2));
    if ~isempty(empty_rows)
        if all(p.F_struc(p.K.f+empty_rows,1)>=0)            
            p.F_struc(p.K.f+empty_rows,:)=[];
            p.K.l = p.K.l - length(empty_rows);         
        else
            feasible = 0;
        end
    end
end

% ********************************
% Tighten bounds at root
% ********************************
if p.options.bmibnb.roottight & feasible
    lowersolver = eval(['@' p.solver.lowercall]);
    c = p.c;
    Q = p.Q;
    mt = p.monomtable;
    p.monomtable = eye(length(c));
    i = 1;
    while i<=length(p.linears) & feasible   
        j = p.linears(i);  
        p.c = eyev(length(p.c),j);
        output = feval(lowersolver,p);  
        if (output.problem == 0) & (output.Primal(j)>p.lb(j))
            p.lb(j) = output.Primal(j);  
            p = updateonenonlinearbound(p,j);
            p = clean_bounds(p);
        end
        if output.problem == 1
            feasible = 0;
        else
            % p = updatenonlinearbounds(p,0,1);
            p.c = -eyev(length(p.c),j);
            output = feval(lowersolver,p);
            if (output.problem == 0) & (output.Primal(j) < p.ub(j))
                p.ub(j) = output.Primal(j);
                p = updateonenonlinearbound(p,j);
                p = clean_bounds(p);
            end    
            if output.problem == 1
                feasible = 0;
            end             
            i = i+1;
        end
    end
    p.lb(p.lb<-1e10) = -inf;
    p.ub(p.ub>1e10) = inf;
    p.c = c;
    p.Q = Q;
    p.monomtable = mt;
end

if feasible
    
    % *******************************
    % Bounded domain?
    % *******************************
    involved = unique(p.nonlins(:,2:3));
    if isinf(p.lb(involved)) | isinf(p.ub(involved))    
        error('You have to bound all complicating variables explicitely (I cannot deduce bounds on all variables)')
        output.Primal = [];
        output.problem = -1;
    end
        
    % *******************************
    % We don't need to save node data
    % *******************************
    p.options.savesolverinput  = 0;
    p.options.savesolveroutput = 0;
    
    % *******************************
    % RUN BRANCH & BOUND
    % *******************************
    [x_min,solved_nodes,lower,upper] = branch_and_bound(p);
    
    % **********************************
    % 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
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)

% ***************************************
% LPs ARE USED IN  BOX-REDUCTION 
% (this is essentially a cutting plane pool)
% ***************************************
p.lpcuts = p.F_struc(1+p.K.f:1:p.K.l,:);

% ***************************************
% Create function handles to solvers
% ***************************************
try
    lowersolver = eval(['@' p.solver.lowercall]); % Local LMI solver
    uppersolver = eval(['@' p.solver.uppercall]); % Local BMI solver
    lpsolver    = eval(['@' p.solver.lpcall]);    % 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;
p.options.saveduals = 0;
options = p.options;

% ************************************************
% ORIGINAL PROBLEM (used in LOCAL BMI solver)
% ************************************************
p_upper = p;

% ************************************************
% Remove linear cutting planes from problem
% ************************************************
p_upper.F_struc(p_upper.K.f+p_upper.KCut.l,:)=[];
p_upper.K.l = p_upper.K.l - length(p_upper.KCut.l);

% ************************************************
% Remove sdp cutting planes from 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

% ************************************************
% INITILIZATION
% ************************************************
p.depth = 0;
p.dpos = 0; 
p.lower = NaN;
upper   = inf;
lower   = NaN;
gap = inf;
x_min   = zeros(length(p.c),1);
stack   = [];
solved_nodes = 0;
solved_lower = 0;
solved_upper = 0;
solved_lp = 0;

if isempty(p.x0)
    p.x0    = zeros(length(p.c),1);
end

x0 = evaluate_nonlinear(p,p.x0);
upper_residual = resids(p,x0);
x0_feasible = all(upper_residual(1:p.K.f)>=-options.bmibnb.eqtol) & all(upper_residual(1+p.K.f:end)>=options.bmibnb.pdtol);
if p.options.usex0 & x0_feasible
    x_min = x0;
    upper = p.f+p.c'*x0+x0'*Q*x0;
end

% ************************************************
% Branch & bound loop
% ************************************************
if options.bmibnb.verbose>0
    fprintf('******************************************************************************************************************\n')
    fprintf('#node          Was'' up                         gap      upper      node    lower  dpth  stk     Memory Vol-red\n')
    fprintf('******************************************************************************************************************\n')
end

doplot = 0;
if doplot
    close all;
    hold on;
end

t_start = cputime;
go_on  = 1;
while go_on
    
    if doplot;ellipplot(diag([200 50]),1,'y',[p.dpos;-p.depth]);drawnow;end;
    % ********************************************
    % ASSUME THAT WE WON'T FATHOME
    % ********************************************
    keep_digging = 1;
    % ********************************************
    % REDUCE BOX
    % ******************************************** 
    if ~options.bmibnb.lpreduce
        % [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,[]);  
        vol_reduction = 1;
        feasible = 1;
    else     
        [p,feasible,vol_reduction] =  boxreduce(p,upper,lower,lpsolver,options);  
    end
    
    % ********************************************
    % SOLVE LOWER AND UPPER
    % ********************************************
    if feasible               
        
        output = solvelower(p,options,lowersolver); 
        
        info_text = '';
        switch output.problem
            case 1
                if doplot;ellipplot(diag([200 25]),1,'r',[p.dpos;-p.depth]);drawnow;end;
                info_text = 'Infeasible node';
                keep_digging = 0;
                cost = inf;
                feasible = 0;
                
            case 2
                cost = -inf;
                
            case {0,3,4}
                
                x = output.Primal; 
                
                cost = f+c'*x+x'*Q*x;