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

danieldu

看看看看那

madio

1. function [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)       / B* t& ]% P$ @, b6 F* o; |$ M
   
2. % [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)
   1 h) t  d% v4 l+ j
3. % 自由始端和终端的动态规划,求指标函数最小值的逆序算法递归8 I+ d; R6 J$ x9 W
   
4. % 计算程序。x是状态变量，一列代表一个阶段状态；M-函数8 `+ P! A0 y9 `5 x5 k* ?( k
   
5. % DecisFun(k,x)由阶段k的状态变量x求出相应的允许决策变量;
   1 L& Z; ]6 z' f, k( t
6. % M-函数ObjFun(k,x,u)是阶段指标函数，M-函数TransFun(k,x,u)  a, V5 i% o: q; S: ]
   
7. % 是状态转移函数,其中x是阶段k的某状态变量，u是相应的决策变量；
   % H% H  z7 o4 |( O* G& K
8. % 输出p_opt由4列构成，p_opt=[序号组;最优策略组;最优轨线组;
   % S  s& o# X. W0 L
9. % 指标函数值组]；fval是一个列向量，各元素分别表示p_opt各
   3 _$ c0 g$ l; e
10. % 最优策略组对应始端状态x的最优函数值；
    ; C) Y' P5 D( `\" z9 R
11. %. d' Q  l  b: N% q' K( l1 _  g
    
12. %例（参看胡良剑等编《数学实验--使用MATLAB》P180
    & z$ G! \3 p6 {* Z* d4 r1 k
13. %先写3个函数
    4 h7 w3 {& U\" c* w8 i
14. %                 eg13f1_2.m' _1 |5 {; G4 A) t\" u
    
15. %    function u=DecisF_1(k,x)8 Z* |/ \8 z' j) c2 s\" n5 e
    
16. %    在阶段k由状态变量x的值求出其相应的决策变量所有的取值
    0 I9 b2 j: ]1 f: b
17. %    c=[70,72,80,76];q=10*[6,7,12,6];' e9 S; U7 Z6 o% f5 T0 }
    
18. %    if q(k)-x<0,u=0:100;       %决策变量不能取为负值
    - d& l\" m6 S7 a; `# Z. \' M5 P  L
19. %    else,u=q(k)-x:100;end;     %产量满足需求且不超过100
    ' K/ I, }& M3 J6 j4 R$ F- ?
20. %    u=u(:);
    & ^8 i# ]( T, {( c* M: f
21. %                 eg13f2_2.m9 g$ J4 G) X\" q: e# J, \, a% n( n
    
22. %    function v=ObjF_1(k,x,u)
    7 @: {) c) D. c8 G& z: y# h; p
23. %    阶段k的指标函数
    ) u: h+ G+ o8 r6 \
24. %    c=[70,72,80,76];v=c(k)*u+2*x;
    ; l( H5 D8 n7 d! }2 I
25. %                 eg13f3_2.m
    ' X# Y& ^0 X4 m; A
26. %    function y=TransF_1(k,x,u)7 B0 p$ b' h( o2 U; n9 H
    
27. %     状态转移方程\" X/ T9 k! }\" p' z+ A, h
    
28. %     q=10*[6,7,12,6];y=x+u-q(k);7 ^7 P3 G+ S2 d1 b  u( P
    
29. %调用DynProg.m计算如下：
    . ?$ Z5 e7 @4 U4 J: r
30. %    clear;x=nan*ones(14,4);% x是10的倍数，最大范围0≤x≤130,
    8 l. E- x) W! W7 ^
31. %       %因此x=0,1,...13，所以x初始化取14行，nan表示无意义元素
    % k8 D- M( B; y: u/ b, \+ L9 |: F
32. %    x(1:7,1)=10*(0:6)';     % 按月定义x的可能取值
      z: w5 d, t+ k! S
33. %    x(1:11,2)=10*(0:10)';x(1:12,3)=10*(2:13)';/ w. ~- y) A$ E; }
    
34. %    x(1:7,4)=10*(0:6)';) P. L4 _6 f$ V* w/ O  `& W$ c: _% @8 b
    
35. %    [p,f]=dynprog(x,'eg13f1_2','eg13f2_2','eg13f3_2')3 l& e: b- H\" J2 P6 X
    
36. : p# s/ d3 V  n
    
37. % By X.D. Ding June 2000. w( u6 P% k! o
    
38. 
    - H, M9 }2 z% O( z9 C* n) q
39. k=length(x(1,:));f_opt=nan*ones(size(x));d_opt=f_opt;8 N+ @$ ]. _1 I5 W- r- h
    
40. t_vubm=inf*ones(size(x));x_isnan=~isnan(x);t_vub=inf;' @1 L+ S# j! H. W
    
