例1：最短路径

clear

clc

% user input data

% 'a':每个阶段的站点数目

% 'b':所有站点的距离

% num:所有站点的个数总和

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

b=[5 3 -1 -1;

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

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;

% c:最小路径

% c1:第一个决策点所产生的四个距离

% c2:其他决策点所产生的四个距离，并将较小值放入c1中

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（生产-存贮问题）

模型建立：这是一个多阶段决策问题，每一月份为一个阶段。

（1）状态变量Sk（k=1,2,3,4，5）表示第k个月份初具有的产品数，实际是上一月份满足订货需求后的剩余产品数，S1=0，S5表示第四月份末剩余的产品数；

（2）决策变量uk（k=1，2,3,4），表示第k个月份需要生产的产品数；

（3）状态转移方程，由状态Sk采取决策uk后的状态转移方程为：

Sk+1=Sk+uk-Ak

其中Ak表示第k月份的订货需求量；

（4）阶段效益，表示该月份的所有成本费用之和，由题意，有：

dk（Sk，uk）=3+uk+0.5Sk

（5）基本方程，用fk(Sk)表示从状态Sk出发采用最优策略到第四月份结束的最小费用，则有下面的动态规划模型

fk（Sk）=min（uk≥Ak-Sk）{dk（Sk,uk）+f（k+1）（Sk+1）}

=min（uk≥Ak-Sk）（3+uk+0.5Sk+f（k+1）（Sk+1））

f5（S5）=0，k=4,3,2,1

模型求解：按照动态规划法的基本思想（递推算法）求解以上模型。

K=4时，求：

f4（S4）=min（u4≥4-S4）{3+u4+0.5Sk}

显然应取u4=4-S4,于是得：

f4（S4）=7-0.5*S4

K=3时，此时S4=S3+u3-2，

u3=2-S3求

f3(S3)=min(u3≥2-S3){3+u3+0.5*S3+7-0.5*S4}=12.5-0.5*S3

K=2时，注意到S3=S2+u2-3，

u2=3-S2求

f2（S2）=min(u2≥3-S2){3+u2+0.5*S2+12.5-0.5*S3}=18.5-0.5*S2

K=1时，注意到S2=S1+u1-2，

显然u1=2-S1,求：

f1（S1）=min(u1≥2-S1){3+u1+0.5*S1+18.5-0.5*S2}=23.5-0.5*S1

由于S1=0,于是u1=2，得S2=0，u2=3，S3=0，u3=2，S4=0，u4=4

Matlab代码及运行结果：

function dynamic_li2

%计算各状态变量可能取值，第k列代表第k个状态变量的可能取值，没有的使用NaN代替

x1=0:4;

s=nan*ones(5,1);

s(1)=0;

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

%直接调用dynprog函数

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

% 输出p_opt由4列构成，p_opt=[序号组;最优策略组;最优轨线组;

% 指标函数值组]；fval是一个列向量，各元素分别表示p_opt各

% 最优策略组对应始端状态x的最优函数值；

function u=DecisFun(k,x)    %决策函数

d=[2 3 2 4];

m=6;

if k==4

   u=d(k)-x;

else

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

end

function f=SubObjFun(k,x,u)    %阶段效益函数

d=[2 3 2 4];

if u==0

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

else if u>6

       f=10^6;

   else

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

   end

end

function s=TransFun(k,x,u)    %状态转移函数

d=[2 3 2 4];

s=x+u-d(k);

运行结果：

>> dynamic_li2

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（生产-库存管理系统）

3、生产——库存管理系统

function dynamic_li3

x1=100*(0:4);

s=nan*ones(5,1);

s(1)=0;

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

[p_opt,fval]=dynprog(x,@DecisFun_3,@SubObjFun_3,@TransFun_3)

function V=SubObjFun_3(k,x,u)

V=x+0.005*u^2;

end

function u=DecisFun_3(k,x)

d=100*[6 7 5 12];

if k==4

   u=d(k)-x;

else

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

end

end

function s_next=TransFun_3(k,x,u)

d=[600 700 500 1200];

s_next=x+u-d(k);

end

运行结果：

dynamic_li3

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

5、设备更新问题

function dynamic_li5

x1=(1:5);

s=nan*ones(5,1);

s(1)=5;

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

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

function V=ObjFun_5(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

end

function s_next=TransFun_5(k,x,u)

if u==1

   s_next=x+1;

else

   s_next=1;

end

end

function u=DecisFun_5(k,x)

u=[1 0];

end

dynamic_li5

p_opt =

   1.0000    5.0000         0    2.0000

   2.0000    1.0000    1.0000    4.5000

   3.0000    2.0000    1.0000    3.5000

   4.0000    3.0000    1.0000    2.5000

   5.0000    4.0000    1.0000    1.7500

fval =

  14.2500

      NaN

      NaN

      NaN

      NaN