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

北航最优化方法大作业参考-标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII
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×13
)T，
根据上述四个约束条件，分别求得四个情况下的最优决策变量x=((x
12,x
13
,…,x
75
)
1×
13
)。

图 1 网络拓扑和流量需求
1.27节点算例求解
1.2.1算例1（b1=[4;-4;0;0;0;0;0]T）
转化为线性规划问题：
Minimize c T x1
Subject to Ax1=b1
x1>=0利用Matlab编写对偶单纯形法程序，可求得:
最优解为x1*=[4 0 0 0 0 0 0 0 0 0 0 0 0]T
对应的最优值c T x1=20
1.2.2算例2（b2=[4;0;-4;0;0;0;0]T）
Minimize c T x2
Subject to Ax2=b2
X2>=0 利用Matlab编写对偶单纯形法程序，可求得:
最优解为x2*=[0 4 0 0 0 0 0 0 0 0 0 0 0]T
对应的最优值c T x2=20
1.2.3算例3（b3=[0;-4;4;0;0;0;0]T）
Minimize c T x3
Subject to Ax3=b3
X3>=0 利用Matlab编写对偶单纯形法程序，可求得:
最优解为x3*=[4 0 0 0 4 0 0 0 0 0 0 0 0]T
对应的最优值c T x3=40
1.2.4算例4（b4=[4;0;0;0;0;0;-4]T）
Minimize c T x4
Subject to Ax4=b4
X4>=0
利用Matlab编写对偶单纯形法程序，可求得:
最优解为x4*=[4 0 0 4 0 0 0 0 0 4 0 0 0]T
对应的最优值c T x4=60
1.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'()97
f x x =-
- 2
4
"()(7)f x x =
-
故，牛顿迭代法的确切公式为：
2
4974(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.19
N=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 flag
kk=kk+1;
if A(:,nA)>=0
flag=0;
sol=zeros(1,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&A(i,1:nA-1)>=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
temp=-inf;
for i=1:nA-1
if sita(i)<0&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
end
end
3.2最速降线法求Rosenbrock函数最小值matlab程序如下：
function rb = rbfun(x,y)
rb=100*(y-x^2)^2+(1-x)^2
end
clear
clc
syms x y g G
g=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=1
while(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+1
end
3.3牛顿法求Rosenbrock函数最小值matlab程序如下：
function rb = rbfun(x,y)
rb=100*(y-x^2)^2+(1-x)^2
end
clear
clc
syms x y g G
g=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=1
while (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+1
end
3.4共轭梯度法求解习题5.19程序如下：
clear
clc
K=40
G=zeros(K,K)
for m = 1: K
for n = 1: K
G(m,n)=1/(m+n-1)
end
end
X(1,:)=zeros(1,K)
b=ones(1,K)
A(1,:)=X(1,:)*G-b
P(1,:)=-A(1,:)
i=1
while (norm(A(i,:),4)>10^(-6)) %收敛条件
d=P(i,:)*G
fz(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)*d
beta(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+1
end。