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

动态规划方法的matlab实现及其应用
（龙京鹏，张华庆，罗明良，刘水林） (南昌航空大学，数学与信息科学学院，江西，南昌)
摘要：本文运用matlab语言实现了动态规划的逆序算法，根据状态变量的维数，编写了指标函数最小值的逆序算法递归计算程序。两个实例的应用检验了该程序的有效性，同时也表明了该算法程序对众多类典型的动态规划应用问题尤其是确定离散型的应用问题的通用性，提供了求解各种动态规划问题的有效工具。关键词：动态规划 基本方程的逆序算法 MATLAB实现
MATLAB Achieve For Dynamic Programming and Its Application
(JingpengLong，HuaqingZhang，MingliangLuo，ShuilinLiu)
（School of Mathematics and Information Science,Nanchang Hangkong
University,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 achievement

8
动态规划是一类解决多阶段决策问题的数学方法, 在工程技术、科学管理、工农业生产及军事等领域都有广泛的应用。在理论上,动态规划是求解这类问题全局最优解的一种有效方法,特别是对于实际中某些非线性规划问题可能是最优解的唯一方法。然而,动态规划仅仅决多阶段决策问题的一种方法,或者说是考查问题的一种途径,而不是一种具体的算法。就目前而言,动态规划没有统一的标准模型,其解法也没有标准算法,在实际应用中,需要具体问题具体分析。动态规划模型的求解问题是影响动态规划理论和方法应用的关键所在,而子问题的求解和大量结果的存储、调用更是一个难点所在。然而, 随着计算机技术的快速发展,特别是内存容量和计算速度的增加,使求解较小规模的动态规划问题成为可能,从而使得动态规划的理论和方法在实际中的应用范围迅速增加。
目前,在计算机上实现动态规划的一般求解方法并不多见,尤其是用来解决较复杂的具体问题的成果甚少。本文从实际出发,利用数学工具软件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( , ) 。
自由始端和终端的动态规划，求指标函数最小值的
逆序算法递归计算程序：
function
[p_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:tmp2
u=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 end
end for ii=k-1:-1:1
% 从后往前面递推求出f_opt以及d_opt tmp10=find(x_isnan(:,ii));tmp20=length(t
mp10);
for i=1:tmp20
u=feval(DecisFun,ii,x(tmp10(i),ii));
tmp30=length(u); for j=1:tmp30
tmp00=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 nargin<5
tmp00=tmp00+f_opt(tmp60(1),ii+1); else
tmp00=feval(ObjFun,tmp00,f_opt(tmp60(1), ii+1)); end
if tmp00<=t_vubm(i,ii)
f_opt(tmp10(i),ii)=tmp00;d_opt(i,ii)=u(j)
;
t_vubm(tmp10(i),ii)=tmp00; end
end
end
end end
fval=f_opt(find(x_isnan(:,1)),1);
% fval即为最有函数值矩阵
p_opt=[];tmpx=[];tmpd=[];tmpf=[]; tmp0=find(x_isnan(:,1));tmp01=length(tmp
0); for i=1:tmp01 tmpd(i)=d_opt(tmp0(i),1);
% 求出第一阶段的决策值
tmpx(i)=x(tmp0(i),1);
% 求出第一阶段的状态值
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 % 按顺序求出各
阶段的决策值、状态值以及指标函数值
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,[1 2 3
4])=[ii,tmpx(i),tmpd(i),tmpf(i)]; end
end
（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 t
k 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，如下：
function
[p_opt,fval]=dynprog1(x1,x2,DecisFun,Sub
ObjFun,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(t
mp21); for i=1:tmp12 for t=1:tmp22
[u1,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:tmp14
tmp=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; end
end
end
end
end
for 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:tmp022
[u1,u2]=feval(DecisFun,ii,x1
% 若决策变量为一维，那么在定义DecisFun函数时，就
(tmp011(i),ii),x2(tmp021(t),ii)); 令u2恒为1
tmp013=length(u1); tmp014=length(u2); for j=1:tmp013
for l=1:tmp014
tmp000=feval(SubObjFun,ii,x1(tmp011(i),i
i),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) if nargin<6
tmp000=tmp000+f_opt(tmp400(1),tmp500(1),
ii+1); else tmp000=feval(ObjFun,tmp000, f_opt(tmp400(1),mp500(1),ii+1));
end if tmp000<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; end end end end end end end
fval=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(t mp11);
tmp12=find(x2_isnan(:,1));tmp02=length(t
mp12); for i=1:tmp01
for j=1:tmp02 tmpd1(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,[1 2 3 4 5 6])
=[1,tmpx1(i),tmpx2(j), tmpd1(i),tmpd2(j),tmpf(i,j)]; for ii=2:k u=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); end tmpf(i,j)=feval(SubObjFun,ii,
tmpx1(i),tmpx2(j),mpd1(i),tmpd2(j)); p_opt(t+ii,[1 2 3 4 5 6])
=[ii,tmpx1(i),tmpx2(j), tmpd1(i),tmpd2(j),tmpf(i,j)];
end
end p=k*tmp02;q=i*p;
end
3 实例应用
3.1 生产与存储问题
某工厂要对一种产品制定今后四个时期的生产计划，据估计在今后四个时期内，市场对于该产品的需求量如下表所示：
时期（k)
1
2
3
3
需求量（dk)
2
3
2
4
假定该厂生产每批产品的固定成本为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 ⎫ c xk k( ) =⎨3+1*xk.......当 xk =1,2,...6⎪⎬ ⎪
⎪⎩∞..................当 xk >6 ⎪⎭
在第 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 = −4 4
其中xk ≥ max(0,d vk − k),这是因为一方面每阶段生产的下限为0；另一方面由于要保证供应，故第k阶段末的库存量vk+1 必须非负，即v x dk + −k k ≥ 0 ，所以 x d vk ≥ −k k。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.m
function f=SubObjFun(k,x,u) d=[2 3 2 4]; if u==0 f=0.5*(x+u-d(k));
else if u>6
f=10^6;
else f=3+u+0.5*(x+u-d(k)); end
End % 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','SubObjF
un','TransFun')

