📄 bmibnb.m
字号:
used = find(sum(([usedinquadratic;abs(p.c)';abs(p.F_struc(:,2:end))]),1));
usedinquadratic(p.bilinears(find(ismember(p.bilinears(:,1),used)),3)) = 1;
usedinquadratic(p.bilinears(find(ismember(p.bilinears(:,1),used)),2)) = 1;
not_used = find(~sum(([usedinquadratic;abs(p.c)';abs(p.F_struc(:,2:end))]),1));
%not_used = [];
p.actuallyused = setdiff(1:n_start,not_used);
dd = [];
i = 1;
while i <=length(not_used)
if ~ismember(not_used(i),p.bilinears(:,2)) & ~ismember(not_used(i),p.bilinears(:,3))
dd = [dd not_used(i)];
p.c(not_used(i)) = [];
p.Q(:,not_used(i)) = [];
p.Q(not_used(i),:) = [];
p.F_struc(:,1+not_used(i)) = [];
if find(p.bilinears(:,1)==not_used(i))
p.bilinears(find(p.bilinears(:,1)==not_used(i)),:)=[];
end
p.bilinears(p.bilinears(:,1)>not_used(i),1) = p.bilinears(p.bilinears(:,1)>not_used(i),1) - 1;
p.bilinears(p.bilinears(:,2)>not_used(i),2) = p.bilinears(p.bilinears(:,2)>not_used(i),2) - 1;
p.bilinears(p.bilinears(:,3)>not_used(i),3) = p.bilinears(p.bilinears(:,3)>not_used(i),3) - 1;
p.linears(p.linears == not_used(i)) = [];
p.linears(p.linears > not_used(i)) = p.linears(p.linears > not_used(i)) - 1;
not_used(not_used>not_used(i)) = not_used(not_used>not_used(i)) - 1;
p.monomtable(not_used(i),:)=[];
p.monomtable(:,not_used(i))=[];
end
i = i + 1;
end
else
p.actuallyused = 1:n_start;
end
%p.lb = p.lb(p.actuallyused);
%p.ub = p.ub(p.actuallyused);
%p.x0 = [];
p.branch_variables = [];
function pout = propagatequadratics(p,upper,lower)
pout = p;
if p.bilinears~=0
F_struc = p.F_struc;
p.F_struc = [-p.F_struc(1:p.K.f,:);p.F_struc];
p.K.f=2*p.K.f;
constraintState = [ones(p.K.f,1);p.constraintState];
if upper < inf
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;
constraintState = [1;constraintState];
end
if ~isnan(lower)
p.F_struc = [-(p.lower-abs(p.lower)*0.01)+p.f p.c';p.F_struc];
constraintState = [1;constraintState];
p.K.l = p.K.l + 1;
end
% for i = 1:p.K.l+p.K.f
% A = p.F_struc(i,2:end);
% b = p.F_struc(i,1);
%
% LB = 0;
% UB = 0;
% if any(ismember(p.bilinears(:,1),find(A))) & constraintState(i)==1
% for variable= p.linears'
% usedin = find((p.bilinears(:,2)==variable) | (p.bilinears(:,3)==variable));
% if ~isempty(usedin)
% if any(A(p.bilinears(usedin,1)))
% LB = b;
% UB = b;
% if (variable == 4) & (i==6)
%
% end
% GAINUP = A(variable);
% GAINDOWN = A(variable);
% for j = setdiff(1:size(A,2),variable)
% if isempty(find(p.bilinears(:,1)==j))
% % This is a linear variable
% if j~=variable
% if A(j)>=0
% LB = LB+A(j)*p.lb(j);
% UB = UB+A(j)*p.ub(j);
% else
% LB = LB+A(j)*p.ub(j);
% UB = UB+A(j)*p.lb(j);
% end
% end
% else
% s = find(p.bilinears(:,1)==j);
% if p.bilinears(s,2)==variable
% if A(j)>0
% GAINUP = GAINUP + A(j)*p.ub(p.bilinears(s,3));
% GAINDOWN = GAINDOWN + A(j)*p.lb(p.bilinears(s,3));
% else
% GAINUP = GAINUP + A(j)*p.lb(p.bilinears(s,3));
% GAINDOWN = GAINDOWN + A(j)*p.ub(p.bilinears(s,3));
% end
% elseif p.bilinears(s,3)==variable
% if A(j)>0
% GAINUP = GAINUP + A(j)*p.ub(p.bilinears(s,2));
% GAINDOWN = GAINDOWN + A(j)*p.lb(p.bilinears(s,2));
% else
% GAINUP = GAINUP + A(j)*p.lb(p.bilinears(s,2));
% GAINDOWN = GAINDOWN + A(j)*p.ub(p.bilinears(s,2));
% end
% else
% if A(j)>=0
% LB = LB+A(j)*p.lb(j);
% UB = UB+A(j)*p.ub(j);
% else
% LB = LB+A(j)*p.ub(j);
% UB = UB+A(j)*p.lb(j);
% end
% end
% end
% end
%
% %[GAINDOWN GAINUP LB UB]
% if GAINDOWN > 0
% if -UB/GAINDOWN > p.lb(variable)
% p.lb(variable) = max(p.lb(variable),-UB/GAINDOWN);
% end
% end
% if GAINUP < 0
% if -UB/GAINUP < p.ub(variable)
%
% p.ub(variable) = min(p.ub(variable),-UB/GAINUP);
% end
% end
% end
% end
% end
% end
% end
%[p.lb p.ub]
%
%[p.lb p.ub]
% sdpvar x1 x2 x3 x4 x5 x6 x7;
%
% e1 = - ((9 + (-6*x1) - 16*x2 - 15*x3)*x4 + (15 + (-6*x1) - 16*x2 - 15*x3)*x5) + x6 - 5*x7;
%
% e2= x3*x4 + x3*x5 - 50;
%
% e3= x4 + x6 - 100;
%
% e4= x5 + x7 - 200;
%
% e5= (3*x1 + x2 + x3 - 2.5)*x4 - 0.5*x6 - 0;
%
% e6= (3*x1 + x2 + x3 - 1.5)*x5 + 0.5*x7 - 0;
%
% e7= x1 + x2 + x3-1;
%
%
% F = set([e7] == 0) + set([e2 e3 e4 e5 e6] < 0);
% %F = F + set(0 < [x1 x2 x3 x4 x5 x6 x7] < [1 1 1 100 200 100 200]);
% F = F + set(0 < [x1 x2 x3 x4 x5 x6 x7] < [0.25 1 1 100 200 100 150]);
% %F = F + set((3*x1 + x2 + x3 - 1.5) < 0)% + set(3*x1 + x2 + x3 - 2.5>0);
% ops = sdpsettings('solver','bmibnb'); % Global solver
% ops = sdpsettings(ops,'bmibnb.lowersolver','glpk'); % Lower solver
% ops = sdpsettings(ops,'bmibnb.uppersolver','fmincon'); % Local solver
% ops = sdpsettings(ops,'bmibnb.lpsolver','glpk'); % LP solver
% solvesdp(F,e1,ops)
if p.K.l+p.K.f>0
quadratic_variables = find(p.bilinears(:,2) == p.bilinears(:,3));
if ~isempty(quadratic_variables)
quadratic_variables = p.bilinears(quadratic_variables,1);
for i = 1:length(quadratic_variables)
k = quadratic_variables(i);
if p.lb(k) >= 0
x = p.bilinears(p.bilinears(:,1)==k,2);% x^2
candidates = find((constraintState==1) & p.F_struc(1:p.K.f+p.K.l,1+k))';
for j = candidates
a = p.F_struc(j,2:end);
aij = a(k);
if aij > 0
indNEG = find(a < 0);
indPOS = find(a > 0);
LB = p.lb;
UB = p.ub;
LB(k) = 0;
UB(k) = 0;
a(k) = 0;
newLB = (-p.F_struc(j,1)-a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/aij;
p.lb(k) = max(p.lb(k),newLB);
if p.lb(x)>0
p.lb(x) = max(p.lb(x),sqrt(max(0,newLB)));
elseif p.ub(x)<0
p.ub(x) = min(p.ub(x),-sqrt(max(0,newLB)));
end
elseif aij < 0
indNEG = find(a < 0);
indPOS = find(a > 0);
LB = p.lb;
UB = p.ub;
LB(k) = 0;
UB(k) = 0;
a(k) = 0;
newUB = (p.F_struc(j,1)+a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/(-aij);
p.ub(k) = min(p.ub(k),newUB);
p.ub(x) = min(p.ub(x),sqrt(max(0,newUB)));
end
end
elseif p.lb(k)<0
[p.lb(k) p.ub(k)]
end
end
end
end
if p.K.l+p.K.f>0
bilinear_variables = find(p.bilinears(:,2) ~= p.bilinears(:,3));
if ~isempty(bilinear_variables)
bilinear_variables = p.bilinears(bilinear_variables,1);
for i = 1:length(bilinear_variables)
k = bilinear_variables(i);
if p.lb(k) >= -5000000000
x = p.bilinears(p.bilinears(:,1)==k,2);% x^2
y = p.bilinears(p.bilinears(:,1)==k,3);% x^2
candidates = find((constraintState==1) & p.F_struc(1:p.K.f+p.K.l,1+k))';
for j = candidates
a = p.F_struc(j,2:end);
aij = a(k);
if aij > 0
indNEG = find(a < 0);
indPOS = find(a > 0);
LB = p.lb;
UB = p.ub;
LB(k) = 0;
UB(k) = 0;
a(k) = 0;
newLB = (-p.F_struc(j,1)-a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/aij;
p.lb(k) = max(p.lb(k),newLB);
if p.lb(y) > 0
p.lb(x) = max(p.lb(x),newLB/p.ub(y));
end
if p.lb(x) > 0
p.lb(y) = max(p.lb(y),newLB/p.ub(x));
end
elseif aij < 0
indNEG = find(a < 0);
indPOS = find(a > 0);
LB = p.lb;
UB = p.ub;
LB(k) = 0;
UB(k) = 0;
a(k) = 0;
newUB = (p.F_struc(j,1)+a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/(-aij);
p.ub(k) = min(p.ub(k),newUB);
if p.lb(y) > 0
p.ub(x) = min(p.ub(x),newUB/p.lb(y));
end
if p.lb(x) > 0
p.ub(y) = min(p.ub(y),newUB/p.lb(x));
end
end
end
elseif p.lb(k)<0
[p.lb(k) p.ub(k)]
end
end
end
end
end
% if p.K.l+p.K.f>0
% linear_variables = setdiff(1:length(p.c),p.bilinears(:,1));
% if ~isempty(linear_variables)
% for i = 1:length(linear_variables)
% k = linear_variables(i);
% candidates = find(constraintState & p.F_struc(1:p.K.f+p.K.l,1+k))';
% for j = candidates
% a = p.F_struc(j,2:end);
% aij = a(k);
% if aij > 0
% indNEG = find(a < 0);
% indPOS = find(a > 0);
% LB = p.lb;
% UB = p.ub;
% LB(k) = 0;
% UB(k) = 0;
% a(k) = 0;
% newLB = (-p.F_struc(j,1)-a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/aij;
% if newLB > p.lb(k)
% [newLB p.lb(k)];
% end
% p.lb(k) = max(p.lb(k),newLB);
% if any(p.lb>p.ub)
% [p.lb p.ub];
% end
% elseif aij < 0
% indNEG = find(a < 0);
% indPOS = find(a > 0);
% LB = p.lb;
% UB = p.ub;
% LB(k) = 0;
% UB(k) = 0;
% a(k) = 0;
% newUB = (p.F_struc(j,1)+a([indPOS(:);indNEG(:)])*[UB(indPOS);LB(indNEG)])/(-aij);
% if newUB < p.ub(k)
% [newUB p.ub(k)];
% end
% p.ub(k) = min(p.ub(k),newUB);
% if any(p.lb>p.ub)
% [p.lb p.ub];
% end
% end
% end
% end
% end
% end
pout.lb = p.lb;
pout.ub = p.ub;
function [p,x_min,upper] = initializesolution(p);
x_min = zeros(length(p.c),1);
upper = inf;
if p.options.usex0
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
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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -