bmibnb.m

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

need_to_add = find(violation < -1e-4);
if ~isempty(need_to_add)
    %    need_to_add
    %[i,j] = sort(violation(need_to_add));
    % p.constraintState(inactiveCuts(need_to_add(j(1:10)))) = 1;
    p.constraintState(inactiveCuts(need_to_add)) = 1;
end


% *************************************************************************
% Strategy for deciding which variable to branch on
% *************************************************************************
function spliton = branchvariable(p,options,x)
% Split if box is to narrow
width = abs(p.ub(p.branch_variables)-p.lb(p.branch_variables));
if (min(width)/max(width) < 0.1) | (size(p.bilinears,1)==0) % 
    [i,j] = max(width);
    spliton = p.branch_variables(j);
else
    res = x(p.bilinears(:,1))-x(p.bilinears(:,2)).*x(p.bilinears(:,3));
    [ii,jj] = sort(abs(res));
    v1 = p.bilinears(jj(end),2);
    v2 = p.bilinears(jj(end),3);
    
    acc_res1 = sum(abs(res(find((p.bilinears(:,2)==v1) |  p.bilinears(:,3)==v1))));
    acc_res2 = sum(abs(res(find((p.bilinears(:,2)==v2) |  p.bilinears(:,3)==v2))));
    
    if ~ismember(v2,p.branch_variables) | (acc_res1>acc_res2)
        spliton = v1;
    else
        spliton = v2;
    end
end

% *************************************************************************
% Strategy for diving the search space
% *************************************************************************
function bounds = partition(p,options,spliton,x_min)
x = x_min;
if isinf(p.lb(spliton))
    p.lb(spliton) = -1e6;
end
if isinf(p.ub(spliton))
    p.ub(spliton) = 1e6;
end

switch options.bmibnb.branchrule
case 'omega'
    if ~isempty(x_min)
        U = p.ub(spliton);
        L = p.lb(spliton);
        x = x(spliton);
        bounds = [p.lb(spliton) 0.5*max(p.lb(spliton),min(x_min(spliton),p.ub(spliton)))+0.5*(p.lb(spliton)+p.ub(spliton))/2 p.ub(spliton)];
    else
        bounds = [p.lb(spliton) (p.lb(spliton)+p.ub(spliton))/2 p.ub(spliton)];
    end
case 'bisect'
    bounds = [p.lb(spliton) (p.lb(spliton)+p.ub(spliton))/2 p.ub(spliton)];
otherwise
    bounds = [p.lb(spliton) (p.lb(spliton)+p.ub(spliton))/2 p.ub(spliton)];
end
if isnan(bounds(2)) %FIX
    if isinf(p.lb(spliton))
        p.lb(spliton) = -1e6;
    end
    if isinf(p.ub(spliton))
        p.ub(spliton) = 1e6;
    end
    bounds(2) = (p.lb(spliton)+p.ub(spliton))/2;
end

function [p,stack] = selectbranch(p,options,stack,x_min,upper);
switch options.bmibnb.branchmethod
case 'maxvol'
    [node,stack] = pull(stack,'maxvol',x_min,upper,p.branch_variables);
case 'best'
    [node,stack] = pull(stack,'best',x_min,upper);
otherwise
    [node,stack] = pull(stack,'best',x_min,upper);
end
% Copy node data to p
if isempty(node)
    p = [];
else
    p.depth = node.depth;
    p.dpos = node.dpos;
    p.lb = node.lb;
    p.ub = node.ub;
    p.lower = node.lower;
    p.lpcuts = node.lpcuts;
    p.x0 = node.x0;
    p.constraintState = node.constraintState;
    p.cutState = node.cutState;
end



function output = solvelower(p,options,lowersolver)

if 0%p.options.do
    removeThese = find(p.constraintState~=1);
    p.F_struc(p.K.f + removeThese,:) = [];
    p.K.l = p.K.l - length(removeThese);
    p_with_bilinear_cuts = p;
else
    removeThese = find(p.constraintState==inf);
    p.F_struc(p.K.f + removeThese,:) = [];
    p.K.l = p.K.l - length(removeThese);
    p_with_bilinear_cuts = p;
