北航最优化方法大作业参考.docx

1、北航最优化方法大作业参考1流量工程问题1.1问题重述定义一个有向网络G=(N,E)，其中N是节点集，E是弧集。令A是网络G的点弧关联矩阵，即NE阶矩阵，且第l列与弧里(I,j)对应，仅第i行元素为1，第j行元素为-1，其余元素为0。再令bm=(bm1,bmN)T，fm=(fm1,fmE)T，则可将等式约束表示成：Afm=bm本算例为一经典TE算例。算例网络有7个节点和13条弧，每条弧的容量是5个单位。此外有四个需求量均为4个单位的源一目的对，具体的源节点、目的节点信息如图所示。这里为了简单，省区了未用到的弧。此外，弧上的数字表示弧的编号。此时，c=(5,5,5)113)T，根据上述四个约束条件

2、，分别求得四个情况下的最优决策变量x=(x12,x13,x75)113)。图 1 网络拓扑和流量需求1.27节点算例求解1.2.1算例1（b1=4;-4;0;0;0;0;0T）转化为线性规划问题：Minimize cTx1Subject to Ax1=b1 x1=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;0T）Minimize cTx2Subject to Ax2=b2 X2=0利用Matlab编写对偶单纯形法程序，可求得:最优解为x2*=0 4

3、 0 0 0 0 0 0 0 0 0 0 0T对应的最优值cTx2=201.2.3算例3（b3=0;-4;4;0;0;0;0T）Minimize cTx3Subject to Ax3=b3 X3=0利用Matlab编写对偶单纯形法程序，可求得:最优解为x3*=4 0 0 0 4 0 0 0 0 0 0 0 0T对应的最优值cTx3=401.2.4算例4（b4=4;0;0;0;0;0;-4T）Minimize cTx4Subject to Ax4=b4 X4=0利用Matlab编写对偶单纯形法程序，可求得:最优解为x4*=4 0 0 4 0 0 0 0 0 4 0 0 0T对应的最优值cTx4=

4、601.3计算结果及结果说明1.3.1算例1（b1=4;-4;0;0;0;0;0T）算例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算例2（b2=4;0;-4;0;0;0;0T）算例2中，由b2可知，节点3为需求节点，节点1为供给节点，由节点1将信息传输至节点2的最短路径为弧2。图 3 算例2最优传输示意图求得的

5、最优解为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;0T）算例3中，由b3可知，节点2为需求节点，节点3为供给节点，由节点3将信息传输至节点2的最短路径为弧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。经分析，计算

6、结果合理可信。1.3.4算例4（b4=4;0;0;0;0;0;-4T）算例4中，由b4可知，节点7为需求节点，节点1为供给节点，由节点1将信息传输至节点7的最短路径为弧1-弧4-弧10。图 5 算例3最优传输示意图求得的最优解为x4*=4 0 0 4 0 0 0 0 0 4 0 0 0T，即经过弧1运输4个单位流量至节点2，再经弧4运输4个单位流量至节点5，最后经弧5运输4个单位流量至节点7，其余弧无流量。又因为，每条弧的费用均为5，所以最小费用为60。经分析，计算结果合理可信。2重要算法编写与观察2.1习题5.6（a）初值为（0,0）时本算法令g的2范数在10-4时，停止迭代，经过86次迭代

7、收敛。收敛因子（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 收敛因子截

8、图图 10 自变量（x1,x2）截图总结：最速降线法的收敛因子随着初值的不同而变化，对于个别初值（如本习题初值取（11,0）时），算法可迅速收敛。因此，初值的选取对于最速降线法的收敛速度有较大影响。2.2习题5.7（a）由可得：故，牛顿迭代法的确切公式为：（b）从以下五个初值开始迭代（1）x(0)=7.40表 1 初值1牛顿法迭代结果表迭代次数x值梯度值17.4-127.44-0.0909090937.4444-0.0009000947.4444444-9.00E-0857.4444444-9.00E-08（2）x(0)=7.20表 2 初值2牛顿法迭代结果表迭代次数x值梯度值17.2-112

9、7.31-3.637.403775-0.747.440723-7.60E-0257.444413-0.000631068（3）x(0)=7.01表 3 初值3牛顿法迭代结果表迭代次数x值梯度值17.01-39127.019775-193.275600537.03867-94.47.073976-4.51E+0157.135638-20.（4）x(0)=7.80表 4 初值4牛顿法迭代结果表迭代次数x值梯度值17.8427.16-1637.2624-6.947.369879-1.81E+0057.431934-0.3（5）x(0)=7.88表 5 初值5牛顿法迭代结果表迭代次数x值梯度值17.8

10、84.527.0176-218.272727337.034503-106.931813547.066328-5.13E+0157.122757-23.（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时，（

11、未完成）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范数在=0 flag=0; sol=zeros(1,

12、nA); for i=1:mA-1 sol(N(i)=A(i,nA); end val=sol*(B(mA,:); else for i=1:mA-1 if A(i,nA)=0 disp(have infinite solution!); flag=0; break; end end if flag temp=0; for i=1:mA-1 if A(i,nA)temp temp=A(i,nA); outb=i; end end sita=zeros(1,nA-1); for i=1:nA-1 if A(outb,i)0 sita(i)=A(mA,i)/A(outb,i); end end t

13、emp=-inf; for i=1:nA-1 if sita(i)temp temp=sita(i); inb=i; end end for i=1:mA-1 if i=outb N(i)=inb; end end A(outb,:)=A(outb,:)/A(outb,inb); for i=1:mA if i=outb A(i,:)=A(i,:)-A(outb,:)*A(i,inb); A(mA,nA)=0; end end end endend3.2最速降线法求Rosenbrock函数最小值matlab程序如下：function rb = rbfun(x,y)rb=100*(y-x2)2+

14、(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,:) %精确搜索法

15、步长的分母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) %定义梯度向量G=hessian(rb(x,y),x y) %定义海森阵X(1,:)=-1.4 1; %定义初值

16、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 %得到新的增量si=i+1end3.4共轭梯度法求解习

17、题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-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