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

1、.動態(tài)規(guī)劃方法的matlab實現及其應用（龍京鵬，張華慶，羅明良，劉水林） (南昌航空大學，數學與信息科學學院，江西，南昌)摘要：本文運用matlab語言實現了動態(tài)規(guī)劃的逆序算法，根據狀態(tài)變量的維數，編寫了指標函數最小值的逆序算法遞歸計算程序。兩個實例的應用檢驗了該程序的有效性，同時也表明了該算法程序對眾多類典型的動態(tài)規(guī)劃應用問題尤其是確定離散型的應用問題的通用性，提供了求解各種動態(tài)規(guī)劃問題的有效工具。關鍵詞：動態(tài)規(guī)劃 基本方程的逆序算法 matlab實現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精品.動態(tài)規(guī)劃是一類解決多階段決策問題的數學方法, 在工程技術、科學管理、工農業(yè)生產及軍事等領域都有廣泛的應用。在理論上,動態(tài)規(guī)劃是求解這類問題全局最優(yōu)解的一種有效方法,特別是對于實際中某些非線性規(guī)劃問題可能是最優(yōu)解的唯一方法。然而,動態(tài)規(guī)劃僅僅決多階段決策問題的一種方法,或者說是考查問題的一種途徑,而不是一種具體的算法。就目前而言,動態(tài)規(guī)劃沒有統(tǒng)一的標

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

6、 和“背包問題”1進行了應用與檢驗,實際證明結果是令人滿意的。1 動態(tài)規(guī)劃的基本模型實際中,要構造一個標準的動態(tài)規(guī)劃模型,通常需要采用以下幾個步驟:劃分階段 按照問題的時間或空間特征,把問題分為若干個階段。這些階段必須是有序的或者是可排序的(即無后向性) ,否則,應用無效。選擇狀態(tài) 將問題發(fā)展到各個階段時所處的各種客觀情況用不同的狀態(tài)表示,即稱為狀態(tài)。狀態(tài)的選擇要滿足無后效性和可知性,即狀態(tài)不僅依賴于狀態(tài)的轉移規(guī)律,還依賴于允許決策集合和指標函數結構。確定決策變量與狀態(tài)轉移方程 當過程處于某一階段的某個狀態(tài)時,可以做出不同的決策,描述決策的變量稱為決策變量。在決策過程中,由一個狀態(tài)到另一個狀態(tài)

7、的演變過程稱為狀態(tài)轉移。狀態(tài)轉移就是根據上一階段的狀態(tài)和決策來導出本階段的狀態(tài)。寫出動態(tài)規(guī)劃的基本方程 動態(tài)規(guī)劃的基本方程一般根據實際問題可分為兩種形式,逆序形式和順序形式。這里只考慮逆序形式。動態(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)轉移方程為sk+1 =t s xk k k( , ), 態(tài)處于k 階段的允

8、許決策集合記為d sk k() , v s xk k k(,) 為指標函數。當求解時,由邊界條件從k n= 開始, 由后向前逆推,逐階段求出最優(yōu)決策和過程的最優(yōu)值, 直到最后求出 f s1( 1) 即得到問題的最優(yōu)解。動態(tài)規(guī)劃逆序解法計算程序框圖如下：精品.精品.2 基本方程逆序算法的matlab程序（1）動態(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(,) 。自由始端和終端的動態(tài)

9、規(guī)劃，求指標函數最小值的逆序算法遞歸計算程序：functionp_opt,fval=dynprog(x,decisfun,subobjfu n,transfun,objfun)% x為狀態(tài)變量，一列代表一個階段的狀態(tài)%m_函數decisfun（k,x)表示由階段k的狀態(tài)值x求出相應的允許決策集合% m_函數subobjfun（k,x,u)表示階段k的指標函數% m_函數transfun（k,x,u)是狀態(tài)轉移函數，其中x 是階段k的狀態(tài)值，u是其決策集合% m_函數objfun（v,f)是第k階段到最后階段的指標函數，當objfun（v,f)=v+f時，輸入objfun（v,f)可以省略% 輸

