bnb.m

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

        node0.ub = p0.ub;
        node0.depth = p0.depth;
        node0.lower = p0.lower;
        node0.x0 = p0.x0;
        node0.binary_variables = p0.binary_variables;

        stack = push(stack,node1);
        stack = push(stack,node0);
    end

    % Lowest cost in any open node
    if ~isempty(stack)
        lower = min([stack.lower]);
        if can_use_ceil_lower
            lower = ceil(lower);
        end
    end
    
    % **********************************
    % Get a new node to solve
    % **********************************
    [node,stack] = pull(stack,p.options.bnb.method,x_min,upper);
    if ~isempty(node)
        p.lb = node.lb;
        p.ub = node.ub;
        p.depth = node.depth;
        p.lower = node.lower;
        p.x0 = node.x0;
        p.binary_variables = node.binary_variables;
    end
    gap = abs((upper-lower)/(1e-3+abs(upper)+abs(lower)));
    if isnan(gap)
        gap = inf;
    end

    if p.options.bnb.verbose;fprintf(' %4.0f : %12.3E  %7.2f   %12.3E  %2.0f \n',solved_nodes,upper,100*gap,lower,length(stack)+length(node));end

end
if p.options.bnb.verbose;showprogress([num2str2(solved_nodes,3)  ' Finishing.  Cost: ' num2str(upper) ],p.options.bnb.verbose);end


%%
function stack = push(stackin,p)
if ~isempty(stackin)
    stack = [p;stackin];
else
    stack(1)=p;
end

%%
function [p,stack] = pull(stack,method,x_min,upper);

if ~isempty(stack)
    switch method
        case {'depth','depthfirst','depthbreadth','depthproject','depthbest'}
            [i,j]=max([stack.depth]);
            p=stack(j);
            stack = stack([1:1:j-1 j+1:1:end]);

        case 'breadth'
            [i,j]=min([stack.depth]);
            p=stack(j);
            stack = stack([1:1:j-1 j+1:1:end]);

        case 'best'
            [i,j]=min([stack.lower]);
            p=stack(j);
            stack = stack([1:1:j-1 j+1:1:end]);

        otherwise
    end
else
    p = [];
end

% **********************************
%% BRANCH VARIABLE
% **********************************
function [index,whatsplit] = branchvariable(x,integer_variables,binary_variables,options,x_min,Weight)
all_variables = [integer_variables(:);binary_variables(:)];
switch options.bnb.branchrule
    case 'first'
        index = min(find(abs(x(all_variables)-round(x(all_variables)))>options.bnb.inttol));
    case 'last'
        index = max(find(abs(x(all_variables)-round(x(all_variables)))>options.bnb.inttol));
    case 'min'
        nint = find(abs(x(all_variables)-round(x(all_variables)))>options.bnb.inttol);
        [val,index] = min(abs(x(nint)));
        index = nint(index);
    case 'max'
        [val,index] = max(abs(x(all_variables)-round(x(all_variables))));
        %   index=all_variables(index);
    otherwise
        error
end
if index<=length(integer_variables)
    whatsplit = 'integer';
else
    index = index-length(integer_variables);
    whatsplit = 'binary';
end

% **********************************
% SPLIT PROBLEM
% **********************************
function [p0,p1,variable] = binarysplit(p,x,index,lower,options)
p0 = p;
p1 = p;
%variable = index;
variable = p.binary_variables(index);

p0.ub(variable)=0;
p0.lb(variable)=0;

p0.lower = lower;
p0.depth = p.depth+1;
p0.binary_variables = setdiff1D(p0.binary_variables,variable);

p1.ub(variable)=1;
p1.lb(variable)=1;

p1.binary_variables = setdiff1D(p1.binary_variables,variable);
p1.lower = lower;
p1.depth = p.depth+1;

% % *****************************
% % PROCESS MOST PROMISING FIRST
% % (p0 in top of stack)
% % *****************************
if x(variable)>0.5
    pt=p1;
    p1=p0;
    p0=pt;
end

function [p0,p1] = integersplit(p,x,index,lower,options,x_min)

variable = p.integer_variables(index);
current = x(p.integer_variables(index));
lb = floor(current)+1;
ub = floor(current);

% xi<ub
p0 = p;
p0.lower = lower;
p0.depth = p.depth+1;
p0.x0(variable) = ub;
p0.ub(variable)=min(p0.ub(variable),ub);

% xi>lb
p1 = p;
p1.lower = lower;
p1.depth = p.depth+1;
p1.x0(variable) = lb;
p1.lb(variable)=max(p1.lb(variable),lb);

% *****************************
% PROCESS MOST PROMISING FIRST
% *****************************
if lb-current<0.5
    pt=p1;
    p1=p0;
    p0=pt;
end


function s = num2str2(x,d,c);
if nargin==3
    s = num2str(x,c);
else
    s = num2str(x);
end
s = [repmat(' ',1,d-length(s)) s];


function [stack,lower] = prune(stack,upper,options,solved_nodes,p)
% *********************************
% PRUNE STACK W.R.T NEW UPPER BOUND
% *********************************
if ~isempty(stack)
    toolarge = find([stack.lower]>upper*(1-1e-4));
    if ~isempty(toolarge)
        stack(toolarge)=[];
    end
