北航最优化方法大作业参考1 流量工程问题1.1 问题重述定义一个有向网络G=(N,E)，其中N是节点集，E是弧集。

令A是网络G的点弧关联矩阵，即N×E阶矩阵，且第l列与弧里(I,j)对应，仅第i行元素为1，第j行元素为-1，其余元素为0。

再令b m=(b m1,…,b mN)T，f m=(f m1,…,f mE)T，则可将等式约束表示成：Af m=b m本算例为一经典TE算例。

算例网络有7个节点和13条弧，每条弧的容量是5个单位。

此外有四个需求量均为4个单位的源一目的对，具体的源节点、目的节点信息如图所示。

这里为了简单，省区了未用到的弧。

此外，弧上的数字表示弧的编号。

此时，c=((5,5…,5)1 )T，×13根据上述四个约束条件，分别求得四个情况下的最优决策变量x=((x12,x13,…,x75)1×)。

13图 1 网络拓扑和流量需求1.2 7节点算例求解1.2.1 算例1（b1=[4;-4;0;0;0;0;0]T）转化为线性规划问题：Minimize c T x1Subject to Ax1=b1x1>=0 利用Matlab编写对偶单纯形法程序，可求得:最优解为x1*=[4 0 0 0 0 0 0 0 0 0 0 0 0]T对应的最优值c T x1=201.2.2 算例2（b2=[4;0;-4;0;0;0;0]T）Minimize c T x2Subject to Ax2=b2X2>=0 利用Matlab编写对偶单纯形法程序，可求得:最优解为x2*=[0 4 0 0 0 0 0 0 0 0 0 0 0]T对应的最优值c T x2=201.2.3 算例3（b3=[0;-4;4;0;0;0;0]T）Minimize c T x3Subject to Ax3=b3X3>=0 利用Matlab编写对偶单纯形法程序，可求得:最优解为x3*=[4 0 0 0 4 0 0 0 0 0 0 0 0]T对应的最优值c T x3=401.2.4 算例4（b4=[4;0;0;0;0;0;-4]T）Minimize c T x4Subject to Ax4=b4X4>=0利用Matlab编写对偶单纯形法程序，可求得:最优解为x4*=[4 0 0 4 0 0 0 0 0 4 0 0 0]T对应的最优值c T x4=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 0]T，即只经过弧1运输4个单位流量，其余弧无流量。

又因为，每条弧的费用均为5，所以最小费用为20。

经分析，计算结果合理可信。

1.3.2 算例2（b2=[4;0;-4;0;0;0;0]T）算例2中，由b2可知，节点3为需求节点，节点1为供给节点，由节点1将信息传输至节点2的最短路径为弧2。

图 3 算例2最优传输示意图求得的最优解为x2*=[0 4 0 0 0 0 0 0 0 0 0 0 0]T，即只经过弧2运输4个单位流量，其余弧无流量。

又因为，每条弧的费用均为5，所以最小费用为20。

经分析，计算结果合理可信。

1.3.3 算例3（b3=[0;-4;4;0;0;0;0]T）算例3中，由b3可知，节点2为需求节点，节点3为供给节点，由节点3将信息传输至节点2的最短路径为弧5->弧1。

图 4 算例3最优传输示意图求得的最优解为x3*=[4 0 0 0 4 0 0 0 0 0 0 0 0]T，即经过弧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 算例3最优传输示意图求得的最优解为x4*=[4 0 0 4 0 0 0 0 0 4 0 0 0]T，即经过弧1运输4个单位流量至节点2，再经弧4运输4个单位流量至节点5，最后经弧5运输4个单位流量至节点7，其余弧无流量。

又因为，每条弧的费用均为5，所以最小费用为60。

经分析，计算结果合理可信。

2 重要算法编写与观察2.1 习题5.6（a）初值为（0,0）时本算法令g的2范数在<10-4时，停止迭代，经过86次迭代收敛。

收敛因子（f(k+1)-f*）/（f(k)-f*）=0.7623图 6 收敛因子截图（b）初值为（-0.4,0）时本算法令g的2范数在<10-4时，停止迭代，经过112次迭代收敛。

收敛因子（f(k+1)-f*）/（f(k)-f*）=0.81图7 收敛因子截图（c）初值为（10,0）时本算法令g的2范数在<10-4时，停止迭代，经过5次迭代收敛。

