动态规划方法matlab实现应用

1、- . z.动态规划方法的matlab实现及其应用龙京鹏，华庆，罗明良，水林 (航空大学，数学与信息科学学院，)摘要：本文运用matlab语言实现了动态规划的逆序算法，根据状态变量的维数，编写了指标函数最小值的逆序算法递归计算程序。两个实例的应用检验了该程序的有效性，同时也说明了该算法程序对众多类典型的动态规划应用问题尤其是确定离散型的应用问题的通用性，提供了求解各种动态规划问题的有效工具。关键词：动态规划 根本方程的逆序算法 MATLAB实现MATLAB Achieve For Dynamic Programming and Its Application(JingpengLong，Huaq

2、ingZhang，MingliangLuo，ShuilinLiu)School of Mathematics and Information Science,Nanchang HangkongUniversity,Nanchang,ChinaAbstract:This article achieves the reverse algorithm of dynamic programming by using the matlab language，and prepares the recursive calculation program of reverse algorithm which

3、thetargetfunctionvalueisthesmallest.Theapplicationoftwoe*amplesshowthattheprogram 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 solutio

4、n of a variety of dynamic programming problems. Key words:dynamic programming；reverse algorithm；Matlab achievement-. z动态规划是一类解决多阶段决策问题的数学方法, 在工程技术、科学管理、工农业生产及军事等领域都有广泛的应用。在理论上,动态规划是求解这类问题全局最优解的一种有效方法,特别是对于实际中*些非线性规划问题可能是最优解的唯一方法。然而,动态规划仅仅决多阶段决策问题的一种方法,或者说是考察问题的一种途径,而不是一种具体的算法。就目前而言,动态规划没有统一的标准模型,其解法

5、也没有标准算法,在实际应用中,需要具体问题具体分析。动态规划模型的求解问题是影响动态规划理论和方法应用的关键所在,而子问题的求解和大量结果的存储、调用更是一个难点所在。然而, 随着计算机技术的快速开展,特别是存容量和计算速度的增加,使求解较小规模的动态规划问题成为可能,从而使得动态规划的理论和方法在实际中的应用围迅速增加。目前,在计算机上实现动态规划的一般求解方法并不多见,尤其是用来解决较复杂的具体问题的成果甚少。本文从实际出发,利用数学工具软件matlab 的强大功能, 对动态规划模型的求解方法做了尝试,编写出了动态规划逆序算法的matlab程序，并结合生产与存储问题1 和背包问题1进展了应

6、用与检验,实际证明结果是令人满意的。1 动态规划的根本模型实际中,要构造一个标准的动态规划模型,通常需要采用以下几个步骤:划分阶段 按照问题的时间或空间特征,把问题分为假设干个阶段。这些阶段必须是有序的或者是可排序的(即无后向性) ,否则,应用无效。选择状态 将问题开展到各个阶段时所处的各种客观情况用不同的状态表示,即称为状态。状态的选择要满足无后效性和可知性,即状态不仅依赖于状态的转移规律,还依赖于允许决策集合和指标函数构造。确定决策变量与状态转移方程 当过程处于*一阶段的*个状态时,可以做出不同的决策,描述决策的变量称为决策变量。在决策过程中,由一个状态到另一个状态的演变过程称为状态转移。

7、状态转移就是根据上一阶段的状态和决策来导出本阶段的状态。写出动态规划的根本方程 动态规划的根本方程一般根据实际问题可分为两种形式,逆序形式和顺序形式。这里只考虑逆序形式。动态规划根本方程的逆序形式为f sk k() = opt gv s * ( k k k(,)+f sk+1( k+1) * D sk k k()k nn= , 1,1边界条件f sn+1( n+1) = 0或 f s v s *n n() = n n n(,)其中第k 阶段的状态为sk,其决策变量*k表示状 sk的决策,状态转移方程为sk+1 =T s *k k k( , ), 态处于k 阶段的允许决策集合记为D sk k()

8、 , v s *k k k(,) 为指标函数。当求解时,由边界条件从k n= 开场, 由后向前逆推,逐阶段求出最优决策和过程的最优值, 直到最后求出 f s1( 1) 即得到问题的最优解。动态规划逆序解法计算程序框图如下：2 根本方程逆序算法的matlab程序1动态规划逆序求最小值的根本方程：f sk k() = * D skmin k k() gv s *( k k k(,) +f sk+1( k+1)k nn= , 1,1边界条件 f s v s *n n() = n n n(,) sk+1 =T s *k k k(,) 。自由始端和终端的动态规划，求指标函数最小值的逆序算法递归计算程序：

9、functionp_opt,fval=dynprog(*,DecisFun,SubObjFu n,TransFun,ObjFun)% *为状态变量，一列代表一个阶段的状态%M_函数DecisFunk,*)表示由阶段k的状态值*求出相应的允许决策集合% M_函数SubObjFunk,*,u)表示阶段k的指标函数% M_函数TransFunk,*,u)是状态转移函数，其中* 是阶段k的状态值，u是其决策集合% M_函数ObjFunv,f)是第k阶段到最后阶段的指标函数，当ObjFunv,f)=v+f时，输入ObjFunv,f)可以省略% 输出p_opt由4列组成，p_opt=序号组，最优轨线组，最