end

if ~isempty(p.bilinears)
    p_with_bilinear_cuts.F_struc(1:p.K.f,:)=[];
    p_with_bilinear_cuts = addmcgormick(p_with_bilinear_cuts);
    p_with_bilinear_cuts.F_struc = [p.F_struc(1:p.K.f,:);p_with_bilinear_cuts.F_struc];
end

% **************************************
% SOLVE NODE PROBLEM
% **************************************
if any(p_with_bilinear_cuts.ub<p_with_bilinear_cuts.lb)
    output.problem=1;
else
    % We are solving relaxed problem (penbmi might be local solver)
    p_with_bilinear_cuts.monomtable = eye(length(p_with_bilinear_cuts.c));
    
    output = feval(lowersolver,p_with_bilinear_cuts);
    
    % Minor clean-up
    output.Primal(output.Primal<p.lb) = p.lb(output.Primal<p.lb);
    output.Primal(output.Primal>p.ub) = p.ub(output.Primal>p.ub);
end




function output = solveupper(p,p_original,x,options,uppersolver)

p_upper = p_original;

% Pick an initial point (this can be a bit tricky...)
% Use relaxed point, shifted towards center of box
if all(x<=p.ub) & all(x>=p.lb)
    p_upper.x0 = 0.1*x + 0.9*(p.lb+p.ub)/2;
else
    p_upper.x0 = (p.lb+p.ub)/2;
end
% Shift towards interior for variables with unbounded lower or upper
lbinfbounds = find(isinf(p.lb));
ubinfbounds = find(isinf(p.ub));
p_upper.x0(ubinfbounds) = x(ubinfbounds)+0.01;
p_upper.x0(lbinfbounds) = x(lbinfbounds)-0.01;
ublbinfbounds = find(isinf(p.lb) & isinf(p.ub));
p_upper.x0(ublbinfbounds) = x(ublbinfbounds);

% ...expand the current node just slightly
p_upper.lb = p.lb;
p_upper.ub = p.ub;
p_upper.lb(~isinf(p_original.lb)) = 0.99*p.lb(~isinf(p_original.lb))+p_original.lb(~isinf(p_original.lb))*0.01;
p_upper.ub(~isinf(p_original.ub)) = 0.99*p.ub(~isinf(p_original.ub))+p_original.ub(~isinf(p_original.ub))*0.01;
p_upper.lb(isinf(p_original.lb)) = p_upper.lb(isinf(p_original.lb)) - 0.001;
p_upper.ub(isinf(p_original.ub)) = p_upper.ub(isinf(p_original.ub)) + 0.001;
p_upper.options.saveduals = 0;

ub = p_upper.ub ;
lb = p_upper.lb ;
p_upper.x0 = p_upper.x0(p_upper.actuallyused);
p_upper.ub = p_upper.ub(p_upper.actuallyused);
p_upper.lb = p_upper.lb(p_upper.actuallyused);

% Solve upper bounding problem
p_upper.options.usex0 = 1;
output = feval(uppersolver,p_upper);
% Project into the box (numerical issue)
output.Primal(output.Primal<p_upper.lb) = p_upper.lb(output.Primal<p_upper.lb);
output.Primal(output.Primal>p_upper.ub) = p_upper.ub(output.Primal>p_upper.ub);


function vars = decide_branch_variables(p)

if size(p.bilinears,1)==0
    if p.K.s(1)>0
        if any(p.K.s>p.K.rank)
            vars = p.linears;
            return
        end
    end
end

if p.options.bmibnb.lowrank==0
    nonlinear = find(~(sum(p.monomtable~=0,2)==1 & sum(p.monomtable,2)==1));
    vars      =  find(sum(abs(full(p.monomtable(nonlinear,:))),1));
