例1（最短路线问题）

程序：

a=[1 2 4 3 3 2 1];

b=[5 3 -1 -1;

1 3 6 -1;-1 8 7 6;

6 8 -1 -1;3 5 -1 -1;-1 3 3 -1;-1 8 4 -1;

2 2 -1 -1;-1 1 2 -1;-1 3 3 -1;

3 5 -1 -1;5 2 -1 -1;6 6 -1 -1;

4 -1 -1 -1;3 -1 -1 -1];

num=16;

max_col=4;

a_length=size(a);

route=zeros(a_length(2),max_col);

row=2;

c=b(1,:);

c1=zeros(1,max_col);c2=zeros(1,max_col);

for i=2:a_length(2)-1

   n=a(i);

   for k=1:max_col

       if b(row,k)>0

           c1(k)=b(row,k)+c(1);

           route(i,k)=1;

       end

   end

   if n>1

      for j=1:n-1

          for k=1:max_col

            if b(row+j,k)>0

                c2(k)=b(row+j,k)+c(j+1);

                if c2(k)<c1(k)||c1(k)==0

                    c1(k)=c2(k);

                    route(i,k)=j+1;

                end

            end

          end

      end

   end

   row=row+n;

   c=c1;

   c1=zeros(1,max_col);c2=zeros(1,max_col);

end

correct_col=route(end-1,1);

l_route=route(:,correct_col);

l_route(1)=1;

l_route(end-1)=correct_col;

l_route(end)=1;

disp('最小路径为');disp(c(1));disp('路径规划为');disp(l_route');

结果：

最小路径为

    18

路径规划为

     1     1     2     1     2     2     1

例2（生产—存贮问题）

程序：

function [p_opt,fval]=dynprog(x,DecisFun,ObjFun,TransFun)

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=DecisFun(k,x(i,k));

    tmp3=length(u);

    for j=1:tmp3

        tmp=ObjFun(k,x(tmp1(i),k),u(j));

        if tmp<=t_vubm(i,k)

            f_opt(i,k)=tmp;

            d_opt(i,k)=u(j);

            t_vubm(i,k)=tmp;

        end

    end

end

%%  逆推计算各阶段的递归调用程序

for ii=k-1:-1:1

    tmp10=find(x_isnan(:,ii));

    tmp20=length(tmp10);

    for i=1:tmp20      % 求出当前状态下，最有可能的决策

        u=DecisFun(ii,x(i,ii));

        tmp30=length(u);

        for j=1:tmp30  % 求出当前状态下所有决策的最小性能值

            tmp00=ObjFun(ii,x(tmp10(i),ii),u(j));% 单步性能指标

            tmp40=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);

%%  记录最优决策、最优轨线和相应指标函数值

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)=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)=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)=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

function V=ObjFun(k,x,u)

d = [ 2 3 2 4 ];

if u == 0

    V = 0.5 * (x + u -d(k));

else

    if u>6

        V = 10^6;   

    else

        V = 3 + u +0.5*(x + u -d(k));

    end

end

function s_next=TransFun(k,x,u)

d = [ 2 3 2 4 ];

s_next = x + u - d(k);

function u=DecisFun(k,x)

m=6;

d = [ 2 3 2 4 ];

if k == 4

    u = d(k) - x;

else

    u = max (0,d(k)-x):m;

end

x1 = 0:4;

s=nan*ones(5,1);

s(1) = 0;

x = [s x1' x1' x1'];

[p_opt,fval]=dynprog(x,@DecisFun,@ObjFun,@TransFun)

结果：

p_opt =

    1.0000         0    5.0000    9.5000

    2.0000    3.0000         0         0

    3.0000         0    6.0000   11.0000

    4.0000    4.0000         0         0

fval =

   20.5000

       NaN

       NaN

       NaN

       NaN

例3（生产—库存管理系统）

程序：

function u=DecisFun(k,x)

d =[ 600 700 500 1200];

if k == 4

    u = d(k) - x;

else

    u = max (0,d(k)-x):1200;

end

function V=ObjFun(k,x,u)

V = x+ 0.005*u^2;

function s_next=TransFun(k,x,u)

d = [ 600 700 500 1200];

s_next = x + u - d(k);

x1 =100*(0:4);

s = nan*ones(5,1);

s(1)=0;

x = [s x1' x1' x1'];

[p_opt,fval]=dynprog(x,@DecisFun,@ObjFun,@TransFun)

结果：

p_opt =

           1           0         600        1800

           2           0         700        2450

           3           0         800        3200

           4         300         900        4350

fval =

       11800

         NaN

         NaN

         NaN

         NaN

例4（企业生产管理问题）

程序：

function V=ObjFun(k,x,u)

V =8*u  + 5*(x-u);

End

function s_next=TransFun(k,x,u)

a = 0.7;

b = 0.9;

s_next =(a - b)*u + b*x;

function u=DecisFun(k,x)

u = 0 : x;

x1 = (500:999);

s=nan*ones(500,1);

s(1) = 1000;

x = [s x1' x1' x1' x1'];

[p_opt,fval]=dynprog(x,@DecisFun,@ObjFun,@TransFun)

结果：

p_opt =

           1        1000        1000        8000

           2         700          10        3530

           3         628           6        3158

           4         564           3        2829

           5         507           0        2535

fval =

       20052

例5（设备更新问题）

function V=ObjFun(k,x,u)

R = [ 5 4.5 4 3.75 3 2.5]; % 产生的经济效益

U = [ 0.5 1 1.5 2 2.5 3]; % 维修费用

C = [ 0.5 1 1.5 2 2.5 3]; % 换购一台新的设备所需要的费用

if u == 1

    V = R(x) - U(x);

else

    V = R(1) - U(1) - C(x);

End

function s_next=TransFun(k,x,u)

if u == 1

    s_next = x + 1;

else

    s_next = 1;

end

function u=DecisFun(k,x)

u = [1 0];

x1 = (1:6);

s=nan*ones(6,1);

s(1) = 1;

x = [s x1' x1' x1' x1'];

[p_opt,fval]=dynprog(x,@DecisFun,@ObjFun,@TransFun)

结果：

p_opt =

    1.0000    1.0000    1.0000    4.5000

    2.0000    2.0000    1.0000    3.5000

    3.0000    3.0000    1.0000    2.5000

    4.0000    4.0000    1.0000    1.7500

    5.0000    5.0000    1.0000    0.5000

fval =

   12.7500

       NaN

       NaN

       NaN

       NaN

       NaN