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

动态规划方法的matlab实现及其应用（龙京鹏，张华庆，罗明良，刘水林） (南昌航空大学，数学与信息科学学院，江西，南昌)摘要：本文运用matlab语言实现了动态规划的逆序算法，根据状态变量的维数，编写了指标函数最小值的逆序算法递归计算程序。两个实例的应用检验了该程序的有效性，同时也表明了该算法程序对众多类典型的动态规划应用问题尤其是确定离散型的应用问题的通用性，提供了求解各种动态规划问题的有效工具。关键词：动态规划 基本方程的逆序算法 MATLAB实现MATLAB Achieve For Dynamic Programming and Its Application(JingpengLong，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 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 solution of a variety of dynamic programming problems. Key words:dynamic programming；reverse algorithm；Matlab achievement8动态规划是一类解决多阶段决策问题的数学方法, 在工程技术、科学管理、工农业生产及军事等领域都有广泛的应用。在理论上,动态规划是求解这类问题全局最优解的一种有效方法,特别是对于实际中某些非线性规划问题可能是最优解的唯一方法。然而,动态规划仅仅决多阶段决策问题的一种方法,或者说是考查问题的一种途径,而不是一种具体的算法。就目前而言,动态规划没有统一的标准模型,其解法也没有标准算法,在实际应用中,需要具体问题具体分析。动态规划模型的求解问题是影响动态规划理论和方法应用的关键所在,而子问题的求解和大量结果的存储、调用更是一个难点所在。然而, 随着计算机技术的快速发展,特别是内存容量和计算速度的增加,使求解较小规模的动态规划问题成为可能,从而使得动态规划的理论和方法在实际中的应用范围迅速增加。目前,在计算机上实现动态规划的一般求解方法并不多见,尤其是用来解决较复杂的具体问题的成果甚少。本文从实际出发,利用数学工具软件matlab 的强大功能, 对动态规划模型的求解方法做了尝试,编写出了动态规划逆序算法的matlab程序，并结合“生产与存储问题”1 和“背包问题”1进行了应用与检验,实际证明结果是令人满意的。1 动态规划的基本模型实际中,要构造一个标准的动态规划模型,通常需要采用以下几个步骤:划分阶段 按照问题的时间或空间特征,把问题分为若干个阶段。这些阶段必须是有序的或者是可排序的(即无后向性) ,否则,应用无效。选择状态 将问题发展到各个阶段时所处的各种客观情况用不同的状态表示,即称为状态。状态的选择要满足无后效性和可知性,即状态不仅依赖于状态的转移规律,还依赖于允许决策集合和指标函数结构。确定决策变量与状态转移方程 当过程处于某一阶段的某个状态时,可以做出不同的决策,描述决策的变量称为决策变量。在决策过程中,由一个状态到另一个状态的演变过程称为状态转移。状态转移就是根据上一阶段的状态和决策来导出本阶段的状态。写出动态规划的基本方程 动态规划的基本方程一般根据实际问题可分为两种形式,逆序形式和顺序形式。这里只考虑逆序形式。动态规划基本方程的逆序形式为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 阶段的状态为sk,其决策变量xk表示状 sk的决策,状态转移方程为sk+1 =T s xk k k( , ), 态处于k 阶段的允许决策集合记为D sk k() , v s xk k k(,) 为指标函数。当求解时,由边界条件从k n= 开始, 由后向前逆推,逐阶段求出最优决策和过程的最优值, 直到最后求出 f s1( 1) 即得到问题的最优解。动态规划逆序解法计算程序框图如下：2 基本方程逆序算法的matlab程序（1）动态规划逆序求最小值的基本方程：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(,) 。自由始端和终端的动态规划，求指标函数最小值的逆序算法递归计算程序：functionp_opt,fval=dynprog(x,DecisFun,SubObjFu n,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); for i=1:tmp2u=feval(DecisFun,k,x(tmp1(i),k);% 求出相应的允许决策向量 tmp3=length(u);for j=1:tmp3% 该for语句是为了求出相应的最有函数值以及最优决策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 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);% 由该状态值及相应的决策值求出下一阶段的状态值 tmp50=x(:,ii+1)-tmp40;tmp60=find(tmp50=0);% 找出下一阶段的状态值在x(:,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=;tmpf=; 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);% 求出第一阶段的状态值tmpf(i)=feval(SubObjFun,1,tmpx(i),tmpd(i);% 求出第一阶段的指标函数值p_opt(k*(i-1)+1,12 3 4)=1,tmpx(i),tmpd(i),tmpf(i);for ii=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);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)当状态变量是二维时，也即有两个状态变量，此时动态规划逆序求最小值的基本方程：f sk k() = k k k k k t (min), ( ) (gv s x t u f sk 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基础上进行拓展，我们得到状态变量为二维情况下的指标函数最小值的逆序算法递归计算程序: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,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(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=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=length(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),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);ifisempty(tmp400)&isempty(tmp500) if nargin6tmp000=tmp000+f_opt(tmp400(1),tmp500(1),ii+1); else tmp000=feval(ObjFun,tmp000, 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()最优值函数 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 ( 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=为求出该问题的最优值，利用上面的计算程序 dynprog.m。根据上面所述的阶段指标函数、状态转移函数和递推关系式，编写出下面3个M_函数，以备主程序调用。% 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 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(x,DecisFun,SubObjFun,TransFun)运行结果如下：p_opt =1.000005.00009.50002.00003.0000003.000006.000011.00004.0000 fval =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=3 4 5;v=8 6 4;c=4 5 6; 解：用动态规划方法求解，按物品的种类数将该问题分为3各个阶段；si表示用于装入第i种物品到第3种物品的总重量；ti表示用于装入第i种物品到第3种物品的总体积； ui表示装入第i种物品的件数；可得状态转移方程：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最优值函数 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=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,S ubObjFun1,TransFun1)运行结果如下： p_opt =1 102221-82 4611-53 00010 fval =-13从 上面 的结 果可 知： 最优 装入 方案 为： u1=2,u2=1,u3=0;也即各种物品分别装入2件，1件，0件，此时产生的最大效益值为13。此程序得出的结果与事实相符合，说明此程序是可行的。4 程序使用的几点说明程序dynprog.m只能使用于具体问题中只有一个状态变量的情形，程序dynprog1.m适用于状态变量为二维的情形。这两个程序都要求各阶段状态变量的取值是离散的。要使用好这两个程序，关键要做到以下三点：一、 要掌握动态规划的基本原理与基本概念，能对具体问题写出基本方程的逆序形式，要认真读懂这两个程序，理解程序中每个变量所代表的含义，只有理解了程序，才能更好地使用程序，才能对运行出的结果进行分析；二、 对具体问题要作具体的分析，要用动态规划的方法求解问题，要能够写出状态转移方程、基本方程以及允许决策集合，并要根据这些方程在matlabruan软件中建立四个M_函数，以备主程序调用；三、 每个阶段的状态变量的取值一定要合理地离散化。5 结束语本文中运用matlab语言实现了动态规划（包括状态变量二维情形）的逆序算法，两个实例的应用检验了该程序的有效性，同时也表明了该算法程序对众多类典型的动态规划应用问题尤其是确定离散型的应用问题的通用性。本文拓展了文2中有关动态规划逆序算法的 matlab程序实现，由原来的一维情形扩展到二维情形，这是一个进步。当问题的阶段数和各阶段的状态数较小时，这两个程序能够运行出结果，但当问题的阶段数较和各阶段的状态数较大时，这两个程序运行时就要花费较长的时间，有时甚至是运行不出结果来，因为花费的时间太长。参考文献：1运筹学教材编写组，运筹学M 北京:清华大学出版社，2005.62 胡良剑 丁晓东 孙晓君，数学实验-使用MATLABM 上海:上海科技出版社，20013 张志涌，精通MATLAB6.5M 北京:北京航空航天大学出版社，20034 刘保柱 苏彦华 张宏林， MATLAB7.0从入门到精通（修订版） 北京：人民邮电出版社，2010 5771001803090012095 5790368228596330825771001803090012386 5761373997357606965771001803090013594 5780775799025155125771001803090012387 5771649826018180515771001803090012138 5721311921589183265771001803090012359 5790368223610760535771001803090012356 5761352861437917425771001803090012355 57508786970469327917088100343355274 10122994432583337917088100343355275 10186673293883200817088100343356107 10158115250150052217088100343356108 10100018005987173217088100343354295 10107419414268701717088100343356184 10187866086962880217088100343356185 10177583117408667417088100343356109 10108601437357284617088100343356110 10115220721601491617088100343355237 10102704160570270917088100343355238 10122936486142541417088100343356169 10186220440263571817088100343354928 101760654089788804