整数规划matlab程序,运筹学整数规划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命令窗口参数输入，执行函数，得到最优解

标签：AA,end,BB,max,整数,flag,matlab,BIndex,运筹学

来源： https://blog.csdn.net/Gin_ger/article/details/106797901