实验2_解无约束非线性规划

实 验 报 告
解无约束非线性规划实验（运筹学与最优化方法，4学时）
一 实验目的
1．掌握迭代算法的思想。

2．掌握0.618法。

3．掌握最速下降法、共轭梯度法和拟牛顿法。

二 实验内容
1．用0.618法解下列问题：
2min ()21f x x x =--
初始区间为：11[,][1,1],0.16a b ε=-=。

2．用最速下降法、共轭梯度法和拟牛顿法分别解下列问题：
221212
min (,)4f x x x x =+ 取(0)(0)12(,)(1,1)x x =
三 实验步骤（算法）与结果
1. 解:首先建立一个黄金分割计算最优值的函数文件并保存为HJFG ，
可供调用：
x(1)=input('a=');
y(1)=input('b=');
k=input('k=');
n=1;m=1;
p(1)=x(1)+0.382*(y(1)-x(1));
q(1)=x(1)+0.618*(y(1)-x(1));
while y(n)-x(n)>=k
if g(p(n))>g(q(n))
n=n+1;
x(n)=p(n-1);y(n)=y(n-1);p(n)=q(n-1);
q(n)=x(n)+0.618*(y(n)-x(n));
elseif g(p(n))<=g(q(n))
n=n+1;
x(n)=x(n-1);y(n)=q(n-1);q(n)=p(n-1);
p(n)=x(n)+0.382*(y(n)-x(n));
end
end
x(n),y(n),min=1/2*(x(n)+y(n))
参数说明：x(1)：计算区间的左极限；
y(1)：计算区间的右极限；
k：计算要求达到的精度；
p(n)，q(n)：黄金分割的计算公式；
g(x)：输入的计算函数；
min=1/2*(x(n)+y(n))：满足条件的最优解.
实际使用：
首先建立函数文件并保存：
function G=g(x)
G=2*x^2-x-1;
然后调用上面的函数
输出结果为：
ans =0.1672
ans =0.2787
min =0.2229
也就是满足条件的解的存在区间之一是[0.1672,0.2787]，取平均值0.2229作为近似最优解.
2.解：
先建立所求目标的函数文件；
function y=fun2(a,b)
y=a.^2+4*b.^2;
以及他们的偏导数：
function dy=dx1_fun(a)
dy=2*a;
function dy=dx2_fun(b)
dy=8*b;
然后是主程序：
（1）最速下降法：
function f=zs(x1,x2,eps)
x(1,1)=x1;x(1,2)=x2;
k=1;
while 1
c(1)=dx1_fun(x(k,1));c(2)=dx2_fun(x(k,2));
if norm(c,2)>=eps
d=-c;
a=sym('a');
g=diff(fun2(x(k,1)+a*d(1),x(k,2)+a*d(2)),a); a=eval(solve(g,'a'));
x(k+1,:)=x(k,:)+a*d;
k=k+1;
else
break
end
end
x
运行：
zs(1,1,0.01)
结果是：
x =
1.0000 1.0000
0.7385 -0.0462
0.1108 0.1108
0.0818 -0.0051
0.0123 0.0123
0.0091 -0.0006
0.0014 0.0014
0.0010 -0.0001
（2）共轭梯度法：
function f2=GE(x1,x2,eps)
x(1,1)=x1;x(1,2)=x2;
k=1;
d=0;
q=[1,0;0,4];
c(1)=dx1_fun(x(k,1));c(2)=dx2_fun(x(k,2));
if norm(c,2)>=eps
d=-c;
a=sym('a');
f=diff(fun2(x(k,1)+a*d(1),x(k,2)+a*d(2)),a); a=eval(solve(f,'a'));
x(k+1,:)=x(k,:)+a*d;
end
while 1
k=k+1;
c(1)=dx1_fun(x(k,1));c(2)=dx2_fun(x(k,2));
if norm(c,2)>eps
r=(d*q*c')/(d*q*d');
d=-c+r*d;
a=sym('a');
f=diff(fun2(x(k,1)+a*d(1),x(k,2)+a*d(2)),a); a=eval(solve(f,'a'));
x(k+1,:)=x(k,:)+a*d;
else
break
end
end
x
运行：
GE(1,1,0.01)
结果为：
x =
1.0000 1.0000
0.7385 -0.0462
0 -0.0000
(3) 拟牛顿法：
syms x1x2a;
y1=diff(f(x1,x2),x1);
y2=diff(f(x1,x2),x2);
x1=input('x1=');x2=input('x2=');
t1(1)=x1;
t2(1)=x2;
A=[y1;y2];
F(:,1)=eval(A);
p(:,1)=-F(:,1);
Q=input('Q=');
H=eye(2);
for k=1:1000
if F(:,k)==0
break;
else
if k==1
d(:,k)=-F(:,k);
b=a*d(:,k);
alpha=solve(diff(f(t1(k)+b(1),t2(k)+b(2)),a),a); B(:,k+1)=[t1(k);t2(k)]+alpha*d(:,k);
t1(k+1)=B(1,k+1);
t2(k+1)=B(2,k+1);
x1=t1(k+1);
x2=t2(k+1);
F(:,k+1)=eval(A);
else
B(:,1)=[t1(1);t2(1)];
C=B(:,k)-B(:,k-1);
R=F(:,k)-F(:,k-1);
H=H+(C*C')/(C'*R)-(H*R*R'*H)/(R'*H*R);
d(:,k)=-H*F(:,k);
b=a*d(:,k);
alpha=solve(diff(f(t1(k)+b(1),t2(k)+b(2)),a),a); B(:,k+1)=[t1(k);t2(k)]+alpha*d(:,k);
t1(k+1)=B(1,k+1);
t2(k+1)=B(2,k+1);
x1=t1(k+1);
x2=t2(k+1);
F(:,k+1)=eval(A);
end
end
end
x1
x2
运行时首先编写函数文件NNT.m：function y=f(x1,x2)
y=x1.^2+4.*x2.^2;
在调用：
NNT
x1=1
x2=1
Q=[1 0;0 4]
运行结果为：
x1 =
x2 =
四实验收获与教师评语
收获：这次实验相对第一次难度低一点点，复习了几个迭代方法顺便认识到这几个方法的共同之处，方便了第二题的编写。

但这种方法使用范围较为有限，结果和预期差不太多，可以接受。