1、動(dòng)態(tài)規(guī)劃方法的matlab實(shí)現(xiàn)及其應(yīng)用(龍京鵬，張華慶，羅明良，劉水林)(南昌航空大學(xué)，數(shù)學(xué)與信息科學(xué)學(xué)院，江西，南昌)摘要：本文運(yùn)用matlab語(yǔ)言實(shí)現(xiàn)了動(dòng)態(tài)規(guī)劃的逆序算法，根據(jù)狀態(tài)變量的維數(shù)，編寫了指標(biāo)函數(shù)最小值的逆序算法遞歸計(jì)算程序。兩個(gè)實(shí)例的應(yīng)用檢驗(yàn)了該程序的有效性，同時(shí)也表明了該算法程序?qū)Ρ姸囝惖湫偷膭?dòng)態(tài)規(guī)劃應(yīng)用問題尤其是確定離散型的應(yīng)用問題的通用性，提供了求解各種動(dòng)態(tài)規(guī)劃問題的有效工具。關(guān)鍵詞：動(dòng)態(tài)規(guī)劃基本方程的逆序算法MATLAB實(shí)現(xiàn)MATLABAchieveForDynamicProgrammingandItsApplication(JingpengLong,HuaqingZh

2、ang,MingliangLuo,ShuilinLiu)(SchoolofMathematicsandInformationScience,NanchangHangkongUniversity,Nanchang,China)Abstract:Thisarticleachievesthereversealgorithmofdynamicprogrammingbyusingthematlablanguage，andpreparestherecursivecalculationprogramofreversealgorithmwhichthetargetfunctionvalueisthesmall

3、est.Theapplicationoftwoexamplesshowthattheprogramiseffective，andthisalgorithmprogramisgeneraltomanytypicalapplicationofdynamicprogramming,especiallytheapplicationofdeterministicdiscrete.Thisalgorithmprogramprovidesaeffectivetooltothesolutionofavarietyofdynamicprogrammingproblems.Keywords:dynamicprog

4、ramming;reversealgorithm;Matlabachievement動(dòng)態(tài)規(guī)劃是一類解決多階段決策問題的數(shù)學(xué)方法，在工程技術(shù)、科學(xué)管理、工農(nóng)業(yè)生產(chǎn)及軍事等領(lǐng)域都有廣泛的應(yīng)用。在理論上，動(dòng)態(tài)規(guī)劃是求解這類問題全局最優(yōu)解的一種有效方法，特別是對(duì)于實(shí)際中某些非線性規(guī)劃問題可能是最優(yōu)解的唯一方法。然而，動(dòng)態(tài)規(guī)劃僅僅決多階段決策問題的一種方法，或者說(shuō)是考查問題的一種途徑，而不是一種具體的算法。就目前而言，動(dòng)態(tài)規(guī)劃沒有統(tǒng)一的標(biāo)準(zhǔn)模型，其解法也沒有標(biāo)準(zhǔn)算法，在實(shí)際應(yīng)用中，需要具體問題具體分析。動(dòng)態(tài)規(guī)劃模型的求解問題是影響動(dòng)態(tài)規(guī)劃理論和方法應(yīng)用的關(guān)鍵所在，而子問題的求解和大量結(jié)果的存儲(chǔ)、調(diào)用更是

5、一個(gè)難點(diǎn)所在。然而，隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，特別是內(nèi)存容量和計(jì)算速度的增加，使求解較小規(guī)模的動(dòng)態(tài)規(guī)劃問題成為可能，從而使得動(dòng)態(tài)規(guī)劃的理論和方法在實(shí)際中的應(yīng)用范圍迅速增加。目前，在計(jì)算機(jī)上實(shí)現(xiàn)動(dòng)態(tài)規(guī)劃的一般求解方法并不多見，尤其是用來(lái)解決較復(fù)雜的具體問題的成果甚少。本文從實(shí)際出發(fā)，利用數(shù)學(xué)工具軟件matlab的強(qiáng)大功能，對(duì)動(dòng)態(tài)規(guī)劃模型的求解方法做了嘗試，編寫出了動(dòng)態(tài)規(guī)劃逆序算法的matlab程序，并結(jié)合“生產(chǎn)與存儲(chǔ)問題”1和“背包問題”1進(jìn)行了應(yīng)用與檢驗(yàn)，實(shí)際證明結(jié)果是令人滿意的。1動(dòng)態(tài)規(guī)劃的基本模型實(shí)際中，要構(gòu)造一個(gè)標(biāo)準(zhǔn)的動(dòng)態(tài)規(guī)劃模型，通常需要采用以下幾個(gè)步驟：劃分階段按照問題的時(shí)間或空間特

