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

1、動(dòng)態(tài)規(guī)劃方法的matlab實(shí)現(xiàn)及其應(yīng)用（龍京鵬，張華慶，羅明良，劉水林） (南昌航空大學(xué)，數(shù)學(xué)與信息科學(xué)學(xué)院，江西，南昌)摘要：本文運(yùn)用matlab語言實(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)MATLAB Achieve For Dynamic Programming and Its Application(JingpengLong

2、，HuaqingZhang，MingliangLuo，ShuilinLiu)（School of Mathematics and Information Science,Nanchang HangkongUniversity,Nanchang,China）Abstract:This article achieves the reverse algorithm of dynamic programming by using the matlab language，and prepares the recursive calculation program of reverse algorithm

3、 which thetargetfunctionvalueisthesmallest.Theapplicationoftwoexamplesshowthattheprogram is effective，and this algorithm program is general to many typical application of dynamic programming,especially the application of deterministic discrete.This algorithm program provides a effective tool to the

4、solution of a variety of dynamic programming problems. Key words:dynamic programming；reverse algorithm；Matlab achievement動(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ī)劃僅僅決多階段決策問題的一種方法,或者說是考查問題的一種途徑,而不是一種具體的算法。就目前而言,動(dòng)態(tài)規(guī)劃沒有統(tǒng)一的標(biāo)準(zhǔn)模型,

5、其解法也沒有標(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)用更是一個(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ī)劃的一般求解方法并不多見,尤其是用來解決較復(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 和“背

6、包問題”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í)間或空間特征,把問題分為若干個(gè)階段。這些階段必須是有序的或者是可排序的(即無后向性) ,否則,應(yīng)用無效。選擇狀態(tài) 將問題發(fā)展到各個(gè)階段時(shí)所處的各種客觀情況用不同的狀態(tài)表示,即稱為狀態(tài)。狀態(tài)的選擇要滿足無后效性和可知性,即狀態(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)的演變過

7、程稱為狀態(tài)轉(zhuǎn)移。狀態(tài)轉(zhuǎn)移就是根據(jù)上一階段的狀態(tài)和決策來導(dǎo)出本階段的狀態(tài)。寫出動(dòng)態(tài)規(guī)劃的基本方程 動(dòng)態(tài)規(guī)劃的基本方程一般根據(jù)實(shí)際問題可分為兩種形式,逆序形式和順序形式。這里只考慮逆序形式。動(dòng)態(tài)規(guī)劃基本方程的逆序形式為f sk k() = opt gv s x ( k k k(,)+f sk+1( k+1) x D sk k k()k nn= , 1,1邊界條件f sn+1( n+1) = 0或 f s v s xn n() = n n n(,)其中第k 階段的狀態(tài)為sk,其決策變量xk表示狀 sk的決策,狀態(tài)轉(zhuǎn)移方程為sk+1 =T s xk k k( , ), 態(tài)處于k 階段的允許決策集合記為

8、D sk k() , v s xk k k(,) 為指標(biāo)函數(shù)。當(dāng)求解時(shí),由邊界條件從k n= 開始, 由后向前逆推,逐階段求出最優(yōu)決策和過程的最優(yōu)值, 直到最后求出 f s1( 1) 即得到問題的最優(yōu)解。動(dòng)態(tài)規(guī)劃逆序解法計(jì)算程序框圖如下：2 基本方程逆序算法的matlab程序（1）動(dòng)態(tài)規(guī)劃逆序求最小值的基本方程：f sk k() = x D skmin k k() gv s x( k k k(,) +f sk+1( k+1)k nn= , 1,1邊界條件 f s v s xn n() = n n n(,) sk+1 =T s xk k k(,) 。自由始端和終端的動(dòng)態(tài)規(guī)劃，求指標(biāo)函數(shù)最小值的逆

9、序算法遞歸計(jì)算程序：functionp_opt,fval=dynprog(x,DecisFun,SubObjFu n,TransFun,ObjFun)% x為狀態(tài)變量，一列代表一個(gè)階段的狀態(tài)%M_函數(shù)DecisFun（k,x)表示由階段k的狀態(tài)值x求出相應(yīng)的允許決策集合% M_函數(shù)SubObjFun（k,x,u)表示階段k的指標(biāo)函數(shù)% M_函數(shù)TransFun（k,x,u)是狀態(tài)轉(zhuǎn)移函數(shù)，其中x 是階段k的狀態(tài)值，u是其決策集合% M_函數(shù)ObjFun（v,f)是第k階段到最后階段的指標(biāo)函數(shù)，當(dāng)ObjFun（v,f)=v+f時(shí)，輸入ObjFun（v,f)可以省略% 輸出p_opt由4列組成，p