end

if ~isempty(stack)
    lower = min([stack.lower]);
else
    lower = upper;
end

function p = adaptivestrategy(p,upper,solved_nodes)
% **********************************'
% SWITCH NODE SELECTION STRATEGY?
% **********************************'
if strcmp(p.options.bnb.method,'depthproject') & (upper<inf)
    p.options.bnb.method = 'project';
end
if strcmp(p.options.bnb.method,'depthbest') & (upper<inf)
    p.options.bnb.method = 'best';
end
if strcmp(p.options.bnb.method,'depthbreadth') & (upper<inf)
    p.options.bnb.method = 'breadth';
end
if strcmp(p.options.bnb.method,'depthest') & (upper<inf)
    p.options.bnb.method = 'est';    
end

function  output = solvelower(lowersolver,relaxed_p,upper,lower)

if all(relaxed_p.lb==relaxed_p.ub)
    x = relaxed_p.lb;
    res = resids(relaxed_p,x);
    if all(res>-1e-7)
        output.problem = 0;
    else
        output.problem = 1;
    end
    output.Primal = x;
    return
end

p = relaxed_p;
removethese = p.lb==p.ub;
if nnz(removethese)>0
    
    if ~isinf(upper) & nnz(p.Q)==0
        p.F_struc = [p.F_struc(1:p.K.f,:);upper-p.f -p.c';p.F_struc(1+p.K.f:end,:)];
        p.K.l=p.K.l+1;
    end

%     if ~isempty(p.cuts)
%         p.F_struc = [p.F_struc(1:p.K.f,:);p.cuts;p.F_struc(1+p.K.f:end,:)];
%         p.K.l=p.K.l+size(p.cuts,1);
%     end
% 
    if ~isempty(p.F_struc)
        p.F_struc(:,1)=p.F_struc(:,1)+p.F_struc(:,1+find(removethese))*p.lb(removethese);
        p.F_struc(:,1+find(removethese))=[];
    end

    p.c(removethese)=[];   
    if nnz(p.Q)>0
        p.c = p.c + 2*p.Q(find(~removethese),find(removethese))*p.lb(removethese);
        p.Q(:,find(removethese))=[];
        p.Q(find(removethese),:)=[];
    else
        p.Q = spalloc(length(p.c),length(p.c),0);
    end
    p.lb(removethese)=[];
    p.ub(removethese)=[];
    p.x0(removethese)=[];
    p.monomtable(:,find(removethese))=[];
    p.monomtable(find(removethese),:)=[];

    output = feval(lowersolver,p);

    x=relaxed_p.c*0;
    x(removethese)=relaxed_p.lb(removethese);
    x(~removethese)=output.Primal;
    output.Primal=x;
else
   output = feval(lowersolver,p); 
end


function res = resids(p,x)
res= [];
if p.K.f>0
    res = -abs(p.F_struc(1:p.K.f,:)*[1;x]);
end
if p.K.l>0
    res = [res;p.F_struc(p.K.f+1:p.K.f+p.K.l,:)*[1;x]];
end
if (length(p.K.s)>1) | p.K.s>0
    top = 1+p.K.f+p.K.l;
    for i = 1:length(p.K.s)
        n = p.K.s(i);
        X = p.F_struc(top:top+n^2-1,:)*[1;x];top = top+n^2;
        X = reshape(X,n,n);
        res = [res;min(eig(X))];
    end
end
res = [res;min([p.ub-x;x-p.lb])];


function cuts = addFlowCover(p,x);
X=reshape(p.F_struc(1+p.K.l:end,:)*[1;x],p.K.s(1),p.K.s(1));
[d,v] = eig(X);

for m = 1:1:1
m = 1;
gc = 0;

%     newcuts=1;
%     for j = 1:length(x)+1
%         newF(newcuts,j)= d(:,m)'*reshape(p.F_struc(1+p.K.l:end,j),p.K.s(1),p.K.s(1))*d(:,m);
%     end
% 
%      a = -newF(1,2:end);
%      b = newF(1);
%      
%      f=b-floor(b);
%      fj=a-floor(a);
%      
%      fjf=fj-f;
%      fjf(fjf<0)=0;
%      a=floor(a)+fjf/(1-f);
%      
%      floor(b)-a*x
%       p.F_struc = [p.F_struc(1:p.K.f,:);floor(b) -a;p.F_struc(1+p.K.f:end,:)];
%       p.K.l=p.K.l+1;

%     neg = find(a<0);
%     pos = find(a>0);
% 
%     b = floor(b+sum(abs(a(neg))))-sum(floor(abs(a(neg))));
% 
%     a(pos) = floor(a(pos));
%     a(neg) = -floor(abs(a(neg)));
%     b-a*x
% if 0%1%b-a*x < 0
%         p.F_struc = [p.F_struc(1:p.K.f,:);b -a;p.F_struc(1+p.K.f:end,:)];
%         p.K.l=p.K.l+1;
%     end


    %     newF = p.F_struc(p.K.l,:)*(1+randn(1)^2);
    %     a = -newF(1,2:end);
    %     b = newF(1);
    %     neg = find(a<0);
    %     pos = find(a>0);
    %
    %     a(pos) = floor(a(pos));
    %     a(neg) = -floor(abs(a(neg)));
    %
    %     b = floor(b+sum(abs(a(neg))))-sum(floor(abs(a(neg))));
    %     if b-a*x < 0
    %         p.F_struc = [p.F_struc(1:p.K.f,:);b+sqrt(1e-3) -a;p.F_struc(1+p.K.f:end,:)];
    %         p.K.l=p.K.l+1;
    %     end
    %     try
    %         b=newF(1);
    %         a=-newF(2:end);
    %         aposind=find(a >  0);
    %         anegind=find(a <= 0);
    %         apos=a(aposind);
    %         aneg=a(anegind);
    %         aabs=abs(a);
    %         b=b-sum(aneg);
    %         [i,j]=sort(abs(a));
    %         %       [i,j]=sort(atil);
    %         j=fliplr(j);
    %         i=fliplr(i);
    %         C=min(find(cumsum(aabs(j))>b));
    %         if length(C)>0
    %             used = j(1:C);
    %             pos=used(a(used)>0);
    %             neg=used(a(used)<0);
    %             newF=newF*0;
    %             newF(1)=C-1-length(neg);
    %             newF(1,1+pos)=-1;
    %             newF(1,1+neg)=1;
    %             p.F_struc = [p.F_struc(1:p.K.f,:);newF;p.F_struc(1+p.K.f:end,:)];
    %             Sminus = length(find(newF(2:end)>0));
    %             if 1-Sminus == newF(1)
    %                 Sminus
    %             end
    %             p.K.l=p.K.l+1;
    %         end
    %     catch
    %         error
    %     end
end
cuts = [];

while m <= 1%p.K.s(1) & gc<=10

%     newcuts=1;
%     for j = 1:length(x)+1
%         newF(newcuts,j)= d(:,m)'*reshape(p.F_struc(1+p.K.l:end,j),p.K.s(1),p.K.s(1))*d(:,m);
%     end

% %     j = 1;
%     while j <=15 & gc<=10
%         lambda = 1+rand(1)*10;
%         a = -newF(1,2:end)*lambda;
%         b = newF(1)*lambda;
%         zerothese = find([p.lb==p.ub]);
%         if length(zerothese)>0
%             b = b-a(zerothese)*p.lb(zerothese);
%             a(zerothese) = 0;
%         end
% 
%         if nnz(a)>0
%             neg = find(a<0);
%             pos = find(a>0);
% 
%             b = floor(b+sum(abs(a(neg))))-sum(floor(abs(a(neg))));
% 
%             a(pos) = floor(a(pos));
%             a(neg) = -floor(abs(a(neg)));
% 
%             if (b-a*x < 0) & nnz(a)>0
%                 cuts = [cuts;b -a];
%                 gc  = gc + 1;
%             end
%             j = j+1;
%         end
%     end
% 
%      m = m+1;
% % 
%     b=newF(1);
%     a=-newF(2:end);
%     aposind=find(a >  0);
%     anegind=find(a <= 0);
%     apos=a(aposind);
%     aneg=a(anegind);
%     aabs=abs(a);
%     b=b-sum(aneg);
%     [i,j]=sort(abs(a));
%     %       [i,j]=sort(atil);
%     j=fliplr(j);
%     i=fliplr(i);
%     C=min(find(cumsum(aabs(j))>b));
%     if length(C)>0
%         used = j(1:C);
%         pos=used(a(used)>0);
%         neg=used(a(used)<0);
%         newF=newF*0;
%         newF(1)=C-1-length(neg);
%         newF(1,1+pos)=-1;
%         newF(1,1+neg)=1;
%         cuts = [cuts;newF];
%     end
end

function p = Updatecostbound(p,upper);
if p.simplecost
    ind = find(p.c);
    if p.c(ind)>0
        p.ub(ind) = min(p.ub(ind),(upper-p.f)/p.c(ind));
    else
        p.lb(ind) = max(p.lb(ind),(p.f-upper)/abs(p.c(ind)));
    end
end

function [x_min,upper] = initializesolution(p);

x_min = zeros(length(p.c),1);
upper = inf;
if p.options.usex0
    z = p.x0;
    residual = resids(p,z);
    relaxed_feasible = all(residual(1:p.K.f)>=-1e-12) & all(residual(1+p.K.f:end)>=-1e-6);
    if relaxed_feasible
        upper = p.f+p.c'*z+z'*p.Q*z;
        x_min = z;
    end
else
    p.x0 = zeros(length(p.c),1);
    x = p.x0;
    z = evaluate_nonlinear(p,x);
    residual = resids(p,z);
    relaxed_feasible = all(residual(1:p.K.f)>=-p.options.bmibnb.eqtol) & all(residual(1+p.K.f:end)>=p.options.bmibnb.pdtol);
    if relaxed_feasible
        upper = p.f+p.c'*z+z'*p.Q*z;
        x_min = x;
    end   
end