动态规划方法的matlab实现及其应用

1、动态规划方法的matlab实现及其应用(龙京鹏，张华庆，罗明良，刘水林)(南昌航空大学，数学与信息科学学院，江西，南昌)摘要：本文运用matlab语言实现了动态规划的逆序算法，根据状态变量的维数，编写了指标函数最小值的逆序算法递归计算程序。两个实例的应用检验了该程序的有效性，同时也表明了该算法程序对众多类典型的动态规划应用问题尤其是确定离散型的应用问题的通用性，提供了求解各种动态规划问题的有效工具。关键词：动态规划基本方程的逆序算法MATLAB实现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动态规划是一类解决多阶段决策问题的数学方法，在工程技术、科学管理、工农业生产及军事等领域都有广泛的应用。在理论上，动态规划是求解这类问题全局最优解的一种有效方法，特别是对于实际中某些非线性规划问题可能是最优解的唯一方法。然而，动态规划仅仅决多阶段决策问题的一种方法，或者说是考查问题的一种途径，而不是一种具体的算法。就目前而言，动态规划没有统一的标准模型，其解法也没有标准算法，在实际应用中，需要具体问题具体分析。动态规划模型的求解问题是影响动态规划理论和方法应用的关键所在，而子问题的求解和大量结果的存储、调用更是

5、一个难点所在。然而，随着计算机技术的快速发展，特别是内存容量和计算速度的增加，使求解较小规模的动态规划问题成为可能，从而使得动态规划的理论和方法在实际中的应用范围迅速增加。目前，在计算机上实现动态规划的一般求解方法并不多见，尤其是用来解决较复杂的具体问题的成果甚少。本文从实际出发，利用数学工具软件matlab的强大功能，对动态规划模型的求解方法做了尝试，编写出了动态规划逆序算法的matlab程序，并结合“生产与存储问题”1和“背包问题”1进行了应用与检验，实际证明结果是令人满意的。1动态规划的基本模型实际中，要构造一个标准的动态规划模型，通常需要采用以下几个步骤：划分阶段按照问题的时间或空间特

6、征，把问题分为若干个阶段。这些阶段必须是有序的或者是可排序的(即无后向性)，否则，应用无效。选择状态将问题发展到各个阶段时所处的各种客观情况用不同的状态表示，即称为状态。状态的选择要满足无后效性和可知性，即状态不仅依赖于状态的转移规律，还依赖于允许决策集合和指标函数结构。确定决策变量与状态转移方程当过程处于某一阶段的某个状态时，可以做出不同的决策，描述决策的变量称为决策变量。在决策过程中，由一个状态到另一个状态的演变过程称为状态转移。状态转移就是根据上一阶段的状态和决策来导出本阶段的状态。写出动态规划的基本方程动态规划的基本方程一般根据实际问题可分为两种形式，逆序形式和顺序形式。这里只考虑逆序

