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

前面无约束优化(constrained optimization)问题的解法在这： 大连理工大学 2021年最优化方法大作业（1）_JiangTesla的博客-CSDN博客

大连理工大学 2021年最优化方法大作业（2）_JiangTesla的博客-CSDN博客

这是我们第二题，约束优化(constrained optimization)方法中的乘子法。 这个题目既包含等式约束(equality constraint)也包含不等式约束(inequality constraint)(equality constraint)，所以需要使用Rockafellar的乘子法，绘制框图如下：

从框图可以看出来乘子法的逻辑思路很简单，就是先给出初始乘子，然后利用乘子可以形成一个无约束优化(constrained optimization)问题，从给出的初始点来求解这个问题，得到的解就是下次迭代的点。换句话说，乘子法的迭代方式就是不停地求解新的无约束优化(constrained optimization)问题。

下面补充上面没写的公式

1.无约束优化(constrained optimization)问题公式

2.结束条件的公式

3.乘子迭代的公式

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

3.2下面是不等式约束(inequality constraint)(equality constraint)的乘子迭代

求解无约束优化(constrained optimization)方法之前有提到，我选用的是最好的BFGS方法，下面直接上代码讲解

%这三个设定初始值
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
%这里写出给出的约束条件，一个等式约束(equality constraint)，三个不等式约束(inequality constraint)(equality constraint)
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
%通过原方程和约束方程确定的拉格朗日增广函数，也就是要求的无约束优化(constrained optimization)问题
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
%这是求增广函数的梯度(gradient)                
function g = grad_fai(x,lamda,miu,c)
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);%这个就是调用无约束优化(constrained optimization)方法
[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算法，大家应该能看出来和上一篇一样，只是改了点参数来嵌入(Embedding)乘子法
function x1 = get_xk(x0,lamda,miu,c)
n=2;
g0 = grad_fai(x0,lamda,miu,c);
h0 = eye(2,2);
s0 = -h0*g0.';
k = 0;
count = 0;
lambda = wolfe_powell(x0,s0,lamda,miu,c);
x1 = x0 +lambda*s0;
g1 = grad_fai(x1,lamda,miu,c);
eps = 1e-6;
while (norm(g1) > eps)
    if k<n-1
        detax = x1 - x0;
        detag = g1.' - g0.';
        h1 = get_hk(h0,detax,detag);
        s1 = -h1*g1.';
        k = k+1;
        x0 = x1;
        g0 = grad_fai(x0,lamda,miu,c);
        s0 = s1;
        h0 = h1;
        lambda = wolfe_powell(x0,s0,lamda,miu,c);
        x1 = x0 +lambda*s0;
        g1 = grad_fai(x1,lamda,miu,c);
    else 
      x0 = x1;
      g0 = grad_fai(x0,lamda,miu,c);
      h0 = eye(2,2);
      s0 = -h0*g0.';
      lambda = wolfe_powell(x0,s0,lamda,miu,c);
      x1 = x0 +lambda*s0;
      g1 = grad_fai(x1,lamda,miu,c);
      k = 0;
    end
    count=count+1;
end
end
%一维搜索
function lamda1 = wolfe_powell(xk,dk,lamda,miu,c)
c1 = 0.1;c2=0.5;
a = 0; b =Inf;
lamda1 = 1;
while(1)
    if ~(Lagrange(xk+lamda1*dk,lamda,miu,c)-Lagrange(xk,lamda,miu,c) <= c1*lamda1*grad_fai(xk,lamda,miu,c)*dk)
        b = lamda1;
        lamda1 = (lamda1 + a)/2;
    continue;
    end
    if ~(grad_fai(xk+lamda1*dk,lamda,miu,c)*dk >= c2*grad_fai(xk,lamda,miu,c)*dk)
        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

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