動態(tài)規(guī)劃方法的matlab實現(xiàn)及應(yīng)用

...wd......wd......wd...動態(tài)規(guī)劃方法的matlab實現(xiàn)及其應(yīng)用〔龍京鵬，張華慶，羅明良，劉水林〕(南昌航空大學，數(shù)學與信息科學學院，江西，南昌)摘要：本文運用matlab語言實現(xiàn)了動態(tài)規(guī)劃的逆序算法，根據(jù)狀態(tài)變量的維數(shù)，編寫了指標函數(shù)最小值的逆序算法遞歸計算程序。兩個實例的應(yīng)用檢驗了該程序的有效性，同時也說明了該算法程序?qū)Ρ姸囝惖湫偷膭討B(tài)規(guī)劃應(yīng)用問題尤其是確定離散型的應(yīng)用問題的通用性，提供了求解各種動態(tài)規(guī)劃問題的有效工具。關(guān)鍵詞：動態(tài)規(guī)劃根本方程的逆序算法MATLAB實現(xiàn)MATLABAchieveForDynamicProgrammingandItsApplication(JingpengLong，HuaqingZhang，MingliangLuo，ShuilinLiu)〔SchoolofMathematicsandInformationScience,NanchangHangkongUniversity,Nanchang,China〕Abstract:Thisarticleachievesthereversealgorithmofdynamicprogrammingbyusingthematlablanguage，andpreparestherecursivecalculationprogramofreversealgorithmwhichthetargetfunctionvalueisthesmallest.Theapplicationoftwoexamplesshowthattheprogramiseffective，andthisalgorithmprogramisgeneraltomanytypicalapplicationofdynamicprogramming,especiallytheapplicationofdeterministicdiscrete.Thisalgorithmprogramprovidesaeffectivetooltothesolutionofavarietyofdynamicprogrammingproblems.Keywords:dynamicprogramming；reversealgorithm；Matlabachievement動態(tài)規(guī)劃是一類解決多階段決策問題的數(shù)學方法,在工程技術(shù)、科學管理、工農(nóng)業(yè)生產(chǎn)及軍事等領(lǐng)域都有廣泛的應(yīng)用。在理論上,動態(tài)規(guī)劃是求解這類問題全局最優(yōu)解的一種有效方法,特別是對于實際中某些非線性規(guī)劃問題可能是最優(yōu)解的唯一方法。然而,動態(tài)規(guī)劃僅僅決多階段決策問題的一種方法,或者說是考察問題的一種途徑,而不是一種具體的算法。就目前而言,動態(tài)規(guī)劃沒有統(tǒng)一的標準模型,其解法也沒有標準算法,在實際應(yīng)用中,需要具體問題具體分析。動態(tài)規(guī)劃模型的求解問題是影響動態(tài)規(guī)劃理論和方法應(yīng)用的關(guān)鍵所在,而子問題的求解和大量結(jié)果的存儲、調(diào)用更是一個難點所在。然而,隨著計算機技術(shù)的快速開展,特別是內(nèi)存容量和計算速度的增加,使求解較小規(guī)模的動態(tài)規(guī)劃問題成為可能,從而使得動態(tài)規(guī)劃的理論和方法在實際中的應(yīng)用范圍迅速增加。目前,在計算機上實現(xiàn)動態(tài)規(guī)劃的一般求解方法并不多見,尤其是用來解決較復(fù)雜的具體問題的成果甚少。本文從實際出發(fā),利用數(shù)學工具軟件matlab的強大功能,對動態(tài)規(guī)劃模型的求解方法做了嘗試,編寫出了動態(tài)規(guī)劃逆序算法的matlab程序，并結(jié)合“生產(chǎn)與存儲問題〞[1]和“背包問題〞[1]進展了應(yīng)用與檢驗,實際證明結(jié)果是令人滿意的。1動態(tài)規(guī)劃的根本模型實際中,要構(gòu)造一個標準的動態(tài)規(guī)劃模型,通常需要采用以下幾個步驟:①劃分階段按照問題的時間或空間特征,把問題分為假設(shè)干個階段。這些階段必須是有序的或者是可排序的(即無后向性),否則,應(yīng)用無效。②選擇狀態(tài)將問題開展到各個階段時所處的各種客觀情況用不同的狀態(tài)表示,即稱為狀態(tài)。狀態(tài)的選擇要滿足無后效性和可知性,即狀態(tài)不僅依賴于狀態(tài)的轉(zhuǎn)移規(guī)律,還依賴于允許決策集合和指標函數(shù)構(gòu)造。③確定決策變量與狀態(tài)轉(zhuǎn)移方程當過程處于某一階段的某個狀態(tài)時,可以做出不同的決策,描述決策的變量稱為決策變量。在決策過程中,由一個狀態(tài)到另一個狀態(tài)的演變過程稱為狀態(tài)轉(zhuǎn)移。狀態(tài)轉(zhuǎn)移就是根據(jù)上一階段的狀態(tài)和決策來導(dǎo)出本階段的狀態(tài)。④寫出動態(tài)規(guī)劃的根本方程動態(tài)規(guī)劃的根本方程一般根據(jù)實際問題可分為兩種形式,逆序形式和順序形式。這里只考慮逆序形式。動態(tài)規(guī)劃根本方程的逆序形式為fskk( )=optgvsx{(kkk( , )+fsk+1(k+1))}xDsk∈kk( )knn=,?1,?,1邊界條件fsn+1(n+1)=0或fsvsxnn( )=nnn( , )其中第k階段的狀態(tài)為sk,其決策變量xk表示狀sk的決策,狀態(tài)轉(zhuǎn)移方程為sk+1=Tsxkkk(,),態(tài)處于k階段的允許決策集合記為Dskk( ),vsxkkk( , )為指標函數(shù)。當求解時,由邊界條件從kn=開場,由后向前逆推,逐階段求出最優(yōu)決策和過程的最優(yōu)值,直到最后求出fs1(1)即得到問題的最優(yōu)解。動態(tài)規(guī)劃逆序解法計算程序框圖如下：2根本方程逆序算法的matlab程序〔1〕動態(tài)規(guī)劃逆序求最小值的根本方程：fskk( )=xDskmin{∈kk( )gvsx(kkk( , )+fsk+1(k+1))}knn=,?1,?,1邊界條件fsvsxnn( )=nnn( , )sk+1=Tsxkkk( , )。自由始端和終端的動態(tài)規(guī)劃，求指標函數(shù)最小值的逆序算法遞歸計算程序：function[p_opt,fval]=dynprog(x,DecisFun,SubObjFun,TransFun,ObjFun)%x為狀態(tài)變量，一列代表一個階段的狀態(tài)%M_函數(shù)DecisFun〔k,x)表示由階段k的狀態(tài)值x求出相應(yīng)的允許決策集合%M_函數(shù)SubObjFun〔k,x,u)表示階段k的指標函數(shù)%M_函數(shù)TransFun〔k,x,u)是狀態(tài)轉(zhuǎn)移函數(shù)，其中x是階段k的狀態(tài)值，u是其決策集合%M_函數(shù)ObjFun〔v,f)是第k階段到最后階段的指標函數(shù)，當ObjFun〔v,f)=v+f時，輸入ObjFun〔v,f)可以省略%輸出p_opt由4列組成，p_opt=[序號組，最優(yōu)軌線組，最優(yōu)策略組，指標函數(shù)值組];%輸出fval是列向量，各元素分別表示p_opt各最優(yōu)策略組對應(yīng)始端狀態(tài)x的最優(yōu)函數(shù)值k=length(x(1,:)); %k為階段數(shù)x_isnan=~isnan(x);t_vubm=inf*ones(size(x)); %t_vubm為指標函數(shù)值的上限f_opt=nan*ones(size(x));%f_opt為不同階段、狀態(tài)下的最優(yōu)值矩陣，初值為非數(shù)d_opt=f_opt;%d_opt為不同階段不同狀態(tài)下的決策矩陣，初值為非數(shù)tmp1=find(x_isnan(:,k)); %找出第k階段狀態(tài)值〔不是非數(shù)〕的下標tmp2=length(tmp1);fori=1:tmp2u=feval(DecisFun,k,x(tmp1(i),k));%求出相應(yīng)的允許決策向量tmp3=length(u);forj=1:tmp3 %該for語句是為了求出相應(yīng)的最有函數(shù)值以及最優(yōu)決策tmp=feval(SubObjFun,k,x(tmp1(i),k),u(j));iftmp<=t_vubm(i,k)f_opt(tmp1(i),k)=tmp;d_opt(tmp1(i),k)=u(j);t_vubm(i,k)=tmp;endendendforii=k-1:-1:1%從后往前面遞推求出f_opt以及d_opttmp10=find(x_isnan(:,ii));tmp20=length(tmp10);fori=1:tmp20u=feval(DecisFun,ii,x(tmp10(i),ii));tmp30=length(u);forj=1:tmp30tmp00=feval(SubObjFun,ii,x(tmp10(i),ii),u(j));tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j));%由該狀態(tài)值及相應(yīng)的決策值求出下一階段的狀態(tài)值tmp50=x(:,ii+1)-tmp40;tmp60=find(tmp50==0);%找出下一階段的狀態(tài)值在x(:,ii+1)的下標if~isempty(tmp60)ifnargin<5tmp00=tmp00+f_opt(tmp60(1),ii+1);elsetmp00=feval(ObjFun,tmp00,f_opt(tmp60(1),ii+1));endiftmp00<=t_vubm(i,ii)f_opt(tmp10(i),ii)=tmp00;d_opt(i,ii)=u(j);t_vubm(tmp10(i),ii)=tmp00;endendendendendfval=f_opt(find(x_isnan(:,1)),1);%fval即為最有函數(shù)值矩陣p_opt=[];tmpx=[];tmpd=[];tmpf=[];tmp0=find(x_isnan(:,1));tmp01=length(tmp0);fori=1:tmp01tmpd(i)=d_opt(tmp0(i),1);%求出第一階段的決策值tmpx(i)=x(tmp0(i),1);%求出第一階段的狀態(tài)值tmpf(i)=feval(SubObjFun,1,tmpx(i),tmpd(i));%求出第一階段的指標函數(shù)值p_opt(k*(i-1)+1,[1234])=[1,tmpx(i),tmpd(i),tmpf(i)];forii=2:k %按順序求出各階段的決策值、狀態(tài)值以及指標函數(shù)值tmpx(i)=feval(TransFun,ii-1,tmpx(i),tmpd(i));tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1==0);if~isempty(tmp2)tmpd(i)=d_opt(tmp2(1),ii);endtmpf(i)=feval(SubObjFun,ii,tmpx(i),tmpd(i));p_opt(k*(i-1)+ii,[1 2 34])=[ii,tmpx(i),tmpd(i),tmpf(i)];endend〔2)當狀態(tài)變量是二維時，也即有兩個狀態(tài)變量，此時動態(tài)規(guī)劃逆序求最小值的根本方程：fskk( )=kkkkkt∈(min),∈(){(gvsxtufskkkkk( , , , )+k+1(k+1))}xDsuUtknn=,?1,?,1邊界條件fsvsxtunn( )=nnnnn( , , , ),sk+1=Tsxkkk( , ),tk+1=Ptukkk( , )此時上面的程序就無能為力了，為此在程序dynprog.m根基上進展拓展，我們得到狀態(tài)變量為二維情況下的指標函數(shù)最小值的逆序算法遞歸計算程序:dynprog1.m，如下：function[p_opt,fval]=dynprog1(x1,x2,DecisFun,SubObjFun,TransFun,ObjFun)%(x1,x2)為二維狀態(tài)變量，其中x1,x2的取值是相互獨立的，這里矩陣x1與x2的列數(shù)應(yīng)一樣，該程序考慮決策變量也是二維%DecisFun(k,x1,x2),SubObjFun(k,x1,x2,u1,u2),TransFun(k,x1,x2,u1,u2)等M_函數(shù)的含義與一維的情形一樣，只是它們的參數(shù)相應(yīng)的增加了，ObjFun函數(shù)的含義及參數(shù)保持不變%[p_opt,fval]的含義與一位情形一樣，只是它們的維數(shù)增加了%下面程序的思路與算法同一維根本一樣，只是相應(yīng)矩陣的維數(shù)增加了[k1,k]=size(x1);[k2,k]=size(x2);x1_isnan=~isnan(x1);x2_isnan=~isnan(x2);t_vubm=inf*ones(k1,k2,k);f_opt=nan*ones(k1,k2,k);d_opt1=f_opt;d_opt2=f_opt;tmp11=find(x1_isnan(:,k));tmp12=length(tmp11);tmp21=find(x2_isnan(:,k));tmp22=length(tmp21);fori=1:tmp12fort=1:tmp22[u1,u2]=feval(DecisFun,k,x1(tmp11(i),k),x2(tmp21(t),k));tmp13=length(u1);tmp14=length(u2);forj=1:tmp13forl=1:tmp14tmp=feval(SubObjFun,k,x1(tmp11(i),k),x2(tmp21(t),k),u1(j),u2(l));iftmp<=t_vubm(i,t,k)f_opt(tmp11(i),tmp21(t),k)=tmp;d_opt1(tmp11(i),tmp21(t),k)=u1(j);d_opt2(tmp11(i),tmp21(t),k)=u2(l);t_vubm(i,t,k)=tmp;endendendendendforii=k-1:-1:1tmp011=find(x1_isnan(:,ii));tmp021=find(x2_isnan(:,ii));tmp012=length(tmp011);tmp022=length(tmp021);fori=1:tmp012fort=1:tmp022[u1,u2]=feval(DecisFun,ii,x1%假設(shè)決策變量為一維，那么在定義DecisFun函數(shù)時，就(tmp011(i),ii),x2(tmp021(t),ii));令u2恒為1tmp013=length(u1);tmp014=length(u2);forj=1:tmp013forl=1:tmp014tmp000=feval(SubObjFun,ii,x1(tmp011(i),ii),x2(tmp021(t),ii),u1(j),u2(l));tmp100=feval(TransFun,ii,x1(tmp011(i),ii),x2(tmp021(t),ii),u1(j),u2(l));tmp200=x1(:,ii+1)-tmp100(1);tmp300=x2(:,ii+1)-tmp100(2);tmp400=find(tmp200==0);tmp500=find(tmp300==0);if~isempty(tmp400)&~isempty(tmp500)ifnargin<6tmp000=tmp000+f_opt(tmp400(1),tmp500(1),ii+1);elsetmp000=feval(ObjFun,tmp000,f_opt(tmp400(1),mp500(1),ii+1));endiftmp000<t_vubm(i,t,ii)f_opt(tmp011(i),tmp021(t),ii)=tmp000;d_opt1(tmp011(i),tmp021(t),ii)=u1(j);d_opt2(tmp011(i),tmp021(t),ii)=u2(l);t_vubm(i,t,ii)=tmp000;endendendendendendendfval=f_opt(x1_isnan(:,1),x2_isnan(:,1),1);p_opt=[];tmpx1=[];tmpx2=[];tmpd1=[];tmpf=[];tmpd2=[]; q=0;tmp11=find(x1_isnan(:,1));tmp01=length(tmp11);tmp12=find(x2_isnan(:,1));tmp02=length(tmp12);fori=1:tmp01forj=1:tmp02tmpd1(i)=d_opt1(tmp11(i),tmp12(j),1);tmpd2(j)=d_opt2(tmp11(i),tmp12(j),1);tmpx1(i)=x1(tmp11(i),1);tmpx2(j)=x2(tmp12(j),1);tmpf(i,j)=feval(SubObjFun,1,tmpx1(i),tmpx2(j),mpd1(i),tmpd2(j));t=k*(j-1);t=q+t;p_opt(t+1,[123456])=[1,tmpx1(i),tmpx2(j),tmpd1(i),tmpd2(j),tmpf(i,j)];forii=2:ku=feval(TransFun,ii-1,tmpx1(i),tmpx2(j),tmpd1(i),tmpd2(j));tmpx1(i)=u(1);tmpx2(j)=u(2);tmp1=x1(:,ii)-tmpx1(i);tmp2=x2(:,ii)-tmpx2(j);tmp3=find(tmp1==0);tmp4=find(tmp2==0);if~isempty(tmp3)&~isempty(tmp4)tmpd1(i)=d_opt1(tmp3(1),tmp4(1),ii);tmpd2(j)=d_opt2(tmp3(1),tmp4(1),ii);endtmpf(i,j)=feval(SubObjFun,ii,tmpx1(i),tmpx2(j),mpd1(i),tmpd2(j));p_opt(t+ii,[123456])=[ii,tmpx1(i),tmpx2(j),tmpd1(i),tmpd2(j),tmpf(i,j)];endendp=k*tmp02;q=i*p;end3實例應(yīng)用3.1生產(chǎn)與存儲問題某工廠要對一種產(chǎn)品制定今后四個時期的生產(chǎn)方案，據(jù)估計在今后四個時期內(nèi)，市場對于該產(chǎn)品的需求量如下表所示：時期〔k)1233需求量〔dk)2324假定該廠生產(chǎn)每批產(chǎn)品的固定本錢為3千克，，假設(shè)不生產(chǎn)就為0；每單位產(chǎn)品本錢為1千元；每個時期生產(chǎn)能力所允許的最大生產(chǎn)量不超過6各單位；每個時期末未能售出的產(chǎn)品，每單位需付存儲費0.5千元。還假定在第一個時期的初始庫存量為0，第四個是期末的庫存量也為0。試問該廠應(yīng)假設(shè)何安排各個時期的生產(chǎn)與庫存，才能在滿足市場需要的條件下，使總本錢最小。解：用動態(tài)規(guī)劃方法求解，按四個時期將該問題分為四個階段；記Vk為狀態(tài)變量，它表示第k階段開場時的庫存量；記Xk為決策變量，它表示第k階段的生產(chǎn)量；可得狀態(tài)轉(zhuǎn)移方程：Vk+1=Vk+Xk-dk,k=1,2,3,4由題意知，在第k時期內(nèi)的生產(chǎn)本錢為:?0...................當xk=0?cxkk()=?3+1*xk.......當xk=1,2,...6?????∞..................當xk>6 ??在第k時期末庫存量為vk+1時的存儲費為:hvkk()=0.5*(vxdk+?kk)故第k時期內(nèi)的總本錢為:cxhvkk()+kk()則階段指標函數(shù)為:Vvcxhvkk( )=kk( )+kk( )最優(yōu)值函數(shù)fVkk( )表示從第k階段初始庫存量為Vk時到第四階段末庫存量為0時的最小總費用。則有遞推關(guān)系式:??fvkk( )=max(0,dkvkxkmin?)≤≤6(Vvkk( )+fvk+1(k+1))???fv5(5)=0,xdv4=?4 4其中xk≥max(0,dvk?k),這是因為一方面每階段生產(chǎn)的下限為0；另一方面由于要保證供應(yīng)，故第k階段末的庫存量vk+1必須非負，即vxdk+?k k≥0，所以xdvk≥?k k。vk的取值范圍為[0,min(∑dmdj,?k)],其中v1=0,v5=0。jk=為求出該問題的最優(yōu)值，利用上面的計算程序dynprog.m。根據(jù)上面所述的階段指標函數(shù)、狀態(tài)轉(zhuǎn)移函數(shù)和遞推關(guān)系式，編寫出下面3個M_函數(shù)，以備主程序調(diào)用。%DecisFun.mfunctionu=DecisFun(k,x)d=[2324];m=6;ifk==4u=d(k)-x;elseu=max(0,d(k)-x):m;End%SubObjFun.mfunctionf=SubObjFun(k,x,u)d=[2324];ifu==0f=0.5*(x+u-d(k));elseifu>6f=10^6;elsef=3+u+0.5*(x+u-d(k));endEnd%TransFun.mfunctions=TransFun(k,x,u)d=[2324];s=x+u-d(k);在matlab命令空間輸入：x1=0:4;s=nan*ones(5,1);s(1)=0;x=[sx1'x1'x1'];[p_opt,fval]=dynprog(x,'DecisFun','SubObjFun','TransFun')運行結(jié)果如下：p_opt=1.000005.00009.50002.00003.0000003.000006.000011.00004.0000fval=20.50004.000000從上面的結(jié)果可知，每個時期的最優(yōu)決策為：X1=5,x2=0,x3=6,x4=0。其相應(yīng)的最小總本錢為20.5千元。從上面的結(jié)果中還可以看出，各個時期初的庫存量分別為：V1=0,v2=2,v3=0,v4=4這里的結(jié)果與文[1]的結(jié)果完全符合，這說明該程序是可行的。3.2二維背包問題有一個人帶一個背包上山，其可攜帶物品重量的限度為10公斤，背包體積限制為22立方米。假設(shè)有3種物品可供他選擇裝入背包。第i種物品每件重量為w(i)公斤，體積為v(i)立方米，攜帶該物品xi件產(chǎn)生的效益值為c(i)*xi。問此人該假設(shè)何選擇攜帶物品，才能使產(chǎn)生的效益值最大。其中w=[345];v=[864];c=[456];解：用動態(tài)規(guī)劃方法求解，按物品的種類數(shù)將該問題分為3各個階段；si表示用于裝入第i種物品到第3種物品的總重量；ti表示用于裝入第i種物品到第3種物品的總體積；ui表示裝入第i種物品的件數(shù)；可得狀態(tài)轉(zhuǎn)移方程：sk+1=?sckutk ()*kk,+1=?tckuk ()*k允許決策集合為：stk,k)}Dsu(kk, )={uk|0≤≤uk min(wvk k最優(yōu)值函數(shù)fstkkk( , )表示當總重量不超過sk公斤，總體積不超過tk立方米背包裝入第t種物品到第3種物品產(chǎn)生的最大效益值。可得根本方程：??fstkkk( , )=uDstk∈maxkkk( , )(()*ckufstk+k+1(k+1,k+1)),???fvt4(4,4)=0k=3,2,1下面同樣用計算程序dynprog1.m求解：在使用此程序先要建設(shè)下面3個M_函數(shù)：%DecisFun1.mfunction[u1,u2]=DecisFun1(k,x1,x2)w=[345];v=[864];u1=0:1:min(x1/w(k),x2/v(k));u2=1;%因為這里只有一個決策變量，故令u2恒為1，這樣是程序的需要，%也可減少計算量，此時u2就沒有任何意義，只是一個虛擬的量%SubObjFun1.mfunctionf=SubObjFun1(k,x1,x2,u1,u2)c=[456];f=-c(k)*u1;%求最大值轉(zhuǎn)化為求最小值%TransFun1.mfunctions=TransFun1(k,x1,x2,u