6、征，把問題分為若干個(gè)階段。這些階段必須是有序的或者是可排序的(即無(wú)后向性)，否則，應(yīng)用無(wú)效。選擇狀態(tài)將問題發(fā)展到各個(gè)階段時(shí)所處的各種客觀情況用不同的狀態(tài)表示，即稱為狀態(tài)。狀態(tài)的選擇要滿足無(wú)后效性和可知性，即狀態(tài)不僅依賴于狀態(tài)的轉(zhuǎn)移規(guī)律，還依賴于允許決策集合和指標(biāo)函數(shù)結(jié)構(gòu)。確定決策變量與狀態(tài)轉(zhuǎn)移方程當(dāng)過程處于某一階段的某個(gè)狀態(tài)時(shí)，可以做出不同的決策，描述決策的變量稱為決策變量。在決策過程中，由一個(gè)狀態(tài)到另一個(gè)狀態(tài)的演變過程稱為狀態(tài)轉(zhuǎn)移。狀態(tài)轉(zhuǎn)移就是根據(jù)上一階段的狀態(tài)和決策來(lái)導(dǎo)出本階段的狀態(tài)。寫出動(dòng)態(tài)規(guī)劃的基本方程動(dòng)態(tài)規(guī)劃的基本方程一般根據(jù)實(shí)際問題可分為兩種形式，逆序形式和順序形式。這里只考慮逆序

