求动态规划的MATLAB程序，要通用性强的

danieldu

看看看看那

madio

1. function [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)      
   9 S. w' h! _8 B2 s+ N
2. % [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)
   ; M' O1 V/ D0 W
3. % 自由始端和终端的动态规划,求指标函数最小值的逆序算法递归
   4 b5 |. X  ^- E7 i7 ?
4. % 计算程序。x是状态变量，一列代表一个阶段状态；M-函数. D; R$ |/ F$ F/ f
   
5. % DecisFun(k,x)由阶段k的状态变量x求出相应的允许决策变量;3 Q$ U9 }, i$ L0 J8 g1 `8 x
   
6. % M-函数ObjFun(k,x,u)是阶段指标函数，M-函数TransFun(k,x,u); l8 C9 C( }! {' F' s
   
7. % 是状态转移函数,其中x是阶段k的某状态变量，u是相应的决策变量；
   \" G: D7 R. w3 ]+ h
8. % 输出p_opt由4列构成，p_opt=[序号组;最优策略组;最优轨线组;# F# O) V# k0 g8 `& j6 z: @
   
9. % 指标函数值组]；fval是一个列向量，各元素分别表示p_opt各+ t8 ^8 a. J  P' q$ C\" a
   
10. % 最优策略组对应始端状态x的最优函数值；
    2 o9 m3 ]% g! z4 E) E/ l
11. %
    , v6 q. E+ }2 N# K0 \
12. %例（参看胡良剑等编《数学实验--使用MATLAB》P180
      X/ F( |5 ?  v, R* [, _' v
13. %先写3个函数; h, u) H/ I/ t, b
    
14. %                 eg13f1_2.m! j\" R\" B! n: }* y9 j) C( K% `7 p
    
15. %    function u=DecisF_1(k,x)4 w) v- X\" |5 u- f, ?- m! |8 b
    
16. %    在阶段k由状态变量x的值求出其相应的决策变量所有的取值
    . @: J, }8 `9 z6 l* u/ j5 G
17. %    c=[70,72,80,76];q=10*[6,7,12,6];( l0 Y0 R$ H( U: k
    
18. %    if q(k)-x<0,u=0:100;       %决策变量不能取为负值
    7 J3 v/ S, u  ~+ E6 a
19. %    else,u=q(k)-x:100;end;     %产量满足需求且不超过1007 z$ ]- j  C. Z$ A+ b& p
    
20. %    u=u(:);  L) A' _0 P, a% Y- S
    
21. %                 eg13f2_2.m
    / Z0 O8 Q2 z. l! t
22. %    function v=ObjF_1(k,x,u)2 e; [2 ~; S' y7 {' ?9 J# u
    
23. %    阶段k的指标函数$ ~0 s0 J7 _. v* e' {* E9 Y8 Y/ A
    
24. %    c=[70,72,80,76];v=c(k)*u+2*x;4 K2 s$ L& R$ L0 I
    
25. %                 eg13f3_2.m
    : [# J4 D. z  K3 i9 k7 l
26. %    function y=TransF_1(k,x,u)
    0 W8 b: v7 T' K0 l- S
27. %     状态转移方程1 |/ T3 g: b  t, w
    
28. %     q=10*[6,7,12,6];y=x+u-q(k);9 s+ W  G& a5 k  i2 y4 U' E1 c5 \
    
29. %调用DynProg.m计算如下：: e  V\" i& w5 L/ c
    
30. %    clear;x=nan*ones(14,4);% x是10的倍数，最大范围0≤x≤130,
    5 v% Z- j1 v/ W/ d& C0 r- C5 O
31. %       %因此x=0,1,...13，所以x初始化取14行，nan表示无意义元素# _8 i3 z: d. m2 R4 K( M
    
32. %    x(1:7,1)=10*(0:6)';     % 按月定义x的可能取值: [- E: m. m( v; v\" V* R( S
    
33. %    x(1:11,2)=10*(0:10)';x(1:12,3)=10*(2:13)';: e; c. W3 H5 O. ~$ W7 y8 r
    
34. %    x(1:7,4)=10*(0:6)';1 [% H/ X3 m: r: x\" d\" ~
    
35. %    [p,f]=dynprog(x,'eg13f1_2','eg13f2_2','eg13f3_2')) K* v' a8 e2 e0 i+ y8 d
    
36. 
    , j3 h! v! |& P( p
37. % By X.D. Ding June 2000: Q& n( O4 W/ T5 k4 H. P
    
38. 
    ; Z* u. ^3 H, o) s3 q5 O# b
39. k=length(x(1,:));f_opt=nan*ones(size(x));d_opt=f_opt;9 t9 N$ z0 j/ C3 }
    
40. t_vubm=inf*ones(size(x));x_isnan=~isnan(x);t_vub=inf;6 l/ s; {9 d( J( h( G
    
41. % 计算终端相关值9 w, C3 V9 |1 d7 J( p
    
42. tmp1=find(x_isnan(:,k));tmp2=length(tmp1);( W. a# K1 |+ t( U6 o
    
43. for i=1:tmp2
    4 ]1 b  N8 ^( f( u6 }- p
44.    u=feval(DecisFun,k,x(i,k));tmp3=length(u);
    ) m# [1 h9 c$ c6 F; [. a1 D
45.    for j=1:tmp39 Q  M) {9 w& a% R/ b\" K
    
46.          tmp=feval(ObjFun,k,x(tmp1(i),k),u(j));, U2 K1 E& L6 \  P5 ]\" l
    
47.          if tmp<=t_vub, $ N9 h4 F2 a+ A; E+ C, }
    
48.             f_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vub=tmp;   \- O0 n% ^8 P) p) a
    
49. end;end;end! E' }5 G& j6 n* r# w
    
50. % 逆推计算各阶段的递归调用程序
    % `& {/ O* s2 U& P8 \1 J
51. for ii=k-1:-1:1) R: J$ ^3 r! y$ \* Q6 I. h
    
