求动态规划的MATLAB程序，要通用性强的 - 数学建模社区-数学中国

function [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun) 5 b; @" N) U4 D, |9 A4 U P6 Y
% [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)( ]: f5 d) y z. M2 I$ q( ^# r
% 自由始端和终端的动态规划,求指标函数最小值的逆序算法递归7 r, }4 d. m% x9 H6 R
% 计算程序。x是状态变量，一列代表一个阶段状态；M-函数6 V# z" ]4 s$ @' I W) ~. Z! {
% DecisFun(k,x)由阶段k的状态变量x求出相应的允许决策变量;
9 C* j5 ^) N; e9 O3 S8 _/ M
% M-函数ObjFun(k,x,u)是阶段指标函数，M-函数TransFun(k,x,u)6 F" \/ B) b' Y
% 是状态转移函数,其中x是阶段k的某状态变量，u是相应的决策变量；/ c6 t9 |# `1 M& h4 v+ ]6 X |
% 输出p_opt由4列构成，p_opt=[序号组;最优策略组;最优轨线组;8 l8 m$ q0 G. \: x
% 指标函数值组]；fval是一个列向量，各元素分别表示p_opt各
/ \5 B! ^5 S# m
% 最优策略组对应始端状态x的最优函数值；; z: q" [; b9 K; Y- t
%1 y9 v: k M! ]9 G
%例（参看胡良剑等编《数学实验--使用MATLAB》P180
8 f1 y9 d) V' q' B
%先写3个函数
4 l7 |0 u$ ?/ X% j* f# h
% eg13f1_2.m4 P- a% F& V Y# S6 D$ b
% function u=DecisF_1(k,x)
9 c2 w8 P0 G7 S8 l+ Z2 X
% 在阶段k由状态变量x的值求出其相应的决策变量所有的取值
. ~/ o+ U6 Q+ u/ C) w7 |
% c=[70,72,80,76];q=10*[6,7,12,6];7 v7 k7 I/ Z0 e; {
% if q(k)-x<0,u=0:100; %决策变量不能取为负值
3 R& l% C. Z' E: M+ _! ]7 z( E
% else,u=q(k)-x:100;end; %产量满足需求且不超过100; | n; M/ c4 e0 m/ e8 S# E' m
% u=u(:);
7 ]- ]4 T- F. h. }- [. r$ J
% eg13f2_2.m; M* B7 @9 }) i( G
% function v=ObjF_1(k,x,u)& c1 H! a1 Q+ F
% 阶段k的指标函数
8 K, }0 [: s! Y- ^7 @$ ]* J: \
% c=[70,72,80,76];v=c(k)*u+2*x;
+ Q/ Z" A" y3 W- @, O! U" @& d
% eg13f3_2.m% N" H T. D! F( d
% function y=TransF_1(k,x,u)
6 ^: @1 d* }4 @0 U7 w. ~* R$ ~
% 状态转移方程9 y# S6 ^2 U1 \' d, L+ Q' s
% q=10*[6,7,12,6];y=x+u-q(k);( }0 o8 E+ L; s7 W$ E* M4 j
%调用DynProg.m计算如下：
5 ]# {" d0 p1 @) h9 W& [2 P
% clear;x=nan*ones(14,4);% x是10的倍数，最大范围0≤x≤130,9 Q8 [) G {. e; O
% %因此x=0,1,...13，所以x初始化取14行，nan表示无意义元素
( k% A7 V; X4 M" W. C, L; t
% x(1:7,1)=10*(0:6)'; % 按月定义x的可能取值
$ I0 C6 w8 R& o% e9 E n3 T
% x(1:11,2)=10*(0:10)';x(1:12,3)=10*(2:13)';' H( D4 i$ k6 y& F2 l- O, R9 O
% x(1:7,4)=10*(0:6)';
m: n) E6 z; Y# i! Y, Q; L
% [p,f]=dynprog(x,'eg13f1_2','eg13f2_2','eg13f3_2'); K' {- j9 I& a
$ h1 g6 V3 {7 x' J3 q! R0 e+ ~9 k
% By X.D. Ding June 2000
1 u# z; W, f* K, R0 h# }, T
" O3 U U) t3 }. r) \5 D
k=length(x(1,:));f_opt=nan*ones(size(x));d_opt=f_opt; j) a( h F7 j3 n# i
t_vubm=inf*ones(size(x));x_isnan=~isnan(x);t_vub=inf;8 N4 g; n/ Z+ R+ K e
% 计算终端相关值
4 s' x m& B# k' M
tmp1=find(x_isnan(:,k));tmp2=length(tmp1);9 P. [* @: g( L8 d. g' Q, P5 U
for i=1:tmp2 P1 j5 @! I2 C2 `! h+ W. K! n: j
u=feval(DecisFun,k,x(i,k));tmp3=length(u);. f7 T) y3 v0 `+ z1 V4 {6 l
for j=1:tmp3 b Q3 t- T, T1 D! U9 ^4 M _2 P
tmp=feval(ObjFun,k,x(tmp1(i),k),u(j));
V M' a. h) k& y" h
if tmp<=t_vub, S9 m0 f& }7 V# j
f_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vub=tmp;
) n6 R: t" ]8 c4 k- q
end;end;end- r. ^4 H+ V+ A/ r4 o0 L
% 逆推计算各阶段的递归调用程序% o2 B) g2 j, o
for ii=k-1:-1:1
" s- ^8 N# j. P( {
tmp10=find(x_isnan(:,ii));tmp20=length(tmp10);
7 y; x/ v+ @% e" N/ Z
for i=1:tmp20
) Q) x* H! N- f* `# n, v7 O
u=feval(DecisFun,ii,x(i,ii));tmp30=length(u);; d: c: T Y! M1 B# c- ]
for j=1:tmp30
. s: C+ A1 r# n- _4 L
tmp00=feval(ObjFun,ii,x(tmp10(i),ii),u(j));
1 |( w+ k; e! C* C! G
tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j));
9 O; O* S( w* j: K4 v* i2 D
tmp50=x(:,ii+1)-tmp40;
, U; m$ S" n8 m; b3 Y3 t$ g, p
tmp60=find(tmp50==0);+ Q& e, |1 V$ |; t
if ~isempty(tmp60),
. h5 x' X# l/ W4 h5 o2 l8 L
tmp00=tmp00+f_opt(tmp60(1),ii+1);
2 c5 L4 l* e+ B
if tmp00<=t_vubm(i,ii)6 Q! F: f! T5 N3 P0 r7 C0 W; ]
f_opt(i,ii)=tmp00;d_opt(i,ii)=u(j);2 n$ Q$ P; u [: h% p8 v
t_vubm(i,ii)=tmp00;
: T' D( Q. J E6 j. C4 Y
end;end;end;end;end;
! g7 g: b4 V$ Q2 ^3 u4 n- G, f
fval=f_opt(tmp1,1);2 t" q, u" q R) a& P+ {
% 记录最优决策、最优轨线和相应指标函数值5 E. V8 m- i0 S
p_opt=[];tmpx=[];tmpd=[];tmpf=[];5 j) b) H4 o/ Q/ V+ j9 ^
tmp0=find(x_isnan(:,1));tmp01=length(tmp0);4 C8 ?* u" S* x: C5 B- X
for i=1:tmp01,2 Q" K# u+ C* \) I
tmpd(i)=d_opt(tmp0(i),1);
' \% N' k9 n7 _. ~, V6 y, I# M% F1 F
tmpx(i)=x(tmp0(i),1);
+ U* j0 [4 U# ^3 O+ M0 b
tmpf(i)=feval(ObjFun,1,tmpx(i),tmpd(i));
: J6 U1 e0 e& _
p_opt(k*(i-1)+1,[1,2,3,4])=[1,tmpx(i),...
; F+ B" S1 c. [3 |* q8 f
tmpd(i),tmpf(i)];& S: F+ q; O) `; D5 d: [7 v A/ e8 o
for ii=2:k& w/ G0 ~ d' h
tmpx(i)=feval(TransFun,ii-1,tmpx(i),tmpd(i));( c: o7 D* d: ~6 M
tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1==0);5 R8 L* l: }5 G3 k& ^) K
if ~isempty(tmp2)& m4 s% p6 |/ [% ?
tmpd(i)=d_opt(tmp2(1),ii);
" ]2 @: x2 W# Q; _( m# y- |
end;
& @- W8 e* `7 ~7 D5 E5 ]
tmpf(i)=feval(ObjFun,ii,tmpx(i),tmpd(i));5 f& |* r* b# z- F, |( l/ U
p_opt(k*(i-1)+ii,[1,2,3,4])=[ii,tmpx(i),...
! V7 T, M5 T* r3 @$ X3 \
tmpd(i),tmpf(i)];
; u; ~. A( `, N: X
end;end;3 F& C! g9 N( h* o' S/ ^

复制代码