置換流水車間調(diào)問題的MATLAB求解

1、物流運(yùn)籌實(shí)務(wù)課程設(shè)計(jì)題目：置換流水車間調(diào)度問題的MATLAB求解置換流水車間調(diào)度問題的MATLAB求解目錄一、 前言 5二、 問題描述 6三、 算法設(shè)計(jì) 7四、 實(shí)驗(yàn)結(jié)果 15摘要自從Johnson 1954年發(fā)表第一篇關(guān)于流水車間調(diào)度問題的文章以來流水車間調(diào)度問題引起了許多學(xué)者的關(guān)注。安排合理有效的生產(chǎn)調(diào)度是生產(chǎn)活動能井然有序開展，生產(chǎn)資源得到最佳配置，運(yùn)作過程簡明流暢的有力保證。流水車間調(diào)度問題是許多實(shí)際流水線生產(chǎn)調(diào)度問題的簡化模型。它無論是在離散制造工業(yè)還是在流程工業(yè)中都具有廣泛的應(yīng)用。因此，對進(jìn)行研究具有重要的理論意義和工程價(jià)值。流水線調(diào)度問題中一個非常典型的問題，而置換流水線調(diào)度問題

2、作為FSP問題的子問題，是一個著名的組合優(yōu)化問題。該問題是一個典型的NP難問題，也是生產(chǎn)管理的核心內(nèi)容。隨著生產(chǎn)規(guī)模的擴(kuò)大，流水線調(diào)度問題的優(yōu)化對提高資源利用率的作用越來越大，因此對其研究具有重要的理論和現(xiàn)實(shí)意義。關(guān)鍵字：流水車間，單件小批量生產(chǎn)，jsp模型，Matlab前言企業(yè)資源的合理配置和優(yōu)化利用很大程度上體現(xiàn)在車間一層的生產(chǎn)活動中,所以加強(qiáng)車間層的生產(chǎn)計(jì)劃與控制一直在企業(yè)生產(chǎn)經(jīng)營活動中占有十分重要的地位。車間 生產(chǎn)計(jì)劃與控制的核心理論是調(diào)度理論。車間調(diào)度問題是一類重要的組合優(yōu)化問題。為適應(yīng)訂貨式、多品種、小批量生產(chǎn)的需要,引進(jìn)了置換流水車間調(diào)度概念。在置換流水車間調(diào)度優(yōu)化后,可以避免或

3、大大減少流程工作時間、提高生產(chǎn)效率。因此,研究成組技術(shù)下車間調(diào)度問題是很有必要的。生產(chǎn)調(diào)度，即對生產(chǎn)過程進(jìn)行作業(yè)計(jì)劃，是整個個先進(jìn)生產(chǎn)制造系統(tǒng)實(shí)現(xiàn)管理技術(shù)、優(yōu)化技術(shù)、白動化與計(jì)算機(jī)技術(shù)發(fā)展的核心。置換流水車間調(diào)度問題是許多實(shí)際生產(chǎn)調(diào)度問題的簡化模型。生產(chǎn)計(jì)劃與調(diào)度直接關(guān)系著企業(yè)的產(chǎn)出效率和生產(chǎn)成本，有效的計(jì)劃與調(diào)度算法能最大限度地提高企業(yè)的效益。調(diào)度問題是組合優(yōu)化問題，屬于NP問題，難以用常規(guī)力一法求解。隨著制造業(yè)的快速發(fā)展，大規(guī)模定制生產(chǎn)、全球化制造等思想的提出，使車間調(diào)度問題呈現(xiàn)出以下的新特點(diǎn):約束條件多，時間復(fù)雜度高，空問復(fù)雜度高。這將導(dǎo)致在許多情況下，求解所建立的數(shù)學(xué)模型的快速性無法滿

