bmibnb.m

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

    else        
        info_text = 'Infeasible';
        keep_digging = 0;
        cost = inf;
        feasible = 0;
    end
    solved_nodes = solved_nodes+1;
    
    % ************************************************
    % PRUNE SUBOPTIMAL REGIONS BASED
    % ON NEW UPPER BOUND
    % ************************************************
    if ~isempty(stack)
        [stack,lower] = prune(stack,upper,options,solved_nodes,p);
    end
    lower = min(lower,cost);
    
    % ************************************************
    % CONTINUE SPLITTING?
    % ************************************************
    if keep_digging & max(p.ub(p.branch_variables)-p.lb(p.branch_variables))>options.bmibnb.vartol
        spliton = branchvariable(p,options,x);
        bounds  = partition(p,options,spliton,x);
        for i = 1:length(bounds)-1
            node = savetonode(p,spliton,bounds(i),bounds(i+1),-1,x,cost,p.constraintState,p.cutState);
            node.bilinears = p.bilinears;
            node = updateonenonlinearbound(node,spliton);
            stack = push(stack,node);
        end
        lower = min([stack.lower]);
    end
    
    % ************************************************
    %  Pick and create a suitable node
    % ************************************************
    [p,stack] = selectbranch(p,options,stack,x_min,upper);
    
    if isempty(p)
        if ~isinf(upper)
            relgap = 0;
        end
        if isinf(upper) & isinf(lower)
            relgap = inf;
        end
        depth = 0;
    else
        relgap = 100*(upper-lower)/(1+abs(upper));
        depth = p.depth;
    end
    if options.bmibnb.verbose>0
        fprintf(' %4.0f : %12.3E  %7.2f   %12.3E  %2.0f  %s  \n',solved_nodes,upper,relgap,lower,length(stack)+length(p),info_text);
    end
    
    absgap = upper-lower;
    % ************************************************
    % Continue?
    % ************************************************
    time_ok = cputime-t_start < options.bmibnb.maxtime;
    iter_ok = solved_nodes < options.bmibnb.maxiter;
    any_nodes = ~isempty(p);
    relgap_too_big = (isinf(lower) | isnan(relgap) | relgap>100*options.bmibnb.relgaptol);
    absgap_too_big = (isinf(lower) | isnan(absgap) | absgap>options.bmibnb.absgaptol);
    taget_not_met = upper>options.bmibnb.target;
    go_on = taget_not_met & time_ok & any_nodes & iter_ok & relgap_too_big & absgap_too_big ;
end
if options.bmibnb.verbose>0
    if options.bmibnb.verbose;showprogress([num2str(solved_nodes)  ' Finishing.  Cost: ' num2str(upper) ' Gap: ' num2str(relgap) '%'],options.bnb.verbose);end
end

% *************************************************************************
% Stack functionality
% *************************************************************************
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,branch_variables);
if ~isempty(stack)
    switch method
    case 'maxvol'
        for i = 1:length(stack)
            vol(i) = sum(stack(i).ub(branch_variables)-stack(i).lb(branch_variables));
        end
        [i,j] = max(vol);
        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

function [stack,lower] = prune(stack,upper,options,solved_nodes,p)
if ~isempty(stack)
    toolarge = find([stack.lower]>upper*(1+1e-4));
    if ~isempty(toolarge)        
        stack(toolarge)=[];
    end
    if ~isempty(stack)
        indPOS = find(p.c>0);
        indNEG = find(p.c<0);
        LB = [stack.lb];
        UB = [stack.ub];
        %LOWER = p.c(indPOS)'*LB(indPOS,:)+p.c(indNEG)'*UB(indNEG,:);
        LOWER =  p.c([indPOS(:);indNEG(:)])'*[LB(indPOS,:);UB(indNEG,:)];
        toolarge = find(LOWER > upper*(1+1e-4));
        stack(toolarge)=[];
    end           
end
if ~isempty(stack)
    lower = min([stack.lower]);    
else
    lower = upper;
end

function node = savetonode(p,spliton,bounds1,bounds2,direction,x,cost,constraintState,cutState);
node.lb = p.lb;
node.ub = p.ub;
node.lb(spliton) = bounds1;
node.ub(spliton) = bounds2;
if direction == -1
    node.dpos = p.dpos-1/(2^sqrt(p.depth));
