1、- 1 -1 流量工程问题1.1 问题重述定义一个有向网络 G=(N,E)，其中 N 是节点集，E 是弧集。令 A 是网络 G 的点弧关联矩阵，即 NE 阶矩阵，且第 l 列与弧里(I,j)对应，仅第 i 行元素为 1，第 j 行元素为-1，其余元素为 0。再令 bm=(bm1,bmN)T，f m=(fm1,fmE)T，则可将等式约束表示成：Afm=bm本算例为一经典 TE 算例。算例网络有 7 个节点和 13 条弧，每条弧的容量是 5 个单位。此外有四个需求量均为 4 个单位的源一目的对，具体的源节点、目的节点信息如图所示。这里为了简单，省区了未用到的弧。此外，弧上的数字表示弧的编号。此时，
2、c=(5,5,5)113)T，根据上述四个约束条件，分别求得四个情况下的最优决策变量 x=(x12,x13,x75)113)。图 1 网络拓扑和流量需求- 2 -1.2 7 节点算例求解1.2.1 算例 1（b1=4;-4;0;0;0;0;0 T）转化为线性规划问题：Minimize cTx1Subject to Ax1=b1x1=0利用 Matlab 编写对偶单纯形法程序，可求得:最优解为 x1*=4 0 0 0 0 0 0 0 0 0 0 0 0T对应的最优值 cTx1=201.2.2 算例 2（b2=4;0;-4;0;0;0;0 T）Minimize cTx2Subject to Ax2
3、=b2X2=0利用 Matlab 编写对偶单纯形法程序，可求得:最优解为 x2*=0 4 0 0 0 0 0 0 0 0 0 0 0T对应的最优值 cTx2=201.2.3 算例 3（b3=0;-4;4;0;0;0;0 T）Minimize cTx3Subject to Ax3=b3X3=0利用 Matlab 编写对偶单纯形法程序，可求得:最优解为 x3*=4 0 0 0 4 0 0 0 0 0 0 0 0T对应的最优值 cTx3=40- 3 -1.2.4 算例 4（b4=4;0;0;0;0;0;-4 T）Minimize cTx4Subject to Ax4=b4X4=0利用 Matlab
4、编写对偶单纯形法程序，可求得:最优解为 x4*=4 0 0 4 0 0 0 0 0 4 0 0 0T对应的最优值 cTx4=601.3 计算结果及结果说明1.3.1 算例 1（b1=4;-4;0;0;0;0;0 T）算例 1 中，由 b1 可知，节点 2 为需求节点，节点 1 为供给节点，由节点 1 将信息传输至节点 2 的最短路径为弧 1。图 2 算例 1 最优传输示意图求得的最优解为 x1*=4 0 0 0 0 0 0 0 0 0 0 0 0T，即只经过弧 1 运输 4 个单位流量，其余弧无流量。又因为，每条弧的费用均为 5，所以最小费用为 20。经分析，计算结果合理可信。1.3.2 算例
5、 2（b2=4;0;-4;0;0;0;0 T）算例 2 中，由 b2 可知，节点 3 为需求节点，节点 1 为供给节点，由节点 1 将信息传输至节点 2 的最短路径为弧 2。- 4 -图 3 算例 2 最优传输示意图求得的最优解为 x2*=0 4 0 0 0 0 0 0 0 0 0 0 0T，即只经过弧 2 运输 4 个单位流量，其余弧无流量。又因为，每条弧的费用均为 5，所以最小费用为 20。经分析，计算结果合理可信。1.3.3 算例 3（b3=0;-4;4;0;0;0;0 T）算例 3 中，由 b3 可知，节点 2 为需求节点，节点 3 为供给节点，由节点 3 将信息传输至节点 2 的最短
6、路径为弧 5-弧 1。图 4 算例 3 最优传输示意图求得的最优解为 x3*=4 0 0 0 4 0 0 0 0 0 0 0 0T，即经过弧 5 运输 4 个单位流量至节点 1，再经弧 1 运输 4 个单位流量至节点 2，其余弧无流量。又因为，每条弧的费用均为5，所以最小费用为 40。经分析，计算结果合理可信。1.3.4 算例 4（b4=4;0;0;0;0;0;-4 T）算例 4 中，由 b4 可知，节点 7 为需求节点，节点 1 为供给节点，由节点 1 将信息传输至节点 7 的最短路径为弧 1-弧 4-弧 10。- 5 -图 5 算例 3 最优传输示意图求得的最优解为 x4*=4 0 0 4
7、 0 0 0 0 0 4 0 0 0T，即经过弧 1 运输 4 个单位流量至节点 2，再经弧 4 运输 4 个单位流量至节点 5，最后经弧 5 运输 4 个单位流量至节点 7，其余弧无流量。又因为，每条弧的费用均为 5，所以最小费用为 60。经分析，计算结果合理可信。- 6 -2 重要算法编写与观察2.1 习题 5.6（a） 初值为（0,0）时本算法令 g 的 2 范数在=0 flag=0;sol=zeros(1,nA); for i=1:mA-1sol(N(i)=A(i,nA);endval=sol*(B(mA,:);elsefor i=1:mA-1if A(i,nA)=0 disp(hav
8、e infinite solution!);flag=0;break;endendif flag temp=0;for i=1:mA-1if A(i,nA)temptemp=sita(i);inb=i; endendfor i=1:mA-1if i=outbN(i)=inb;endendA(outb,:)=A(outb,:)/A(outb,inb); for i=1:mAif i=outb A(i,:)=A(i,:)-A(outb,:)*A(i,inb);A(mA,nA)=0; endendendendend3.2 最速降线法求 Rosenbrock 函数最小值 matlab 程序如下：fun
9、ction rb = rbfun(x,y)rb=100*(y-x2)2+(1-x)2end- 17 -clearclcsyms x y g Gg=gradient(rb(x,y),x y) %定义梯度向量G=hessian(rb(x,y),x y) %定义海森阵X(1,:)=-1.4 1; %定义初始点x=X(1,1);y=X(1,4); A(1,:)=subs(g) %给梯度赋初值i=1while(norm(A(i,:),4)10(-4) %收敛条件f(i)=rb(x,y) %记录函数值P(i,:)=-A(i,:) %得到迭代方向fz(i)=-A(i,:)*P(i,:) %-gT*p %精确
10、搜索法步长的分子fm(i)=P(i,:)*subs(G)*P(i,:) %精确搜索法步长的分母a(i)=fz(i)/fm(i) %精确搜索法步长X(i+1,:) = X(i,:)+a(i)*P(i,:) %产生新的点x=X(i+1,1);y=X(i+1,4)A(i+1,:)=subs(g) %产生新的梯度i=i+1end3.3 牛顿法求 Rosenbrock 函数最小值 matlab 程序如下：function rb = rbfun(x,y)rb=100*(y-x2)2+(1-x)2endclearclcsyms x y g Gg=gradient(rb(x,y),x y) %定义梯度向量-
11、18 -G=hessian(rb(x,y),x y) %定义海森阵X(1,:)=-1.4 1; %定义初值x=X(1,1);y=X(1,4);A(1,:)=subs(g) %给梯度赋初值H=subs(inv(G) %得到海森阵初值S(1,:)=-A(1,:)*H %得到 s 初值i=1while (norm(S(i,:),4)10(-6) %收敛条件f(i)=rb(x,y) %定义函数值X(i+1,:) = X(i,:)+S(i,:) %得到下一迭代点x=X(i+1,1);y=X(i+1,4) %给 x,y 分别赋值A(i+1,:)=subs(g) %得到新的梯度值H=subs(inv(G)
12、%得到新的海森阵S(i+1,:)=-A(i+1,:)*H %得到新的增量 si=i+1end3.4 共轭梯度法求解习题 5.19 程序如下：clearclcK=40G=zeros(K,K)for m = 1: Kfor n = 1: KG(m,n)=1/(m+n-1)end endX(1,:)=zeros(1,K)b=ones(1,K)A(1,:)=X(1,:)*G-b- 19 -P(1,:)=-A(1,:) i=1while (norm(A(i,:),4)10(-6) %收敛条件d=P(i,:)*Gfz(i)=A(i,:)*A(i,:) %精确搜索法步长的分子fm(i)=P(i,:)*d %精确搜索法步长的分母a(i)=fz(i)/fm(i) %精确搜索法步长X(i+1,:) = X(i,:)+a(i)*P(i,:) %产生新的点A(i+1,:)=A(i,:)+a(i)*dbeta(i+1)=(A(i+1,:)*A(i+1,:)/(A(i,:)*A(i,:)P(i+1,:)=-A(i+1,:)+beta(i+1)*P(i,:)i=i+1end
