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
endend
endfor 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
endend
if flag_zhengshu==1
disp('已找到整數最優解:')
flag_xunhuan=0;%改變flag_xunhuan,不再進行迴圈
if d==0
z=-z;
endaaa=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;
endend
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)));
endbb(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);%對增加約束條件後的問題用單純形法求最優解
endendend
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);
endfor i=1:n
z=z+(c(i)*x(i));
endaa=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;
endend
[~,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
ende(i,:)=e(i,:)-a(i,h2)*e(h1,:);%將e進行初等行變換
endb=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;
ende(i,:)= e(i,:)-a(i,k1)*e(k2,:);%將e進行初等行變換
endb=e(1:m,1); %將e拆開為新的a和b
a=e(1:m,2:n+1);
end %至此,有關最優性的判別與修正演算法已結束
endend
舉個例子
matlab命令視窗引數輸入,執行函式,得到最優解
運籌學 5 整數規劃
1 整數線性規劃問題定義 2 0 1變數定義限制條件的表達 使用binary variables 0 1變數 來表示邏輯性條件限制。a.總的只能有某幾個變數被選擇 和小於某個數 b.選擇指定某幾個變數 這幾個變數的和等於某個數 c.幾個之中至少有乙個被選擇,和等於一 d.有順序,只有當某個變數x1被...
南郵運籌學實驗3 01整數規劃
題目 某公司計畫在市區的東 西 南 北四區建立銷售中心,擬議中有10個位置 aj j 1,2,3,10 可供選擇,考慮到各地區居民的消費水平及居民居住密集度,規定 在東區由a1 a2 a3 三個點至多選擇兩個 在西區由a4 a5 兩個點中至少選乙個 在南區由a6 a7 兩個點中至少選乙個 在北區由a...
運籌學 0 1整數規劃(隱列舉法)
來整理一下這個學期運籌學的知識點,路過的朋友鼓勵一下 0 1整數規劃的大體的思路可以用樹形圖來說明,直接上圖 目錄 運籌學 0 1整數規劃 1.轉化為標準型 2.直接令所有變數等於0,看是否滿足所有的約束條件。滿足則結束,否則轉下一步。3.令某個變數為0或1 固定變數 其他變數為自由變數 0或1 看...