运筹学整数规划matlab程序

Gomory函数

function [X,Z,AAA,BBB] = Gomory(A,B,C,D)
% 割平面法的实现
% X: 目标函数的最优解
% Z: 目标函数的极小值
% AAA:满足整数条件的最终表中的系数矩阵
% BBB:满足整数条件的最终表中的常数列向量
% A: 约束函数的系数矩阵，输入前需保证每个约束条件均为＜或≤号
% B: 约束函数的常数列向量，可以有负值
% C: 目标函数的系数向量
% D:求max为1，min为0
if D==0
    C=-C;
end
[c,d]=size(A);
BIndex=d+1:1:d+c;%给原始的问题添加c个松弛变量，并用BIndex来记录基向量的下标
for i=1:c
    for j=d+1:d+c
        if i==j-d
            A(i,j)=1;
        else
            A(i,j)=0;  %更新系数矩阵A
        end
    end
end
for i=d+1:1:d+c
    C(1,i)=0;   %更新价值系数C
end
[X,Z,AA,BB,BIndex] = DCXF(A,B,C,BIndex);
flag_xunhuan=1;%用作循环开始的条件
while flag_xunhuan
    [a,b]=size(X);
    flag_zhengshu=1;%是否满足整数条件的标志
    for i=1:a
        if abs(round(X(i))-X(i))>=1e-3
            flag_zhengshu=0;%若不满足整数条件，flag_zhengshu跳变为0；若满足整数条件，flag_zhengshu保持为1
            break
        end
    end
    if flag_zhengshu==1
        disp('已找到整数最优解：')
        flag_xunhuan=0;%改变flag_xunhuan，不再进行循环
        if D==0
            Z=-Z;
        end
        AAA=AA;
        BBB=BB;
    else                %若解不为整数，寻找割平面。
         [m,n]=size(AA);
         max_chazhi=0;
         for i=1:m %循环遍历b矩阵，找出和整数相差最大的那个数和它所在的位置
             if BB(i)-floor(BB(i))>max_chazhi
                 max_chazhi=BB(i)-floor(BB(i));
                 max_chazhi_weizhi=i;
             end
         end
         AA=cat(1,AA,zeros(1,n));%以下操作是对矩阵进行扩展，AA矩阵末尾增加一行一列，BB矩阵末尾增加一行
         AA=cat(2,AA,zeros(m+1,1));
         AA(m+1,n+1)=1;
         BB=cat(1,BB,zeros(1,1)); 
         C=cat(2,C,zeros(1,1));
         BIndex=cat(2,BIndex,zeros(1,1));
         for i=1:n  %以下操作是增加新的约束条件（即寻找的割平面）
             AA(m+1,i)=-(AA(max_chazhi_weizhi,i)-floor(AA(max_chazhi_weizhi,i)));
         end
         BB(m+1,1)=-(BB(max_chazhi_weizhi)-floor(BB(max_chazhi_weizhi)));
         BIndex(m+1)=n+1;
         [X,Z,AA,BB,BIndex] = DCXF(AA,BB,C,BIndex);%对增加约束条件后的问题用单纯形法求最优解
    end
end
end

Gomory中调用了DCXF函数

function [X,Z,AA,BB,BIndex] = DCXF(A,B,C,BIndex)
% 单纯形法的实现(包括对偶单纯形法)，需自行添加好松弛变量
% X: 目标函数的最优解
% Z: 目标函数的极小值
% AA:最终表中的系数矩阵
% BB:最终表中的常数列向量
% A: 约束函数的系数矩阵
% B: 约束函数的常数列向量
% C: 目标函数的系数向量
% BIndex: 记录基变量的下标
flag=1;
Z=0;
[m,n]=size(A);
if m>n
    disp('系数矩阵不符合要求,无法进行计算。')
    X=zeros(1,n);
    return 
end
    while flag
        Cb = C(BIndex);             % 基矩阵对应的目标值b
        Zj = (Cb)*A;
        Sigmaj =C-Zj;                   %判别数σ
        X=zeros(n,1);
        if max(Sigmaj)<=0           %若满足最优性
            if min(B)>=0            %若满足可行性
                for i=1:m
                    X(BIndex(i))=B(i);
                end
                for i=1:n
                    Z=Z+(C(i)*X(i));
                end
                AA=A;
                BB=B;
                flag = 0;   %同时满足最优性与可行性的情况下才会进行到这一步，使得flag=0，不再进行下次循环。
            else          %若已满足最优性，但不满足可行性，使用对偶单纯形法
                [~,h1]=min(B);%获得换出基变量的位置
                for i=1:n
                    if A(h1,i)<0
                        P(i)=Sigmaj(i)/A(h1,i);
                    else
                        P(i)=inf;
                    end
                end
                [~,h2]=min(P);%确定换入变量的位置
                BIndex(h1)=h2;%新的基变量
                E=[B,A];           %B,A合成一个新的矩阵
                E(h1,:)=E(h1,:)/A(h1,h2);%将E中h1行，除以主元素（即A中h1行，h2列的那个数）；
                for i=1:m                   %更新E
                    if(i==h1)
                        continue
                    end
                    E(i,:)=E(i,:)-A(i,h2)*E(h1,:);%将E进行初等行变换
                end
                B=E(1:m,1); %将E拆开为新的A和B
                A=E(1:m,2:n+1);
            end %至此，有关可行性的判别与修正算法已结束
            else                    %若不满足最优性，继续进行迭代
                 [~, k1] = max(Sigmaj); %获得换入基变量的位置
                 for i=1:m
                     if A(i,k1)>0
                        Q(i)=B(i)/A(i,k1);
                     else 
                        Q(i)=inf;%
                     end   
                 end
                 [~, k2]=min(Q);    %获得换出基变量的位置
                 BIndex(k2)=k1;     %新的基变量
                 E=[B,A];           %B,A合成一个新的矩阵
                 E(k2,:)=E(k2,:)/A(k2,k1); %将E中k2行，除以主元素（即A中k2行，k1列的那个数）；
                 for i=1:m                   %更新E
                     if(i==k2)
                         continue;
                     end
                     E(i,:)= E(i,:)-A(i,k1)*E(k2,:);%将E进行初等行变换
                 end
                 B=E(1:m,1); %将E拆开为新的A和B
                 A=E(1:m,2:n+1);
        end %至此，有关最优性的判别与修正算法已结束
    end
end

举个例子

matlab命令窗口参数输入，执行函数，得到最优解