動(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語(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)用問(wèn)題尤其是確定離散型的應(yīng)用問(wèn)題的通用性，提供了求解各種動(dòng)態(tài)規(guī)劃問(wèn)題的有效工具。關(guān)鍵詞：動(dòng)態(tài)規(guī)劃 基本方程的逆序算法 matlab實(shí)現(xiàn)matlab achieve for dynamic programming and its application(jingpenglon

2、g，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 algorith

3、m 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ī)劃是一類解決多階段決策問(wèn)題的數(shù)學(xué)方法, 在工程技術(shù)、科學(xué)管理、工農(nóng)業(yè)生產(chǎn)及軍事等領(lǐng)域都有廣泛的應(yīng)用。在理論上,動(dòng)態(tài)規(guī)劃是求解這類問(wèn)題全局最優(yōu)解的一種有效方法,特別是對(duì)于實(shí)際中某些非線性規(guī)劃問(wèn)題可能是最優(yōu)解的唯一方法。然而,動(dòng)態(tài)規(guī)劃僅僅決多階段決策問(wèn)題的一種方法,或者說(shuō)是考查問(wèn)題的一種途徑,而不是一種具體的算法。就目前而言,動(dòng)態(tài)規(guī)劃沒(méi)有統(tǒng)一的標(biāo)

5、準(zhǔn)模型,其解法也沒(méi)有標(biāo)準(zhǔn)算法,在實(shí)際應(yīng)用中,需要具體問(wèn)題具體分析。動(dòng)態(tài)規(guī)劃模型的求解問(wèn)題是影響動(dòng)態(tài)規(guī)劃理論和方法應(yīng)用的關(guān)鍵所在,而子問(wèn)題的求解和大量結(jié)果的存儲(chǔ)、調(diào)用更是一個(gè)難點(diǎn)所在。然而, 隨著計(jì)算機(jī)技術(shù)的快速發(fā)展,特別是內(nèi)存容量和計(jì)算速度的增加,使求解較小規(guī)模的動(dòng)態(tài)規(guī)劃問(wèn)題成為可能,從而使得動(dòng)態(tài)規(guī)劃的理論和方法在實(shí)際中的應(yīng)用范圍迅速增加。目前,在計(jì)算機(jī)上實(shí)現(xiàn)動(dòng)態(tài)規(guī)劃的一般求解方法并不多見(jiàn),尤其是用來(lái)解決較復(fù)雜的具體問(wèn)題的成果甚少。本文從實(shí)際出發(fā),利用數(shù)學(xué)工具軟件matlab 的強(qiáng)大功能, 對(duì)動(dòng)態(tài)規(guī)劃模型的求解方法做了嘗試,編寫出了動(dòng)態(tài)規(guī)劃逆序算法的matlab程序，并結(jié)合“生產(chǎn)與存儲(chǔ)問(wèn)題”1

6、 和“背包問(wèn)題”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è)步驟:劃分階段 按照問(wèn)題的時(shí)間或空間特征,把問(wèn)題分為若干個(gè)階段。這些階段必須是有序的或者是可排序的(即無(wú)后向性) ,否則,應(yīng)用無(wú)效。選擇狀態(tài) 將問(wèn)題發(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)過(guò)程處于某一階段的某個(gè)狀態(tài)時(shí),可以做出不同的決策,描述決策的變量稱為決策變量。在決策過(guò)程中,由一個(gè)狀態(tài)到另一個(gè)狀態(tài)

7、的演變過(guò)程稱為狀態(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í)際問(wèn)題可分為兩種形式,逆序形式和順序形式。這里只考慮逆序形式。動(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= 開(kāi)始, 由后向前逆推,逐階段求出最優(yōu)決策和過(guò)程的最優(yōu)值, 直到最后求出 f s1( 1) 即得到問(wèn)題的最優(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)

9、規(guī)劃，求指標(biāo)函數(shù)最小值的逆序算法遞歸計(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)可以省略% 輸

10、出p_opt由4列組成，p_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ù)）

