基于Matlab的动态规划程序实现

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

u=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_vub

f_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vubm(i,k)=tmp;

end;end;end

for ii=k-1:-1:1

tmp10=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 value

else,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:tmp01

tmpd(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:k

tmpx(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;