else
    % Try to find a bi-partite structure
    pool1 = p.bilinears(1,2);
    pool2 = p.bilinears(1,3);
    
    for i = 2:size(p.bilinears,1)
        v1 = p.bilinears(i,2);
        v2 = p.bilinears(i,3);
        if v1==v2
            % We are fucked
            pool1 = [pool1 v1];
            pool2 = [pool2 v2];
        else
            if ismember(v1,pool1)
                pool2 = [pool2 v2];
            elseif ismember(v1,pool2)
                pool1 = [pool1 v2];
            elseif ismember(v2,pool1)
                pool2 = [pool2 v1];
            elseif ismember(v2,pool2)
                pool1 = [pool1 v1];
            else
                % No member yet
                pool1 = [pool1 v1];
                pool2 = [pool2 v2];
            end
        end
    end
    pool1 = unique(pool1);
    pool2 = unique(pool2);
    if isempty(intersect(pool1,pool2))
        if length(pool1)<=length(pool2)
            vars = pool1;
        else
            vars = pool2;
        end
    else
        nonlinear = find(~(sum(p.monomtable~=0,2)==1 & sum(p.monomtable,2)==1));
        vars =  find(sum(abs(full(p.monomtable(nonlinear,:))),1));
    end
end


function x = evaluate_nonlinear(p,x)
if ~isempty(p.bilinears)
    x(p.bilinears(:,1)) = x(p.bilinears(:,2)).*x(p.bilinears(:,3));

    badones = find(sum(p.monomtable(p.bilinears(:,1),:),2)>2);
    if ~isempty(badones)
        redo = p.bilinears(badones,1);
        for i = 1:length(redo)
            these =  find(p.monomtable(redo(i),:));
            x(redo(i)) = prod(x(these).^((p.monomtable(redo(i),these))'));
        end
    end
end

function p = clean_bounds(p)
close = 1e-6>abs(p.ub - round(p.ub));
p.ub(close) = round(p.ub(close));
close = 1e-6>abs(p.lb - round(p.lb));
p.lb(close) = round(p.lb(close));
p.ub(p.binary_variables) = floor(p.ub(p.binary_variables) + 1e-2);
p.lb(p.lb<-1e12) = -inf;
p.ub(p.ub>1e12) = inf;


% *************************************************************************
% Some pre-calc of indicies for bilinear variables
% *************************************************************************
function [linears,bilinears,nonlinears] = compile_nonlinear_table(p)
linears = find(sum(p.monomtable,2)==1);
nonlinears = find(~(sum(p.monomtable~=0,2)==1 & sum(p.monomtable,2)==1));
bilinears   = [];
for i = 1:length(nonlinears)
    z = find(p.monomtable(nonlinears(i),:));
    if length(z)==1
        bilinears = [bilinears;nonlinears(i) z z];
    else
        bilinears = [bilinears;nonlinears(i) z(1) z(2)];
    end
end


% *************************************************************************
% Extract bounds on variable from inequalities
% *************************************************************************
function p = preprocess_bounds(p)

if ~isempty(p.bilinears)
    quadratics = find(p.bilinears(:,2)  == p.bilinears(:,3));
    for i = quadratics(:)'
        if ismember(p.bilinears(i,2),p.binary_variables)
            p.binary_variables = [p.binary_variables p.bilinears(i,1)];
        elseif ismember(p.bilinears(i,2),p.integer_variables)
            p.integer_variables = [p.integer_variables p.bilinears(i,1)];
        end
    end
end

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+p.nonlinears),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.lb(p.binary_variables) = max(0,p.lb(p.binary_variables));
p.ub(p.binary_variables) = min(1,p.ub(p.binary_variables));
p = clean_bounds(p);
if ~isempty(p.bilinears)
    p = updatenonlinearbounds(p);
end


% *************************************************************************
% Pre-solve (nothing coded yet)
% *************************************************************************
function [p,feasible] = simple_presolve(p)
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
            p.feasible = 0;
        end
    end
end


% *************************************************************************
% Tighten bounds at root
% *************************************************************************
function p = root_node_tighten(p,upper);
p.feasible = all(p.lb<=p.ub) & p.feasible;
if p.options.bmibnb.roottight & p.feasible
    lowersolver = eval(['@' p.solver.lowercall]);
    c = p.c;
    Q = p.Q;
    mt = p.monomtable;
    p.monomtable = eye(length(c));