逆序法代码

1 dpReverse函数

function [opt,fval]=dpReverse(S,DecisFun,SubObjFun,TransFun,ObjFun)
%动态规划逆序法
%变量或函数名       类型    　释义
%--------------------------------------------------------------------------
% S               | 矩阵 | 状态变量，每一列代表一个阶段的状态，空缺处用NAN补齐
% DecisFun(k,s)   | 函数 | 代入阶段k及状态变量s求出允许决策集合
% SubObjFun(k,s,u)| 函数 | 阶段指标函数 V(sk,uk)
% TransFun(k,s,u) | 函数 | Tk(sk,uk)，其中s是阶段k的状态，u是相应的决策变量
% ObjFun(v,f)     | 函数 | 第k阶段直到最后阶段的指标函数
% opt             | 矩阵 | [序号,最优策略,最优轨线,指标]
% fval            | 数值 | 最优解数值%第一阶段：变量初始化======================================================
N=sum(sum(~isnan(S))); %节点总个数
opt_F=ones(1,N).*inf;  %各节点到达目标点最短路长
opt_U=ones(N,N).*inf;  %各节点通向目标点最短路的下一个节点
opt_V=ones(1,N).*inf;  %各节点到通向目标点最短路的下一个节点路长
ePnt=S(:,size(S,2));   %目标点
ePnt(isnan(ePnt))=[];
opt_F(ePnt)=0;
opt_U(1,ePnt)=ePnt;
opt_V(ePnt)=0;%第二阶段：动态规划逆序求解================================================
%注：程序中以大写字母开头的变量表示集合，小写字母开头变量表示元素
for k=size(S,2)-1:-1:1Sk=S(:,k);Sk(isnan(Sk))=[];          %找到Skfor sk=Sk'                           Uk=feval(DecisFun,k,sk);         %找到各sk对应的UkUk=feval(TransFun,k,sk,Uk);      %使用TransFun找到连线节点F=ones(1,length(Uk)).*inf;       %记录各路径长度for i=1:length(Uk)uk=Uk(i);v=feval(SubObjFun,k,sk,uk); f=feval(ObjFun,v,opt_F(uk)); F(i)=f;end%找到最短路径并更新各部分数据opt_pos=find(F==min(F));opt_F(sk)=min(F);opt_U(1:length(opt_pos),sk)=Uk(opt_pos)';opt_V(sk)=feval(SubObjFun,k,sk,opt_U(1,sk));end
end
%第三阶段：路径还原========================================================
sPnt=S(:,1);%定义路径搜索起始点
sPnt(isnan(sPnt))=[];
Path=zeros(1,size(S,2)+1);      
Path(1)=sPnt;
for i=1:size(S,2)NewPath=zeros(1,size(S,2)+1);         NewPath(:,:)=[];for j=1:size(Path,1)%当改点到达目标点最小路径不止一条时%自动生成分支路径branch=opt_U(:,Path(j,i));branch(isinf(branch))=[];for m=branch'tempPath=Path(j,:);tempPath(1,i+1)=m;NewPath=[NewPath;tempPath];endend  Path=NewPath;
end
%temp_n:序号 | temp_s:最优策略 | temp_e:最优轨线 | temp_v:指标
temp_n=reshape((1:size(S,2))'*ones(1,size(Path,1)),size(S,2)*size(Path,1),1); 
temp_s=reshape(Path(:,1:size(S,2))',size(S,2)*size(Path,1),1);
temp_e=reshape(Path(:,2:size(S,2)+1)',size(S,2)*size(Path,1),1);
temp_v=opt_V(temp_s)';
opt=[temp_n,temp_s,temp_e,temp_v];
opt(find(diff(opt(:,2))==0)+1,:)=[];
fval=sum(opt(1:size(S,2),4));
end

2 应用实例

计算由V1到达V10的最短路径。

---

dpReverse函数为了应对动态规划的不同实例，需要将
状态变量，允许决策集合，阶段指标函数 ，第k阶段直到最后阶段的指标函数等部分作为参数传入dpReverse函数进行计算即：
dpReverse(S,DecisFun,SubObjFun,TransFun,ObjFun)

---

2.1 状态变量矩阵定义

即函数中的参数S,S以矩阵形式被表示,空缺部分使用nan补齐
就该图片而言，我们将其分成5个阶段，每个阶段含有的状态变量分别为{V1}，{V2，V3，V4}，{V5，V6，V7}，{V8，V9}，{V10}
转化为matlab内矩阵形式即为：

这个过程可以通过编写如下代码来实现：

S=nan.*ones(3,5);
S(1,1)=1;
S(1:3,2)=[2;3;4];
S(1:3,3)=[5;6;7];
S(1:2,4)=[8;9];
S(1,5)=10;

---

2.2 允许决策集合函数

代入阶段k及状态变量s求出允许决策集合函数
可以简单的使用switch构造：
如图中case 1,u=[2;3;4];表明节点V1的下一个节点可以为V2，V3，V4，由于该实例较为简单，这里实际上只使用的了s一个参数，k参数的引入可以刻画更加复杂的问题，这里不展开讨论。

function u=dpMinDis_DecisFun(k,s)switch scase 1,u=[2;3;4];case 2,u=[5;6];case 3,u=[5;6;7];case 4,u=[6;7];case 5,u=8;   case 6,u=[8;9]; case 7,u=[8;9];     case 8,u=10;     case 9,u=10; case 10,u=10;     end
end

---

2.3 阶段指标函数 V(sk,uk)

在该题中指的便是两点之间的路长，例如case (s==6&u =8),v=2;即为V6和V8点之间距离为2。同样由于实例较为简单并未使用参数k。

function v=dpMinDis_SubObjFun(k,s,u)v=inf;switch 1case (s==6&u==8),v=2;case (s==1&u==2)|(s==4&u==7)|(s==6&u==9)|(s==8&u==10),v=3;case (s==3&u==5)|(s==5&u==8)|(s==9&u==10),v=4;case (s==1&u==4)|(s==3&u==6)|(s==4&u==6),v=5;case (s==1&u==3)|(s==3&u==7)|(s==7&u==8)|(s==7&u==9),v=6;case (s==2&u==6),v=7;case (s==2&u==5),v=8;end
end

---

2.4 状态转移方程Tk(sk,uk)

$S_{k+1}=T_k(S_k,U_k)$
对于该题来说，实际上$S_{k+1}$即为$U_k$

function s=dpMinDis_TransFun(k,s,u)s=u;
end

---

2.5 第k阶段直到最后阶段的指标函数

当前问题下f没有劳损或增值等因素，故直接采用f=v+f

function f=dpMinDis_ObjFun(v,f)f=v+f;
end

---

2.6 程序调用

通过如下方式进行程序调用：

[opt,fval]=dpReverse(S,@dpMinDis_DecisFun,@dpMinDis_SubObjFun,...@dpMinDis_TransFun,@dpMinDis_ObjFun)

2.7 计算结果

计算结果如下：

其中opt为最优路径，fval为最短距离，我们详细来看opt结果，opt最左列是用来区分路径的序号，它由1增长到5之后再次由1增长到5，说明最优结果包含两条路径，且都经过了五个点，

opt中
第二列是每一段路的起点
第三列是每一段路的终点
第四列是每一段路的路长
如图所示，下图第一行表示最短路首先由点1到达点2且经过距离为3，而整幅图说明最短路为：
1—>2—>6—>8—>10
且最短路长为3+7+2+3+0=15

因此我们从计算结果opt可以看出，该题有两个最短路，分别为:
1—>2—>6—>8—>10
1—>4—>6—>8—>10
最短路长为15