帶時間窗物流配送車輛路徑問題摘要本題是一個帶有時間窗的車輛路徑安排問題（VRPTW問題）。根據(jù)題目條件，本文建立了一個求解最小派送費用的VRPTW優(yōu)化模型，采用遺傳算法，給出了該模型的求解方法。然后，對一個實際問題進行求解，給出了一個比較好的路線安排方式。模型一（見5.1.2）針對問題一，在需求量、接貨時間段、各種費用消耗已知的情況下，決定采用規(guī)劃模型，引入0-1變量，建立各個約束條件，包括車輛的容量限制，到達每個客戶的車輛和離開每個客戶的車輛均為1的限制，總車輛數(shù)的限制，目標函數(shù)為費用的最小化，費用包括車輛的行駛費用，車輛早到或晚到造成的損失。模型一的求解采用遺傳算法（見5.1.3），對題目給出的實際問題進行求解，得到3條行車路線，總路線長度為910公里，具體結果如下：車輛編號所執(zhí)行的任務路線到達各點的時間路線長度貨運量10-3-1-2-00—1.5—3.3—5.6—8.875+40+65+60=2404.5+2+1.5=820-8-5-7-00-1.6-3.9-7.7-13.980+75+90+160=4053+1.5+2.5=730-6-4-00-2-6-10.8100+75+90=2654+3=7考慮在車輛返回時選擇最短路線，我們采用Dijkstra算法求出中心倉庫到各個客戶的最短距離，將結果改進為885公里,具體結果如下：車輛編號所執(zhí)行的任務路線到達各點的時間路線長度貨運量10-3-1-2-00—1.5—3.3—5.6—8.875+40+65+60=2404.5+2+1.5二820-8-5-7-2-00-1.6-3.9-7.7-13.980+75+90+135=3803+1.5+2.5二730-6-4-00-2-6-10.8100+75+90=2654+3=7模型二考慮需求量隨機變化時的安排方案，假設客戶需求量遵循正態(tài)分布，首先按照需求期望根據(jù)模型一得到一個比較好的方案，然后按照這一方案進行送貨，在送貨過程中，如果出現(xiàn)需求量過大的情況，允許車輛返回倉庫進行補充。模型一的思路清晰，考慮條件全面。但最優(yōu)解解決起來困難，遺傳算法只是一種相對好的解決方法，可以找出最優(yōu)解的近似解。模型二的想法比較合理，易于實施，但還有待改進。關鍵詞：規(guī)劃時間窗物流車輛路徑遺傳算法一、 問題重述一個中心倉庫，擁有一定數(shù)量容量為Q的車輛，負責對N個客戶進行貨物派送工作，客戶i的貨物需求量為q，且q<Q，車輛必須在一定的時間范圍ii［a,b］內(nèi)到達，早于a到達將產(chǎn)生等待損失，遲于b到達將處以一定的懲罰，請iiii解決如下問題：（1）給出使派送費用最小的車輛行駛路徑問題的數(shù)學模型及其求解算法。并具體求解以下算例：客戶總數(shù)N=8,每輛車的容量Q=8（噸/輛），各項任務的貨運量q（單位：i噸）、裝貨（或卸貨）時間s（單位：小時）以及要求每項任務開始執(zhí)行的時間i范圍［a,b］由附錄1給出，車場0與各任務點以及各任務點間的距離（單位：公ii里）由附件二給出，這里假設車輛的行駛時間與距離成正比，每輛車的平均行駛速度為50公里/小時，問如何安排車輛的行駛路線使總運行距離最短；（2）進一步請討論當客戶i的貨物需求量q為隨機參數(shù)時的數(shù)學模型及處理方i法。二、 問題分析本題主要在兩種不同情況下，研究使派送費用最小的車輛行駛路徑問題。車輛行駛派送的費用主要包括運輸成本、車輛在客戶要求到達時間之前到達產(chǎn)生的等待損失和車輛在客戶要求到達時間之后到達所受懲罰等等。為滿足派送費用最小的需求，即要使所選行車路徑產(chǎn)生的總費用最小，從而確定出最佳的車輛派送當客戶i的貨物需求量q固定時，首先，我們根據(jù)題意，取若干輛車進行送i貨，然后，主要考慮每輛車各負責哪些客戶的送貨任務，我們可以給出滿足題中限制條件的很多參考方案供選用，并考慮以所選行車路徑產(chǎn)生的總費用最小為目標的情況下，建立最優(yōu)化模型確定最佳的車輛派送方案。進一步討論，當客戶i的貨物需求量q為隨機參數(shù)時，我們首先可以簡化隨i機模型，根據(jù)客戶i的貨物需求量的期望與方差，確定每天應該運送給客戶i的貨物量，即q，再根據(jù)第一題，確定最佳的車輛派送方案。i但考慮到客戶的儲存能力有限及貨物在客戶處的儲存費用，客戶不需要將一天的貨物一次性接收完，只要滿足缺貨的情況出現(xiàn)的概率很低，客戶可以讓配送中心一天幾次送貨，這樣可以得到很多滿足約束的方案，考慮以單位時間的儲存費用最小為目標，建立最優(yōu)化模型，確定配送中心給每位客戶每次的配送量、配送周期與最有車輛行駛路徑。三、 模型假設（1） 每個客戶的需求只能由一輛配送車滿足；（2） 每輛車送貨時行駛的路程不超過它所能行駛的最遠路程；（3） 中心倉庫的車輛總數(shù)大于或等于當派送費用最小時所需的車輛數(shù)；（4） 從配送中心到各個用戶、各個用戶之間的運輸距離已知；（5） 配送中心有足夠的資源以供配送。四、 符號說明Q:運貨車的容量N:該配送中心服務的客戶總數(shù)m:派送費用最小時所需的車輛數(shù)q:第i位客戶所需貨物ix:車k由i駛向jy：點i的貨運任務友s車完成isa:第i位客戶最早允許接貨時間ib:第i位客戶最晚允許接貨時間iD:車輛在第i位客戶處等待時間iX:車輛在第i位客戶處遲到時間it:第i位客戶處車輛到達時間it..:從第i位客戶到第j位客戶所需時1]間s:第i位客戶處裝貨（或卸貨）所需i時間c：第i位客戶與第j位客戶之間的距離c:車輛行駛單位距離的運輸成本d:車輛早到單位時間產(chǎn)生的等待損失e:車輛遲到單位時間應承擔的懲罰Z:派送貨物產(chǎn)生的總損失A：運輸成本B：車輛早到所產(chǎn)生總的等待損失C：車輛遲到所受的總懲罰五、 模型的建立和求解問題一模型的建立及求解：5.1.1問題的分析中心倉庫為了給N個客戶派送貨物，供發(fā)出m輛車，為了派貨的節(jié)約和方便，每輛車載著適量的貨物出發(fā)，可以給某一片的若干個滿足約束條件的客戶派送貨物，見圖一：