11、的下標(biāo)tmp2=length(tmp1); for i=1:tmp2u=feval(decisfun,k,x(tmp1(i),k);% 求出相應(yīng)的允許決策向量 tmp3=length(u);for j=1:tmp3% 該for語(yǔ)句是為了求出相應(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

12、_opt tmp10=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(:

13、,ii+1)的下標(biāo) if isempty(tmp60)if nargin5tmp00=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=;tm

14、pf=; 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),tm精品.pd(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,i

15、i-1,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

16、k k 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í)上面的程序就無(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的取

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*on

18、es(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(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

19、=feval(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),i精品.i),x2(tmp021(t),ii),u1(j),u2(l)

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

22、 f_opt(tmp400(1),mp500(1),ii+1);end if tmp0006在第 k時(shí)期末庫(kù)存量為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階段初始庫(kù)存量為vk時(shí)到第四階段末庫(kù)存量為0時(shí)的最小總費(fèi)用。則有遞推關(guān)系式:f vk k() = max(0,dk vk xkmin ) 6(v vk k() +f vk+1( k+1)f v5 (

23、5 ) = 0, x d v4 = 44其中xk max(0,d vk k),這是因?yàn)橐环矫婷侩A段生產(chǎn)的下限為0；另一方面由于要保證供應(yīng)，故第k階段末的庫(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=為求出該問(wèn)題的最優(yōu)值，利用上面的計(jì)算程序 dynprog.m。根據(jù)上面所述的階段指標(biāo)函數(shù)、狀態(tài)轉(zhuǎn)移函數(shù)和遞推關(guān)系式，編寫出下面3個(gè)m_函數(shù)，以備主程序調(diào)用。% decisfun.m function u=decisfun(k,x) d=2 3 2 4;m=6; if

24、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 u6f=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命令空間輸入：x1=0:4;s=nan*ones(5,1);s(1)=0;x=s x1 x1 x1;p_opt,fval=dynprog(

25、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é)果中還可以看出，各個(gè)時(shí)期初的庫(kù)存量分別為：v1=0,v2=2,v3=0,v4=4這里的結(jié)果與文1的結(jié)果完全符合，這說(shuō)明該程序是可行的。3.2 二維背包問(wèn)題有一個(gè)人帶一個(gè)背包上山，其可攜帶物品重量的限度為10公斤，背包體積限制為

26、22立方米。假設(shè)有3種物品可供他選擇裝入背包。已知第i種物品每件重量為w(i)公斤，體積為v(i)立方米，攜帶該物品xi件產(chǎn)生的效益值為c(i)*xi。問(wèn)此人該如何選擇攜帶物品，才能使產(chǎn)生的效益值最大。其中 w=3 4 5;v=8 6 4;c=4 5 6; 解：用動(dòng)態(tài)規(guī)劃方法求解，按物品的種類數(shù)將該問(wèn)題分為3各個(gè)階段；si表示用于裝入第i種物品到第3種物品的總重量；ti表示用于裝入第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(

27、 k k,) = uk | 0 ukmin(w vkk精品.最優(yōu)值函數(shù) f s tk k k(,) 表示當(dāng)總重量不超過(guò)sk公斤，總體積不超過(guò)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 ) = 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

28、=8 6 4;u1=0:1:min(x1/w(k),x2/v(k);u2=1;% 因?yàn)檫@里只有一個(gè)決策變量，故令u2恒為1，這樣是程序的需要，% 也可減少計(jì)算量，此時(shí)u2就沒(méi)有任何意義，只是一個(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命令空間輸入：a1=0:10;b1=0:22;s1=nan*ones(11,1);s1(1)=10; s2=nan*ones(23,1);s2(1)=22;x1=s1 a1 a1;x2=s2 b1 b1;p_opt,fval=dynprog1(x1,x2,decisfun1