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

function [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)
8 B i' J$ W+ V/ ~1 {! Z
% [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun). A& M' q2 m$ I2 ]/ V! a: L; e
% 自由始端和终端的动态规划,求指标函数最小值的逆序算法递归' N; J+ d6 d, R0 {7 x# K$ c
% 计算程序。x是状态变量，一列代表一个阶段状态；M-函数$ J9 t% A0 m5 \- W1 F- e' y3 ^
% DecisFun(k,x)由阶段k的状态变量x求出相应的允许决策变量;8 U" }; K( t+ r) Z, R5 e- X+ n
% M-函数ObjFun(k,x,u)是阶段指标函数，M-函数TransFun(k,x,u)% ^0 B- `# D( O
% 是状态转移函数,其中x是阶段k的某状态变量，u是相应的决策变量；
( g* ~# e2 f* K
% 输出p_opt由4列构成，p_opt=[序号组;最优策略组;最优轨线组;
% v( ~! Y9 ^6 N; ~( b
% 指标函数值组]；fval是一个列向量，各元素分别表示p_opt各4 D6 ~$ E2 u- T3 \# y# m
% 最优策略组对应始端状态x的最优函数值；, t3 Q% Z, u [( c" J' p
%
2 A* w+ k2 Q& j& V
%例（参看胡良剑等编《数学实验--使用MATLAB》P180
) I3 f. V9 ]' |* i- c$ X
%先写3个函数
" g! l( k9 o4 e. a1 r5 X7 u; p
% eg13f1_2.m
* W) z) f7 m3 U- e
% function u=DecisF_1(k,x)
8 n# `( G1 F( E* E3 Y8 v
% 在阶段k由状态变量x的值求出其相应的决策变量所有的取值
( U6 e0 x( ?3 ^8 k6 ?# Y
% c=[70,72,80,76];q=10*[6,7,12,6];
& T, s' r& f6 y
% if q(k)-x<0,u=0:100; %决策变量不能取为负值
3 A* H- @7 h. Q
% else,u=q(k)-x:100;end; %产量满足需求且不超过100* E' V6 { Q3 e8 |. z
% u=u(:);- \& ?3 U) z3 }* f3 J
% eg13f2_2.m
+ w0 U8 b2 y* v
% function v=ObjF_1(k,x,u)( Q* R5 B1 W1 E$ k
% 阶段k的指标函数
' j9 B8 {1 p% n* R7 X3 q0 T
% c=[70,72,80,76];v=c(k)*u+2*x;
: |% _, ]& ?0 I9 s* H8 D2 X( G
% eg13f3_2.m7 L& w" y. `5 ?) X% v
% function y=TransF_1(k,x,u)
$ h- [ o" T- o. `8 E/ ]
% 状态转移方程4 i, C: G0 T$ S
% q=10*[6,7,12,6];y=x+u-q(k);
0 b/ g1 _# Z2 v& i1 M
%调用DynProg.m计算如下：- E# S; m% O' j. }
% clear;x=nan*ones(14,4);% x是10的倍数，最大范围0≤x≤130,6 O' R K% N) P" H8 _
% %因此x=0,1,...13，所以x初始化取14行，nan表示无意义元素# h h0 ~; u1 f4 p& W4 [
% x(1:7,1)=10*(0:6)'; % 按月定义x的可能取值! ?' W' T+ d i! w
% x(1:11,2)=10*(0:10)';x(1:12,3)=10*(2:13)';3 G) G- P3 J0 Z( T, ^7 l+ s
% x(1:7,4)=10*(0:6)';
- }, Z @4 r8 U3 g, H# Z
% [p,f]=dynprog(x,'eg13f1_2','eg13f2_2','eg13f3_2')
! r8 I# g; c4 Y4 l u: O
# p% e- m( q' o2 ^0 F( t8 c$ I) S
% By X.D. Ding June 2000
7 E" d" w Q( g" T3 c. t* L3 Z. d
- ?" z4 D C. O6 o& }
k=length(x(1,:));f_opt=nan*ones(size(x));d_opt=f_opt;
# c+ X3 g! Q; S) j0 R5 Y! b% K
t_vubm=inf*ones(size(x));x_isnan=~isnan(x);t_vub=inf;4 `$ @* z8 p$ ~, A
% 计算终端相关值" ?' m# K' I* D4 ]9 j! L, r- C/ @( j
tmp1=find(x_isnan(:,k));tmp2=length(tmp1);
, c* Y6 ]+ B. R& U9 o
for i=1:tmp2. n a' X. \" w: k& y
u=feval(DecisFun,k,x(i,k));tmp3=length(u);/ z, ^) p# I; ]! T. |- R
for j=1:tmp3" W0 \, }, M [& Z
tmp=feval(ObjFun,k,x(tmp1(i),k),u(j));
! G* J0 P f( ~
if tmp<=t_vub,
. H+ E1 S" Y" O
f_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vub=tmp; $ L$ U9 J$ Z j" E
end;end;end$ d# e' a" P0 M ?
% 逆推计算各阶段的递归调用程序
5 t; s1 S0 _5 u% W
for ii=k-1:-1:17 {5 R: f% i% K4 D
tmp10=find(x_isnan(:,ii));tmp20=length(tmp10);
5 ]" v& G ?2 J
for i=1:tmp20
; ] r2 _0 g* e' |, B5 M. k
u=feval(DecisFun,ii,x(i,ii));tmp30=length(u);. ?( _0 C9 P4 v9 I
for j=1:tmp30
1 [4 B. c8 t$ d1 J8 S B5 }, t( {7 [
tmp00=feval(ObjFun,ii,x(tmp10(i),ii),u(j));
) H/ m4 l, ]! q9 A7 C
tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j));
8 [' v7 _! W- H) ^' }
tmp50=x(:,ii+1)-tmp40;* C. V, R @0 [
tmp60=find(tmp50==0);
6 ?$ d/ ]$ R2 N" I
if ~isempty(tmp60),1 c. n& l& `' v4 {# w
tmp00=tmp00+f_opt(tmp60(1),ii+1);
i! Q, I1 N4 {
if tmp00<=t_vubm(i,ii)
& q0 d- L; W7 Z: {5 _5 ?6 c( f
f_opt(i,ii)=tmp00;d_opt(i,ii)=u(j);$ r. f$ L/ Y& ~* L* p
t_vubm(i,ii)=tmp00;
1 q6 j+ B! U" y/ c2 H3 |
end;end;end;end;end;
& Z* m- m+ R5 h u, w
fval=f_opt(tmp1,1);
: k9 F# {9 v2 @0 D, S4 r% i
% 记录最优决策、最优轨线和相应指标函数值8 j% v' s. j6 o- j0 W
p_opt=[];tmpx=[];tmpd=[];tmpf=[];# d8 x! R& R9 \3 {
tmp0=find(x_isnan(:,1));tmp01=length(tmp0);
0 m. M0 J5 |' F
for i=1:tmp01,. y) m9 q/ |: S; f
tmpd(i)=d_opt(tmp0(i),1); : [0 [! E; F* ?2 F V- o
tmpx(i)=x(tmp0(i),1);- k; J" ^# Q/ w: B( P: O. B4 ^
tmpf(i)=feval(ObjFun,1,tmpx(i),tmpd(i));
9 ~; C/ d# |7 ]7 J* l
p_opt(k*(i-1)+1,[1,2,3,4])=[1,tmpx(i),..." c f3 ^& o) u. r, P: T
tmpd(i),tmpf(i)];7 T3 u/ Z# _% }4 U; u/ `1 n) E
for ii=2:k
/ l0 q4 Y6 N2 v5 n$ u
tmpx(i)=feval(TransFun,ii-1,tmpx(i),tmpd(i));) X: M5 q" a8 R% t F) E: J. |* D
tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1==0);+ r' Q! Y! D* E1 X$ R5 K
if ~isempty(tmp2)
8 y4 ?& L- ^7 ~
tmpd(i)=d_opt(tmp2(1),ii);0 |1 A6 y. ~1 o3 W, ^
end;
. o7 p; I$ F# E, y3 C7 f9 [
tmpf(i)=feval(ObjFun,ii,tmpx(i),tmpd(i));
& s* s' l9 M# k2 B5 N: W
p_opt(k*(i-1)+ii,[1,2,3,4])=[ii,tmpx(i),...
i) `1 L8 t( t# `0 b
tmpd(i),tmpf(i)];4 p0 b* O. t0 }9 J/ G8 a! s
end;end;# ` M9 q& ~' @0 _, |

复制代码