52.    tmp10=find(x_isnan(:,ii));tmp20=length(tmp10);
    ( r9 R, }. X\" d: F+ M) M
53.    for i=1:tmp20/ L9 G( r( c; K\" U# N% v: R
    
54.       u=feval(DecisFun,ii,x(i,ii));tmp30=length(u);
    $ D' ?: ]6 e% \8 ^/ A9 V
55.       for j=1:tmp30' c6 m; B  {  Y
    
56.          tmp00=feval(ObjFun,ii,x(tmp10(i),ii),u(j));3 z* q2 z; n3 Y* v
    
57.          tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j));6 X5 X! q: q0 D\" k% \, ?
    
58.          tmp50=x(:,ii+1)-tmp40;
    $ a3 u) P, z4 T5 N
59.          tmp60=find(tmp50==0);$ n; `, a' _  x2 t
    
60.          if ~isempty(tmp60),\" @; p, S\" q5 G4 I
    
61.             tmp00=tmp00+f_opt(tmp60(1),ii+1);   
    ; A8 \7 J2 W% l\" X- Y/ @
62.             if tmp00<=t_vubm(i,ii)) U8 t) d# W1 j
    
63.                f_opt(i,ii)=tmp00;d_opt(i,ii)=u(j);3 u& V: Y9 v/ t9 P3 B
    
64.                t_vubm(i,ii)=tmp00;
    0 v. Q+ ~2 n9 z2 @& @5 R9 ^7 ]0 F
65. end;end;end;end;end;9 q( _5 G& t4 I4 b/ w& x
    
66. fval=f_opt(tmp1,1);
    2 d( V4 F0 a- X/ n\" n
67. % 记录最优决策、最优轨线和相应指标函数值$ E* n\" r/ n& E$ m+ o\" v\" f
    
68. p_opt=[];tmpx=[];tmpd=[];tmpf=[];5 w8 H4 ?# f\" [% `, |# w  @# ]
    
69. tmp0=find(x_isnan(:,1));tmp01=length(tmp0);
    7 o9 F9 N' D6 p$ @/ Q  Q/ H; ?
70. for i=1:tmp01,/ L0 R8 h. r; l7 d) e
    
71.   tmpd(i)=d_opt(tmp0(i),1);
    , u5 X: W# S6 `
72.   tmpx(i)=x(tmp0(i),1);- K8 F- f* h) |5 T2 X- e
    
73.   tmpf(i)=feval(ObjFun,1,tmpx(i),tmpd(i));
    8 R# J* D9 _% _- U; M
74.   p_opt(k*(i-1)+1,[1,2,3,4])=[1,tmpx(i),.... u+ A  E( A+ i9 M+ }. q+ B
    
75. tmpd(i),tmpf(i)];# U4 p' _' S9 o/ n2 A; C/ O
    
76.   for ii=2:k
    ) a) S( N: l% Q: N& q3 [
77.      tmpx(i)=feval(TransFun,ii-1,tmpx(i),tmpd(i));# {2 F( \, i6 c9 w3 U6 R' k* D
    
78.      tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1==0);$ b) e2 x+ |5 L7 M3 K
    
79.      if ~isempty(tmp2)7 {6 J7 y( w+ K% E1 k( U& ]* }\" k
    
80.         tmpd(i)=d_opt(tmp2(1),ii);
    & C  I* v; ^  w, K# {
81.      end;
    \" D' ?. G0 q2 @2 @. W: a0 t9 Y
82.      tmpf(i)=feval(ObjFun,ii,tmpx(i),tmpd(i));% y9 D6 n: T0 P& `; ~4 u7 H* y8 s
    
83.      p_opt(k*(i-1)+ii,[1,2,3,4])=[ii,tmpx(i),...
    4 @5 ^/ q' |2 A# I+ N5 n
84. tmpd(i),tmpf(i)];# u. W, m! H  b5 [0 l% z
    
85. end;end;/ o* p3 c# K) Y2 S& N9 r1 t
    

复制代码

m存在

张坤14

hugoczx 发表于 2009-12-21 23:20
有是有                              
1 q( C, W. j, J  g8 l' H) }不过不能白给！

/ X! ?$ E- Y) t7 E0 P" c5 N同求动态规划的源程序，现在其他的算法基本上都可以集成，就是动态规划没有了。希望大神赐教，我的qq邮箱1240258213@qq.com，谢谢。关于其他算法的集成， 可以与大神共享。

无语的昭昭

额额额，其他地方很多啊