收敛因子（f(k+1)-f*）/（f(k)-f*）=3.9022e-4图8 收敛因子截图（d）初值为（11,0）时本算法令g的2范数在<10-4时，停止迭代，经过2次迭代收敛。

收敛因子（f(k+1)-f*）/（f(k)-f*）= 0图9 收敛因子截图图10 自变量（x1,x2）截图总结：最速降线法的收敛因子随着初值的不同而变化，对于个别初值（如本习题初值取（11,0）时），算法可迅速收敛。

因此，初值的选取对于最速降线法的收敛速度有较大影响。

2.2 习题5.7（a ） 由()94ln(7)f x x x =--可得：4'()97f x x =-- 24"()(7)f x x =- 故，牛顿迭代法的确切公式为：24974(7)x s x --=-- （b ）从以下五个初值开始迭代 （1）x(0)=7.40（2）x(0)=7.20（3）x(0)=7.01（4）x(0)=7.80（5）x(0)=7.88（c）本问题的最优值为7.4444444。

由上述五个初值点的前五步迭代可以看出：当初值点在区间（7.4444444，7.8888）内时，第二次迭代点将落在（7，7.4444444）之间，随后逐渐增加，直至逼近最优值。

当初值点在区间（7，7.4444444）内时，则迭代点逐渐增加，逼近最优值。

当取初值不在（7,7.8888）内时，牛顿法不收敛。

2.3 习题5.8（a）没有线搜索的牛顿法μ=0.1时，μ=1时，（b）具有线搜索的牛顿法μ=0.1时，μ=1时，（未完成）2.4 习题5.9（a）初值选（1.2,1.2）时，◆最速降线法：本算法令g的2范数在<10-2时，停止迭代，经过3262次迭代得到以下结果。

图11 最速降线法初值为（1.2,1.2）的等值线图及迭代轨迹◆牛顿法：本算法令s的4范数在<10-6时，停止迭代，经过4次迭代得到以下结果。

图12 牛顿法初值为（1.2,1.2）的等值线图及迭代轨迹（b）初值选（-1.2,1）时，最速降线法：本算法令g的4范数在<10-2时，停止迭代，经过6835次迭代得到以下结果。

图13 最速降线法初值为（-1.2,1）的等值线图及迭代轨迹牛顿法：本算法令s的4范数在<10-6时，停止迭代，经过6次迭代得到以下结果。

图14 牛顿法初值为（-1.2,1）的等值线图及迭代轨迹2.5 习题5.19N=5迭代6次后，满足收敛条件。

N=8迭代19次后，满足收敛条件。

N=14迭代49次后，满足收敛条件。

（表略）N=40迭代74次后，满足收敛条件。

（表略）2.6 习题5.27调用MATLAB自带的lsqnonlin.m函数，计算可得对应的x(1)、x(2)和标准差如下表所示。

由上可知，标准差值较为恒定，随初值变化不十分显著；x1和x2值随初值选取的不同而不同。

2.7 习题6.4（未完成）3 附录3.1 对偶单纯形法函数MATLAB程序function [sol,val,kk]=duioudanchun(A,N)B=A;[mA,nA]=size(A);kk=0;flag=1;while flagkk=kk+1;if A(:,nA)>=0flag=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&A(i,1:nA-1)>=0disp('have infinite solution!');flag=0;break;endendif flagtemp=0;for i=1:mA-1if A(i,nA)<temptemp=A(i,nA);outb=i;endendsita=zeros(1,nA-1);for i=1:nA-1if A(outb,i)<0sita(i)=A(mA,i)/A(outb,i);endendtemp=-inf;for i=1:nA-1if sita(i)<0&sita(i)>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~=outbA(i,:)=A(i,:)-A(outb,:)*A(i,inb);A(mA,nA)=0;endendendendend3.2 最速降线法求Rosenbrock函数最小值matlab程序如下：function rb = rbfun(x,y)rb=100*(y-x^2)^2+(1-x)^2endclearclcsyms 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 %精确搜索法步长的分子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-x^2)^2+(1-x)^2endclearclcsyms 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) %给梯度赋初值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)) %得到新的海森阵S(i+1,:)=-A(i+1,:)*H %得到新的增量s i=i+1end3.4 共轭梯度法求解习题5.19程序如下：clearclcK=40G=zeros(K,K)for m = 1: Kfor n = 1: KG(m,n)=1/(m+n-1)endendX(1,:)=zeros(1,K)b=ones(1,K)A(1,:)=X(1,:)*G-bP(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。