运行结果如下：
p_opt =

1.0000
0
5.0000
9.5000
2.0000
3.0000
0
0
3.0000
0
6.0000
11.0000
4.0000 fval =
20.5000
4.0000
0
0
从上面的结果可知，每个时期的最优决策为：
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 ≤ ≤uk min(
w vk k
最优值函数 f s tk k k( , ) 表示当总重量不超过sk公斤，总体积不超过tk立方米背包装入第t种物品到第3种物品产生的最大效益值。
可得基本方程：
⎧⎪f s tk k k( , ) =u D s tk∈maxk 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 function
[u1,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 function
f=SubObjFun1(k,x1,x2,u1,u2) c=[4 5 6];
f=-c(k)*u1; % 求最大值转化为求最小值
% TransFun1.m function
s=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 10 22 2 1 -8
2 4 6 1 1 -5
3 0 0 0 1 0 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.6
[2] 胡良剑 丁晓东 孙晓君，数学实验-使用MATLAB[M] 上海:上海科技出版社，2001
[3] 张志涌，精通MATLAB6.5[M] 北京:北京航空航天大学出版社，2003
[4] 刘保柱 苏彦华 张宏林， MATLAB7.0从入门到精通
（修订版） 北京：人民邮电出版社，2010
5771001803090012095 579036822859633082
5771001803090012386 576137399735760696
5771001803090013594 578077579902515512
5771001803090012387 577164982601818051
5771001803090012138 572131192158918326
5771001803090012359 579036822361076053
5771001803090012356 576135286143791742
5771001803090012355 575087869704693279

17088100343355274 101229944325833379
17088100343355275 101866732938832008
17088100343356107 101581152501500522
17088100343356108 101000180059871732
17088100343354295 101074194142687017
17088100343356184 101878660869628802
17088100343356185 101775831174086674
17088100343356109 101086014373572846
17088100343356110 101152207216014916
17088100343355237 101027041605702709
17088100343355238 101229364861425414
17088100343356169 101862204402635718
17088100343354928 101760654089788804