最优化实验报告
无约束最优化上机实验报告
数学与应用数学091
梁启凡
090080
实验一
实验目的：
编程并测试算法A.1 Gaussian Elimination with Row Partial Pivoting 和A.2 Cholesky Factorization
实验过程
一、A.1 Gaussian Elimination with Row Partial Pivoting
根据书上附录提供的伪代码，用MATLAB写出程序。用随机方阵测试，将测试结果与MATLAB自带的lu()函数产生的结果对比。
程序代码：
function [L,U,P]=GaussElimination(A)
%判断A是否奇异
if det(A)==0
error('matrix A is singular!')
stop
end
n=length(A);%求方阵A的维数
L=zeros(n);
P=eye(n);
for i=1:1:n
[m,t]=max(abs(A(i:n,i)));%找到A下三角阵中每列的最大值
j=t+i-1;%返回该最大值的行坐标
if(i~=j)
%交换A的i,j列
temp=A(i,1:n);
A(i,1:n)=A(j,1:n);
A(j,1:n)=temp;
%交换P的i,j列
temp=P(i,1:n);
P(i,1:n)=P(j,1:n);
P(j,1:n)=temp;
%交换L的i,j列
temp=L(i,1:n);
L(i,1:n)=L(j,1:n);
L(j,1:n)=temp;
end
L(i,i)=1;
for k=i+1:n
L(k,i)=A(k,i)/A(i,i);
for l=i+1:n
A(k,l)=A(k,l)-L(k,i)*A(i,l);
end
end
end
U=triu(A);%U为变化后的A的上三角部分
测试：
生成5维随机矩阵A：
>> A=rand(5)
A =
0.8147 0.0975 0.1576 0.1419 0.6557
0.9058 0.2785 0.9706 0.4218 0.0357
0.1270 0.5469 0.9572 0.9157 0.8491
0.9134 0.9575 0.4854 0.7922 0.9340
0.6324 0.9649 0.8003 0.9595 0.6787
对A进行分解：
>> [L,U,P]=GaussElimination(A)
L =
1.0000 0 0 0 0
0.8920 1.0000 0 0 0
0.1390 -0.5469 1.0000 0 0
0.9917 0.8870 0.9924 1.0000 0
0.6923 -0.3991 0.4794 0.1476 1.0000
U =
0.9134 0.9575 0.4854 0.7922 0.9340
0 -0.7565 -0.2753 -0.5648 -0.1774
0 0 0.7391 0.4967 0.6223
0 0 0 -0.3559 -1.3507
0 0 0 0 -0.1376
P =
0 0 0 1 0
1 0 0 0 0
0 0 1 0 0
0 1 0 0 0
0 0 0 0 1
与MATLAB自带的lu函数结果完全相同：
二、A.2 Cholesky Factorization
根据书上附录提供的伪代码，用MATLAB写出程序。用随机方阵测试，将测试结果与MATLAB自带的chol()函数产生的结果对比。
程序代码：
function [L]=cholesky(A)
n=length(A);
for i=1:n
L(i,i)=sqrt(A(i,i));
for j=i+1:n
L(j,i)=A(j,i)/L(i,i);
for k=i+1:j
A(j,k)=A(j,k)-L(j,i)*L(k,i);
end
end
end
测试：
由于Cholesky分解要求被分解方阵是对称正定的，这里选取5阶帕斯卡矩阵进行测试。
即A。
用MATLAB自带函数chol()
与自己编写的分解函数所得结果互为转置。
实验分析：
以上两个程序均得到了正确结果，变换不同的维数的矩阵，也皆能得到正确的结果，说明程序是正确的。
对第二个程序，由于没有对称正定检测，所以哪怕输如不满足条件的矩阵也能得出结果。所以应当在输入A时确保A是对称正定的。
无约束最优化测试问题：
, 
, 
, 
, 
实验二
实验目的：
使用最速下降法，编程并求解
。比较不同的终止参数和步长的选取。
实验过程及方案：
问题一的解决
终止参数选择：梯度的范数小于目标精确值
步长选择：固定步长 ，Armijo搜索
A：代码（固定步长）
function SDP1
syms x1 x2;
f=100*(x2-x1^2)^2+(1-x1)^2;
x0=[-1.2 1]';
eps=10e-3;
d1=jacobian(f);
d1=d1';
d2=jacobian(d1);
a=subs(d1,{x1,x2},{x0(1),x0(2)});
n=0;
while(norm(a)>eps)
p=-a;
step=0.0005;
x0=x0+step*p
a=subs(d1,{x1,x2},{x0(1),x0(2)});
n=n+1;
end
x0
n %迭代步数
结果：
x0 =
0.9889
0.9779
n =
18179
迭代步长为0.0005，目标精度值为
，得
，共迭代了18179步，非常慢。
B:代码（Armijo搜索）
定义两个函数
%目标函数
function f=fun(x)
f=100*(x(1)^2-x(2))^2+(x(1)-1)^2+(x(3)-1)^2+90*(x(3)^2-x(4))^2+10.1*((x(2)-1)^2+(x(4)-1)^2)+19.8*(x(2)-1)*(x(4)-1);
%目标函数的梯度，用jacobian（）函数求得
function gf=gfun(x)
gf=[ 2*x(1) - 400*x(1)*(- x(1)^2 + x(2)) - 2, - 200*x(1)^2 + (1101*x(2))/5 + (99*x(4))/5 - 40, 2*x(3) - 360*x(3)*(- x(3)^2 + x(4)) - 2, - 180*x(3)^2 + (99*x(2))/5 + (1001*x(4))/5 - 40]';
function SDP1(fun,gfun)
syms x
x0=[-1.2 1]';
maxk=5000000;%最大迭代步数
rho=0.5;
sigma=0.4;
k=0;
eps=10e-3;
while(k<maxk)
g=feval(gfun,x0);计算梯度值
d=-g;
if(norm(d)<eps),break;end
m=0;mk=0;
while(m<20)%Armoji搜索
if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)
mk=m;break;
end
m=m+1;
end
x0=x0+rho^mk*d;
k=k+1;
end
x0 %最优点
k %迭代次数
结果：
x0 =
0.9895
0.9791
k =
1202
达到同样的精度速度是很快的。
将精度改为10e-5
x0 =
1.0000
1.0000
k =
1435
精度相差10e2，搜索次数仅仅增加了200余次，效果拔群。
问题二的解决
终止参数选择：梯度的范数小于目标精确值
步长选择：Armijo搜索，wolfe线搜索条件
A：代码（Armijo搜索）
定义两个函数
%目标函数
function f=fun(x)
f=100*(x(1)^2-x(2))^2+(x(1)-1)^2+(x(3)-1)^2+90*(x(3)^2-x(4))^2+10.1*((x(2)-1)^2+(x(4)-1)^2)+19.8*(x(2)-1)*(x(4)-1);
%目标函数的梯度，用jacobian（）函数求得
function gf=gfun(x)
gf=[ 2*x(1) - 400*x(1)*(- x(1)^2 + x(2)) - 2, - 200*x(1)^2 + (1101*x(2))/5 + (99*x(4))/5 - 40, 2*x(3) - 360*x(3)*(- x(3)^2 + x(4)) - 2, - 180*x(3)^2 + (99*x(2))/5 + (1001*x(4))/5 - 40]';
function SDP2(fun,gfun)
syms x
x0=[-3 -1 -3 -1]';
maxk=5000000;%最大迭代步数
rho=0.5;
sigma=0.4;
k=0;
eps=10e-5;
while(k<maxk)
g=feval(gfun,x0);%计算梯度值
d=-g;
if(norm(d)<eps),break;end
m=0;mk=0;
while(m<20)%Armoji搜索
if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)
mk=m;break;
end
m=m+1;
end
x0=x0+rho^mk*d;
k=k+1;
end
x0 %最优点
k %迭代次数
结果：
x0 =
1.0000
1.0001
1.0000
0.9999
k =
7156
非常精确。
B:代码（wolfe线搜索条件）fun函数和gfun函数不变
function SDP2(fun,gfun)
syms x
x0=[-3 -1 -3 -1]';
maxk=5000000;%最大迭代步数
c1=10e-4;
c2=0.9;
k=0;
eps=10e-5;
while(k<maxk)
g=feval(gfun,x0);%计算梯度值
d=-g;
if(norm(d)<eps),break;end
m=0;mk=0;
step=1;
while(step>0.0005)%wolfe搜索÷
if(feval(fun,x0+step*d)<feval(fun,x0)+c1*step*feval(gfun,x0)'*d && feval(gfun,x0)'*d>c2*feval(gfun,x0)'*d)
step=step;break;
end
step=step*0.8;
end
x0=x0+step*d
k=k+1
end
x0 %最优点
k %迭代次数
结果：
x0 =
1.0000
1.0001
1.0000
0.9999
k =
21471
相同精度要求下，没有Armoji搜索快。
问题三的解决
终止参数选择：梯度的范数小于目标精确值
步长选择：Armijo搜索，wolfe线搜索条件
A：代码（Armijo搜索）
定义两个函数
%目标函数
function f=fun(x)
f=(x(1)+10*x(2))^2+5*(x(3)-x(4))^2+(x(2)-2*x(3))^4+10*(x(1)-x(4))^4;
%目标函数的梯度，用jacobian（）函数求得
function gf=gfun(x)
gf=[ 2*x(1) + 20*x(2) + 40*(x(1) - x(4))^3, 20*x(1) + 200*x(2) + 4*(x(2) - 2*x(3))^3, 10*x(3) - 10*x(4) - 8*(x(2) - 2*x(3))^3, 10*x(4) - 10*x(3) - 40*(x(1) - x(4))^3]’;
function SDP3(fun,gfun)
syms x
x0=[3 -1 0 1]';
maxk=5000000;%最大迭代步数
rho=0.5;
sigma=0.4;
k=0;
eps=10e-5;
while(k<maxk)
g=feval(gfun,x0);计算梯度值
d=-g;
if(norm(d)<eps),break;end
m=0;mk=0;
while(m<20)%Armoji搜索
if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)
mk=m;break;
end
m=m+1;
end
x0=x0+rho^mk*d;
k=k+1;
end
x0 %最优点
k %迭代次数
结果：
x0 =
0.0237
-0.0024
0.0118
0.0118
k =
13624
将精度eps改为10e-6
x0 =
0.0110
-0.0011
0.0055
0.0055
k =
62671
越来越接近实际值，但是搜索速度非常慢。
B:代码（wolfe线搜索条件）fun函数和gfun函数不变
代码同问题二的B方案一致。将代码中初始迭代点x0改为x0=[3 -1 0 1]'
结果：
速度非常慢，在进行到105416步的时候停止了计算，得到
x0 =
0.0468
-0.0047
0.0233
0.0234 这个结果离精确值差距还很大。
实验分析：
问题一用了2种步长，如果固定步长会使得搜索过程很慢，而逐一选择步长可以加快搜索过程和精度。问题二和问题三分别用Armoji搜索和wolfe搜索算步长的方法，Armoji搜索的速度要比wolfe快很多
实验三
实验目的：
使用牛顿法，编程并求解。
实验过程及方案：
问题一的解决
代码：
function NTP1
syms x1 x2;
f=100*(x2-x1^2)^2+(1-x1)^2;
x0=[-1.2 1]';
eps=10e-3;
d1=jacobian(f);
d2=jacobian(d1);
d1=d1';
a=subs(d1,{x1,x2},{x0(1),x0(2)});
b=subs(d2,{x1,x2},{x0(1),x0(2)});
n=0;
while(norm(a)>eps)
x0=x0-inv(b)*a
a=subs(d1,{x1,x2},{x0(1),x0(2)});
b=subs(d2,{x1,x2},{x0(1),x0(2)});
n=n+1;
end
x0
n
结果：
x0 =
1.0000
1.0000
n =
5
非常快速，只用5步就得到精度为10e-3的解。
问题二的解决
代码：
function nt
syms x1 x2 x3 x4;
f=100*(x1^2-x2)^2+(x1-1)^2+(x3-1)^2+90*(x3^2-x4)^2+10.1*((x2-1)^2+(x4-1)^2)+19.8*(x2-1)*(x4-1);
x0=[-3 -1 -3 -1]';
eps=10e-5;
d1=jacobian(f);
d2=jacobian(d1);
d1=d1';
a=subs(d1,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)});
b=subs(d2,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)});
n=0;
while(norm(a)>eps)
x0=x0-inv(b)*a
a=subs(d1,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)});
b=subs(d2,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)});
n=n+1;
end
x0
n
结果：
x0 =
-0.9680
0.9471
-0.9695
0.9512
n =
13
得到精度为10e-5的解只用了13步。
问题三的解决
代码：
function nt
syms x1 x2 x3 x4;
f=(x1+10*x2)^2+5*(x3-x4)^2+(x2-2*x3)^4+10*(x1-x4)^4;;
x0=[3 -1 0 1]';
eps=10e-20;
d1=jacobian(f);
d2=jacobian(d1);
d1=d1';
a=subs(d1,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)});
b=subs(d2,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)});
n=0;
while(norm(a)>eps)
x0=x0-inv(b)*a
a=subs(d1,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)});
b=subs(d2,{x1,x2,x3,x4},{x0(1),x0(2),x0(3),x0(4)});
n=n+1;
end
x0
n
结果：
x0 =
1.0e-006 *
0.1435
-0.0144
0.0229
0.0229
n =
41
经过41步，得到精确度为10e-20的解。
实验分析：
牛顿法的速度较最速下降法有非常明显的改善，速度很快
实验四
实验目的：
调用MATLAB自带函数求解。
实验过程及方案：
分别用fminsearch函数和fminunc函数测试，对比两个函数的结果。
P1
1、使用fminsearch函数
>> f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2','x');
>> x0=[-1.2; 1]; ff=optimset; ff.Display='iter';
>> x=fminsearch(f,x0,ff)
Iteration Func-count min f(x) Procedure
0 1 24.2
1 3 20.05 initial simplex
2 5 5.1618 expand
3 7 4.4978 reflect
4 9 4.4978 contract outside
………………………………………………………………………………
81 151 1.12293e-009 contract inside
82 153 1.12293e-009 contract outside
83 155 1.12293e-009 contract inside
84 157 1.10755e-009 contract outside
85 159 8.17766e-010 contract inside
Optimization terminated:
the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004
and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-004
x =
1.0000
1.0000
总共用了85步得到结果。
2、使用fminunc函数
在刚才基础上输入
>>x=fminunc(f,[-1.2;.1],ff)
First-order
Iteration Func-count f(x) Step-size optimality
0 3 24.2 216
1 9 4.28049 0.000861873 15.2
2 12 4.12869 1 3
3 15 4.12069 1 1.5
……………………………………………………………………………………
31 123 0.00108975 0.5 0.877
32 126 9.14432e-005 1 0.143
33 129 7.27182e-006 1 0.0807
34 132 4.44847e-007 1 0.022
35 135 1.47502e-008 1 0.00272
36 138 2.83357e-011 1 1.91e-005
Local minimum found.
Optimization completed because the size of the gradient is less than
the default value of the function tolerance.
<stopping criteria details>
x =
1.0000
1.0000
只用了36步便得到正确结果，比fminsearch函数快很多。
P2
1、使用fminsearch函数
>>f=inline('100*(x(1)^2-x(2))^2+(x(1)-1)^2+(x(3)-1)^2+90*(x(3)^2-x(4))^2+10.1*((x(2)-1)^2+(x(4)-1)^2)+19.8*(x(2)-1)*(x(4)-1)','x');
>> x0=[-3;-1;-3;-1]; ff=optimset; ff.Display='iter';
>> x=fminsearch(f,x0,ff)
314 527 1.94483e-009 contract inside
用了314步得到结果
x =
1.0000
1.0000
1.0000
1.0000
2、使用fminunc函数
>>x=fminunc(f,[-3;-1;-3;-1],ff)
35 215 1.84511e-007 1 0.00339
用35步得到结果
x =
1.0002
1.0004
0.9998
0.9996
结果没有fminsearch函数精确。
P3
1、使用fminsearch函数
>>f=inline('(x(1)+10*x(2))^2+5*(x(3)-x(4))^2+(x(2)-2*x(3))^4+10*(x(1)-x(4))^4','x');
>> x0=[3;-1;0;1]; ff=optimset; ff.Display='iter';
>> x=fminsearch(f,x0,ff)
185 305 1.39059e-006 reflect
用185步得到结果
x =
0.0094
-0.0009
0.0166
0.0166
与真实值
还有一些差距。
norm(x-[0 0 0 0]') =0.0253
2、使用fminunc函数
>>x=fminunc(f,[3;-1;0; 1],ff)
33 170 6.38904e-009 1 0.000217
用33步得到结果
x =
0.0015
-0.0002
-0.0031
-0.0031
norm(x-[0 0 0 0]')=0.0047 比fminsearch函数得到的结果与精确值的差的范数要小。
实验分析：
fminunc函数比fminisearch函数速度更快，快很多，不过精确性上是不分伯仲的。fminisearch与fminunc这两个函数都是最优化工具箱中的，如果需要速度快些，就用fminiunc函数。
- 产品质量检验报告
-
解读《血液化验报告》
解读血液化验报告AKPALPALBALBGAMS淀粉酶AMS淀粉酶APOAAPOB碱性42磷酸12酶8白蛋35白5515白2球比5…
-
如何看化验报告单
如何看化验报告单检验报告单是检查所得出的客观数据记录检验项目很多这里只从定性和定量两个方面概略地作一介绍定性检验是看送检的标本中有…
- 化 验 报 告 单
-
化验报告
徐州市泉山区泉山社区卫生服务中心检验报告单样本号3姓名李曾慧病人类型床号标本类型血清性别女病人编号病区备注年龄28岁科室内科送检医…
-
RF优化总结报告
富川县城现场RF优化报告3月10日3月11日对富川县城进行了现场RF优化并对县城区主要道路进行了DT复测重点对基站远点处信号质量覆…
-
网络优化开题报告
毕业设计论文开题报告系别电子信息工程系专业通信工程班级学生姓名学号指导教师报告日期3月22日毕业设计论文开题报告表注内容用小四宋体
-
网站优化分析报告
旅途旅游网优化分析报告报告分析人邓然分析网站旅途旅游网网站URL域名年龄3年9个月创建于20xx年10月26日优化关键词无锡旅行社…
-
网站优化报告
网站优化分析报告一整体测试与分析利用网络工具分析网站域名得到结果上图为综合查询分析结果以下我会一一进行分析1百度排名在这里可以清楚…
-
搜索引擎优化实验报告
华北科技学院管理系实验报告册实验报告实验时间20xx年12月20日关键词新闻凤凰网凤凰网是一个集图文资讯视频点播专题报道虚拟社区免…