fd.m

这是一个用于优化计算的算法,名为分之定界算法.能够优化约束问题.

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