else
    node.dpos = p.dpos+1/(2^sqrt(p.depth));
end
node.depth = p.depth+1;
node.x0 = x;
node.lpcuts = p.lpcuts;
node.lower = cost;
node.constraintState = constraintState;
node.cutState = cutState;


% *************************************************************************
% Code for setting the numerical values of nonlinear terms
% *************************************************************************
function p = updateonenonlinearbound(p,changed_var)
if ~isempty(p.bilinears)
    impactedVariables = find((p.bilinears(:,2) == changed_var) | (p.bilinears(:,3) == changed_var));
    x = p.bilinears(impactedVariables,2);
    y = p.bilinears(impactedVariables,3);
    z = p.bilinears(impactedVariables,1);
    x_lb = p.lb(x);
    x_ub = p.ub(x);
    y_lb = p.lb(y);
    y_ub = p.ub(y);
    bounds = [x_lb.*y_lb x_lb.*y_ub x_ub.*y_lb x_ub.*y_ub];
    p.lb(z) = max([p.lb(z) min(bounds,[],2)],[],2);
    p.ub(z) = min([p.ub(z) max(bounds,[],2)],[],2)';
    p.lb(impactedVariables(x==y)<0) = 0;
end

function p = updatenonlinearbounds(p,changed_var,keepbest);
if ~isempty(p.bilinears)
    x = p.bilinears(:,2);
    y = p.bilinears(:,3);
    z = p.bilinears(:,1);
    x_lb = p.lb(x);
    x_ub = p.ub(x);
    y_lb = p.lb(y);
    y_ub = p.ub(y);
    bounds = [x_lb.*y_lb x_lb.*y_ub x_ub.*y_lb x_ub.*y_ub];
    p.lb(z) = max([p.lb(z) min(bounds,[],2)],[],2);
    p.ub(z) = min([p.ub(z) max(bounds,[],2)],[],2);
    quadratic_variables = p.bilinears(x==y,1);
    %p.lb(p.bilinears((x==y)<0,1)) = 0;
    p.lb(quadratic_variables(p.lb(quadratic_variables)<0)) = 0;
end

% *************************************************************************
% Compute residual in constraints (smallest eigenvalue etc)
% *************************************************************************
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])];


% *************************************************************************
% Define constraints for the convex hull of bilinear terms
% *************************************************************************
function pcut = addmcgormick(p)

pcut = p;

z = p.bilinears(:,1);
x = p.bilinears(:,2);
y = p.bilinears(:,3);
x_lb = p.lb(x);
x_ub = p.ub(x);
y_lb = p.lb(y);
y_ub = p.ub(y);
m = length(x);
one = ones(m,1);
general_vals =[x_lb.*y_lb one -y_lb -x_lb,x_ub.*y_ub one -y_ub -x_ub,-x_ub.*y_lb -one y_lb x_ub,-x_lb.*y_ub -one y_ub x_lb]';  %3
general_cols = [one z+1 x+1 y+1 one z+1  x+1 y+1 one z+1 x+1 y+1 one z+1 x+1 y+1]';
general_row = [1;1;1;1;2;2;2;2;3;3;3;3;4;4;4;4];
quadratic_row = [1;1;1;2;2 ;2; 3; 3; 3];

quadratic_cols =  [one  z+1 x+1 one z+1 x+1 one  z+1 x+1]';
quadratic_vals = [-x_ub.*x_lb -one x_lb+x_ub x_lb.*y_lb one -y_lb-x_lb x_ub.*y_ub one -y_ub-x_ub]';

m = 1+length(p.c);
rows = [];
cols = [];
vals = [];
nrow = 0;
dummy = ismembc(p.bilinears(:,1),p.r);
for i = 1:size(p.bilinears,1)
    x = p.bilinears(i,2);
    y = p.bilinears(i,3);
    if x~=y
        if dummy(i)
            here = find(p.bilinears(i,1) == p.r);
            if p.s(here) > 0
                rows = [rows;general_row(1:8,:)+nrow];
                vals = [vals;general_vals(1:8,i)];
                cols = [cols;general_cols(1:8,i)];
                nrow = nrow + 2;
            elseif p.s(here)<0
                rows = [rows;general_row(1:8,:)+nrow];
                vals = [vals;general_vals(9:end,i)];
                cols = [cols;general_cols(9:end,i)];
                nrow = nrow + 2;
            end
        else
            rows = [rows;general_row+nrow];
            vals = [vals;general_vals(:,i)];
            cols = [cols;general_cols(:,i)];
            nrow = nrow + 4;
        end
    else
        if dummy(i)
            %col = quadratic_cols(:,i);
            %val = quadratic_vals(:,i);
            here = find(p.bilinears(i,1) == p.r);
            if p.s(here) > 0
                rows = [rows;quadratic_row(1:6,:)+nrow];
                vals = [vals;quadratic_vals(4:end,i)];
                cols = [cols;quadratic_cols(4:end,i)];
                nrow = nrow + 2;
            elseif p.s(here)<0
                rows = [rows;quadratic_row(1:3,:)+nrow];
                vals = [vals;quadratic_vals(1:3,i)];
                cols = [cols;quadratic_cols(1:3,i)];
                nrow = nrow + 1;
                %rows = [rows;general_row+nrow];
                %vals = [vals;val];
                %cols = [cols;col];
                %nrow = nrow + 4;
            end
        else
            col = quadratic_cols(:,i);
            val = quadratic_vals(:,i);
            rows = [rows;quadratic_row+nrow];
            vals = [vals;val];
            cols = [cols;col];
            nrow = nrow + 3;
        end
    end
end

F_temp = sparse(rows,cols,vals,nrow,m);
keep = find(~isinf(F_temp(:,1)));
F_temp = F_temp(keep,:);
pcut.F_struc = [F_temp;pcut.F_struc];
pcut.K.l = pcut.K.l+size(F_temp,1);


% *************************************
% DERIVE LINEAR CUTS FROM SDPs
% *************************************
function p = addsdpcut(p,x)
if p.K.s > 0
    top = p.K.f+p.K.l+1;
    newcuts = 1;
    newF = [];
    for i = 1:length(p.K.s)
        n = p.K.s(i);
        X = p.F_struc(top:top+n^2-1,:)*[1;x];
        X = reshape(X,n,n);
        [d,v] = eig(X);
        for m = 1:length(v)
            if v(m,m)<0
                for j = 1:length(x)+1;
                    newF(newcuts,j)= d(:,m)'*reshape(p.F_struc(top:top+n^2-1,j),n,n)*d(:,m);
                end
                % max(abs(newF(:,2:end)),[],2)
                newF(newcuts,1)=newF(newcuts,1)+1e-6;
                newcuts = newcuts + 1;
                if size(p.lpcuts,1)>0
                    dist = p.lpcuts*newF(newcuts-1,:)'/(newF(newcuts-1,:)*newF(newcuts-1,:)');
                    if any(abs(dist-1)<1e-3)
                        newF = newF(1:end-1,:);
                        newcuts = newcuts - 1;
                    end
                end
            end
        end
        top = top+n^2;
    end
    
    if ~isempty(newF)
        % Don't keep all
        m = size(newF,2);
        %  size(p.lpcuts)
        p.lpcuts = [newF;p.lpcuts];
        p.cutState = [ones(size(newF,1),1);p.cutState];
        violations = p.lpcuts*[1;x];
        p.lpcuts = p.lpcuts(violations<0.1,:);
        
        if size(p.lpcuts,1)>15*m
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!');
            violations = p.lpcuts*[1;x];
            [i,j] = sort(violations);
            %p.lpcuts = p.lpcuts(j(1:15*m),:);
            %p.cutState = lpcuts = p.lpcuts(j(1:15*m),:);
            %p.lpcuts = p.lpcuts(end-15*m+1:end,:);
        end
    end
end

function p = addlpcuts(p,z)
inactiveCuts = find(~p.cutState);
violation = p.lpcuts(inactiveCuts,:)*[1;z];
need_to_add = find(violation < -1e-4);
if ~isempty(need_to_add)
    %    need_to_add
    %[i,j] = sort(violation(need_to_add));
    p.cutState(inactiveCuts(need_to_add)) = 1;
    %p.cutState(inactiveCuts(need_to_add(j(1:10)))) = 1;
end
inactiveCuts = find(p.constraintState == 0 );
violation = p.F_struc(p.K.f+inactiveCuts,:)*[1;z];