动态规划方法的Matlab 实现与应用动态规划(Dynamic Programming)是求解决策过程最优化的有效数学方法，它是根据“最优决策的任何截断仍是最优的”这最优性原理，通过将多阶段决策过程转化为一系列单段决策问题，然后从最后一段状态开始逆向递推到初始状态为止的一套最优化求解方法。
1．动态规划基本组成(1) 阶段 整个问题的解决可分为若干个阶段依次进行，描述阶段的变量称为阶段变量，记为k(2) 状态 状态表示每个阶段开始所处的自然状况或客观条件，它描述了研究问题过程的状况。
各阶段状态通常用状态变量描述，用k x 表示第k 阶段状态变量，n 个阶段决策过程有n+ 1个状态。
(3) 决策 从一确定的状态作出各种选择从而演变到下一阶段某一状态，这种选择手段称为决策。
描述决策的变量称为决策变量，决策变量限制的取值范围称为允许决策集合。
用()k k u x 表示第k 阶段处于状态k x 时的决策变量，它是k x 的函数。
用()k k D x Dk(xk)表示k x 的允许决策的集合。
(4) 策略 每个阶段的决策按顺序组成的集合称为策略。
由第k 阶段的状态k x 开始到终止状态的后部子过程的策略记为{}11(),(),,()k k k k n n u x u x u x ++ 。
可供选择的策略的范围称为允许策略集合，允许策略集合中达到最优效果的策略称为最优策略。
从初始状态*11()x x =出发，过程按照最优策略和状态转移方程演变所经历的状态序列{}****121,,,,n n x x x x + 称为最优轨线。
(5) 状态转移方程 如果第k 个阶段状态变量为k x ，作出的决策为k u ，那么第k+ 1阶段的状态变量1k x +也被完全确定。
用状态转移方程表示这种演变规律，记为1(,)k k k x T x u +=。
(6) 指标函数 指标函数是系统执行某一策略所产生结果的数量表示，是衡量策略优劣的数量指标，它定义在全过程和所有后部子过程上，用()k k f x 表示。
过程在某阶段j 的阶段指标函数是衡量该阶段决策优劣数量指标，取决于状态j x 和决策j u ，用(,)j j j v x u 表示。
2．动态规划基本方程(){}11()min ,,(),()k k k k k k k k k k f x g v x u f x u D x ++=∈⎡⎤⎣⎦Matlab 实现 (dynprog.m 文件)function [p_opt,fval]=dynprog (x,DecisFun,SubObjFun,TransFun,ObjFun)% x 是状态变量，一列代表一个阶段的所有状态；% M-函数DecisFun(k,x) 由阶段k 的状态变量x 求出相应的允许决策变量; % M-函数SubObjFun(k,x,u) 是阶段指标函数,% M-函数ObjFun(v,f) 是第k 阶段至最后阶段的总指标函数% M-函数TransFun(k,x,u) 是状态转移函数, 其中x 是阶段k 的某状态变量, u 是相应的决策变量; %输出 p_opt 由4列构成，p_opt=[序号组;最优策略组;最优轨线组;指标函数值组]；%输出 fval 是一个列向量，各元素分别表示p_opt 各最优策略组对应始端状态x 的最优函数值。
k=length(x(1,:)); % 判断决策级数x_isnan=~isnan(x); % 非空状态矩阵t_vubm=inf*ones(size(x)); % 性能指标中间矩阵f_opt=nan*ones(size(x)); % 总性能指标矩阵d_opt=f_opt; % 每步决策矩阵tmp1=find(x_isnan(:,k)); % 最后一步状态向量tmp2=length(tmp1); % 最后一步状态个数for i=1:tmp2u=feval(DecisFun,k,x(tmp1(i),k)); tmp3=length(u); % 决策变量for j=1:tmp3 % 求出当前状态下所有决策的最小性能指标tmp=feval(SubObjFun,k,x(tmp1(i),k),u(j));if tmp <= t_vubm(i,k) %t_vubf_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vubm(i,k)=tmp;end;end;endfor ii=k-1:-1:1tmp10=find(x_isnan(:,ii));tmp20=length(tmp10);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; % 找出下一状态在 x 矩阵的位置tmp60=find(tmp50==0);if~isempty(tmp60),if nargin<5,tmp00=tmp00+f_opt(tmp60(1),ii+1); % set the default object valueelse,tmp00=feval(ObjFun,tmp00,f_opt(tmp60(1),ii+1)); end %当前状态的性能指标 if tmp00<=t_vubm(i,ii) f_opt(i,ii)=tmp00;d_opt(i,ii)=u(j);t_vubm(i,ii)=tmp00; end; end;end; end;end;fval=f_opt(:,1);tmp0 = find(~isnan(fval));fval=fval(tmp0,1);p_opt=[];tmpx=[];tmpd=[];tmpf=[];tmp01=length(tmp0);for i=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,[1,2,3,4])=[1,tmpx(i),tmpd(i),tmpf(i)];for ii=2:ktmpx(i)=feval(TransFun,ii,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;3．动态规划方法应用例 某电子设备由5种元件1,2,3,4,5组成，其可靠性分别为0.9,0.8,0.5,0.7,0.6 为保证电子设备系统的可靠性，同种元件可并联多个。
现允许设备使用元件的总数为15个，问如何设计使设备可靠性最大。
将该问题看成一个5阶段动态规划问题，每个元件的配置看成一个阶段。
记k x 为配置第k 个元件时可用元件的总数(状态变量)；k u 为第k 个元件并联的数目(决策变量);k c 为第k 个元件的可靠性，阶段指标函数为 (,)1(1)k uk k k k v x u c =−−; 状态转移方程为1k k k x x u +=−基本方程为(){}{}111()min ,,()min k k k k k k k k k f x g v x u f x v f +++==⎡⎤⎣⎦i根据上述的允许决策、阶段指标函数、状态转移方程和基本方程写出下面的4个M-函数%DecisF1.mfunction u=DecisF1(k,x) %在阶段k 由状态变量x 的值求出其相应的决策变量所有的取值 if k==5,u=x; else u=1:x-1;end%SubObjF1.mfunction v=SubObjF1(k,x,u) %阶段k 的指标函数 c=[0.9,0.8,0.5,0.5,0.4]; v=1-(1-c(k)).^ u; v= -v; %将求max 转换为求min%TransF1.mfunction y=TransF1(k,x,u) %状态转移方程 y=x-u;%ObjF1.mfunction y=ObjF1(v,f) %基本方程中的函数g y=v*f; y=-y; %将求max 转换为求min调用DynProg.m 计算的主程序如下(example1.m): clear;n=15; %15个元件 x1=[n;nan*ones(n-1,1)];x2=1:n; x2=x2’; x=[x1,x2,x2,x2,x2]；[p,f]=dynprog(x,’DecisF1’,’SubObjF1’,’TransF1’,’ObjF1’)运行结果： p = 1.0000 15.0000 2.0000 -0.9900 2.0000 13.0000 2.0000 -0.9600 3.0000 11.0000 4.0000 -0.9375 4.0000 7.0000 3.0000 -0.9730 5.0000 4.0000 4.0000 -0.9744 f = -0.8447结果表明1,2,3,4和5号元件分别并联2,2,4,3,4个，系统总可靠性最大为0.84474．实验要求(1) 根据前面2，3部分的说明，编写出dynprog.m ，DecisF1.m ，SubObjF1.m ，TransF1.m ，ObjF1.m 五个文件，并运行example1.m 文件，查看结果是否正确。
(2) 现有4种不同的车床1,2,3,4, 同时加工500件相同的零件。
各车床加工一个零件的时间分别为0.5,0.1,0.2和0.05 . 问如何给4个车床分配加工零件数目,使完工时间最短? 要求参考(1) 编写出DecisF2.m ，SubObjF2.m ，TransF2.m ，ObjF2.m 四个文件，通过调用dynprog 函数，求解本问题。
[提示：可将问题按车床编号分为4个阶段。
设状态变量k x 表示分配给第k 号车床至第4号车床的零件总数。
决策变量k u 表示分配给第k 号车床的零件数, 1,2,,k k u x = 且有44u x =.状态转移方程1k k k x x u +=−, 阶段指标函数()k k v u 表示k u 个零件分配到第k 号车床加工所需的时间()k k k k v u u t =, k t 是k 号车床的加工时间. ()k k f x 表示k x 个零件分配给第k 至第4号车床加工所需的最短时间 (用时最长的车床所需时间为总加工时间)](3) 已知离散系统方程 x(k+1)= x(k)+u(k), x(1)=1, }1,1{)(−∈k u ,4321()()()()()k J x k u k x k u k x k =⎡⎤=++⎣⎦∑, 求控制序列u*(k)和状态x*(k), 使J=min.要求参考(1) 编写出DecisF3.m ，SubObjF3.m ，TransF3.m ，ObjF3.m 四个文件，通过调用dynprog 函数，求解本问题。