动态规划方法的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

动态规划是一类解决多阶段决策问题的数学方法, 在工程技术、科学管理、工农业生产及军事等领域都有广泛的应用。在理论上,动态规划是求解这类问题全局最优解的一种有效方法,特别是对于实际中某些非线性规划问题可能是最优解的唯一方法。然而,动态规划仅仅决多阶段决策问题的一种方法,或者说是考查问题的一种途径,而不是一种具体的算法。就目前而言,动态规划没有统一的标准模型,其解法也没有标准算法,在实际应用中,需要具体问题具体分析。动态规划模型的求解问题是影响动态规划理论和方法应用的关键所在,而子问题的求解和大量结果的存储、调用更是一个难点所在。然而, 随着计算机技术的快速发展,特别是内存容量和计算速度的增加,使求解较小规模的动态规划问题成为可能,从而使得动态规划的理论和方法在实际中的应用范围迅速增加。

目前,在计算机上实现动态规划的一般求解方法并不多见,尤其是用来解决较复杂的具体问题的成果甚少。本文从实际出发,利用数学工具软件matlab 的强大功能, 对动态规划模型的求解方法做了尝试,编写出了动态规划逆序算法的matlab程序，并结合“生产与存储问题”和“背包问题”进行了应用与检验,实际证明结果是令人满意的。

[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

边界条件

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

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))

3

; if tmp<=t_vubm(i,k)

tmpx(i)=feval(TransFun,ii-1,tmpx(i),tmpd f_opt(tmp1(i),k)=tmp;

(i));

d_opt(tmp1(i),k)=u(j); t_vubm(i,k)=tmp; end end tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1==0);

end for ii=k-if ~isempty(tmp2)

1:-1:1

tmpd(i)=d_opt(tmp2(1),ii); % 从后往前面递推求出f_opt以及d_opt

end

tmp10=find(x_isnan(:,ii));tmp20=length(t tmpf(i)=feval(SubObjFun,ii,tmpx(i),tmpd( mp10);

i)); for i=1:tmp20

p_opt(k*(i-1)+ii,[1 2 3 u=feval(DecisFun,ii,x(tmp10(i),ii)); 4])=[ii,tmpx(i),tmpd(i),tmpf(i)];

tmp30=length(u); end for j=1:tmp30 end

tmp00=feval(SubObjFun,ii,x(tmp10(i),ii), （2)当状态变量是二维时，也即有两个状态变量，此时u(j));

动态规划逆序求最小值的基本方程：

tmp40=feval(TransFun,ii,x(tmp10(i),ii),u (j));

% 由该状态值及相应的决策值求出下一阶段的状态值

f sk k( ) = k k k k k t∈ (tmp50=x(:,ii+1)-tmp40;

min), ∈ ( ){ (gv s x t u f sk k k k k(

, , x D s u U ttmp60=find(tmp50==0);

% 找出下一阶段的状态值在x(:,ii+1)的下标 k nn= , ?1,?,1

if ~isempty(tmp60)

if nargin<5 边界条件

tmp00=tmp00+f_opt(tmp60(1),ii+1);

f s v s x t un n() = n n n n n( ,

,

else

, ), tmp00=feval(ObjFun,tmp00,f_opt(tmp60(1), ii+1)); end

if tmp00<=t_vubm(i,ii) sk+1 =T s xk k k( , ), f_opt(tmp10(i),ii)=tmp00;d_opt(i,ii)=u(j)

tk+1 =P t uk k k( , )

;

此时上面的程序就无能为力了，为此在程序

t_vubm(tmp10(i),ii)=tmp00;

dynprog.m基础上进行拓展，我们得到状态变量为二end 维情况下的指标函数最小值的逆序算法递归计算程end 序:dynprog1.m，如下： end function

end [p_opt,fval]=dynprog1(x1,x2,DecisFun,Sub end

ObjFun,TransFun,ObjFun)

fval=f_opt(find(x_isnan(:,1)),1); % (x1,x2)为二维状态变量，其中x1,x2的取值是相% fval即为最有函数值矩阵

互独立的，这里矩阵x1与x2的列数应相同，该程序p_opt=[];tmpx=[];tmpd=[];tmpf=[];

考虑决策变量也是二维

tmp0=find(x_isnan(:,1));tmp01=length(tmp TcisFun(k,x1,x2),SubObjFun(k,x1,x2,u1, 0); for i=1:tmp01

u2),TransFun(k,x1,x2,u1,u2)等M_函数的含义 tmpd(i)=d_opt(tmp0(i),1);

与一维的情形一样，只是它们的参数相应的增加了， % 求出第一阶段的决策值

ObjFun函数的含义及参数保持不变

tmpx(i)=x(tmp0(i),1);

% [p_opt,fval]的含义与一位情形一样，只是它们的% 求出第一阶段的状态值

维数增加了

tmpf(i)=feval(SubObjFun,1,tmpx(i),tm % 下面程序的思路与算法同一维基本相同，只是相应矩pd(i));

阵的维数增加了

% 求出第一阶段的指标函数值 p_opt(k*(i-1)+1,[12 3 4])=

[k1,k]=size(x1);[k2,k]=size(x2);

[1,tmpx(i),tmpd(i),tmpf(i)];

x1_isnan=~isnan(x1);x2_isnan=~isnan(x2);

for ii=2:k

% 按顺序求出各

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);

4

, )+ k+

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

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);

% 若决策变量为一维，那么在定义DecisFun函数时，就

if ~isempty(tmp3)&~isempty(tmp4) (tmp011(i),ii),x2(tmp021(t),ii)); 令

tmpd1(i)=d_opt1(tmp3(1),tmp4(1),ii); u2恒为1

tmpd2(j)=d_opt2(tmp3(1),tmp4(1),ii); end tmp013=length(u1);

tmpf(i,j)=feval(SubObjFun,ii, tmp014=length(u2);

tmpx1(i),tmpx2(j),mpd1(i),tmpd2(j)); for j=1:tmp013

p_opt(t+ii,[1 2 3 4 5 6]) for l=1:tmp014

=[ii,tmpx1(i),tmpx2(j),

tmpd1(i),tmpd2(j),tmpf(i,j)]; tmp000=feval(SubObjFun,ii,x1(tmp011(i),i

end i),x2(tmp021(t),ii),u1(j),u2(l));

end tmp100=feval(TransFun,ii,x1(tmp011(i),ii)

p=k*tmp02;q=i*p; ,x2(tmp021(t),ii),u1(j),u2(l));

end tmp200=x1(:,ii+1)-tmp100(1);

tmp300=x2(:,ii+1)-tmp100(2); tmp400=find(tmp200==0);

3 实例应用 tmp500=find(tmp300==0);

if~isempty(tmp400)&~isempty(tmp500) if nargin<6 3.1 生产与存储问题 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

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);

某工厂要对一种产品制定今后四个时期的生产计划，据估计在今后四个时期内，市场对于该产品的需求量如下表所示： 时期（k) 需求量（d) k1 2 2 3 3 2 3 4 5