实验2 最速下降法和共轭梯度法的程序设计
一、实验目的
1、熟悉无约束优化问题的最速下降算法和共轭梯度法。
2、培养matlab 编程与上机调试能力。
二、实验课时：2个课时
三、实验准备
1、预习无约束优化问题的最速下降算法和共轭梯度法。
2、熟悉matlab 软件的基本操作及程序编写。
四、实验内容
课堂实验演示
根据最速下降法编写程序，求函数
21222121342)(min x x x x x x x f -++-=
的极小值，其中初始点为()01,1T x =
算法步骤如下：
Step1:：给出初始点0x ，和精度1;0k ε<<=；
Step2：计算()k f x ∇，如果()k f x ε∇≤，则停止迭代，输出结果；否则转step3； Step3：令下降方向()k k d f x =-∇，计算步长因子k λ使得0()min ()k k k k k f x d f x d λλλ≥+=+，令1,1k k k k x x d k k λ+=+=+，转step2。
其程序如下：
function [x,iter,val,dval] = Steepest_Descent_Method(x,eps) k = 1;
dy = grad_obj(x);
x_mat(:,1) = x;%存储每一次迭代得到的点x
while norm(dy)>eps
d = -dy; % 搜索方向
lambda = line_search(x,d);%步长
x = x + d*lambda;
k = k + 1;
x_mat(:,k) = x;
dy = grad_obj(x);
end
iter = k - 1;
val = obj(x);%目标函数在极值点处的函数值
dval = grad_obj(x);%目标函数在极值点处的梯度
%-------------------------------------------------------------- x1 = linspace(-1.2,1.2,40);
x2 = linspace(-0.2,1.2,40);
[xx,yy] = meshgrid(x1,x2);
for i = 1:length(x1)
for j = 1:length(x2)
z(i,j) = obj([xx(i,j);yy(i,j)]);
end
end
contour(xx,yy,z,10);%画出目标函数的等高线
hold on
plot(x_mat(1,:),x_mat(2,:),'-o')%画出最速下降法的迭代路径
hold off
function y = obj(x)
%目标函数
y = x(1).^2 - 2*x(1).*x(2) + 4*x(2).^2 + x(1) - 3*x(2);
function dy = grad_obj(x)
%目标函数的梯度
dy = [2*x(1) - 2*x(2) + 1; -2*x(1) + 8*x(2) - 3];
function lambda = line_search(xk,dk)
%作线搜索，求步长
%phi(lambda) = obj( xk + lambda*dk )
%d_phi(lambda) = dk'*grad_obj( xk + lambda*dk )
syms eqn lambda
eqn = dk'*grad_obj(xk+lambda*dk);
lambda = solve(eqn); %用符号计算命令solve 求方程d_phi(\lambda)=0的根 lambda = eval(lambda);%将符号计算的结果转化为数值类型
>> x = [1;1]; eps = 1.0e-6;
>> [x,iter,val,dval] = Steepest_Descent_Method(x,eps) x =[ -0.1667 0.3333]’
val = -0.5833
dval = [0.5280*1.0e-006 -0.1760*1.0e-006]’
iter = 43
共轭梯度法的计算步骤：
Step1: 给出初始点0x ，令)(00x f d -∇=，精度1;0k ε<<=；
Step2：计算()k f x ∇，如果()k f x ε∇≤，则停止迭代，输出结果；否则转step3； Step3：计算k k k k d x x λ+=+1，其中步长因子k λ使得0
()min ()k k k k k f x d f x d λλλ≥+=+，计算下降方向k k k k k d x f x f x f d 22111||
)(||||)(||)(∇∇+-∇=+++； 令1+=k k ，转step2。
其程序如下：
function [x,iter,val,dval] = Conjugate_Gradient_Method(x,eps)
k = 1;
x_mat(:,1) = x;%存储每一次迭代得到的点x
x_old = x;
dy_old = grad_obj(x_old);
d_old = -dy_old;
while norm(dy_old)>eps
lambda = line_search(x_old,d_old);%步长
x_new = x_old + lambda*d_old;
dy_new = grad_obj(x_new);
coef = norm(dy_new)/norm(dy_old);
d_new = -dy_new + coef^2*d_old; % 搜索方向
k = k + 1;
x_old = x_new;
dy_old = dy_new;
d_old = d_new;
x_mat(:,k) = x_new;
%防止死循环
if k > 100
break;
end
end
x = x_new;
iter = k - 1;
val = obj(x_new);%目标函数在极值点处的函数值
dval = grad_obj(x_new);%目标函数在极值点处的梯度
%-------------------------------------------------------------- x1 = linspace(-1.2,1.2,40);
x2 = linspace(-0.2,1.2,40);
[xx,yy] = meshgrid(x1,x2);
for i = 1:length(x1)
for j = 1:length(x2)
z(i,j) = obj([xx(i,j);yy(i,j)]);
end
end
contour(xx,yy,z,10);%画出目标函数的等高线
hold on
plot(x_mat(1,:),x_mat(2,:),'-o')%画出最速下降法的迭代路径
hold off
function y = obj(x)
%目标函数
y = x(1)^2 - 2*x(1)*x(2) + 4*x(2)^2 + x(1) - 3*x(2);
function dy = grad_obj(x)
