文档介绍：[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