function [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)       
% [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)
% 自由始端和終端的動態(tài)規(guī)劃,求指標函數(shù)最小值的逆序算法遞歸
% 計算程序。x是狀態(tài)變量，一列代表一個階段狀態(tài)；M-函數(shù)
% DecisFun(k,x)由階段k的狀態(tài)變量x求出相應的允許決策變量;
% M-函數(shù)ObjFun(k,x,u)是階段指標函數(shù)，M-函數(shù)TransFun(k,x,u)
% 是狀態(tài)轉(zhuǎn)移函數(shù),其中x是階段k的某狀態(tài)變量，u是相應的決策變量；
% 輸出p_opt由4列構(gòu)成，p_opt=[序號組;最優(yōu)策略組;最優(yōu)軌線組;
% 指標函數(shù)值組]；fval是一個列向量，各元素分別表示p_opt各
% 最優(yōu)策略組對應始端狀態(tài)x的最優(yōu)函數(shù)值；
%
%例（參看胡良劍等編《數(shù)學實驗--使用MATLAB》P180
%先寫3個函數(shù)
%                 eg13f1_2.m
%    function u=DecisF_1(k,x)
%    在階段k由狀態(tài)變量x的值求出其相應的決策變量所有的取值
%    c=[70,72,80,76];q=10*[6,7,12,6];
%    if q(k)-x<0,u=0:100;       %決策變量不能取為負值
%    else,u=q(k)-x:100;end;     %產(chǎn)量滿足需求且不超過100
%    u=u(:);
%                 eg13f2_2.m
%    function v=ObjF_1(k,x,u)
%    階段k的指標函數(shù)
%    c=[70,72,80,76];v=c(k)*u+2*x;
%                 eg13f3_2.m
%    function y=TransF_1(k,x,u)
%     狀態(tài)轉(zhuǎn)移方程
%     q=10*[6,7,12,6];y=x+u-q(k);
%調(diào)用DynProg.m計算如下：
%    clear;x=nan*ones(14,4);% x是10的倍數(shù)，最大范圍0≤x≤130,
%       %因此x=0,1,...13，所以x初始化取14行，nan表示無意義元素
%    x(1:7,1)=10*(0:6)';     % 按月定義x的可能取值
%    x(1:11,2)=10*(0:10)';x(1:12,3)=10*(2:13)';
%    x(1:7,4)=10*(0:6)';
%    [p,f]=dynprog(x,'eg13f1_2','eg13f2_2','eg13f3_2')

% By X.D. Ding June 2000

k=length(x(1,:));f_opt=nan*ones(size(x));d_opt=f_opt;
t_vubm=inf*ones(size(x));x_isnan=~isnan(x);t_vub=inf;
% 計算終端相關值
tmp1=find(x_isnan(:,k));tmp2=length(tmp1);
for i=1:tmp2
   u=feval(DecisFun,k,x(i,k));tmp3=length(u);
   for j=1:tmp3
         tmp=feval(ObjFun,k,x(tmp1(i),k),u(j));
         if tmp<=t_vub, 
            f_opt(i,k)=tmp;d_opt(i,k)=u(j);t_vub=tmp; 
end;end;end
% 逆推計算各階段的遞歸調(diào)用程序
for ii=k-1:-1:1
   tmp10=find(x_isnan(:,ii));tmp20=length(tmp10);
   for i=1:tmp20
      u=feval(DecisFun,ii,x(i,ii));tmp30=length(u);
      for j=1:tmp30
         tmp00=feval(ObjFun,ii,x(tmp10(i),ii),u(j));
         tmp40=feval(TransFun,ii,x(tmp10(i),ii),u(j));
         tmp50=x(:,ii+1)-tmp40;
         tmp60=find(tmp50==0);
         if ~isempty(tmp60),
            tmp00=tmp00+f_opt(tmp60(1),ii+1);    
            if tmp00<=t_vubm(i,ii)
               f_opt(i,ii)=tmp00;d_opt(i,ii)=u(j);
               t_vubm(i,ii)=tmp00;
end;end;end;end;end;
fval=f_opt(tmp1,1);
% 記錄最優(yōu)決策、最優(yōu)軌線和相應指標函數(shù)值
p_opt=[];tmpx=[];tmpd=[];tmpf=[];
tmp0=find(x_isnan(:,1));tmp01=length(tmp0);
for i=1:tmp01,
  tmpd(i)=d_opt(tmp0(i),1); 
  tmpx(i)=x(tmp0(i),1);
  tmpf(i)=feval(ObjFun,1,tmpx(i),tmpd(i));
  p_opt(k*(i-1)+1,[1,2,3,4])=[1,tmpx(i),...
tmpd(i),tmpf(i)];
  for ii=2:k
     tmpx(i)=feval(TransFun,ii-1,tmpx(i),tmpd(i));
     tmp1=x(:,ii)-tmpx(i);tmp2=find(tmp1==0);
     if ~isempty(tmp2)
        tmpd(i)=d_opt(tmp2(1),ii);
     end;
     tmpf(i)=feval(ObjFun,ii,tmpx(i),tmpd(i));
     p_opt(k*(i-1)+ii,[1,2,3,4])=[ii,tmpx(i),...
tmpd(i),tmpf(i)];
end;end;