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

...wd......wd......wd...动态规划方法的matlab实现及其应用〔龙京鹏，张华庆，罗明良，刘水林〕(南昌航空大学，数学与信息科学学院，江西，南昌)摘要：本文运用matlab语言实现了动态规划的逆序算法，根据状态变量的维数，编写了指标函数最小值的逆序算法递归计算程序。两个实例的应用检验了该程序的有效性，同时也说明了该算法程序对众多类典型的动态规划应用问题尤其是确定离散型的应用问题的通用性，提供了求解各种动态规划问题的有效工具。关键词：动态规划根本方程的逆序算法MATLAB实现MATLABAchieveForDynamicProgrammingandItsApplication(JingpengLong，HuaqingZhang，MingliangLuo，ShuilinLiu)〔SchoolofMathematicsandInformationScience,NanchangHangkongUniversity,Nanchang,China〕Abstract:Thisarticleachievesthereversealgorithmofdynamicprogrammingbyusingthematlablanguage，andpreparestherecursivecalculationprogramofreversealgorithmwhichthetargetfunctionvalueisthesmallest.Theapplicationoftwoexamplesshowthattheprogramiseffective，andthisalgorithmprogramisgeneraltomanytypicalapplicationofdynamicprogramming,especiallytheapplicationofdeterministicdiscrete.Thisalgorithmprogramprovidesaeffectivetooltothesolutionofavarietyofdynamicprogrammingproblems.Keywords:dynamicprogramming；reversealgorithm；Matlabachievement动态规划是一类解决多阶段决策问题的数学方法,在工程技术、科学管理、工农业生产及军事等领域都有广泛的应用。在理论上,动态规划是求解这类问题全局最优解的一种有效方法,特别是对于实际中某些非线性规划问题可能是最优解的唯一方法。然而,动态规划仅仅决多阶段决策问题的一种方法,或者说是考察问题的一种途径,而不是一种具体的算法。就目前而言,动态规划没有统一的标准模型,其解法也没有标准算法,在实际应用中,需要具体问题具体分析。动态规划模型的求解问题是影响动态规划理论和方法应用的关键所在,而子问题的求解和大量结果的存储、调用更是一个难点所在。然而,随着计算机技术的快速开展,特别是内存容量和计算速度的增加,使求解较小规模的动态规划问题成为可能,从而使得动态规划的理论和方法在实际中的应用范围迅速增加。目前,在计算机上实现动态规划的一般求解方法并不多见,尤其是用来解决较复杂的具体问题的成果甚少。本文从实际出发,利用数学工具软件matlab的强大功能,对动态规划模型的求解方法做了尝试,编写出了动态规划逆序算法的matlab程序，并结合“生产与存储问题〞[1]和“背包问题〞[1]进展了应用与检验,实际证明结果是令人满意的。1动态规划的根本模型实际中,要构造一个标准的动态规划模型,通常需要采用以下几个步骤:①划分阶段按照问题的时间或空间特征,把问题分为假设干个阶段。这些阶段必须是有序的或者是可排序的(即无后向性),否则,应用无效。②选择状态将问题开展到各个阶段时所处的各种客观情况用不同的状态表示,即称为状态。状态的选择要满足无后效性和可知性,即状态不仅依赖于状态的转移规律,还依赖于允许决策集合和指标函数构造。③确定决策变量与状态转移方程当过程处于某一阶段的某个状态时,可以做出不同的决策,描述决策的变量称为决策变量。在决策过程中,由一个状态到另一个状态的演变过程称为状态转移。状态转移就是根据上一阶段的状态和决策来导出本阶段的状态。④写出动态规划的根本方程动态规划的根本方程一般根据实际问题可分为两种形式,逆序形式和顺序形式。这里只考虑逆序形式。动态规划根本方程的逆序形式为fskk( )=optgvsx{(kkk( , )+fsk+1(k+1))}xDsk∈kk( )knn=,−1,⋯,1边界条件fsn+1(n+1)=0或fsvsxnn( )=nnn( , )其中第k阶段的状态为sk,其决策变量xk表示状sk的决策,状态转移方程为sk+1=Tsxkkk(,),态处于k阶段的允许决策集合记为Dskk( ),vsxkkk( , )为指标函数。当求解时,由边界条件从kn=开场,由后向前逆推,逐阶段求出最优决策和过程的最优值,直到最后求出fs1(1)即得到问题的最优解。动态规划逆序解法计算程序框图如下：2根本方程逆序算法的matlab程序〔1〕动态规划逆序求最小值的根本方程：fskk( )=xDskmin{∈kk( )gvsx(kkk( , )+fsk+1(k+1))}knn=,−1,⋯,1边界条件fsvsxnn( )=nnn( , )sk+1=Tsxkkk( , )。自由始端和终端的动态规划，求指标函数最小值的逆序算法递归计算程序：function[p_opt,fval]=dynprog(x,DecisFun,SubObjFun,TransFun,ObjFun)%x为状态变量，一列代表一个阶段的状态%M_函数DecisFun〔k,x)表示由阶段k的状态值x求出相应的允许决策集合%M_函数SubObjFun〔k,x,u)表示阶段k的指标函数%M_函数TransFun〔k,x,u)是状态转移函数，其中x是阶段k的状态值，u是其决策集合%M_函数ObjFun〔v,f)是第k阶段到最后阶段的指标函数，当ObjFun〔v,f)=v+f时，输入ObjFun〔v,f)可以省略%输出p_opt由4列组成，p_opt=[序号组，最优轨线组，最优策略组，指标函数值组];%输出fval是列向量，各元素分别表示p_opt各最优策略组对应始端状态x的最优函数值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为不同阶段、状态下的最优值矩阵，初值为非数d_opt=f_opt;%d_opt为不同阶段不同状态下的决策矩阵，初值为非数tmp1=find(x_isnan(:,k)); %找出第k阶段状态值〔不是非数〕的下标tmp2=length(tmp1);fori=1:tmp2u=feval(DecisFun,k,x(tmp1(i),k));%求出相应的允许决策向量tmp3=length(u);forj=1:tmp3 %该for语句是为了求出相应的最有函数值以及最优决策tmp=feval(SubObjFun,k,x(tmp1(i),k),u(j));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_opttmp10=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~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(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(:,ii)-tmpx(i);tmp2=find(tmp1==0);if~isempty(tmp2)tmpd(i)=d_opt(tmp2(1),ii);endtmpf(i)=feval(SubObjFun,ii,tmpx(i),tmpd(i));p_opt(k*(i-1)+ii,[1 2 34])=[ii,tmpx(i),tmpd(i),tmpf(i)];endend〔2)当状态变量是二维时，也即有两个状态变量，此时动态规划逆序求最小值的根本方程：fskk( )=kkkkkt∈(min),∈(){(gvsxtufskkkkk( , , , )+k+1(k+1))}xDsuUtknn=,−1,⋯,1边界条件fsvsxtunn( )=nnnnn( , , , ),sk+1=Tsxkkk( , ),tk+1=Ptukkk( , )此时上面的程序就无能为力了，为此在程序dynprog.m根基上进展拓展，我们得到状态变量为二维情况下的指标函数最小值的逆序算法递归计算程序:dynprog1.m，如下：function[p_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,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(tmp11);tmp21=find(x2_isnan(:,k));tmp22=length(tmp21);fori=1:tmp12fort=1:tmp22[u1,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;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:tmp022[u1,u2]=feval(DecisFun,ii,x1%假设决策变量为一维，那么在定义DecisFun函数时，就(tmp011(i),ii),x2(tmp021(t),ii));令u2恒为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(tmp300==0);if~isempty(tmp400)&~isempty(tmp500)ifnargin<6tmp000=tmp000+f_opt(tmp400(1),tmp500(1),ii+1);elsetmp000=feval(ObjFun,tmp000,f_opt(tmp400(1),mp500(1),ii+1));endiftmp000<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;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_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(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);tmpx2(j)=u(2);tmp1=x1(:,ii)-tmpx1(i);tmp2=x2(:,ii)-tmpx2(j);tmp3=find(tmp1==0);tmp4=find(tmp2==0);if~isempty(tmp3)&~isempty(tmp4)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_opt(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。试问该厂应假设何安排各个时期的生产与库存，才能在满足市场需要的条件下，使总本钱最小。解：用动态规划方法求解，按四个时期将该问题分为四个阶段；记Vk为状态变量，它表示第k阶段开场时的库存量；记Xk为决策变量，它表示第k阶段的生产量；可得状态转移方程：Vk+1=Vk+Xk-dk,k=1,2,3,4由题意知，在第k时期内的生产本钱为:⎧0...................当xk=0⎫cxkk()=⎨3+1*xk.......当xk=1,2,...6⎪⎬⎪⎪⎩∞..................当xk>6 ⎪⎭在第k时期末库存量为vk+1时的存储费为:hvkk()=0.5*(vxdk+−kk)故第k时期内的总本钱为:cxhvkk()+kk()则阶段指标函数为:Vvcxhvkk( )=kk( )+kk( )最优值函数fVkk( )表示从第k阶段初始库存量为Vk时到第四阶段末库存量为0时的最小总费用。则有递推关系式:⎧⎪fvkk( )=max(0,dkvkxkmin−)≤≤6(Vvkk( )+fvk+1(k+1))⎨⎪⎩fv5(5)=0,xdv4=−4 4其中xk≥max(0,dvk−k),这是因为一方面每阶段生产的下限为0；另一方面由于要保证供应，故第k阶段末的库存量vk+1必须非负，即vxdk+−k k≥0，所以xdvk≥−k k。vk的取值范围为[0,min(∑dmdj,−k)],其中v1=0,v5=0。jk=为求出该问题的最优值，利用上面的计算程序dynprog.m。根据上面所述的阶段指标函数、状态转移函数和递推关系式，编写出下面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=10^6;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: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.0000fval=20.50004.000000从上面的结果可知，每个时期的最优决策为：X1=5,x2=0,x3=6,x4=0。其相应的最小总本钱为20.5千元。从上面的结果中还可以看出，各个时期初的库存量分别为：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表示装入第i种物品的件数；可得状态转移方程：sk+1=−sckutk ()*kk,+1=−tckuk ()*k允许决策集合为：stk,k)}Dsu(kk, )={uk|0≤≤uk min(wvk k最优值函数fstkkk( , )表示当总重量不超过sk公斤，总体积不超过tk立方米背包装入第t种物品到第3种物品产生的最大效益值。可得根本方程：⎧⎪fstkkk( , )=uDstk∈maxkkk( , )(()*ckufstk+k+1(k+1,k+1)),⎨⎪⎩fvt4(4,4)=0k=3,2,1下面同样用计算程序dynprog1.m求解：在使用此程序先要建设下面3个M_函数：%DecisFun1.mfunction[u1,u2]=DecisFun1(k,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,u