10、出p_opt由4列組成，p_opt=序號組，最優(yōu)軌線組，最優(yōu)策略組，指標函數值組; % 輸出fval是列向量，各元素分別表示p_opt各最優(yōu)策略組對應始端狀態(tài)x的最優(yōu)函數值k=length(x(1,:);% k為階段數x_isnan=isnan(x); t_vubm=inf*ones(size(x);%t_vubm為指標函數值的上限f_opt=nan*ones(size(x);% f_opt為不同階段、狀態(tài)下的最優(yōu)值矩陣，初值為非數 d_opt=f_opt; % d_opt為不同階段不同狀態(tài)下的決策矩陣，初值為非數tmp1=find(x_isnan(:,k);% 找出第k階段狀態(tài)值（不是非數）

11、的下標tmp2=length(tmp1); for i=1:tmp2u=feval(decisfun,k,x(tmp1(i),k);% 求出相應的允許決策向量 tmp3=length(u);for j=1:tmp3% 該for語句是為了求出相應的最有函數值以及最優(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)值及相應的決策值求出下一階段的狀態(tài)值 tmp50=x(:,ii+1)-tmp40;tmp60=find(tmp50=0);% 找出下一階段的狀態(tài)值在x(:

13、,ii+1)的下標 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即為最有函數值矩陣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);% 求出第一階段的指標函數值p_opt(k*(i-1)+1,12 3 4)=1,tmpx(i),tmpd(i),tmpf(i);for ii=2:k% 按順序求出各階段的決策值、狀態(tài)值以及指標函數值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)當狀態(tài)變量是二維時，也即有兩個狀態(tài)變量，此時動態(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(,)此時上面的程序就無能為力了，為此在程序 dynprog.m基礎上進行拓展，我們得到狀態(tài)變量為二維情況下的指標函數最小值的逆序算法遞歸計算程序:dynprog1.m，如下：functionp_opt,fval=dynprog1(x1,x2,decisfun,subobjfun,transfun,objfun)% (x1,x2)為二維狀態(tài)變量，其中x1,x2的取

17、值是相互獨立的，這里矩陣x1與x2的列數應相同，該程序考慮決策變量也是二維%decisfun(k,x1,x2),subobjfun(k,x1,x2,u1, u2),transfun(k,x1,x2,u1,u2)等m_函數的含義與一維的情形一樣，只是它們的參數相應的增加了，objfun函數的含義及參數保持不變% p_opt,fval的含義與一位情形一樣，只是它們的維數增加了% 下面程序的思路與算法同一維基本相同，只是相應矩陣的維數增加了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函數時，就(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時期末庫存量為vk+1 時的存儲費為:h vk k( ) = 0.5*(v x dk + k k) 故第k時期內的總成本為: c x h vk k( ) + k k( ) 則階段指標函數為:精品.v v c x h vk k() = k k() + k k()最優(yōu)值函數 f vk k() 表示從第k階段初始庫存量為vk時到第四階段末庫存量為0時的最小總費用。則有遞推關系式: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),這是因為一方面每階段生產的下限為0；另一方面由于要保證供應，故第k階段末的庫存量vk+1 必須非負，即v x dk + kk 0 ，所以 x d vk kk。vk的取值范圍為0, min(d m dj, k) ,其中v1 = 0,v5 = 0 。j k=為求出該問題的最優(yōu)值，利用上面的計算程序 dynprog.m。根據上面所述的階段指標函數、狀態(tài)轉移函數和遞推關系式，編寫出下面3個m_函數，以備主程序調用。% 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)運行結果如下：p_opt =1.000005.00009.50002.00003.0000003.000006.000011.00004.0000 fval =20.50004.000000從上面的結果可知，每個時期的最優(yōu)決策為：x1=5,x2=0,x3=6,x4=0。其相應的最小總成本為20.5千元。從上面的結果中還可以看出，各個時期初的庫存量分別為：v1=0,v2=2,v3=0,v4=4這里的結果與文1的結果完全符合，這說明該程序是可行的。3.2 二維背包問題有一個人帶一個背包上山，其可攜帶物品重量的限度為10公斤，背包體積限制為

26、22立方米。假設有3種物品可供他選擇裝入背包。已知第i種物品每件重量為w(i)公斤，體積為v(i)立方米，攜帶該物品xi件產生的效益值為c(i)*xi。問此人該如何選擇攜帶物品，才能使產生的效益值最大。其中 w=3 4 5;v=8 6 4;c=4 5 6; 解：用動態(tài)規(guī)劃方法求解，按物品的種類數將該問題分為3各個階段；si表示用于裝入第i種物品到第3種物品的總重量；ti表示用于裝入第i種物品到第3種物品的總體積； ui表示裝入第i種物品的件數；可得狀態(tài)轉移方程：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)值函數 f s tk k k(,) 表示當總重量不超過sk公斤，總體積不超過tk立方米背包裝入第t種物品到第3種物品產生的最大效益值。可得基本方程：f 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下面同樣用計算程序dynprog1.m求解：在使用此程序先要建立下面3個m_函數：% 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;% 因為這里只有一個決策變量，故令u2恒為1，這樣是程序的需要，% 也可減少計算量，此時u2就沒有任何意義，只是一個虛擬的量% subobjfun1.m functionf=subobjfun1(k,x1,x2,u1,u2) c=4 5 6;f=-c(k)*u1; % 求最大值轉化為求最小值% 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