# 大连理工大学2021最优化方法大作业（3）

1.无约束优化问题公式

2.结束条件的公式

3.乘子迭代的公式

3.1这是等式约束的乘子迭代，这里的μ对应于框图里的λ

3.2下面是不等式约束的乘子迭代

``````%这三个设定初始值
x = [0;0];
eps = 1e-4;
start_Lagrange(x,eps);

%定义原题的方程式
function f= fun(x)
f=(x(1)-2)^4+(x(1) -2*x(2))^2;
end

%这里写出给出的约束条件，一个等式约束，三个不等式约束
function [h,g] =constrains(x)
h=x(1)^2-x(2);
g=zeros(3,1);
g(1)=4-x(1)^2-x(2)^2;
g(2)=x(1);
g(3)=x(2);
end

%通过原方程和约束方程确定的拉格朗日增广函数，也就是要求的无约束优化问题
function fei=Lagrange(x,lamda,miu,c) %拉格朗日增广函数
[h,gk] = constrains(x);
fei = fun(x) + lamda*h + (c/2)*h^2;
for i = 1:3
fei = fei + (1/(2*c))*((min(0,miu(i) + c*gk(i)))^2 - miu(i)^2);
end
end

%这是求增广函数的梯度
g = zeros(1,2);
[hh,gg] = constrains(x);
g(1) = 2*x(1) - 4*x(2) + 4*(x(1) - 2)^3 + 2*lamda*x(1) - 2*c*x(1)*(- x(1)^2 + x(2));
g(2) = 8*x(2) - 4*x(1) - lamda + (c*(- 2*x(1)^2 + 2*x(2)))/2;
if miu(1) + c*gg(1) < 0
g(1) =g(1)-2*x(1)*(miu(1) - c*(x(1)^2 + x(2)^2 - 4));
g(2) = g(2)-2*x(2)*(miu(1) - c*(x(1)^2 + x(2)^2 - 4));
end
if miu(2) + c*gg(2) < 0
g(1) =g(1)+ miu(2) + c*x(1);
g(2) = g(2);
end
if miu(3) + c*gg(3) < 0
g(1) =g(1);
g(2) = g(2)+miu(3) + c*x(2);
end
end

%乘子法的核心部分
function start_Lagrange(x0,eps)
r=0.25;
arfa = 2;%框图里的常量
c = 4;
a = 2;
k = 0;
lamda = 0;
miu = zeros(3,1);

xk = get_xk(x0,lamda,miu,c);%这个就是调用无约束优化方法
[h0,g0] =constrains(xk);
panduan = h0^2+(min(g0(1),-miu(1)/c))^2+(min(g0(2),-miu(2)/c))^2+(min(g0(3),-miu(3)/c))^2;%判断语句
while panduan > eps^2
[h1,g1]=constrains(xk);
beita = norm(h1)/norm(h0);
if beita > r
c = a*c;
end

for i=1:3
miu(i)= min(0,miu(i)+c*g1(i));
end
lamda = lamda + c*h1;
fprintf('The %d-th iteration, the residual is %f\n',k,panduan);
k = k+1;
x0 = xk;
xk = get_xk(x0,lamda,miu,c);

[h0,g0] =constrains(x0);
panduan = h0^2+(min(g0(1),-miu(1)/c))^2+(min(g0(2),-miu(2)/c))^2+(min(g0(3),-miu(3)/c))^2;

end
fprintf('The %d-th iteration, the residual is %f\n',k,panduan);
fprintf('x=[%f,%f],min(f):%f\n',xk(1),xk(2),fun(xk));
end

%bfgs算法，大家应该能看出来和上一篇一样，只是改了点参数来嵌入乘子法
function x1 = get_xk(x0,lamda,miu,c)
n=2;
h0 = eye(2,2);
s0 = -h0*g0.';
k = 0;
count = 0;
lambda = wolfe_powell(x0,s0,lamda,miu,c);
x1 = x0 +lambda*s0;
eps = 1e-6;
while (norm(g1) > eps)

a = lamda1;
lamda1 = min([2*lamda1,(b+lamda1)/2]);
continue;
end
break;
end
end

function hk = get_hk(h,x,g)%进来的是列向量
miu = 1 + g.'*h*g/(x.'*g);
fenzi = miu*x*x.'-h*g*x.'-x*g.'*h;
hk = h + fenzi/(x.'*g);
end``````

创作不易，还请各位校友感觉有用的话点个赞，嘿嘿嘿~~~~~~~~~~~~