动态规划

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

f_opt=nan*ones(size(x));

% f_opt为不同阶段、状态下的最优值矩阵，初值为非数

d_opt=f_opt; % d_opt为不同阶段不同状态下的决策矩阵，初值为非数 tmp1=find(x_isnan(:,1)); % 找出第1阶段状态值（不是非数）的下标 tmp2=length(tmp1); for i=1:tmp2

u=feval(DecisFun,1,x(tmp1(i),1)); % 求出相应的允许决策向量 tmp3=length(u); t_vubm=inf;

for j=1:tmp3 % 该for语句是为了求出相应的最有函数值以及最优决策 tmp=feval(SubObjFun,1,x(tmp1(i),1),u(j)); if tmp<=t_vubm

f_opt(tmp1(i),1)=tmp; d_opt(tmp1(i),1)=u(j); t_vubm=tmp; end end end

for ii=2:k

% 从前往后面递推求出f_opt以及d_opt tmp10=find(x_isnan(:,ii)); tmp20=length(tmp10); for i=1:tmp20

u=feval(DecisFun,ii,x(tmp10(i),ii)); tmp30=length(u); t_vubm=inf; 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);

1

% 找出上一阶段的状态值在x(:,ii-1)的下标 if ~isempty(tmp60) if nargin<5

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

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

if tmp00<=t_vubm

f_opt(tmp10(i),ii)=tmp00;d_opt(i,ii)=u(j); t_vubm=tmp00; end end end end end

fval=f_opt(find(x_isnan(:,k)),k); % fval即为最有函数值矩阵 p_opt=[];tmpx=[];tmpd=[];tmpf=[];

tmp0=find(x_isnan(:,k));tmp01=length(tmp0); for i=1:tmp01

tmpd(i)=d_opt(tmp0(i),k); % 求出第k阶段的决策值 tmpx(i)=x(tmp0(i),k); % 求出第k阶段的状态值

tmpf(i)=feval(SubObjFun,k,tmpx(i),tmpd(i)); % 求出第k阶段的指标函数值 p_opt(k*(i-1)+k,[1 2 3 4])=[k,tmpx(i),tmpd(i),tmpf(i)];

for ii=k-1:-1:1 % 按逆序求出各阶段的决策值、状态值以及指标函数值

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

function u=DecisFun(k,x) d=[2 3 2 4];m=6; if k>1

u=0:min(x+d(k),m); else

u=x+d(1); end

2

function s=TransFun(k,x,u) d=[2 3 2 4]; s=x-u+d(k);

function f=SubObjFun(k,x,u) d=[2 3 2 4]; if u==0

f=0.5*x; else if u>6

f=10^6; else

f=3+u+0.5*x; end end

x1=0:4;s=nan*ones(5,1);s(1)=0; x=[x1' x1' x1' s];

[p_opt,fval]=dynprog(x,'DecisFun','SubObjFun','TransFun')

3

二维情形

function

[p_opt,fval]=mindynprog1(x1,x2,DecisFun,StageObjFun,Stage_TransFun,ObjFun)

% (x1,x2)为二维状态变量,这里矩阵x1与x2的列数应相同,该程序考虑的决策变量也是二维

% M_函数[u1,u2]=DecisFun(k,x1,x2)表示由阶段k的状态值(x1,x2)求出相应的允许决策集合

% M_函数v=StageObjFun(k,x1,x2,u1,u2)表示阶段k的指标函数 % M_函数[s1,s2]=Stage_TransFun(k,x1,x2,u1,u2)是状态转移函数,其中(x1,x2)是阶段k的状态

%值,(u1,u2)是相应的决策值,(s1,s2)为下一阶段即阶段k+1的状态值 % M_函数ObjFun（v,f)是第k阶段到最后阶段的指标函数,当ObjFun（v,f)=v+f时,输入ObjFun（v,f)可以省略

% 输出p_opt由4列组成，p_opt=[序号组，最优轨线组，最优策略组，指标函数值组];

% 输出fval是一个列向量，各元素分别表示p_opt各最优策略组对应始端状态x的最优函数值

[k1,k]=size(x1);[k2,k]=size(x2); % k为阶段数 x1_isnan=~isnan(x1); x2_isnan=~isnan(x2);

4

f_opt=nan*ones(k1,k2,k); % f_opt(i,j,m)为第m阶段初状态值为(x1(i,m),x2(j,m))下的最优值,初值为非数

d_opt1=f_opt;d_opt2=f_opt; % (d_opt1(i,j,m),d_opt2(i,j,m))为第m阶段初状态值为(x1(i,m),x2(j,m))下的最优决策值,初值为非数

tmp11=find(x1_isnan(:,1));

tmp21=find(x2_isnan(:,1)); % 找出第k阶段状态值(不是非数)的下标 tmp12=length(tmp11); tmp22=length(tmp21); for i=1:tmp12 for t=1:tmp22

[u1,u2]=feval(DecisFun,1,x1(tmp11(i),1),x2(tmp21(t),1)); % 求出相应的允许决策向量, 若决策变量为一维，那么在定义DecisFun函数时，就令u2恒为1 tmp13=length(u1);tmp14=length(u2); t_vubm=inf;

% 下面两个for语句是为了求出第k阶段初状态值为x1(tmp11(i),k),x2(tmp21(t),k)时的最优函数值以及最优决策值 for j=1:tmp13 for l=1:tmp14

5

本文来源：https://www.bwwdw.com/article/h4od.html