function branch_boundary
clc
%format long
tic
A=[-1 0 0;0 -1 0;0 0 -1;1 1 1;1 1 -1/4;-2 -2 1;0 0 1];
[m,n]=size(A);
R=2;
b=[0 0 0 2 1 1 3];
D=[A b'];
[B,L,v0]=ba(A,b); %---------利用约束条件找到一个顶点-----------------
%---------B是对顶点v0的约束指标集,L是B的补集---
P=[];
ff=A(B(1),:);
P=[ff b(B(1))];
for i=2:n
j=B(i);
P=[P;A(j,:) b(j)];
ff=ff+A(j,:);
end
[a,gama0]=linprog(ff,A,b');
gama0=-gama0;
K=P;
P=[P;-ff gama0];
V=[v0];
T=K;
for i=1:n
T(i,:)=P(n+1,:);
vi=T(:,1:n)\T(:,n+1);
V=[V vi];
T=K;
end
tt=fun(V(:,1));
for i=2:n+1
tt=[tt fun(V(:,i))];
end
[alf0,u]=linprog(tt,A*V,b',ones(1,1+n),1,zeros(n+1,1));%-----u最优值下界
x=V*alf0;%-----最优解
r=fun(x);%-----最优值上界
flag=1;
k=1;
while(flag)
fprintf('迭代%d\n',k)
disp('单纯型为：')
V
disp('上界为：')
r
disp('下界为：')
u
if r==u
flag=0;
else
[V1,V2]=fj(V);
%-----计算下界--------
tt=[];
tt=fun(V1(:,1));
for i=2:n+1
tt=[tt fun(V1(:,i))];
end
[alf1,u1]=linprog(tt,A*V1,b',ones(1,1+n),1,zeros(n+1,1));
tt=[];
tt=fun(V2(:,1));
for i=2:n+1
tt=[tt fun(V2(:,i))];
end
[alf2,u2]=linprog(tt,A*V2,b',ones(1,1+n),1,zeros(n+1,1));
%---------更新上界r-------------
x1=V1*alf1;
r1=fun(x1);
if r1<r
x=x1;
r=r1;
end
x2=V2*alf2;
r2=fun(x2);
if r2<r
x=x2;
r=r2;
end
%--------删除不好的下界-----------------
if u1>=r
V1=[];
end
if u2>=r
V2=[];
end
%---------更新下界u---------
if u1<u2
u=u1;
V=V1;
end
if u2<u1
u=u2;
V=V2;
end
if isempty(V)
u=r;
end
end
k=k+1;
end
%------------------输出结果--------
disp('*****************************************')
disp('计算时间为：')
toc
disp('最优解为：')
x
disp('最优值为：')
r
disp('******************************************')
%-------------函数------------------
function result=fun(x)
sol=x;
result=-abs(sol(1)+1/2*sol(2)+2/3*sol(3))^1.5-sol(1)^2;
%------求一个顶点-------------------------------
function [B,L,x]=ba(A,b)
[m,n]=size(A);
B=1:n;
L=1:(m-n);
a=zeros(n,n);
p=zeros(1,n);
a(1:n,:)=A(B(1:n),:);
p(1:n)=b(B(1:n));
for i=1:m-n+1;
if (rank(a)==n)
mid=a\p';
if A*mid<=b'
x=mid;
L=tt(B,m);
return
end
end
for j=1:n;
for k=i+n:m;
B(j)=k;
a(1:n,:)=A(B(1:n),:);
p(1:n)=b(B(1:n));
if (rank(a)==n)
mid=a\p';
if A*mid<=b';
x=mid;
L=tt(B,m);
return
end
end
end
B(j)=j+i-1;
a(1:n,:)=A(B(1:n),:);
p(1:n)=b(B(1:n));
end
if(i~=m-n+1)
B=B+1;
a(1:n,:)=A(B(1:n),:);
p(1:n)=b(B(1:n));
end
end
if(i==m-n+1&rank(a)<n)
B=[];
disp('请用启发式方法求顶点！');
end
%----------------------------这是求顶点集v0的补集---------------------------
function [L]=tt(B,m)
n=length(B);
L=zeros(1,m-n);
s=1;
t=0;
for i=1:m
t=0;
for j=1:n
t=t+1;
if B(j)==i
t=0;
break;
end
if t==n
L(s)=i;
s=s+1;
end
end
end
%---------函数fj用于求解单纯型的二分----------------
function [V1,V2]=fj(V)
n=size(V,1);
V1=V;
V2=V;
dm=JL(V(:,1),V(:,2));
u=1;
v=2;
for i=1:n
for j=i+1:n+1
d=JL(V(:,i),V(:,j));
if d>dm
dm=d;
u=i;
v=j;
end
end
end
w=0.5*V(:,u)+0.5*V(:,v);
V1(:,u)=w;
V2(:,v)=w;
%-----------距离函数--------------
function [d]=JL(x1,x2)
n=size(x1,1);
d=0;
for i=1:n
d=d+(x1(i)-x2(i))^2;
end
