bmibnb.m

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

    i = 1;
    
    % Add an upper bound cut?
    if (upper < inf)
        % p.c'*x+p.f < upper
        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
    
    while i<=length(p.linears) & p.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
            p.feasible = 0;
        else
            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
                p.feasible = 0;
            end
            i = i+1;
        end
    end
    if upper < inf
        p.F_struc(1+p.K.f,:) = [];
        p.K.l = p.K.l - 1;
    end
    p.lb(p.lb<-1e10) = -inf;
    p.ub(p.ub>1e10) = inf;
    p.c = c;
    p.Q = Q;
    p.monomtable = mt;
end

% *************************************************************************
% Bound strengthening, main entrance
% *************************************************************************
function [p,feasible,vol_reduction] = domain_reduction(p,upper,lower,lpsolver);
% This is just too expensive
t1 = p.binary_variables;
t2 = p.integer_variables;
p.binary_variables = [];
p.integer_variables = [];
if ~p.options.bmibnb.lpreduce | (size(p.lpcuts,1)==0)
    vol_reduction = 1;
    feasible = 1;
else
    [p,feasible,vol_reduction] =  boxreduce(p,upper,lower,lpsolver,p.options);
end
p.binary_variables  = t1;
p.integer_variables = t2;

% *************************************************************************
% Bound strengthening
% *************************************************************************
function [p,feasible,vol_reduction] = boxreduce(p,upper,lower,lpsolver,options);

if options.bmibnb.lpreduce
    
    vol_start    = prod(p.ub(p.branch_variables)-p.lb(p.branch_variables));
    diag_before  =  sum(p.ub(p.branch_variables)-p.lb(p.branch_variables));
    
    [pcut,feasible,lower] = lpbmitighten(p,lower,upper,lpsolver);
    diag_after = sum(pcut.ub(p.branch_variables)-pcut.lb(p.branch_variables));
    iterations = 0;
    while (diag_after/(1e-18+diag_before) < 0.75    ) & feasible & iterations<4
        [pcut,feasible,lower] = lpbmitighten(pcut,lower,upper,lpsolver);
        diag_before = diag_after;
        diag_after = sum(pcut.ub(p.branch_variables)-pcut.lb(p.branch_variables));
        iterations = iterations + 1;
    end
    
    % Clean up...
    for i = 1:length(pcut.lb)
        if (pcut.lb(i)>pcut.ub(i)) & (pcut.lb-pcut.ub < 1e-3)
            pcut.lb(i)=pcut.ub(i);
            pcut = updatenonlinearbounds(pcut,i);
        end
    end
    p.lb = pcut.lb;
    p.ub = pcut.ub;
    
    % Metric = (V0/V)^(1/n)
    vol_reduction = max(0,min(1,(prod(p.ub(p.branch_variables)-p.lb(p.branch_variables))/(1e-12+vol_start))^(1/length(p.branch_variables))));
    p.lb(p.lb<-1e12) = -inf;
    p.ub(p.ub>1e12) = inf;
else
    vol_reduction = 1;
    feasible = 1;
end

% *************************************************************************
% Bound strengthening, here is where we actually solve LPs
% *************************************************************************
function [p,feasible,lower] = lpbmitighten(p,lower,upper,lpsolver)

% Construct problem with only linear terms
% and add cuts from lower/ upper bounds

