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

动态规划方法的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动态规划的基本模型实际中,要构造一个标准的动态规划模型,通常需要采用以下几个步骤:①划分阶段按照问题的时间或空间特征,把问题分为若干个阶段。这些阶段必须是有序的或者是可排序的(即无后向性),否则,应用无效。②选择状态将问题发展到各个阶段时所处的各种客观情况用不同的状态表示,即称为状态。状态的选择要满足无后效性和可知性,即状态不仅依赖于状态的转移规律,还依赖于允许决策集合和指标函数结构。③确定决策变量与状态转移方程当过程处于某一阶段的某个状态时,可以做出不同的决策,描述决策的变量称为决策变量。在决策过程中,由一个状态到另一个状态的演变过程称为状态转移。状态转移就是根据上一阶段的状态和决策来导出本阶段的状态。④写出动态规划的基本方程动态规划的基本方程一般根据实际问题可分为两种形式,逆序形式和顺序形式。这里只考虑逆序形式。动态规划基本方程的逆序形式为11()(){((,)())}kkkkkkkkkkxDsfsoptgvsxfs++∈=+,1,,1knn=−⋯边界条件11()0nnfs++=或()(,)nnnnnfsvsx=其中第k阶段的状态为ks,其决策变量kx表示状态处于ks的决策,状态转移方程为1(,)kkkksTsx+=,k阶段的允许决策集合记为()kkDs,(,)kkkvsx为2指标函数。当求解时,由边界条件从kn=开始,由后向前逆推,逐阶段求出最优决策和过程的最优值,直到最后求出11()fs即得到问题的最优解。动态规划逆序解法计算程序框图如下：2基本方程逆序算法的matlab程序（1）动态规划逆序求最小值的基本方程：11()()min{((,)())}kkkkkkkkkkxDsfsgvsxfs++∈=+,1,,1knn=−⋯边界条件()(,)nnnnnfsvsx=1(,)kkkksTsx+=。自由始端和终端的动态规划，求指标函数最小值的逆序算法递归计算程序：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)ifnargin5tmp00=tmp00+f_opt(tmp60(1),ii+1);elsetmp00=feval(ObjFun,tmp00,f_opt(tmp60(1),ii+1));end3iftmp00=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,[1234])=[ii,tmpx(i),tmpd(i),tmpf(i)];endend（2)2)2)2)当状态变量是二维时，也即有两个状态变量，此时动态规划逆序求最小值的基本方程：11(),()()min{((,,,)())}kkkkktkkkkkkkkkxDsuUtfsgvsxtufs++∈∈=+,1,,1knn=−⋯边界条件11()(,,,),(,),(,)nnnnnnnkkkkkkkkfsvsxtusTsxtPtu++===此时上面的程序就无能为力了，为此在程序dynprog.m基础上进行拓展，我们得到状态变量为二维情况下的指标函数最小值的逆序算法递归计算程序:dynprog1dynprog1dynprog1dynprog1.m.m.m.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);t