10、_opt=序號(hào)組，最優(yōu)軌線組，最優(yōu)策略組，指標(biāo)函數(shù)值組; % 輸出fval是列向量，各元素分別表示p_opt各最優(yōu)策略組對(duì)應(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為指標(biāo)函數(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ù)）的下標(biāo)tmp2=lengt

11、h(tmp1); for i=1:tmp2u=feval(DecisFun,k,x(tmp1(i),k);% 求出相應(yīng)的允許決策向量 tmp3=length(u);for j=1:tmp3% 該for語句是為了求出相應(yīng)的最有函數(shù)值以及最優(yōu)決策tmp=feval(SubObjFun,k,x(tmp1(i),k),u(j); if tmp<=t_vubm(i,k)f_opt(tmp1(i),k)=tmp; d_opt(tmp1(i),k)=u(j); t_vubm(i,k)=tmp; end endend for ii=k-1:-1:1% 從后往前面遞推求出f_opt以及d_opt tmp1

12、0=find(x_isnan(:,ii);tmp20=length(tmp10);for i=1:tmp20u=feval(DecisFun,ii,x(tmp10(i),ii);tmp30=length(u); for j=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)

13、 if isempty(tmp60)if nargin<5tmp00=tmp00+f_opt(tmp60(1),ii+1); elsetmp00=feval(ObjFun,tmp00,f_opt(tmp60(1), ii+1); endif tmp00<=t_vubm(i,ii)f_opt(tmp10(i),ii)=tmp00;d_opt(i,ii)=u(j);t_vubm(tmp10(i),ii)=tmp00; endendendend endfval=f_opt(find(x_isnan(:,1),1);% fval即為最有函數(shù)值矩陣p_opt=;tmpx=;tmpd=;tmp

14、f=; tmp0=find(x_isnan(:,1);tmp01=length(tmp0); for i=1:tmp01 tmpd(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,12 3 4)=1,tmpx(i),tmpd(i),tmpf(i);for ii=2:k% 按順序求出各階段的決策值、狀態(tài)值以及指標(biāo)函數(shù)值tmpx(i)=feval(TransFun,ii-1,

15、tmpx(i),tmpd (i);tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1=0); if isempty(tmp2) tmpd(i)=d_opt(tmp2(1),ii);end tmpf(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ī)劃逆序求最小值的基本方程：f sk k() = k k k k k t (min), ( ) (gv s x t u f sk k k

16、k k(,)+ k+1( k+1)x D s u U tk nn= , 1,1邊界條件f s v s x t un n() = n n n n n(,),sk+1 =T s xk k k(,),tk+1 =P t uk k k(,)此時(shí)上面的程序就無能為力了，為此在程序 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的取值是相互

17、獨(dú)立的，這里矩陣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(k

18、1,k2,k);f_opt=nan*ones( k1,k2,k);d_opt1=f_opt;d_opt2=f_opt;tmp11=find(x1_isnan(:,k);tmp12=length(t mp11);tmp21=find(x2_isnan(:,k);tmp22=length(tmp21); for i=1:tmp12 for t=1:tmp22u1,u2=feval(DecisFun,k,x1(tmp11(i),k),x2(tmp21(t),k); tmp13=length(u1);tmp14=length(u2);for j=1:tmp13 for l=1:tmp14tmp=fev

19、al(SubObjFun,k,x1(tmp11(i),k),x2(tmp21(t),k),u1(j),u2(l); if tmp<=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; endendendendendfor ii=k-1:-1:1 tmp011=find(x1_isnan(:,ii); tmp021=find(x2_isnan(:,ii); tmp012=le

20、ngth(tmp011); tmp022=length(tmp021); for i=1:tmp012 for t=1:tmp022u1,u2=feval(DecisFun,ii,x1% 若決策變量為一維，那么在定義DecisFun函數(shù)時(shí)，就(tmp011(i),ii),x2(tmp021(t),ii); 令u2恒為1tmp013=length(u1); tmp014=length(u2); for j=1:tmp013for l=1:tmp014tmp000=feval(SubObjFun,ii,x1(tmp011(i),ii),x2(tmp021(t),ii),u1(j),u2(l); t

21、mp100=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);ifisempty(tmp400)&isempty(tmp500) if nargin<6tmp000=tmp000+f_opt(tmp400(1),tmp500(1),ii+1); else tmp000=feval(ObjFun,tm

22、p000, f_opt(tmp400(1),mp500(1),ii+1);end if tmp000<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; end end end end end end endfval=f_opt(x1_isnan(:,1),x2_isnan(:,1),1);p_opt=;tmpx1=;tmpx2=;tmpd

23、1=;tmpf =;tmpd2=;q=0; tmp11=find(x1_isnan(:,1);tmp01=length(t mp11);tmp12=find(x2_isnan(:,1);tmp02=length(tmp12); for i=1:tmp01for j=1:tmp02 tmpd1(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(

24、i),tmpx2(j), mpd1(i),tmpd2(j);t=k*(j-1);t=q+t;p_opt(t+1,1 2 3 4 5 6)=1,tmpx1(i),tmpx2(j), tmpd1(i),tmpd2(j),tmpf(i,j); for ii=2:k u=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

25、=0);if isempty(tmp3)&isempty(tmp4) tmpd1(i)=d_opt1(tmp3(1),tmp4(1),ii); tmpd2(j)=d_opt2(tmp3(1),tmp4(1),ii); end tmpf(i,j)=feval(SubObjFun,ii,tmpx1(i),tmpx2(j),mpd1(i),tmpd2(j); p_opt(t+ii,1 2 3 4 5 6)=ii,tmpx1(i),tmpx2(j), tmpd1(i),tmpd2(j),tmpf(i,j);endend p=k*tmp02;q=i*p;end3 實(shí)例應(yīng)用3.1 生產(chǎn)與存儲(chǔ)問題某

26、工廠要對(duì)一種產(chǎn)品制定今后四個(gè)時(shí)期的生產(chǎn)計(jì)劃，據(jù)估計(jì)在今后四個(gè)時(shí)期內(nèi)，市場(chǎng)對(duì)于該產(chǎn)品的需求量如下表所示：時(shí)期（k)1233需求量（dk)2324假定該廠生產(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í)期的初始庫存量為0，第四個(gè)是期末的庫存量也為0。試問該廠應(yīng)如何安排各個(gè)時(shí)期的生產(chǎn)與庫存，才能在滿足市場(chǎng)需要的條件下，使總成本最小。解：用動(dòng)態(tài)規(guī)劃方法求解，按四個(gè)時(shí)期將該問題分為四個(gè)階段；記Vk為狀態(tài)變量，它表示第k階段開始時(shí)的庫存量；記Xk為決策變量，它表示

27、第k階段的生產(chǎn)量；可得狀態(tài)轉(zhuǎn)移方程：Vk+1=Vk+Xk-dk,k=1,2,3,4 由題意知，在第k時(shí)期內(nèi)的生產(chǎn)成本為:0.當(dāng) xk =0 c xk k( ) =3+1*xk.當(dāng) xk =1,2,.6 .當(dāng) xk >6在第 k時(shí)期末庫存量為vk+1 時(shí)的存儲(chǔ)費(fèi)為:h vk k( ) = 0.5*(v x dk + k k) 故第k時(shí)期內(nèi)的總成本為: c x h vk k( ) + k k( ) 則階段指標(biāo)函數(shù)為:V v c x h vk k() = k k() + k k()最優(yōu)值函數(shù) f Vk k() 表示從第k階段初始庫存量為Vk時(shí)到第四階段末庫存量為0時(shí)的最小總費(fèi)用。則有遞推關(guān)系式

28、:f vk k() = max(0,dk vk xkmin ) 6(V vk k() +f vk+1( k+1)f v5 ( 5 ) = 0, x d v4 = 44其中xk max(0,d vk k),這是因?yàn)橐环矫婷侩A段生產(chǎn)的下限為0；另一方面由于要保證供應(yīng)，故第k階段末的庫存量vk+1 必須非負(fù)，即v x dk + kk 0 ，所以 x d vk kk。vk的取值范圍為0, min(d m dj, k) ,其中v1 = 0,v5 = 0 。j k=為求出該問題的最優(yōu)值，利用上面的計(jì)算程序 dynprog.m。根據(jù)上面所述的階段指標(biāo)函數(shù)、狀態(tài)轉(zhuǎn)移函數(shù)和遞推關(guān)系式，編寫出下面3個(gè)M_函數(shù)，以

29、備主程序調(diào)用。% DecisFun.m function u=DecisFun(k,x) d=2 3 2 4;m=6; if k=4 u=d(k)-x;else u=max(0,d(k)-x):m; End% SubObjFun.mfunction f=SubObjFun(k,x,u) d=2 3 2 4; if u=0 f=0.5*(x+u-d(k);else if u>6f=106;else f=3+u+0.5*(x+u-d(k); endEnd % TransFun.m function s=TransFun(k,x,u)d=2 3 2 4;s=x+u-d(k);在matlab命令

30、空間輸入：x1=0:4;s=nan*ones(5,1);s(1)=0;x=s x1' x1' x1'p_opt,fval=dynprog(x,'DecisFun','SubObjFun','TransFun')運(yùn)行結(jié)果如下：p_opt =1.000005.00009.50002.00003.0000003.000006.000011.00004.0000 fval =20.50004.000000從上面的結(jié)果可知，每個(gè)時(shí)期的最優(yōu)決策為：X1=5,x2=0,x3=6,x4=0。其相應(yīng)的最小總成本為20.5千元。從上面的結(jié)果中

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

32、于裝入第i種物品到第3種物品的總體積； ui表示裝入第i種物品的件數(shù)；可得狀態(tài)轉(zhuǎn)移方程：sk+1 = s ck u tk( )* k k, +1 = t ck uk( )* k允許決策集合為：s tk , k )Ds u( k k,) = uk | 0 ukmin(w vkk最優(yōu)值函數(shù) f s tk k k(,) 表示當(dāng)總重量不超過sk公斤，總體積不超過tk立方米背包裝入第t種物品到第3種物品產(chǎn)生的最大效益值?？傻没痉匠蹋篺 s tk k k(,) =u D s tkmaxk k k(,)( ( )*ck u f s tk+ k+1( k+1, k+1),f v t4 ( 4 , 4 ) =

33、 0 k= 3,2,1下面同樣用計(jì)算程序dynprog1.m求解：在使用此程序先要建立下面3個(gè)M_函數(shù)：% DecisFun1.m functionu1,u2=DecisFun1(k,x1,x2) w=3 4 5;v=8 6 4;u1=0:1:min(x1/w(k),x2/v(k);u2=1;% 因?yàn)檫@里只有一個(gè)決策變量，故令u2恒為1，這樣是程序的需要，% 也可減少計(jì)算量，此時(shí)u2就沒有任何意義，只是一個(gè)虛擬的量% SubObjFun1.m functionf=SubObjFun1(k,x1,x2,u1,u2) c=4 5 6;f=-c(k)*u1; % 求最大值轉(zhuǎn)化為求最小值% TransFun1.m functions=TransFun1(k,x1,x2,u1,u2) w=3 4 5;v=8 6 4;s(1)=x1-u1*w(k);s(2)=x2-u1*v(k);在matlab