圖一中心倉庫派送貨物圖中心倉如上圖庫派送貨物時，必須滿足約束條件：（1） 各個客戶群的總需求小于或等于運輸車的裝載量；（2） 每個客戶都必須且只能由一輛運輸車運輸所需貨物；（3） 運輸車為每位客戶開始服務的時間必須盡可能在時間窗內(nèi)。根據(jù)如上的約束條件，我們可以得到很多可行解，但考慮到以所選行車路徑產(chǎn)生的總費用最小為目標的情況下，我們可以建立最優(yōu)化模型確定最佳的車輛派送方案，最優(yōu)路徑產(chǎn)生圖如下：

圖二最優(yōu)路徑產(chǎn)生圖5.1.2模型的建立（1）中心倉庫使用車輛數(shù)量的確定設配送中心需要向N個客戶送貨，每個客戶的貨物需求量是gi（i=l,2,…?.N）,每輛配送車的載重量是Q,且gi〈Q。首先為了安排路線需要對要使用的車輛數(shù)有一個估計。在現(xiàn)實情況中,貨物裝（卸）車越復雜,約束條件越多,一輛車的實際載貨量就越小。在本文中使用文獻［1］的公式來確定需要的車輛數(shù)m：m=￡g/aQii=l［］表示取整，a為參數(shù)，0〈a〈1。約束條件越多，貨物裝（卸）越復雜，a值越小。參考文獻［2］,取a為0?85。2）引入0—1變量：1）xjs表示車輛s是否從客戶'行駛到客戶j。定義其為0—1變量’則_J1車輛S從客戶i行駛到客戶j 1Xijs_[0—‘ s=1,,m否則 sm2）兒表示客戶i的任務由車輛s完成。同樣定義其為0—1變量，則is

