大连理工大学优化方法上机大作业

学院：

专业：

班级：

学号：

姓名：

上机大作业1：

1.最速下降法：

function f = fun(x)

f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; end

function g = grad(x)

g = zeros(2,1);

g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);

end

function x_star = steepest(x0,eps) gk = grad(x0);

res = norm(gk);

k = 0;

while res > eps && k<=1000

dk = -gk;

ak =1; f0 = fun(x0);

f1 = fun(x0+ak*dk);

slope = dot(gk,dk);

while f1 > f0 + *ak*slope

ak = ak/4;

xk = x0 + ak*dk;

f1 = fun(xk);

end

k = k+1;

x0 = xk;

gk = grad(xk);

res = norm(gk);

fprintf('--The %d-th iter, the residual is %f\n',k,res); end

x_star = xk;

end

>> clear

>> x0=[0,0]';

>> eps=1e-4;

>> x=steepest(x0,eps)

2.牛顿法：

function f = fun(x)

f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; end

function g = grad2(x)

g = zeros(2,2);

g(1,1)=2+400*(3*x(1)^2-x(2));

g(1,2)=-400*x(1);

g(2,1)=-400*x(1);

g(2,2)=200;

end

function g = grad(x)

g = zeros(2,1);

g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);

end

function x_star = newton(x0,eps)

gk = grad(x0);

bk = [grad2(x0)]^(-1);

res = norm(gk);

k = 0;

while res > eps && k<=1000

dk=-bk*gk;

xk=x0+dk;

k = k+1;

x0 = xk;

gk = grad(xk);

bk = [grad2(xk)]^(-1);

res = norm(gk);

fprintf('--The %d-th iter, the residual is %f\n',k,res); end

x_star = xk;

end

>> clear

>> x0=[0,0]';

>> eps=1e-4;

>> x1=newton(x0,eps)

--The 1-th iter, the residual is

--The 2-th iter, the residual is

x1 =

法：

function f = fun(x)

f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; end

function g = grad(x)

g = zeros(2,1);

g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);

end

function x_star = bfgs(x0,eps) g0 = grad(x0);

gk=g0;

res = norm(gk);

Hk=eye(2);

k = 0;

while res > eps && k<=1000

dk = -Hk*gk;

ak =1; f0 = fun(x0);

f1 = fun(x0+ak*dk);

slope = dot(gk,dk);

while f1 > f0 + *ak*slope

ak = ak/4;

xk = x0 + ak*dk;

f1 = fun(xk);

end

k = k+1;

fa0=xk-x0;

x0 = xk;

go=gk;

gk = grad(xk);

y0=gk-g0;

Hk=((eye(2)-fa0*(y0)')/((fa0)'*(y0)))*((eye(2)-(y0)*(fa0)')/((fa0)'*(y0)))+(fa0*(fa 0)')/((fa0)'*(y0));

res = norm(gk);

fprintf('--The %d-th iter, the residual is %f\n',k,res);

end

x_star = xk;

End

>> clear

>> x0=[0,0]';

>> eps=1e-4;

>> x=bfgs(x0,eps)

4.共轭梯度法：

function f = fun(x)

f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; end

function g = grad(x)

g = zeros(2,1);

g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);

end

function x_star =CG(x0,eps) gk = grad(x0);

res = norm(gk);

k = 0;

dk = -gk;

while res > eps && k<=1000 ak =1; f0 = fun(x0);

f1 = fun(x0+ak*dk);

slope = dot(gk,dk);

while f1 > f0 + *ak*slope ak = ak/4;

xk = x0 + ak*dk;

f1 = fun(xk);

end

k = k+1;

x0 = xk;

g0=gk;

gk = grad(xk);

res = norm(gk);

p=(gk/g0)^2;

dk1=dk;

dk=-gk+p*dk1;

fprintf('--The %d-th iter, the residual is %f\n',k,res); end

x_star = xk;

end

>> clear

>> x0=[0,0]';

>> eps=1e-4;

>> x=CG(x0,eps)

上机大作业2：

function f= obj(x)

f=4*x(1)-x(2)^2-12;

end

function [h,g] =constrains(x) h=x(1)^2+x(2)^2-25;

g=zeros(3,1);

g(1)=-10*x(1)+x(1)^2-10*x(2)+x(2)^2+34;

g(2)=-x(1);

g(3)=-x(2);

end

function f=alobj(x) %拉格朗日增广函数

%N_equ等式约束个数?

%N_inequ不等式约束个数

N_equ=1;

N_inequ=3;

global r_al pena;%全局变量

h_equ=0;

h_inequ=0;

[h,g]=constrains(x);

%等式约束部分?

for i=1:N_equ

h_equ=h_equ+h(i)*r_al(i)+(pena/2)*h(i).^2;

end

%不等式约束部分

for i=1:N_inequ

h_inequ=h_inequ+pena)*(max(0,(r_al(i)+pena*g(i))).^2-r_al(i).^2); end

%拉格朗日增广函数值

f=obj(x)+h_equ+h_inequ;

function f=compare(x)

global r_al pena N_equ N_inequ;

N_equ=1;

N_inequ=3;

h_inequ=zeros(3,1);

[h,g]=constrains(x);

%等式部分

for i=1:1

h_equ=abs(h(i));

end

%不等式部分

for i=1:3

h_inequ=abs(max(g(i),-r_al(i+1)/pena));

end

h1 = max(h_inequ);

f= max(abs(h_equ),h1); %sqrt(h_equ+h_inequ);

function [ x,fmin,k] =almain(x_al)

%本程序为拉格朗日乘子算法示例算法%函数输入：

% x_al:初始迭代点

% r_al：初始拉格朗日乘子N-equ:等式约束个数N_inequ:不等式约束个数?%函数输出% X：最优函数点FVAL:最优函数值

%============================程序开始================================ global r_al pena ; %参数（全局变量）

pena=10; %惩罚系数

r_al=[1,1,1,1];

c_scale=2; %乘法系数乘数

cta=; %下降标准系数

e_al=1e-4; %误差控制范围

max_itera=25;

out_itera=1; %迭代次数

%===========================算法迭代开始============================= while out_itera

x_al0=x_al;

r_al0=r_al;

%判断函数?

compareFlag=compare(x_al0);

%无约束的拟牛顿法BFGS

[X,fmin]=fminunc(@alobj,x_al0);

x_al=X; %得到新迭代点

%判断停止条件?

if compare(x_al)

disp('we get the opt point');

break

end

%c判断函数下降度?

if compare(x_al)

pena=1*pena; %可以根据需要修改惩罚系数变量

else

pena=min(1000,c_scale*pena); %%乘法系数最大1000

disp('pena=2*pena');

end

%%?更新拉格朗日乘子

[h,g]=constrains(x_al);

for i=1:1

%%等式约束部分

r_al(i)= r_al0(i)+pena*h(i);

end

for i=1:3

%%不等式约束部分

r_al(i+1)=max(0,(r_al0(i+1)+pena*g(i)));

end

out_itera=out_itera+1;

end

%+++++++++++++++++++++++++++迭代结束+++++++++++++++++++++++++++++++++ disp('the iteration number');

k=out_itera;

disp('the value of constrains'); compare(x_al)

disp('the opt point');

x=x_al;

fmin=obj(X);

>> clear

>> x_al=[0,0];

>> [x,fmin,k]=almain(x_al)

上机大作业3： 1、