4、足，如果采用適度線形化處理之后求解，將會因簡化太多而使結(jié)果嚴(yán)承失真。所以需選擇功能強(qiáng)大的數(shù)值計(jì)算工具來實(shí)現(xiàn)這一問題的求解。MATLAB恰好提供了這樣的平臺。MATLAB是一個高度集成的系統(tǒng)，集科學(xué)計(jì)算、圖像處理、聲音處理于一體，具有極高的編程效率。典型JSP模型分析與Matlab的應(yīng)用結(jié)合使流水車間調(diào)度問題迎刃而解。最大完工時間是生產(chǎn)調(diào)度中最常用的性能度量指標(biāo)之一，最大完工時間越短，則說明產(chǎn)品總的生產(chǎn)周期越短，生產(chǎn)能力越大；此類調(diào)度問題的優(yōu)化研究有助于提高企業(yè)的生產(chǎn)效率與資源利用率。一、問題描述流水車間調(diào)度問題通?？梢悦枋鰹閭€工件要在臺機(jī)器上加工，每個工件有道工序，每道工序都要在不同的機(jī)器上加

5、工，所有工件的加工順序都相同，問題的目標(biāo)是確定每臺機(jī)器上工件的加工順序及開工時間，使得特定的性能指標(biāo)最優(yōu)。置換流水車間調(diào)度問題PFSP是對流水車間調(diào)度問題的進(jìn)一步約束，即約定每臺機(jī)器上所有工件的加工順序相同，其解空間的規(guī)模為，遠(yuǎn)遠(yuǎn)小于流水車間調(diào)度問題的規(guī)模。本次課程實(shí)驗(yàn)主要研究PFSP中的最小化最大完工時間問題，利用三元組表示法（）求解Carlier (1978)提出的8個算例、以及Reeves (1995)提出的21個算。由于三臺機(jī)器以上的調(diào)度問題被證明是NP難問題，對于大規(guī)模的調(diào)度，至今仍未出現(xiàn)求解最優(yōu)的方法，常常采用啟發(fā)式算法來求解近優(yōu)解。本案例主要采用instance car2進(jìn)行求解

6、。案例：某產(chǎn)品，需要經(jīng)過4道工序?qū)?3個工件進(jìn)行加工，這13個工件的生產(chǎn)流程是一樣的。加工時間表見下：表4-3 某產(chǎn)品加工時間表工件號加工時間12345678910111213tj178963021457321865821420778569653212457tj293021425789653214254786532112412345678tj321475320124752147532145763214257854123tj4320142753214528653514527536214528888999計(jì)算步驟如下：首先確定n/m/F/Cmax的最大完工時間為： k=2,.,m i=2,.,n

7、則 Cmax=二、算法設(shè)計(jì)(一)假設(shè)工件在機(jī)器上的加工順序是相同的，同時假定各工件準(zhǔn)備就緒，機(jī)器一開動就投入生產(chǎn)，開工時間為0，則最大完工時間等于最大流程時間。同時3臺機(jī)器以上的流水車間調(diào)度是NP難問題，所以本文只考慮了2臺、3臺機(jī)器的情況，解決3臺機(jī)器以上的問題方法也可運(yùn)用人工智能算法，解的質(zhì)量更高，但因該類算法需良好的軟件編程能力，故本文不加探究。n個工件在m臺機(jī)器上的加工順序相同。工件在機(jī)器上的加工時間是給定的。問題的目標(biāo)是求n個工件在每合機(jī)器上的最大完工時間等于最大流程時間。這種流水線調(diào)度問題要在滿足以下兩個約束條件的前提下，使得加工完所有的工件所花的時間盡可能地少：1、工件約束每個工

8、件在每臺機(jī)器上恰好加工一次，每個工件在各機(jī)器上加工順序相同。不失一般性，假設(shè)各工件按機(jī)器1至m的順序進(jìn)行加工。各工件在各機(jī)器上的加工時間已知。2、機(jī)器約束每臺機(jī)器在任何時刻至多加工一個工件，每臺機(jī)器加工的各工件的順序相同。置換流水線調(diào)度問題實(shí)質(zhì)是如何調(diào)整加工工件的序列，提高機(jī)器的利用率的問題，即在同一時刻正在加工的機(jī)攫數(shù)越多，機(jī)器利用率越大口根據(jù)該原則，我們根據(jù)下面規(guī)則安排工件的加工順序: (l)在前面機(jī)器加工時間較短、后面機(jī)器加工時間較長的工件，安排在序列前。這樣可以使得后面的機(jī)器盡快參加工作，并且后面的機(jī)器不需要作空等待， (2)機(jī)器加工時間較為平均且加工時間較長的工件，安排在序列的中部。