41. % 计算终端相关值
    * M6 N. w7 g6 O0 `8 P
42. tmp1=find(x_isnan(:,k));tmp2=length(tmp1);
    & y. V$ R7 s% D( B1 h  A8 O
43. for i=1:tmp26 M' z8 p: w9 y8 ~; Q4 @: C
    
44.    u=feval(DecisFun,k,x(i,k));tmp3=length(u);
    9 ~+ E: @. V- M9 ?  j
45.    for j=1:tmp3
    ; n! o5 Z2 u5 P% d3 }
46.          tmp=feval(ObjFun,k,x(tmp1(i),k),u(j));6 I: ^2 j8 a4 }' O# K# N
    
47.          if tmp<=t_vub,
    1 m& c- Q. y. J1 h; s- M
48.             f_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vub=tmp;
    * \- M! x. ]% o# {$ k- i
49. end;end;end
    ; w. Q/ C/ n: M$ V
50. % 逆推计算各阶段的递归调用程序
    / y: R9 @; O6 h, V
51. for ii=k-1:-1:1
    , z. s2 v7 d4 F' H
52.    tmp10=find(x_isnan(:,ii));tmp20=length(tmp10);
    - T+ l3 M\" g: G  U. M* E
53.    for i=1:tmp20, w4 u+ N( K\" `4 @9 J
    
54.       u=feval(DecisFun,ii,x(i,ii));tmp30=length(u);\" G. q( N, l1 M( e7 x; [; b
    
55.       for j=1:tmp30
    - r\" d. J, Z- [
56.          tmp00=feval(ObjFun,ii,x(tmp10(i),ii),u(j));; x- S: g7 [# F
    
57.          tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j));; w% O0 b' U0 T5 A5 x  G
    
58.          tmp50=x(:,ii+1)-tmp40;
    ! L2 P0 i/ @9 T: K
59.          tmp60=find(tmp50==0);
    / ]9 R1 ~' K* u& |
60.          if ~isempty(tmp60),
    $ T! K7 S, S, F8 E- U
61.             tmp00=tmp00+f_opt(tmp60(1),ii+1);   
    0 d7 z6 y& k\" y! R
62.             if tmp00<=t_vubm(i,ii)
    % i  C; d6 C5 j, {2 [\" [
63.                f_opt(i,ii)=tmp00;d_opt(i,ii)=u(j);
    5 z4 ~! E9 p/ L% E
64.                t_vubm(i,ii)=tmp00;
    . m. \7 a; ?6 }1 ]. j5 x
65. end;end;end;end;end;
      I/ n: P9 a* r5 {, n\" O$ U\" Y
66. fval=f_opt(tmp1,1);
    4 b3 j3 Q+ z8 I; K% U0 e: U% t
67. % 记录最优决策、最优轨线和相应指标函数值
    # [2 |- U3 K4 m# c' }, U/ ^! r/ O3 S
68. p_opt=[];tmpx=[];tmpd=[];tmpf=[];
    7 F! p7 X& `4 {/ w( ^) V/ R: X
69. tmp0=find(x_isnan(:,1));tmp01=length(tmp0);( k  P1 [8 Q* k
    
70. for i=1:tmp01,& h7 a4 t0 J0 _$ v# G' g
    
71.   tmpd(i)=d_opt(tmp0(i),1); 0 R, [6 B& j  c- K! Y( }6 K/ P$ q
    
72.   tmpx(i)=x(tmp0(i),1);( B# N3 V3 C; O, [6 L1 N+ _
    
73.   tmpf(i)=feval(ObjFun,1,tmpx(i),tmpd(i));
    ; ]\" C7 ^2 Z; N; K3 o  K
74.   p_opt(k*(i-1)+1,[1,2,3,4])=[1,tmpx(i),...! X/ j- d( @' d5 {  \$ ^0 x
    
75. tmpd(i),tmpf(i)];
    , v- p2 A\" }1 V
76.   for ii=2:k: V\" A- E% i( j
    
77.      tmpx(i)=feval(TransFun,ii-1,tmpx(i),tmpd(i));
    ; c# R: L3 G4 O; Y' M$ }
78.      tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1==0);! g/ M& }, C( ~
    
79.      if ~isempty(tmp2)
    7 C0 [% ~' A% ]# g# i; D. l  u
80.         tmpd(i)=d_opt(tmp2(1),ii);
    9 B/ q2 x* W) x) p
81.      end;: E/ z2 L6 N' v\" L# n
    
82.      tmpf(i)=feval(ObjFun,ii,tmpx(i),tmpd(i));
    1 X' F2 m( J! d' {9 H$ C7 ]
83.      p_opt(k*(i-1)+ii,[1,2,3,4])=[ii,tmpx(i),...
    - g& }- x# B% s( ^& {7 B\" A
84. tmpd(i),tmpf(i)];
    6 T4 t/ C! V* M% y8 D
85. end;end;
    2 ?- M- U2 v\" w# ]9 T

复制代码

m存在

张坤14

hugoczx 发表于 2009-12-21 23:20
有是有                              
# k1 h+ y8 h. J* x不过不能白给！

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

无语的昭昭

额额额，其他地方很多啊
4 [6 [) m0 i3 r