p_test = p;
p_test.F_struc = p.lpcuts;
p_test.K.l = size(p.lpcuts,2);
p_test.K.f = 0;
p_test.K.s = 0;
p_test.K.q = 0;
if ~isnan(lower)
    p_test.F_struc = [-(p.lower-abs(p.lower)*0.01)+p.f p_test.c';p_test.F_struc];
end
if upper < inf
    p_test.F_struc = [upper+abs(upper)*0.01-p.f -p_test.c';p_test.F_struc];
end
if ~isempty(p.bilinears)
    p_test = addmcgormick(p_test);
end
p_test.K.l = size(p_test.F_struc,1);
p_test.F_struc = [p.F_struc(1:1:p.K.f,:);p_test.F_struc];
p_test.K.f = p.K.f;

feasible = 1;

% Perform reduction on most violating approximations at current solution
if p.options.bmibnb.lpreduce ~= 1
    n = ceil(max(p.options.bmibnb.lpreduce*length(p_test.linears),1));
    res = zeros(length(p.lb),1);
    for i = 1:size(p.bilinears,1)
        res(p.bilinears(i,2)) = res(p.bilinears(i,2)) + abs( p.x0(p.bilinears(i,1))-p.x0(p.bilinears(i,2)).*p.x0(p.bilinears(i,3)));
        res(p.bilinears(i,3)) = res(p.bilinears(i,3)) + abs( p.x0(p.bilinears(i,1))-p.x0(p.bilinears(i,2)).*p.x0(p.bilinears(i,3)));
    end
    res = res(p.linears);
    [ii,jj] = sort(abs(res));
    jj = jj(end-n+1:end);
else
    jj=1:length(p_test.linears);
end

j = 1;
while feasible & j<=length(jj)
    i = p_test.linears(jj(j));
    if abs(p.ub(i)-p.lb(i)>0.1)
        p_test.c = eyev(length(p_test.c),i);
        
        output = feval(lpsolver,p_test);
        if output.problem == 0
            if p_test.lb(i) < output.Primal(i)-1e-5
                p_test.lb(i) = output.Primal(i);
                p_test = updateonenonlinearbound(p_test,i);
            end
            p_test.c = -p_test.c;
            output = feval(lpsolver,p_test);
            if output.problem == 0
                if p_test.ub(i) > output.Primal(i)+1e-5
                    p_test.ub(i) = output.Primal(i);
                    p_test = updateonenonlinearbound(p_test,i);
                end
            end
            if output.problem == 1
                feasible = 0;
            end
        end
        if output.problem == 1
            feasible = 0;
        end
    end
    j = j + 1;
end
p_test = clean_bounds(p_test);
p.lb = p_test.lb;
p.ub = p_test.ub;


% *************************************************************************
% Heuristics from relaxed
% Basically nothing coded yet. Just check feasibility...
% *************************************************************************
function [upper,x_min,cost,info_text,numglobals] = heuristics_from_relaxed(p_upper,x,upper,x_min,cost,numglobals)
x(p_upper.binary_variables) = round(x(p_upper.binary_variables));
x(p_upper.integer_variables) = round(x(p_upper.integer_variables));

x = evaluate_nonlinear(p_upper,x);
z = x(p_upper.actuallyused);
%binaryused = find(ismember(p_upper.actuallyused,p_upper.binary_variables));
%z(binaryused) = round(z(binaryused));
p_upper.lb = p_upper.lb(p_upper.actuallyused);
p_upper.ub = p_upper.ub(p_upper.actuallyused);

relaxed_residual = resids(p_upper,z);

eq_ok = all(relaxed_residual(1:p_upper.K.f)>=-p_upper.options.bmibnb.eqtol);
iq_ok = all(relaxed_residual(1+p_upper.K.f:end)>=p_upper.options.bmibnb.pdtol);



% temporary hack for ranks
ranks_ok = check_rank_feasibility(p_upper,z);

relaxed_feasible = eq_ok & iq_ok & ranks_ok;
info_text = '';
if relaxed_feasible
    this_upper = p_upper.f+p_upper.c'*z+z'*p_upper.Q*z;
    if (this_upper < (1-1e-5)*upper) & (this_upper < upper - 1e-5)
        x_min = x;
        upper = this_upper;
        info_text = 'Improved solution';
        cost = cost-1e-10; % Otherwise we'll fathome!
        numglobals = numglobals + 1;
    end
end

% *************************************************************************
% Solve local upper bound problem
% *************************************************************************
function [upper,x_min,info_text,numglobals] = solve_upper_in_node(p,p_upper,x,upper,x_min,uppersolver,info_text,numglobals);
output = solveupper(p,p_upper,x,p.options,uppersolver);
integerused = find(ismember(p_upper.actuallyused,p_upper.integer_variables));
output.Primal(integerused) = round(output.Primal(integerused));
binaryused = find(ismember(p_upper.actuallyused,p_upper.binary_variables));
output.Primal(binaryused) = round(output.Primal(binaryused));

xu = evaluate_nonlinear(p_upper,output.Primal);
%binaryused = find(ismember(p_upper.actuallyused,p_upper.binary_variables));
%xu(binaryused) = round(xu(binaryused));
p_upper.lb = p_upper.lb(p_upper.actuallyused);
p_upper.ub = p_upper.ub(p_upper.actuallyused);

upper_residual = resids(p_upper,xu);
%feasible = ((output.problem == 0) | (output.problem == 5) | (output.problem == 3)  | (output.problem == 5)) & (all(upper_residual(1:p_upper.K.f)>=-p.options.bmibnb.eqtol) & all(upper_residual(1+p_upper.K.f:end)>=p.options.bmibnb.pdtol));
feasible = ~isempty(xu) & ~any(isnan(xu)) & all(upper_residual(1:p_upper.K.f)>=-p.options.bmibnb.eqtol) & all(upper_residual(1+p_upper.K.f:end)>=p.options.bmibnb.pdtol);

feasible = feasible &  check_rank_feasibility(p,xu);

%feasible = ((output.problem == 0) ) & (all(upper_residual(1:p_upper.K.f)>=-p.options.bmibnb.eqtol) & all(upper_residual(1+p_upper.K.f:end)>=p.options.bmibnb.pdtol));
if feasible
    this_upper = p_upper.f+p_upper.c'*xu+xu'*p_upper.Q*xu;
    if (this_upper < (1-1e-5)*upper) & (this_upper < upper - 1e-5)
        x_min = zeros(length(p.c),1);
        x_min(p_upper.actuallyused) = xu;
        upper = this_upper;
        info_text = 'Improved solution';
        numglobals = numglobals + 1;
    end
end


function ranks_ok = check_rank_feasibility(p,x);

ranks_ok = 1;
if p.K.s(1)>0
    if any(p.K.s>p.K.rank)
        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);
            ranks(i) = rank(X,1e-8);
        end
        ranks_ok = all(ranks <= p.K.rank);
    end
end


% *************************************************************************
% Detect redundant constraints
% *************************************************************************
function p = remove_redundant(p);

b = p.F_struc(1+p.K.f:p.K.l+p.K.f,1);
A = -p.F_struc(1+p.K.f:p.K.l+p.K.f,2:end);

redundant = find(((A>0).*A*(p.ub-p.lb) - (b-A*p.lb) <-1e-2));

% [row,col,vals] = find(A);
% pos = vals > 0;
% neg = vals < 0;
% vals(pos) = vals(pos).*p.ub(col(pos));
% vals(neg) = vals(neg).*p.lb(col(neg));
% redundant =find(sum(sparse(row,col,vals),2) < b);

if length(redundant)>1
    p.constraintState(redundant) = inf;
end

if p.options.bmibnb.lpreduce
    b = p.lpcuts(:,1);
    A = -p.lpcuts(:,2:end);
    redundant = find(((A>0).*A*(p.ub-p.lb) - (b-A*p.lb) <-1e-2));
    if length(redundant)>1
        p.lpcuts(redundant,:) = [];
        p.cutState(redundant) = [];
    end
end

% *************************************************************************
% Use equality constraints to derive bound
% *************************************************************************
function p = propagateequalities(p)
if p.K.f >0
    for j = 1:p.K.f
        if p.F_struc(j,1)==0
            [row,col,val] = find(p.F_struc(j,:));
            if length(row) == 2
                if val(1) == -val(2)
                    p.lb(col(1)-1) = max(p.lb(col(1)-1),p.lb(col(2)-1));
                    p.lb(col(2)-1) = max(p.lb(col(1)-1),p.lb(col(2)-1));
                    p.ub(col(1)-1) = min(p.ub(col(1)-1),p.ub(col(2)-1));
                    p.ub(col(2)-1) = min(p.ub(col(1)-1),p.ub(col(2)-1));
                end
            end
        end
    end
end


% *************************************************************************
% Clean the upper bound model
% Remove cut constraints, and
% possibly unused variables not used
% *************************************************************************
function p = cleanuppermodel(p);

% Remove cuts
p.F_struc(p.K.f+p.KCut.l,:)=[];
p.K.l = p.K.l - length(p.KCut.l);
n_start = length(p.c);
% Quadratic mode, and quadratic aware solver?
if ~isempty(p.bilinears)
    used_in_c = find(p.c);
    quadraticterms = used_in_c(find(ismember(used_in_c,p.bilinears(:,1))));
    if ~isempty(quadraticterms) & p.solver.uppersolver.objective.quadratic.nonconvex
        usedinquadratic = zeros(1,length(p.c));
        for i = 1:length(quadraticterms)
            Qij = p.c(quadraticterms(i));
            j = find(p.bilinears(:,1) == quadraticterms(i));
            x = p.bilinears(j,2);
            y = p.bilinears(j,3);
            p.Q(x,y) = Qij/2;
            p.Q(y,x) = p.Q(y,x)+Qij/2;
            p.c(quadraticterms(i)) = 0;
            usedinquadratic(x)=1;
            usedinquadratic(y)=1;
        end
    else
        usedinquadratic = zeros(1,length(p.c));
    end
    
    
    % OK, what is left