】 客戶i的任務由車輛s完成y=<%[0 否則非線性規(guī)劃模型的建立： …a?目標函數(shù)的確定。題目要求所選行車路徑產(chǎn)生的總費用最小，我們確定總費用為目標函數(shù)，記為Z。總費用由運輸成本A、等待損失B和遲到所收懲罰C組成，根據(jù)題意有：A二c為H遲cxijijsi=1j=1s=1B=為d*max{D,0}ii=1C=Xe*max{X,0}ii=1所以，總費用Z最小化為：minZcXNXNminZcXNXNXmcxijijsi=0 j=0 s=1+Xd*max{D,0}+Xe*max{X,0}iii=1i=1b?約束條件的確定。約束1：車輛k的運送總重量應不超過車輛的最大載重，即車輛有一定的運送能力,則可引入約束條件，Xqy<Q (Vkg{1,2,,m})iisi=1約束2：每個客戶只能由一輛車來配送，則可引入約束條件，申 \1 i=1,2,3,???NXy= ] .nis [m i=0s=1約束3：保證到達一個客戶的車輛也離開該客戶，則可引入約束條件，mN(j=1,2,3,,N;)XXx=1(j=1,2,3,,N;)ijss=1i=1

ijks=1j=1(ijks=1j=1(i=1,2,3, ,N;)C?變量之間關系的確定由上可確定等待時間D，超時時間X為iiD=a—t i=1,2,....NiiiX=t—b i=1,2, Niii車輛k從客戶i到客戶j需經(jīng)過兩段時間t為：ijt=c/v i,j=1,2,Nijij設車輛為客戶i運送完貨物后即為客戶j運送，則到達客戶i處時間t和到達i客戶j處時間t之間的關系為：jtjxtjx(t+max{D,0}+1+s)ijsi iijii=0 j=0s=1i=0 j=0s=1i=1i=1為qy<Qiisi=1y=iss=1s=1,2, ,mi=1,2,3,...Ni=0mN審X=1ijss=1i=1江X=1ijks=1j=1D=a—tiiiX=t—biiit=c/vijijj=1,2,3, ,N;i=123,…,N;i=1,2, Ni=1,2,....Nj=1,2,Nd此非線性規(guī)劃模型為:ijijsminZ=正工區(qū)cx+為d*max{D,0}+He*max{X,0}ijijsiit=x(t+max{D,0}+1+s)j ijsi iijii=1s=15.1.3模型的求解我們采用遺傳算法解決上面的問題：1．編碼采用自然數(shù)編碼，即序數(shù)編碼。貨物運輸路線可以編成長度為N+m的染色體(o,i,i,,i,0,i,i,0,,0,i,,i)，其中，i表示第i項任務。01112 1s21 2t m1 mw ik ik表示車場，m表示完成任務所需的車輛數(shù)。出生初始群體初始群體隨機產(chǎn)生，?即產(chǎn)生N項貨物?運輸任務點的全排列，如i，i，，i,12N如果藝q..§Q， >Q，將S至N的數(shù)向后移動一位，將0插入第s位。ij ijj=1 j=1 ?…接著，繼續(xù)上述操作，直到m個0全部插入為止。這樣就構成了一條初始染色體。用這種方法構造一個群體的染色體。如：82576314，該編碼插零之后變成08250763014。它代表著需要三輛車運輸貨物。其中，第1輛車行走路線為08250，即從倉庫出發(fā)到依次到8、2、5商店再回到倉庫。第2輛車行走路線為07630，第3輛車行走路線為0140。適應度函數(shù)bz'適應度函數(shù)取f= ，其中f為染色體v的適應度，b為常數(shù)，z'為初kzkkk始種群中最好的染色體的運輸成本，Z為染色體V對應的運輸成本。kk遺傳算子選取最佳保留的輪盤賭復制法進行染色體的復制。變異算子采用反轉(zhuǎn)變異。交叉算子用最大保留交叉，其操作過程為：若染色體交叉點處的兩個基因都為0，則直接進行順序交叉運算；若染色體交叉點處的基因不全為0，則將交叉點左移(右移)，直到左右兩個交叉處的基因都為0，再進行順序交叉運算。算法的實現(xiàn)步驟Stepl：采用自然數(shù)編碼的方式，構造表示可行車路線的染色體；Step2：設置控制參數(shù)，包括交叉率p=0.7、變異率p=0.1、群體規(guī)模n=10；cmStep3：初始化，令d=0，隨機產(chǎn)生初始群體P(0)，群體中包括n個染色體，每個染色體代表一條行車線路；Step4:令i—1；Step5：將群體p(d)中的第i個染色體譯為線路長度；Step6:計算適應度；Step7:若滿足算法終止條件，則停止，否則繼續(xù)；Step8：i=i+1；Step9：若i—n，回到step5，否則，轉(zhuǎn)steplO；Stepll:進行最大保留交叉、基于位的變異和倒位操作;Step12:d—d+1；Step13:若滿足算法終止條件，則停止，否則轉(zhuǎn)step4。運用matlab軟件編寫程序得到在車輛總行程最短的條件小的派送方案為:車輛編號所執(zhí)行的任務路線到達各點的時間路線長度貨運量10-3-1-2-00—1.5—3.3—5.6—8.875+40+65+60=2404.5+2+1.5二820-8-5-7-00—1.6—3.9—7.7—13.980+75+90+160=4053+1.5+2.5二730-6-4-00-2-6-10.8100+75+90=2654+3=7從上表可以看出，按照上表的行車路線，總路線長度為910公里。5.1.4結果分析從上面解出的派送方案可以看出，上面的每條路線在車輛送完貨物后，直接從最后一個客戶處返回中心倉庫，我們通過求中心倉庫到各個客戶的最短路徑可以發(fā)現(xiàn)，上面的返回路線不是最優(yōu)的，返回路線可以經(jīng)過某些客戶使得所走路線最短。我們用Dijkstra算法求出中心倉庫到各個客戶的最短路徑如下:編號012345678最短距離0406075909010013580根據(jù)上表，我們對5.1.2所求的結果進行改進，得到下表:車輛編號所執(zhí)行的任務路線到達各點的時間路線長度貨運量10—3—1—2—00—1.5—3.3—5.6—8.875+40+65+60=2404.5+2+1.5=820—8—5—7—2—00—1.6—3.9—7.7—13.980+75+90+135=3803+1.5+2.5=730—6—4—00—2—6—10.8100+75+90=2654+3=7根據(jù)上表計算結果，我們在不改變原路線的情況下，使總路線減小為885問題二模型的建立及求解：5.2.1問題的分析與假設在問題一中，根據(jù)已知的各個商店的需求分布，根據(jù)遺傳算法求解，預先確定一條總費用最小的路徑（或者雖不是最優(yōu)路徑，但是此結果能夠接受的）。車輛沿該路徑服務商店，因此服務商店的順序固定不變。而問題二中，在車輛服務商店的過程中，事先并不知道商店的需求量。而這個不確定信息要隨著車輛的服務逐步確定，商店具體的需求只有車輛到達用戶后才確定。這樣第一問的求解方法得到的路徑并不適用于第二問的求解。假設：（1） 商店的需求變化符合正態(tài)分布，第i個商店的需求量的期望和方差分別為卩i和Q2。i（2） 商店的供貨可以分為多次補充但在每次供貨中最大程度滿足用戶需求，即只有出現(xiàn)當前車載余量小于用戶需求量時才出現(xiàn)下一次的供貨。模型的建立基于第一問解決了在已知用戶需求概率情況下，確定一個服務方案，滿足所有n個商店的需求并且使車輛期望行程費用最小這個問題。我們由假設可知，第i個商店的需求量的期望為卩，則由需求量的期望得到一個使車輛期望行程費用i最小的服務方案，稱該方案為A。當用戶的需求未知（當車輛服務商店時需求變成已知量）時，由于用戶的需求量在區(qū)間（卩土3d）的概率是96%,而在區(qū)間外的事件可以看成小概率事件，由ii小概率事件定理可知，在一次試驗中，小概率事件可以看成不可能事件。由此可知，用戶需求量就在區(qū)間（卩土3d）里。ii用戶需求在區(qū)間（卩土3d）里，而需求決定服務方案。由上面可知，服務方ii案在A方案附近變化，而變化的幅度由方差d2決定。當d2越小時（說明需求量ii接近一個常數(shù)期望卩），最優(yōu)方案（或滿意方案）與A方案越接近（即在A方案i上面稍作改動即可）。反之，A方案需作較大的方案。由假設（2）知，車輛只要滿足當前用戶部分需求，就服務該用戶，用戶未滿足的部分以后再服務。在服務用戶后，車輛根據(jù)當前的位置信息、車載余量以及需要服務用戶的信息，決定下一個服務用戶（包括當前還未滿足需求的用戶）或回庫房裝載貨物。先按A方案進行分配車輛路線（若增加車輛數(shù)目雖可以更好滿足條件，但會增加多于開支）。假設m輛車都從倉庫0出發(fā)，按A方案中的預定路線進行服務，當服務完第一個商店時，再判斷剩余的貨物量。此時貨物量為Q'二Q-q（其i中Q為貨車滿載量，q為車輛到達商店i時確定了i個商店的需求量）。然后判斷i該剩余量是否在大于下一個商店需求量（卩+3q），若大于則進行下一個商店服務，若小于則判斷下個商店是否滿足大ii于（卩+3。）這個條件，若滿足，則對其進行服務，否則繼續(xù)下個判斷，直至其ii所有負責區(qū)域遍歷完一遍。當判斷其負責區(qū)域所有的商店不滿足條件時再判斷該剩余量是否大于下一個商店需求量（卩-3q），若大于則進行下一個商店服務，若小于則判斷下個ii商店是否滿足大于（卩-3q）這個條件，若滿足，則對其進行服務，否則繼續(xù)下ii個判斷，直至其所有負責區(qū)域遍歷完一遍。若所有商店的需求量（卩-3a）大于車輛貨物剩余量，則說明此車輛的剩余ii量不能滿足其所負責的區(qū)域，因此該車輛需要回倉庫進行貨物補充。當貨物補充完之后進行判斷剩余未服務商店的時間窗口和路程距離進行判斷（產(chǎn)生方法同于第二問A方案的產(chǎn)生方法），然后再進行服務。服務商店之后再進行前面的判斷,直至其所有負責商店都被服務完回到倉庫。六、 模型的評價和推廣模型的評價由5.1和5.2建立的模型，我們得到了有軟時間限制的物流配送車輛路徑問題的解決辦法。首先我們對問題進行了合理的假設，模型簡化了實際中的復雜問題，考慮了主要的約束條件與目標，思路清晰，結果也基本符合實際。但是，模型對實際進行簡化的同時忽略了實際中很多次要的條件，由累加效果可能會產(chǎn)生很大誤差，同時，我們進行模型的求解時只是簡單使用了遺傳算法，沒有對其進行改進。模型的推廣1．模型中不合實際的假設：（1） 在問題一和問題二總費用的組成上，我們沒有考慮組織送貨的費用，即每組織一輛車進行送貨都需要一定的費用，我們設為S元/（次、輛）；（2） 在問題一中，我們假設每輛車送貨時行駛的路程不超過它所能行駛的最遠路程，但是實際上，考慮到車輛行駛中的耗油等因素，每輛車的行駛最大距離都有限制，我們設車輛行駛的最大距離為L,我們可以得到約束條件：迓》xc<L s=1,2,...,mijsiji=0j=0（3） 在實際中，為了公平，每位送貨車輛司機行駛的路程相差必須控制在一定范圍之內(nèi)，我們?nèi)∵@個差值的最大允許值為M,則有：迓另xc一區(qū)另xc<M s,k=1,2,...,mijsij ijkij

2.模型的推廣廣目標函數(shù)：cxijijs正d*max{D,0}+為e*max{X,0}+mS2.模型的推廣廣目標函數(shù)：cxijijs正d*max{D,0}+為e*max{X,0}+mSiii=0 j=0 s=1i=1i=1約束條件為qyaQiis

i=1m乂yriss=1藝Kx=1ijss=1i=1KmKNx =1ijks=1j=1X=t—biiit=c/vijijs=1,2, ,mi=1,2,3,...Ni=0j=1,2,3, ,N;i=1,2,3,…,N;i=1,2,....Ni=1,2,Ni,j=1,2, Nt=KKx(t+max{D,0}+t+s)ijij ijsi iijii=1s=1s=1,2,...,mKNKNxcaLs=1,2,...,mijsiji=0j=0ii=0j=0i=0j=0xcaMijkiji=0j=0s,k=1,2,...,m根據(jù)以上目標函數(shù)與約束條件，帶入實際數(shù)據(jù)，由遺傳算法即可求解出總費用最小的車輛行駛路徑方案。七、 參考文獻[1]李軍，謝秉磊，郭耀煌，非滿載車輛調(diào)度問題的遺傳算法。系統(tǒng)工程理論與實踐，2000，20（3）：235—239；[2]邢文訓，謝金星編著現(xiàn)代優(yōu)化計算方法。北京：清華大學出版社，1999.140；魏明高成修胡潤洲，一種帶時間窗和容量約束的車輛路線問題及其TabuSearch算法，運籌與管理，Vol.11,No.3：P49-P53,2002；胡大偉朱志強胡勇，車輛路徑問題的模擬退火算法，中國公路學報，Vol.19No.4：P123-P126,2006閻慶鮑遠律，新型遺傳模擬退火算法求解物流配送路徑問題，/grid2008/detail.aspx?filename=HNSZ200804014&dbname=CJFD2008，2009/8/16張志霞邵必林，基于改進蟻群算法的運輸調(diào)度規(guī)劃，公路交通科技，Vol.25No．4：P137-P140，2008；杜文衰慶達周再玲，一類隨機庫存/運輸聯(lián)合優(yōu)化問題求解過程分析，中國公路學報，Vol.17No.1：P114-P118，2004.八、 附錄9.1附錄一表1任務的特征及其要求任務i12345678q(噸)i21.54.531.542.53s（小時）i121322.530.8[a,b]ii[1,4][4,6][1,2][4,7][3,5.5][2,5][5,8][1.5,4]表2點對之間的距離0123456780040607590200100160801400654010050751101002606507510010075757537540750100509090150490100100100010075751005200501005010007090756100757590757007010071601107590759070010088010075150100751001000附錄二Matlab7.0程序functiongatsp(s)sum=0;M=1000000;%無窮大數(shù)inn=100;citynum=8;K=2;gnmax=1000;%最大代數(shù)pm=0.8;%變異概率pc=0.2;c=[0,40,60,75,90,200,100,160,80;40,0,65,40,100,50,75,110,100;60,65,0,75,100,100,75,75,75;75,40,75,0,100,50,90,90,150;90,100,100,100,0,100,75,75,100;200,50,100,50,100,0,70,90,75;100,75,75,90,75,70,0,70,100;160,110,75,90,75,90,70,0,100;80,100,75,150,100,75,100,100,0];q=[21.54.531.542.53];s=[121322.530.8];a=[14143251];b=[46275.5584];% %產(chǎn)生初始種群m=zeros(1,inn);，m=m';KM=citynum+K-1;s=zeros(inn,citynum+K-1);fori=1:1:inns(i,:)=randperm(KM);ends=[ms];fori=1:innforj=1:KM-1ifs(i,j)>citynums(i,j)=0;endendend% %主程序functionga[f,p]=objf(s)gn=1;whilegn<gnmax+1forj=1:2:innseln=sell(s,ps);%選擇操作scro=cross(s,seln,pc);%交叉操作scnew(j,:)=scross(1,:);scnew(j+1,:)=scross(2,:);smnew(j,:)=chang(scnew(j,:),pm);%變異操作smnew(j+1,:)=chang(scnew(j+1,:),pm);ends=smnew;%產(chǎn)生了新的種群[f,p]=objf(s,dislist);%計算新種群的適應度%記錄當前代最好和平均的適應度[fmax,nmax]=max(f);ymean(gn)=1/mean(f);ymax(gn)=1/fmax;%記錄當前代的最佳個體x=s(nmax,:);gn=gn+1;%pause;endgn=gn-1;figure(2);plot(ymax,'r');holdon;plot(ymean,'b');grid;title('搜索過程')；legend('最優(yōu)解'，’平均解')；functionpcc=pro(pc)；test(1:100)=0；l=round(100*pc)；test(1:l)=1；n=round(rand*99)+1；pcc=test(n)；% %計算適應度function[f，p]=objf(s)；y=zeros(citynum+1，citynum+1)；fori=1:inn-1a=s(i，:)；forj=1:KM-1m=a(j)；n=a(j+1)；m=m+1；n=n+1；endy(m，n)=1；y=y'；fori=1:citynumforj=1:citynummubiaob=c(i，j)*y(i，:)；endendxuq1=0；fori=1:citynumforj=1:citynumxuq1=xuq1+s(i)*y(i,:)-q(i);endxuqiu=max((xuq1),0)*M;endendshij1=0;shij2=0;fori=1:citynumforj=1:citynumforl=1:citynumshij1=shij1+t(i)-a(i);shij2=shij12+b(i)-t(i);endshij3=max((shij1),0);shij4=max((shij2),0);shijian=M*shij3+M*shij4;endendf=mubiao+xuqiu+shijian;f=1/f;end% %計算選擇概率fsum=0;fori=1:innfsum=fsum+f(i);endfori=1:innps(i)=f(i)/fsum;end%計算累積概率p(1)=ps(1);fori=2:innp(i)=p(i-1)+ps(i);end，p=p';p% %“選擇”