9、這樣可以使得各個機(jī)器在中期的時候都能得到運(yùn)作。 (3前面加工時間較長，后面加一時間較短的上件女排在序列尾部。這樣使得前面的機(jī)器能“延遲”完工，后面的機(jī)器盡快完工。（二）利用Matlab軟件對上面的案例進(jìn)行求解，編程如下:軟件輸出相應(yīng)的結(jié)果，如下：clc; clear all;temp= 0 456 1 856 2 963 3 696 0 789 1 930 2 21 3 320 0 630 1 214 2 475 3 142 0 214 1 257 2 320 3 753 0 573 1 896 2 124 3 214 0 218 1 532 2 752 3 528 0 653 1 142 2

10、 147 3 653 0 214 1 547 2 532 3 214 0 204 1 865 2 145 3 527 0 785 1 321 2 763 3 536 0 696 1 124 2 214 3 214 0 532 1 12 2 257 3 528 0 12 1 345 2 854 3 888 0 457 1 678 2 123 3 999;T=temp(:,2:2:end); n,m=size(T); %n為工件數(shù)，m為機(jī)器數(shù)txm=; %所有m-1個兩臺虛擬機(jī)器問題的加工時間矩陣Cmax=; %所有m-1個兩臺虛擬機(jī)器問題的總完工時間TXY=; %存放m-1個加工順序ticfor

11、 i=1:m-1 for j=1:n tx(j,1)=sum(T(j,1:i); %第一臺虛擬機(jī)器上的加工時間 tx(j,2)=sum(T(j,m+1-i:m); %第二臺虛擬機(jī)器上的加工時間 end cmax,xy=johnson(tx,T,n,m); %調(diào)用Johnson 算法函數(shù) Cmax = Cmax,cmax; txm=txm,tx; TXY=TXY,xy;endtxm, Cmax, TXYoptim=min(Cmax) %近似最優(yōu)總完工時間ind=find(Cmax =optim);optim_seq=TXY(:,ind)' %近似最優(yōu)加工順序rumTime=tocMatl

12、ab運(yùn)行：txm = 456 696 1312 1659 2275 2515 789 320 1719 341 1740 1271 630 142 844 617 1319 831 214 753 471 1073 791 1330 573 214 1469 338 1593 1234 218 528 750 1280 1502 1812 653 653 795 800 942 942 214 214 761 746 1293 1293 204 527 1069 672 1214 1537 785 536 1106 1299 1869 1620 696 214 820 428 1034 552

13、 532 528 544 785 801 797 12 888 357 1742 1211 2087 457 999 1135 1122 1258 1800Cmax = 8679 8423 8667TXY = 13 13 4 9 4 7 4 12 13 8 6 9 6 7 14 1 10 8 14 1 6 7 14 1 10 8 10 12 9 2 2 3 5 5 11 3 11 2 12 3 5 11optim = 8423optim_seq = 13 4 12 6 7 10 1 14 8 9 3 11 2 5（三）繪制甘特圖編程如下：% 已知nbjobs, nbmachines, P% 通

14、過computeCmax計(jì)算完工時間矩陣，然后畫甘特圖temp= 0 456 1 856 2 963 3 696 0 789 1 930 2 21 3 320 0 630 1 214 2 475 3 142 0 214 1 257 2 320 3 753 0 573 1 896 2 124 3 214 0 218 1 532 2 752 3 528 0 653 1 142 2 147 3 653 0 214 1 547 2 532 3 214 0 204 1 865 2 145 3 527 0 785 1 321 2 763 3 536 0 696 1 124 2 214 3 214 0 53

15、2 1 12 2 257 3 528 0 12 1 345 2 854 3 888 0 457 1 678 2 123 3 999;P=temp(:,2:2:end);nbjobs=14;nbmachines=4;% 給定工件的加工順序，計(jì)算完工時間，進(jìn)而求makespanjobOrder= 13 4 12 6 7 10 1 14 8 9 3 11 2 5;C=zeros(nbjobs,nbmachines); % 完工時間矩陣% 第一道工序，工件的完工時間C(jobOrder(1),1)=P(jobOrder(1),1); % 第一個工件的完工時間for i=2:nbjobs % 在加工順序

16、里面的位置 % 第一道工序上，上一個工件加工完下一個工件才可以加工 C(jobOrder(i),1)=C(jobOrder(i-1),1)+P(jobOrder(i),1);end% 其余工序，工件的完工時間for j=2:nbmachines C(jobOrder(1),j)=C(jobOrder(1),j-1)+P(jobOrder(1),j); for i=2:nbjobs % 在上道工序加工完，且前一個工件加工完，才可以加工 C(jobOrder(i),j)=max(C(jobOrder(i),j-1),C(jobOrder(i-1),j)+P(jobOrder(i),j); ende

17、ndCmax=max(C(:,nbmachines); % Cmaxfigure_handle=figure;hold on;% Create titletitle('置換流水車間調(diào)度甘特圖（Cmax= ',num2str(Cmax),'）',. 'FontWeight','bold','FontSize',12);xlabel('時 間');ymax=nbmachines+1;xmax=Cmax;axis(0 xmax 0 ymax);%-YLabel=cell(ymax,1);y=ymax;fo

18、r j=1:nbmachines YLabely='M' int2str(j); y=y-1;endset(gca,'Ygrid','on','YTickLabel',YLabel,'YTick',0:ymax);%hi=0.7;nb=1; % machineyi=ymax-1;for j=1:nbmachines for i=1:nbjobs xi=C(jobOrder(i), j)-P(jobOrder(i),j); wi=P(jobOrder(i), j); rectangle('Position&#

19、39;,xi,yi,wi,hi); text(xi+0.5*wi,yi+0.5*hi,int2str(jobOrder(i),'HorizontalAlignment',. 'center'); end yi=yi-1; nb=nb+1;endyi=yi-1;hold off;三、實(shí)驗(yàn)結(jié)果1、根據(jù)上面matlab的求解得到以下實(shí)驗(yàn)結(jié)果：最優(yōu)排序?yàn)?3 4 12 6 7 10 1 14 8 9 3 11 25min(Cmax)= 84232、甘特圖如下所示：四、流水線型車間作業(yè)調(diào)度問題遺傳算法MATLAB源碼流水線型車間作業(yè)調(diào)度問題可以描述如下：n個任務(wù)在流水線上

20、進(jìn)行m個階段的加工，每一階段至少有一臺機(jī)器且至少有一個階段存在多臺機(jī)器，并且同一階段上各機(jī)器的處理性能相同，在每一階段各任務(wù)均要完成一道工序，各任務(wù)的每道工序可以在相應(yīng)階段上的任意一臺機(jī)器上加工，已知任務(wù)各道工序的處理時間，要求確定所有任務(wù)的排序以及每一階段上機(jī)器的分配情況，使得調(diào)度指標(biāo)(一般求Makespan)最小。 下面的源碼是求解流水線型車間作業(yè)調(diào)度問題的遺傳算法通用MATLAB源碼function Zp,Y1p,Y2p,Y3p,Xp,LC1,LC2=JSPGA(M,N,Pm,T,P)%  輸入?yún)?shù)列表%  M    

21、   遺傳進(jìn)化迭代次數(shù)%  N       種群規(guī)模(取偶數(shù))%  Pm      變異概率%  T       m×n的矩陣，存儲m個工件n個工序的加工時間%  P       1×n的向量，n個工序中，每一個工序所具有的機(jī)床數(shù)目%  輸出參數(shù)列表%  Zp   &#

22、160;  最優(yōu)的Makespan值%  Y1p     最優(yōu)方案中，各工件各工序的開始時刻，可根據(jù)它繪出甘特圖%  Y2p     最優(yōu)方案中，各工件各工序的結(jié)束時刻，可根據(jù)它繪出甘特圖%  Y3p     最優(yōu)方案中，各工件各工序使用的機(jī)器編號%  Xp      最優(yōu)決策變量的值，決策變量是一個實(shí)數(shù)編碼的m×n矩陣% &#

23、160;LC1     收斂曲線1，各代最優(yōu)個體適應(yīng)值的記錄%  LC2     收斂曲線2，各代群體平均適應(yīng)值的記錄%  最后，程序還將繪出三副圖片：兩條收斂曲線圖和甘特圖（各工件的調(diào)度時序圖）%第一步：變量初始化m,n=size(T);%m是總工件數(shù)，n是總工序數(shù)Xp=zeros(m,n);%最優(yōu)決策變量LC1=zeros(1,M);%收斂曲線1LC2=zeros(1,N);%收斂曲線2%第二步：隨機(jī)產(chǎn)生初始種群farm=cell(1,N);%采用細(xì)胞結(jié)構(gòu)存儲種群for

24、k=1:N    X=zeros(m,n);    for j=1:n        for i=1:m            X(i,j)=1+(P(j)-eps)*rand;        end    end    farmk=X;endcounter=0;%設(shè)置迭代計(jì)數(shù)器while counter<M%停

25、止條件為達(dá)到最大迭代次數(shù)        %第三步：交叉    newfarm=cell(1,N);%交叉產(chǎn)生的新種群存在其中    Ser=randperm(N);    for i=1:2:(N-1)        A=farmSer(i);%父代個體        Manner=unidrnd(2);%隨機(jī)選擇交叉方式      &#

26、160; if Manner=1            cp=unidrnd(m-1);%隨機(jī)選擇交叉點(diǎn)            %雙親雙子單點(diǎn)交叉            a=A(1:cp,:);B(cp+1):m,:);%子代個體            b=B(1:cp

27、,:);A(cp+1):m,:);        else            cp=unidrnd(n-1);%隨機(jī)選擇交叉點(diǎn)            b=B(:,1:cp),A(:,(cp+1):n);        end        newfarmi

28、=a;%交叉后的子代存入newfarm        newfarmi+1=b;    end    %新舊種群合并    FARM=farm,newfarm;        %第四步：選擇復(fù)制    FITNESS=zeros(1,2*N);    fitness=zeros(1,N);    plotif=0;    for i=1:(2*N) 

29、;       X=FARMi;        Z=COST(X,T,P,plotif);%調(diào)用計(jì)算費(fèi)用的子函數(shù)        FITNESS(i)=Z;    end    %選擇復(fù)制采取兩兩隨機(jī)配對競爭的方式，具有保留最優(yōu)個體的能力    Ser=randperm(2*N);    for i=1:N     &#

30、160;  f2=FITNESS(Ser(2*i);        if f1<=f2            farmi=FARMSer(2*i-1);            fitness(i)=FITNESS(Ser(2*i-1);        else   

31、0;        farmi=FARMSer(2*i);        end    end    %記錄最佳個體和收斂曲線    minfitness=min(fitness)    meanfitness=mean(fitness)    LC1(counter+1)=minfitness;%收斂曲線1，各代最優(yōu)個體適應(yīng)值的記錄    LC2(counter+1

32、)=meanfitness;%收斂曲線2，各代群體平均適應(yīng)值的記錄    pos=find(fitness=minfitness);    Xp=farmpos(1);        %第五步：變異    for i=1:N        if Pm>rand;%變異概率為Pm            X=farmi;    

33、;        I=unidrnd(m);            J=unidrnd(n);            X(I,J)=1+(P(J)-eps)*rand;            farmi=X;        end 

34、60; end    farmpos(1)=Xp;        counter=counter+1end%輸出結(jié)果并繪圖figure(1);plotif=1;X=Xp;Zp,Y1p,Y2p,Y3p=COST(X,T,P,plotif);figure(2);plot(LC1);figure(3);plot(LC2);function Zp,Y1p,Y2p,Y3p=COST(X,T,P,plotif)%  JSPGA的內(nèi)聯(lián)子函數(shù)，用于求調(diào)度方案的Makespan值%  輸入?yún)?shù)列表% &

35、#160;X       調(diào)度方案的編碼矩陣，是一個實(shí)數(shù)編碼的m×n矩陣%  T       m×n的矩陣，存儲m個工件n個工序的加工時間%  P       1×n的向量，n個工序中，每一個工序所具有的機(jī)床數(shù)目%  plotif  是否繪甘特圖的控制參數(shù)%  輸出參數(shù)列表%  Zp    

36、60;最優(yōu)的Makespan值%  Y1p     最優(yōu)方案中，各工件各工序的開始時刻%  Y2p     最優(yōu)方案中，各工件各工序的結(jié)束時刻%  Y3p     最優(yōu)方案中，各工件各工序使用的機(jī)器編號%第一步：變量初始化m,n=size(X);Y1p=zeros(m,n);Y2p=zeros(m,n);Y3p=zeros(m,n);%第二步：計(jì)算第一道工序的安排Q1=zeros(m,1);Q2=zeros(m,1)

37、;R=X(:,1);%取出第一道工序Q3=floor(R);%向下取整即得到各工件在第一道工序使用的機(jī)器的編號%下面計(jì)算各工件第一道工序的開始時刻和結(jié)束時刻for i=1:P(1)%取出機(jī)器編號    pos=find(Q3=i);%取出使用編號為i的機(jī)器為其加工的工件的編號    lenpos=length(pos);    if lenpos>=1        Q1(pos(1)=0;        if

38、 lenpos>=2            for j=2:lenpos                Q1(pos(j)=Q2(pos(j-1);                Q2(pos(j)=Q2(pos(j-1)+T(pos(j),1);       

39、    end        end    endendY1p(:,1)=Q1;Y3p(:,1)=Q3;%第三步：計(jì)算剩余工序的安排for k=2:n    R=X(:,k);%取出第k道工序    Q3=floor(R);%向下取整即得到各工件在第k道工序使用的機(jī)器的編號    %下面計(jì)算各工件第k道工序的開始時刻和結(jié)束時刻    for i=1:P(k)%取出機(jī)器編號    

40、;    pos=find(Q3=i);%取出使用編號為i的機(jī)器為其加工的工件的編號        lenpos=length(pos);        if lenpos>=1            EndTime=Y2p(pos,k-1);%取出這些機(jī)器在上一個工序中的結(jié)束時刻         &#

41、160;  POS=zeros(1,lenpos);%上一個工序完成時間由早到晚的排序            for jj=1:lenpos                POS(jj)=ppp(1);                EndTime(ppp(1)=Inf;    

42、        end                        %根據(jù)上一個工序完成時刻的早晚，計(jì)算各工件第k道工序的開始時刻和結(jié)束時刻            Q1(pos(POS(1)=Y2p(pos(POS(1),k-1);         

43、;   Q2(pos(POS(1)=Q1(pos(POS(1)+T(pos(POS(1),k);%前一個工件的結(jié)束時刻            if lenpos>=2                for j=2:lenpos                   

44、 Q1(pos(POS(j)=Y2p(pos(POS(j),k-1);%預(yù)定的開始時刻為上一個工序的結(jié)束時刻                    if Q1(pos(POS(j)<Q2(pos(POS(j-1)%如果比前面的工件的結(jié)束時刻還早                       &#

45、160;Q1(pos(POS(j)=Q2(pos(POS(j-1);                    end                end            end        end  &

46、#160; end    Y1p(:,k)=Q1;    Y2p(:,k)=Q2;    Y3p(:,k)=Q3;end%第四步：計(jì)算最優(yōu)的Makespan值Y2m=Y2p(:,n);Zp=max(Y2m);%第五步：繪甘特圖if plotif    for i=1:m        for j=1:n            mPoint1=Y1p(i,j);            mPoint2=Y2p(i,j);