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

danieldu

看看看看那

madio

1. function [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)      
   / c3 r! g5 F+ S9 T9 T
2. % [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)
   , [* j9 p2 M7 y7 [( s' l1 @& q
3. % 自由始端和终端的动态规划,求指标函数最小值的逆序算法递归* k0 _% i; J- \' j. w2 U
   
4. % 计算程序。x是状态变量，一列代表一个阶段状态；M-函数
   ( g6 X( R4 ^  i
5. % DecisFun(k,x)由阶段k的状态变量x求出相应的允许决策变量;
   3 G) t) E% d8 z\" F& R( T
6. % M-函数ObjFun(k,x,u)是阶段指标函数，M-函数TransFun(k,x,u)6 d' T  p, L: J2 q$ o! p4 A: `$ H
   
7. % 是状态转移函数,其中x是阶段k的某状态变量，u是相应的决策变量；
   9 D$ E1 ^7 F8 Y\" }* s; A
8. % 输出p_opt由4列构成，p_opt=[序号组;最优策略组;最优轨线组;
     a; w7 E/ Y+ j
9. % 指标函数值组]；fval是一个列向量，各元素分别表示p_opt各/ M3 F$ i3 U: A0 k8 s
   
10. % 最优策略组对应始端状态x的最优函数值；) F2 C5 ]2 C0 k8 r
    
11. %) v7 x) O: M9 i\" F$ |5 _9 B8 j! U
    
12. %例（参看胡良剑等编《数学实验--使用MATLAB》P180
    3 L% w2 H6 U. T0 p/ u+ ?
13. %先写3个函数
      A3 h# F8 C+ @0 g* `+ D
14. %                 eg13f1_2.m
    7 |, c0 q0 h* J9 b6 i+ e8 d
15. %    function u=DecisF_1(k,x)
    , g: `4 T1 ^) ~* B  }0 K\" e
16. %    在阶段k由状态变量x的值求出其相应的决策变量所有的取值# |+ I1 }0 Q3 H: U+ ]
    
17. %    c=[70,72,80,76];q=10*[6,7,12,6];: ^0 |6 _& t9 S3 }+ p7 L
    
18. %    if q(k)-x<0,u=0:100;       %决策变量不能取为负值
    ; q! S  J3 o# v$ S
19. %    else,u=q(k)-x:100;end;     %产量满足需求且不超过100/ l0 v3 z% y3 @: H1 C
    
20. %    u=u(:);
    ' D$ ^' y+ o( f1 g3 q
21. %                 eg13f2_2.m
    * E/ q( I, B! p0 p2 F
22. %    function v=ObjF_1(k,x,u)
    6 L; h0 f& v2 J% l3 a
23. %    阶段k的指标函数\" `4 h3 Z  Z) m0 z1 `, ?) ?
    
24. %    c=[70,72,80,76];v=c(k)*u+2*x;/ _( X0 o& k. T\" y! A$ b% {
    
25. %                 eg13f3_2.m
    6 C/ a6 p, T# C4 F8 e
26. %    function y=TransF_1(k,x,u)8 r9 L, _/ [+ X. e- P. m, @& o
    
27. %     状态转移方程
    ! h: L. l; ~& x% ~
28. %     q=10*[6,7,12,6];y=x+u-q(k);
    + O6 N- N1 B3 p, Y9 b6 K) x
29. %调用DynProg.m计算如下：
    , x/ \! g\" x5 n2 S
30. %    clear;x=nan*ones(14,4);% x是10的倍数，最大范围0≤x≤130,- q+ P: Q7 I: G; x
    
31. %       %因此x=0,1,...13，所以x初始化取14行，nan表示无意义元素# d* S) F! h0 ^9 P0 ]
    
32. %    x(1:7,1)=10*(0:6)';     % 按月定义x的可能取值
      [/ u, K  U0 F. S& r
33. %    x(1:11,2)=10*(0:10)';x(1:12,3)=10*(2:13)';\" z0 ]0 }3 \$ \7 Y8 ~& e; u- p
    
34. %    x(1:7,4)=10*(0:6)';$ m: R3 {8 i* `3 Z' S3 ]
    
35. %    [p,f]=dynprog(x,'eg13f1_2','eg13f2_2','eg13f3_2')8 S1 m, b4 ~8 D' Z6 Y$ R( h
    
36. - D1 w: D& A9 K- a* t
    
37. % By X.D. Ding June 2000# _\" ^' w/ m1 {# ]% q/ \* O; p! l9 D$ p
    