7、形式。动态规划基本方程的逆序形式为fskk()=optgvsx(kkk(,)+fsk+1(k+1)xDskkk()knn=,-1,1边界条件fsn+1(n+1)=0或fsvsXn()=nnn(,)其中第k阶段的状态为sk,其决策变量xk表示状Sk的决策，状态转移方程为Sk+1=TSXkk(,),态处于k阶段的允许决策集合记为DSkk(),vSXkk(,)为指标函数。当求解时，由边界条件从kn=开始，由后向前逆推，逐阶段求出最优决策和过程的最优值，直到最后求出fsi(i)即得到问题的最优解。动态规划逆序解法计算程序框图如下：2基木方程逆序算法的matlab程序(1二动态关tks)=O(I划逆序求

8、最小值的基本方程:xdskminekk()gvsx(kkk(sk+1(k+1)knknn=,-1,1边界条1斗fsvsxn()=nnn(,)生+1kk(自由始端和终端的动态规划，求指标函数最小值的逆序算法涕归计算程序:function_,SubObjFu齐缶J二口9t,出工第器熄juf黑4S2IM与珑就盘一段的状由阶思态k的状态值x求相应的允许决策集合x,u)表示阶段k的指标函,u)是状态转移函数，其中%M_函数SubObjFun(k存储耳aj和m茶港端屋段4状态叫最优值矩阵初值l_opt态下的决策矩阵，初值为非数为不同阶段不同状tmp1=fing(x_isnan(:国);状态值布是非数是%找

9、出第k阶段下标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);%求出相应的允许决策向量tmp3=length(u);forj=1:tmp3目应的最肓函数值以及最优决策语句是为了求)+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);%由该状态值及相应的决策值求出下一阶段的状态值tmp50=x(:,ii+1)-tmp40;tmp60=find(tmp50=0);%找出下一阶段的状态值在x(:,ii+1)的下标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即为最有函数值矩阵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);%求出第一阶段的状态值tmpf(i)=feval(SubObjFun,1,tmpx(i),tmpd(i);%求出第一阶段的指标函数值p_opt(k*(i-1)+1,1234)=1,tmpx(i),tmpd(i),tmpf(i);forii=2:k%按顺序求出各阶段的决策值、状态值以及指标函数值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)当状态变量是二维时，也即有两个状态变量，此时动态规划逆序求最小值的基本方程：fSkk()=kkkkkte(min),e()(gvsxtufSkkkk(,xDsuUtknn=,-1,1边界条件fSvSxtLIn()=nnnnn(,，)，sk+1=

14、Tsxkkk(,),tk+1=Ptukkk(,)此时上面的程序就无能为力了，为此在程序dynprog.m基础上进行拓展，我们得到状态变量为二维情况下的指标函数最小值的逆序算法递归计算程序：dynprog1.m,如下：functionp_opt,fval=dynprog1(x1,x2,DecisFun,SubObjFun,TransFun,ObjFun)%(x1,x2)为二维状态变量，其中x1,x2的取值是相互独立的，这里矩阵x1与x2的列数应相同，该程序考虑决策变量也是二维%DecisFun(k,x1,x2),SubObjFun(k,x1,x2,u1,u2),TransFun(k,x1,x2,

15、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*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函数时，就令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实例应用3.1生产与存储问题时期（k）1233需求量（dk）2324某工厂要对一种产品制定今后四个时期的生产计划，据估计在今后四个时期内，市场对于该产品的需求量如下表所示：假定该厂生产每批产品的固定成本为3千克，若不生产就为0;每单位产品成本为1千元；每个时期生产能力所允许的最大生产量不超过6各单位；每个时期末未能售出的产品，每单位需付存储费0.5千元。还假定在第一个时期的初始库存量为0,第四个是期末的库存量也为0。

24、试问该厂应如何安排各个时期的生产与库存，才能在满足市场需要的条件下，使总成本最小。解：用动态规划方法求解，按四个时期将该问题分为四个阶段；记Vk为状态变量，它表示第k阶段开始时的库存量；记Xk为决策变量，它表示第k阶段的生产量；可得状态转移方程：V.尸M+X-dk,k=1,2,3,4由题意知，在第k时期内的生产成本为：?0当Xk=0?cxkk()=?3+1*Xk.当Xk=1,2,.6?oo当Xk>6?在第k时期末库存量为Vk+1时的存储费为：hVkk()=0.5*(vx&+-kk)故第k时期内的总成本为：cXhVkk()+kk()则阶段指标函数为：Vvcxhvkk()=kk()+

25、kk()最优值函数fVkk()表示从第k阶段初始库存量为Vk时到第四阶段末库存量为0时的最小总费用。则有递推关系式：?fVkk()=maX(0,dkvkXkmih-尸6(VVkk()+fVk+1(k+1)?V5(5)=0,XdV4=-44其中Xk>maX(0，4k),这是因为一方面每阶段生产的下限为0；另一方面由于要保证供应，故第k阶段末的库存量Vk+1必须非负，即VXdk+-kk>0,所以XdVk>-kkoVk的取值范围为0,min(Xdmd,-k),其中V1=0,V5=0。jk=为求出该问题的最优值，利用上面的计算程序dynprog.m。根据上面所述的阶段指标函数、状态转

26、移函数和递推关系式，编写出下面3个M画数，以备主程序调用。%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')运行结果如下：p_opt=1.000005.00009.50002.00003.0000003.000006.000011.00004.00004.000000fval=20.5000从上面的结果可知，每个时期的最优决策为：X1=5,x2=0,x3=6,x4=0。其相应的最小总成本为20.5千元。从上面的结果中还可以看出，各个时期初的库存量

28、分别为：V1=0,v2=2,v3=0,v4=4这里的结果与文1的结果完全符合，这说明该程序是可行的。3.2二维背包问题有一个人带一个背包上山，其可携带物品重量的限度为10公斤，背包体积限制为22立方米。假设有3种物品可供他选择装入背包。已知第i种物品每件重量为w(i)公斤，体积为v(i)立方米，携带该物品xi件产生的效益值为c(i)*xi。问此人该如何选择携带物品，才能使产生的效益值最大。其中w=345;v=864;c=456;解：用动态规划方法求解，按物品的种类数将该问题分为3各个阶段；si表示用于装入第i种物品到第3种物品的总重量；ti表示用于装入第i种物品到第3种物品的总体积；ui表示装

29、入第i种物品的件数；可得状态转移方程：Sk+1=-sckutk()*kk,+1=-tckUk()*k允许决策集合为：s告,DDsu(kk,)=Uk|0<ukmin(wvkk最优值函数fstkkk(,)表示当总重量不超过sk公斤，总体积不超过tk立方米背包装入第t种物品到第3种物品产生的最大效益值。可得基本方程：?stkkk(,)=udstkmaxkkk()()*ckufsk+k+1(k+1,k+1),?vt4(4,4)=0k=3,2,1下面同样用计算程序dynprog1.m求解：在使用此程序先要建立下面3个M画数：%DecisFun1.mfunctionu1,u2=DecisFun1(k

30、,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;%求最大值转化为求最小值%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