10、优策略组，指标函数值组; % 输出fval是列向量，各元素分别表示p_opt各最优策略组对应始端状态*的最优函数值k=length(*(1,:);% k为阶段数*_isnan=isnan(*); t_vubm=inf*ones(size(*);%t_vubm为指标函数值的上限f_opt=nan*ones(size(*);% f_opt为不同阶段、状态下的最优值矩阵，初值为非数 d_opt=f_opt; % d_opt为不同阶段不同状态下的决策矩阵，初值为非数tmp1=find(*_isnan(:,k);% 找出第k阶段状态值不是非数的下标tmp2=length(tmp1); for i=1:t

11、mp2u=feval(DecisFun,k,*(tmp1(i),k);% 求出相应的允许决策向量 tmp3=length(u);for j=1:tmp3% 该for语句是为了求出相应的最有函数值以及最优决策tmp=feval(SubObjFun,k,*(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(*_isnan(:,ii);t

12、mp20=length(tmp10);for i=1:tmp20u=feval(DecisFun,ii,*(tmp10(i),ii);tmp30=length(u); for j=1:tmp30tmp00=feval(SubObjFun,ii,*(tmp10(i),ii),u(j);tmp40=feval(TransFun,ii,*(tmp10(i),ii),u(j);% 由该状态值及相应的决策值求出下一阶段的状态值 tmp50=*(:,ii+1)-tmp40;tmp60=find(tmp50=0);% 找出下一阶段的状态值在*(:,ii+1)的下标 if isempty(tmp60)if n

13、argin5tmp00=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(*_isnan(:,1),1);% fval即为最有函数值矩阵p_opt=;tmp*=;tmpd=;tmpf=; tmp0=find(*_isnan(:,1);tmp

14、01=length(tmp0); for i=1:tmp01 tmpd(i)=d_opt(tmp0(i),1);% 求出第一阶段的决策值tmp*(i)=*(tmp0(i),1);% 求出第一阶段的状态值tmpf(i)=feval(SubObjFun,1,tmp*(i),tmpd(i);% 求出第一阶段的指标函数值p_opt(k*(i-1)+1,12 3 4)=1,tmp*(i),tmpd(i),tmpf(i);for ii=2:k% 按顺序求出各阶段的决策值、状态值以及指标函数值tmp*(i)=feval(TransFun,ii-1,tmp*(i),tmpd (i);tmp1=*(:,ii)-

15、tmp*(i);tmp2=find(tmp1=0); if isempty(tmp2) tmpd(i)=d_opt(tmp2(1),ii);end tmpf(i)=feval(SubObjFun,ii,tmp*(i),tmpd(i);p_opt(k*(i-1)+ii,1234)=ii,tmp*(i),tmpd(i),tmpf(i); endend2)当状态变量是二维时，也即有两个状态变量，此时动态规划逆序求最小值的根本方程：f sk k() = k k k k k t (min), ( ) (gv s * t u f sk k k k k(,)+ k+1( k+1)* D s u U tk n

16、n= , 1,1边界条件f s v s * t un n() = n n n n n(,),sk+1 =T s *k k k(,),tk+1 =P t uk k k(,)此时上面的程序就无能为力了，为此在程序 dynprog.m根底上进展拓展，我们得到状态变量为二维情况下的指标函数最小值的逆序算法递归计算程序:dynprog1.m，如下：functionp_opt,fval=dynprog1(*1,*2,DecisFun,SubObjFun,TransFun,ObjFun)% (*1,*2)为二维状态变量，其中*1,*2的取值是相互独立的，这里矩阵*1与*2的列数应一样，该程序考虑决策变量也是

17、二维%DecisFun(k,*1,*2),SubObjFun(k,*1,*2,u1, u2),TransFun(k,*1,*2,u1,u2)等M_函数的含义与一维的情形一样，只是它们的参数相应的增加了，ObjFun函数的含义及参数保持不变% p_opt,fval的含义与一位情形一样，只是它们的维数增加了% 下面程序的思路与算法同一维根本一样，只是相应矩阵的维数增加了k1,k=size(*1);k2,k=size(*2); *1_isnan=isnan(*1);*2_isnan=isnan(*2); t_vubm=inf*ones(k1,k2,k);f_opt=nan*ones( k1,k2,k

18、);d_opt1=f_opt;d_opt2=f_opt;tmp11=find(*1_isnan(:,k);tmp12=length(t mp11);tmp21=find(*2_isnan(:,k);tmp22=length(tmp21); for i=1:tmp12 for t=1:tmp22u1,u2=feval(DecisFun,k,*1(tmp11(i),k),*2(tmp21(t),k); tmp13=length(u1);tmp14=length(u2);for j=1:tmp13 for l=1:tmp14tmp=feval(SubObjFun,k,*1(tmp11(i),k),*

19、2(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(*1_isnan(:,ii); tmp021=find(*2_isnan(:,ii); tmp012=length(tmp011); tmp022=length(tmp021)

20、; for i=1:tmp012 for t=1:tmp022u1,u2=feval(DecisFun,ii,*1% 假设决策变量为一维，则在定义DecisFun函数时，就(tmp011(i),ii),*2(tmp021(t),ii); 令u2恒为1tmp013=length(u1); tmp014=length(u2); for j=1:tmp013for l=1:tmp014tmp000=feval(SubObjFun,ii,*1(tmp011(i),ii),*2(tmp021(t),ii),u1(j),u2(l); tmp100=feval(TransFun,ii,*1(tmp011(i

21、),ii),*2(tmp021(t),ii),u1(j),u2(l); tmp200=*1(:,ii+1)-tmp100(1); tmp300=*2(:,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

22、 tmp0006在第 k时期末库存量为vk+1 时的存储费为:h vk k( ) = 0.5*(v * dk + k k) 故第k时期的总本钱为: c * h vk k( ) + k k( ) 则阶段指标函数为:V v c * h vk k() = k k() + k k()最优值函数 f Vk k() 表示从第k阶段初始库存量为Vk时到第四阶段末库存量为0时的最小总费用。则有递推关系式:f vk k() = ma*(0,dk vk *kmin ) 6(V vk k() +f vk+1( k+1)f v5 ( 5 ) = 0, * d v4 = 44其中*k ma*(0,d vk k),这是因

23、为一方面每阶段生产的下限为0；另一方面由于要保证供给，故第k阶段末的库存量vk+1 必须非负，即v * dk + kk 0 ，所以 * 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,*) d=2 3 2 4;m=6; if k=4 u=d(k)-*;else u=ma*(0,d(k)-*):m; End% S

24、ubObjFun.mfunction f=SubObjFun(k,*,u) d=2 3 2 4; if u=0 f=0.5*(*+u-d(k);else if u6f=106;else f=3+u+0.5*(*+u-d(k); endEnd % TransFun.m function s=TransFun(k,*,u)d=2 3 2 4;s=*+u-d(k);在matlab命令空间输入：*1=0:4;s=nan*ones(5,1);s(1)=0;*=s *1 *1 *1;p_opt,fval=dynprog(*,DecisFun,SubObjFun,TransFun)运行结果如下：p_opt

25、=1.000005.00009.50002.00003.0000003.000006.000011.00004.0000 fval =20.50004.000000从上面的结果可知，每个时期的最优决策为：*1=5,*2=0,*3=6,*4=0。其相应的最小总本钱为20.5千元。从上面的结果中还可以看出，各个时期初的库存量分别为：V1=0,v2=2,v3=0,v4=4这里的结果与文1的结果完全符合，这说明该程序是可行的。3.2 二维背包问题有一个人带一个背包上山，其可携带物品重量的限度为10公斤，背包体积限制为22立方米。假设有3种物品可供他选择装入背包。第i种物品每件重量为w(i)公斤，体积为

26、v(i)立方米，携带该物品*i件产生的效益值为c(i)*i。问此人该如何选择携带物品，才能使产生的效益值最大。其中 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(,)

27、 表示当总重量不超过sk公斤，总体积不超过tk立方米背包装入第t种物品到第3种物品产生的最大效益值。可得根本方程：f s tk k k(,) =u D s tkma*k 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,*1,*2) w=3 4 5;v=8 6 4;u1=0:1:min(*1/w(k),*2/v(k);u2=1;% 因为这里只有一个

28、决策变量，故令u2恒为1，这样是程序的需要，% 也可减少计算量，此时u2就没有任何意义，只是一个虚拟的量% SubObjFun1.m functionf=SubObjFun1(k,*1,*2,u1,u2) c=4 5 6;f=-c(k)*u1; % 求最大值转化为求最小值% TransFun1.m functions=TransFun1(k,*1,*2,u1,u2) w=3 4 5;v=8 6 4;s(1)=*1-u1*w(k);s(2)=*2-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)