38. / m5 X: t( F/ Z% ?
    
39. k=length(x(1,:));f_opt=nan*ones(size(x));d_opt=f_opt;
    5 \6 C+ K- W9 ^! c, `6 M( Y4 k
40. t_vubm=inf*ones(size(x));x_isnan=~isnan(x);t_vub=inf;3 x  e+ `$ H& o; y7 ^
    
41. % 计算终端相关值. p! f3 M/ p8 v; e1 ~
    
42. tmp1=find(x_isnan(:,k));tmp2=length(tmp1);
    ( u7 ~4 d& _+ |( E
43. for i=1:tmp2
    . N0 ^4 h, {! n  ?; n7 H
44.    u=feval(DecisFun,k,x(i,k));tmp3=length(u);- d; z- o\" h# f/ e. M- s
    
45.    for j=1:tmp3
    1 [( m% c# t\" d. H  e6 Z2 a5 s
46.          tmp=feval(ObjFun,k,x(tmp1(i),k),u(j));
    % s( f1 X2 x  D. W
47.          if tmp<=t_vub,
    4 L/ V- ]. i- I9 {! [
48.             f_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vub=tmp;
    \" J- f5 [2 H7 Q\" ?) b; R
49. end;end;end
    % q/ O! b. P) `9 v
50. % 逆推计算各阶段的递归调用程序
    & J* t/ }, L0 _
51. for ii=k-1:-1:1
    ( n+ I7 r! o0 L, C
52.    tmp10=find(x_isnan(:,ii));tmp20=length(tmp10);: o. I/ b7 _\" B! k- R0 J; o/ r
    
53.    for i=1:tmp20, d4 }% ~: ^$ _
    
54.       u=feval(DecisFun,ii,x(i,ii));tmp30=length(u);
    , m\" y3 [  H6 K* t/ J) ~+ B4 f
55.       for j=1:tmp30
    / ~; ^0 t. K  m# }
56.          tmp00=feval(ObjFun,ii,x(tmp10(i),ii),u(j));2 |& e- P& n5 s3 l
    
57.          tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j));
    7 C/ ]4 T0 H, S4 V& C
58.          tmp50=x(:,ii+1)-tmp40;
    $ B+ O, R$ o# s9 ^9 Q7 M% Y
59.          tmp60=find(tmp50==0);8 M3 \2 f$ _\" b8 J) j
    
60.          if ~isempty(tmp60),\" \, X7 r3 a4 B$ S\" B
    
61.             tmp00=tmp00+f_opt(tmp60(1),ii+1);   
    % V! _0 M0 `$ L, ]$ V3 c
62.             if tmp00<=t_vubm(i,ii)
    , y# ~% ^! v2 n# W' Y2 W! e  x
63.                f_opt(i,ii)=tmp00;d_opt(i,ii)=u(j);6 D. _- O% A\" n; m3 z
    
64.                t_vubm(i,ii)=tmp00;8 B: g4 j) a7 `2 R; w7 t% @- P  B
    
65. end;end;end;end;end;0 c; R$ I8 U; U6 N
    
66. fval=f_opt(tmp1,1);  O/ s  J\" k! O) v\" L( _( D
    
67. % 记录最优决策、最优轨线和相应指标函数值
    7 G5 _2 j7 v( T
68. p_opt=[];tmpx=[];tmpd=[];tmpf=[];
    & o' ~# m1 F2 H4 a' _
69. tmp0=find(x_isnan(:,1));tmp01=length(tmp0);
    4 N0 i' Z1 ?; S  z. P3 l( |
70. for i=1:tmp01,
    7 d4 a. n% p& E- j8 t  O) \+ s
71.   tmpd(i)=d_opt(tmp0(i),1);
    8 f+ e# t: X* k  z  X# _/ C0 ?- A
72.   tmpx(i)=x(tmp0(i),1);% q% |9 C0 f5 f! k2 D
    
73.   tmpf(i)=feval(ObjFun,1,tmpx(i),tmpd(i));9 G- f, ~$ I8 G7 _
    
74.   p_opt(k*(i-1)+1,[1,2,3,4])=[1,tmpx(i),...
    - h- x: ]. u( p8 P; z, u$ Y' d
75. tmpd(i),tmpf(i)];
    0 N. _7 o5 \* Q! C2 k6 d
76.   for ii=2:k
    # X9 q! T; A3 e& ]: X0 V, p
77.      tmpx(i)=feval(TransFun,ii-1,tmpx(i),tmpd(i));
    # k9 M% F/ g/ I# |( h5 {
78.      tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1==0);
    ( d# h+ M  |: o# ^7 L! V
79.      if ~isempty(tmp2)' W: v; s6 W& _# b) s
    
80.         tmpd(i)=d_opt(tmp2(1),ii);. \' ]) [' v5 _: X% B\" q/ K
    
81.      end;: D4 F+ d& s3 J9 p! K- v# C
    
82.      tmpf(i)=feval(ObjFun,ii,tmpx(i),tmpd(i));1 }7 q/ `- w1 e+ U9 O
    
83.      p_opt(k*(i-1)+ii,[1,2,3,4])=[ii,tmpx(i),...* O/ n' m1 m. J8 M, [# z+ O/ w9 _  K
    
84. tmpd(i),tmpf(i)];
    ' s- ]5 P# x\" R7 ]8 K
85. end;end;
    5 ?) c; \. P1 a( A8 y+ ^, ^

复制代码

m存在

张坤14

hugoczx 发表于 2009-12-21 23:20
有是有                              
5 u/ g' Z# `7 R2 {不过不能白给！

同求动态规划的源程序，现在其他的算法基本上都可以集成，就是动态规划没有了。希望大神赐教，我的qq邮箱1240258213@qq.com，谢谢。关于其他算法的集成， 可以与大神共享。

无语的昭昭

额额额，其他地方很多啊
' V$ ]4 N5 X. }6 S; u  d) k! ?6 E