7、形式。動(dòng)態(tài)規(guī)劃基本方程的逆序形式為fskk()=optgvsx(kkk(,)+fsk+1(k+1)xDskkk()knn=,-1,1邊界條件fsn+1(n+1)=0或fsvsXn()=nnn(,)其中第k階段的狀態(tài)為sk,其決策變量xk表示狀Sk的決策，狀態(tài)轉(zhuǎn)移方程為Sk+1=TSXkk(,),態(tài)處于k階段的允許決策集合記為DSkk(),vSXkk(,)為指標(biāo)函數(shù)。當(dāng)求解時(shí)，由邊界條件從kn=開始，由后向前逆推，逐階段求出最優(yōu)決策和過程的最優(yōu)值，直到最后求出fsi(i)即得到問題的最優(yōu)解。動(dòng)態(tài)規(guī)劃逆序解法計(jì)算程序框圖如下：2基木方程逆序算法的matlab程序(1二動(dòng)態(tài)關(guān)tks)=O(I劃逆序求

8、最小值的基本方程:xdskminekk()gvsx(kkk(sk+1(k+1)knknn=,-1,1邊界條1斗fsvsxn()=nnn(,)生+1kk(自由始端和終端的動(dòng)態(tài)規(guī)劃，求指標(biāo)函數(shù)最小值的逆序算法涕歸計(jì)算程序:function_,SubObjFu齊缶J二口9t,出工第器熄juf黑4S2IM與瓏就盤一段的狀由階思態(tài)k的狀態(tài)值x求相應(yīng)的允許決策集合x,u)表示階段k的指標(biāo)函,u)是狀態(tài)轉(zhuǎn)移函數(shù)，其中%M_函數(shù)SubObjFun(k存儲(chǔ)耳aj和m茶港端屋段4狀態(tài)叫最優(yōu)值矩陣初值l_opt態(tài)下的決策矩陣，初值為非數(shù)為不同階段不同狀tmp1=fing(x_isnan(:國(guó));狀態(tài)值布是非數(shù)是%找

9、出第k階段下標(biāo)k=l中p2=length(tmp1);fOri=1:tmr>2tmp=feval(SubObjFun,k,x(tmp1(i),k),u(j)u=feval(DecisFun,k,x(tmp1(i),k);%求出相應(yīng)的允許決策向量tmp3=length(u);forj=1:tmp3目應(yīng)的最肓函數(shù)值以及最優(yōu)決策語(yǔ)句是為了求)+k;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_opt

10、tmp10=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)的下標(biāo)if

11、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

12、(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);%求出第一階段的指標(biāo)函數(shù)值p_opt(k*(i-1)+1,1234)=1,tmpx(i),tmpd(i),tmpf(i);forii=2:k%按順序求出各階段的決策值、狀態(tài)值以及指標(biāo)函數(shù)值tmpx(i)=feval(TransFun,ii-1,tmpx(i),tmpd(i);tmp1=x(

13、:,ii)-tmpx(i);tmp2=find(tmp1=0);ifisempty(tmp2)tmpd(i)=d_opt(tmp2(1),ii);endtmpf(i)=feval(SubObjFun,ii,tmpx(i),tmpd(i);p_opt(k*(i-1)+ii,1234)=ii,tmpx(i),tmpd(i),tmpf(i);endend(2)當(dāng)狀態(tài)變量是二維時(shí)，也即有兩個(gè)狀態(tài)變量，此時(shí)動(dòng)態(tài)規(guī)劃逆序求最小值的基本方程：fSkk()=kkkkkte(min),e()(gvsxtufSkkkk(,xDsuUtknn=,-1,1邊界條件fSvSxtLIn()=nnnnn(,，)，sk+1=

14、Tsxkkk(,),tk+1=Ptukkk(,)此時(shí)上面的程序就無(wú)能為力了，為此在程序dynprog.m基礎(chǔ)上進(jìn)行拓展，我們得到狀態(tài)變量為二維情況下的指標(biāo)函數(shù)最小值的逆序算法遞歸計(jì)算程序：dynprog1.m,如下：functionp_opt,fval=dynprog1(x1,x2,DecisFun,SubObjFun,TransFun,ObjFun)%(x1,x2)為二維狀態(tài)變量，其中x1,x2的取值是相互獨(dú)立的，這里矩陣x1與x2的列數(shù)應(yīng)相同，該程序考慮決策變量也是二維%DecisFun(k,x1,x2),SubObjFun(k,x1,x2,u1,u2),TransFun(k,x1,x2,

15、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(tmp1

16、1);tmp21=find(x2_isnan(:,k);tmp22=length(tmp21);fori=1:tmp12fort=1:tmp22u1,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

17、;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:tmp022u1,u2=feval(DecisFun,ii,x1(tmp011(i),ii),x2(tmp021(t),ii);u2

18、恒為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(tm

19、p300=0);ifisempty(tmp400)&isempty(tmp500)ifnargin<6tmp000=tmp000+f_opt(tmp400(1),tmp500(1),ii+1);elsetmpO00=feval(ObjFun,tmp000,f_opt(tmp400(1),mp500(1),ii+1);endiftmp000<t_vubm(i,t,ii)f_opt(tmpO11(i),tmp021(t),ii)=tmp000;d_opt1(tmp011(i),tmp021(t),ii)=u1(j);d_opt2(tmp011(i),tmp021(t),ii)=

20、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_o

21、pt2(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);t

22、mpx2(j)=u(2);tmp1=x1(:,ii)-tmpx1(i);tmp2=x2(:,ii)-tmpx2(j);tmp3=find(tmp1=0);tmp4=find(tmp2=0);ifisempty(tmp3)&isempty(tmp4)%若決策變量為一維，那么在定義DecisFun函數(shù)時(shí)，就令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_op

23、t(t+ii,123456)=ii,tmpx1(i),tmpx2(j),tmpd1(i),tmpd2(j),tmpf(i,j);endendp=k*tmp02;q=i*p;end3實(shí)例應(yīng)用3.1生產(chǎn)與存儲(chǔ)問題時(shí)期（k）1233需求量（dk）2324某工廠要對(duì)一種產(chǎn)品制定今后四個(gè)時(shí)期的生產(chǎn)計(jì)劃，據(jù)估計(jì)在今后四個(gè)時(shí)期內(nèi)，市場(chǎng)對(duì)于該產(chǎn)品的需求量如下表所示：假定該廠生產(chǎn)每批產(chǎn)品的固定成本為3千克，若不生產(chǎn)就為0;每單位產(chǎn)品成本為1千元；每個(gè)時(shí)期生產(chǎn)能力所允許的最大生產(chǎn)量不超過6各單位；每個(gè)時(shí)期末未能售出的產(chǎn)品，每單位需付存儲(chǔ)費(fèi)0.5千元。還假定在第一個(gè)時(shí)期的初始庫(kù)存量為0,第四個(gè)是期末的庫(kù)存量也為0。

24、試問該廠應(yīng)如何安排各個(gè)時(shí)期的生產(chǎn)與庫(kù)存，才能在滿足市場(chǎng)需要的條件下，使總成本最小。解：用動(dòng)態(tài)規(guī)劃方法求解，按四個(gè)時(shí)期將該問題分為四個(gè)階段；記Vk為狀態(tài)變量，它表示第k階段開始時(shí)的庫(kù)存量；記Xk為決策變量，它表示第k階段的生產(chǎn)量；可得狀態(tài)轉(zhuǎn)移方程：V.尸M+X-dk,k=1,2,3,4由題意知，在第k時(shí)期內(nèi)的生產(chǎn)成本為：?0當(dāng)Xk=0?cxkk()=?3+1*Xk.當(dāng)Xk=1,2,.6?oo當(dāng)Xk>6?在第k時(shí)期末庫(kù)存量為Vk+1時(shí)的存儲(chǔ)費(fèi)為：hVkk()=0.5*(vx&+-kk)故第k時(shí)期內(nèi)的總成本為：cXhVkk()+kk()則階段指標(biāo)函數(shù)為：Vvcxhvkk()=kk()+

25、kk()最優(yōu)值函數(shù)fVkk()表示從第k階段初始庫(kù)存量為Vk時(shí)到第四階段末庫(kù)存量為0時(shí)的最小總費(fèi)用。則有遞推關(guān)系式：?fVkk()=maX(0,dkvkXkmih-尸6(VVkk()+fVk+1(k+1)?V5(5)=0,XdV4=-44其中Xk>maX(0，4k),這是因?yàn)橐环矫婷侩A段生產(chǎn)的下限為0；另一方面由于要保證供應(yīng)，故第k階段末的庫(kù)存量Vk+1必須非負(fù)，即VXdk+-kk>0,所以XdVk>-kkoVk的取值范圍為0,min(Xdmd,-k),其中V1=0,V5=0。jk=為求出該問題的最優(yōu)值，利用上面的計(jì)算程序dynprog.m。根據(jù)上面所述的階段指標(biāo)函數(shù)、狀態(tài)轉(zhuǎn)

26、移函數(shù)和遞推關(guān)系式，編寫出下面3個(gè)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=10A6;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

27、:4;s=nan*ones(5,1);s(1)=0;X=sx1'x1'x1'p_opt,fVal=dynprog(X,'DecisFun','SubObjFun','TransFun')運(yùn)行結(jié)果如下：p_opt=1.000005.00009.50002.00003.0000003.000006.000011.00004.00004.000000fval=20.5000從上面的結(jié)果可知，每個(gè)時(shí)期的最優(yōu)決策為：X1=5,x2=0,x3=6,x4=0。其相應(yīng)的最小總成本為20.5千元。從上面的結(jié)果中還可以看出，各個(gè)時(shí)期初的庫(kù)存量

28、分別為：V1=0,v2=2,v3=0,v4=4這里的結(jié)果與文1的結(jié)果完全符合，這說(shuō)明該程序是可行的。3.2二維背包問題有一個(gè)人帶一個(gè)背包上山，其可攜帶物品重量的限度為10公斤，背包體積限制為22立方米。假設(shè)有3種物品可供他選擇裝入背包。已知第i種物品每件重量為w(i)公斤，體積為v(i)立方米，攜帶該物品xi件產(chǎn)生的效益值為c(i)*xi。問此人該如何選擇攜帶物品，才能使產(chǎn)生的效益值最大。其中w=345;v=864;c=456;解：用動(dòng)態(tài)規(guī)劃方法求解，按物品的種類數(shù)將該問題分為3各個(gè)階段；si表示用于裝入第i種物品到第3種物品的總重量；ti表示用于裝入第i種物品到第3種物品的總體積；ui表示裝

29、入第i種物品的件數(shù)；可得狀態(tài)轉(zhuǎn)移方程：Sk+1=-sckutk()*kk,+1=-tckUk()*k允許決策集合為：s告,DDsu(kk,)=Uk|0<ukmin(wvkk最優(yōu)值函數(shù)fstkkk(,)表示當(dāng)總重量不超過sk公斤，總體積不超過tk立方米背包裝入第t種物品到第3種物品產(chǎn)生的最大效益值。可得基本方程：?stkkk(,)=udstkmaxkkk()()*ckufsk+k+1(k+1,k+1),?vt4(4,4)=0k=3,2,1下面同樣用計(jì)算程序dynprog1.m求解：在使用此程序先要建立下面3個(gè)M畫數(shù)：%DecisFun1.mfunctionu1,u2=DecisFun1(k

30、,x1,x2)w=345;v=864;u1=0:1:min(x1/w(k),x2/v(k);u2=1;%因?yàn)檫@里只有一個(gè)決策變量，故令u2恒為1,這樣是程序的需要，%也可減少計(jì)算量，此時(shí)u2就沒有任何意義，只是一個(gè)虛擬的量%SubObjFun1.mfunctionf=SubObjFun1(k,x1,x2,u1,u2)c=456;f=-c(k)*u1;%求最大值轉(zhuǎn)化為求最小值%TransFun1.mfunctions=TransFun1(k,x1,x2,u1,u2)w=345;v=864;s(1)=x1-u1*w(k);s(2)=x2-u1*v(k);在matlab命令空間輸入：a1=0:10;b1=0:22;s1=nan*ones(11,1);s1(1)=10;s2=nan*ones(2