44、MATLAB 中的优化算法实现-优快云博客

本文链接：https://blog.youkuaiyun.com/perl8/article/details/155221019

MATLAB 中的优化算法实现

1. 引言

在工程和科学计算中，优化问题是一个常见且重要的问题。本文将介绍几种在 MATLAB 中实现的优化算法，包括广义简约梯度法（GRG）、序列二次规划法（SQP）、循环坐标法和 Hooke - Jeeves 模式搜索法，并给出相应的代码示例和操作步骤。

2. 广义简约梯度法（GRG）

2.1 函数定义

首先定义了一个名为 delgrgg 的函数，用于计算梯度：

function dg = delgrgg(x)
x1 = x(1); x2 = x(2); x3 = x(3); x4 = x(4); 
dg = [2*x1+1 2*x2-1 3*x3+1 2*x4-1 1 0 0;
       2*x1-1 4*x2   2*x3   4*x4-1 0 1 0;
       4*x1+2 2*x2-1 2*x3     -1   0 0 1];
end

2.2 执行命令

通过以下命令执行 grgopt 函数来求解优化问题：

x0 = [0 0 0 0 7 11 6]; 
xl = [-100 -100 -100 -100 0 0 0]; 
xu = 100*ones(1,7); 
kmax = 1e3; 
crit = 1e-4;
[xopt,fopt,iter] = grgopt(@grgfun,@delgrgf,@delgrgg,x0,xl,xu,kmax,crit)

执行结果如下：

xopt =
    0.0081    0.4922    0.6609    0.2043    6.3066   10.2076    6.0011
fopt =
   86.2196
iter =
     5

3. 序列二次规划法（SQP）

3.1 问题描述

考虑一个有约束的最小化问题：
Minimize (f(x))
Subject to (v_i(x) = 0) ((i = 1, \ldots, p)), (w_k(x) \geq 0) ((k = 1, \ldots, q))
其中 (x = [x_1, x_2, \ldots, x_n]^T) 是 (n) 个实值变量的列向量。

3.2 算法步骤

初始化 ：选择 (x_0) 和 (\mu_0) 使得 (w_k(x_0) \geq 0) ((k = 1, \ldots, q)) 且 (\mu_0 \geq 0)。设置 (k = 0) 和 (Z_0 = I_n)。
计算相关量 ：计算 (g_k)，(A_{ek})，(A_{ik})，(v_k) 和 (w_k)。
求解二次规划问题 ：求解二次规划问题得到 (\delta x)，并计算拉格朗日乘子 (\lambda_{k + 1})，(\hat{\mu}_{k + 1}) 和 (\alpha_k)。
更新点 ：设置 (\delta x = \alpha_k\delta x) 并更新 (x_{k + 1} = x_k + \delta x)。
判断终止条件 ：如果 (|\delta x| \leq \varepsilon)，停止计算。
更新 Hessian 矩阵 ：计算 (\gamma_k)，(\theta)，(\eta_k) 和 (Z_{k + 1})。设置 (k = k + 1) 并返回步骤 2。

3.3 MATLAB 实现

主函数 `sqpopt`

function [xopt,fopt,iter] = sqpopt(fun,dfun,x0,lam0,mu0,crit)
x = x0(:); 
n = length(x); 
lam1 = lam0 + 1; 
Hj = eye(n); 
fv = fun(x); 
aj = fv(2:lam1); 
cj = fv((lam1+1):(lam1+mu0));
Gj = dfun(x); 
gj = Gj(:,1); 
Aej = Gj(:,2:lam1)'; 
Aij = Gj(:,(lam1+1):(lam1+mu0))'; 
iter = 0; 
d = 1;
while d >= crit
    delx = quadpr(Hj,gj,Aej,-aj,Aij,-cj,zeros(n,1),crit);
    ad = Aij*(x+delx) + cj; 
    k = find(ad <= crit); 
    nk = length(k); 
    muj = zeros(mu0,1);
    if nk == 0
        lamj = inv(Aej*Aej')*Aej*(Hj*delx+gj);
    else
        Aaik = Aij(k,:); 
        Aaj = [Aej; Aaik]; 
        mun = inv(Aaj*Aaj')*Aaj*(Hj*delx+gj);
        lamj = mun(1:lam0); 
        mujh = mun(lam1:end); 
        muj(k) = mujh;
    end
    alpha = linsearch(fun,x,delx,lam1,muj,crit); 
    delx = alpha*delx; 
    x = x + delx; 
    grd = dfun(x);
    grd1 = grd(:,1); 
    Agrd = grd(:,2:lam1)'; 
    Am = grd(:,(lam1+1):(lam1+mu0))';
    gamj = (grd1-gj)-(Agrd-Aej)'*lamj-(Am-Aij)'*muj; 
    qj = Hj*delx; 
    dg = delx'*gamj; 
    dq = delx'*qj;
    if dg >= 0.2*dq
        theta = 1;
    else
        theta = 0.8*dq/(dq-dg);
    end
    eta = theta*gamj + (1-theta)*qj; 
    invdq = 1/dq; 
    invde = 1/(delx'*eta);
    Hj = Hj + invde*(eta*eta') - invdq*(qj*qj'); 
    Aej = Agrd; 
    Aij = Am; 
    gj = grd1; 
    fv = fun(x);
    aj = fv(2:lam1); 
    cj = fv((lam1+1):(lam1+mu0)); 
    d = norm(delx); 
    iter = iter + 1;
end
xopt = x; 
fopt = fv(1);
end

线性搜索函数 `linsearch`

function alpha = linsearch(fun,xj,dx,lam,muj,crit)
nmuj = length(muj); 
alrange = 0:0.01:1; 
nint = length(alrange);
hz = zeros(nint,1);
for j = 1:nint
    xdj = xj + alrange(j)*dx; 
    fv = fun(xdj); 
    af = fv(2:lam); 
    cf = fv((lam+1):(lam+nmuj));
    hz(j) = fv(1) + 1e2*sum(af.^2)- muj'*cf;
end
[mval,indhz] = min(hz); 
atemp = alrange(indhz); 
indmu = find(muj <= crit); 
mc = length(indmu);
if mc == 0
    alpha = 0.95*atemp;
else
    dv = zeros(mc,1);
    for k = 1:mc
        for j = 1:nint
            aj = alrange(j); 
            xdj = xj + aj*dx; 
            fcj = fun(xdj);
            cj = fcj((lam+1):(lam+nmuj)); 
            hz(j) = cj(indmu(k));
        end
        indhz = find(hz < 0); 
        hc = length(indhz);
        if hc == 0
            dv(k) = 1;
        else
            dv(k) = alrange(indhz(1)-1);
        end
    end
    mdv = min(dv); 
    alpha = 0.95*min(atemp,mdv);
end
end

二次规划函数 `quadpr`

function xqr = quadpr(Q,c,Aeq,beq,Ane,bne,x0,crit)
nbeq = length(beq); 
nbne = length(bne);
rnbne = nbne + 1.5*sqrt(nbne); 
aone = 1-crit; 
x = x0(:); 
y = Ane*x - bne; 
zeta = zeros(nbeq,1); 
gama = ones(nbne,1); 
ym = y.*gama; 
sumy = sum(ym); 
iter = 0;
while sumy > crit
    tau = sumy/rnbne; 
    resid = -Q*x - c + Aeq'*zeta + Ane'*gama;
    diffr = beq - Aeq*x; 
    numt = tau - y.*gama; 
    yj = gama./y; 
    ysj = numt./y; 
    ydm = diag(yj);
    Gr = inv(Q + Ane'*ydm*Ane); 
    ag = Aeq*Gr*Aeq';
    ayj = resid + Ane'*ysj; 
    diffag = diffr - Aeq*Gr*ayj; 
    delz = inv(ag)*diffag;
    dx = Gr*(ayj + Aeq'*delz); 
    dy = Ane*dx; 
    dgama = (numt - (gama.*dy))./y;
    indny = find(dy < 0); 
    yr = min(y(indny)./(-dy(indny)));
    indny = find(dgama < 0); 
    gamr = min(gama(indny)./(-dgama(indny)));
    aj = aone*min([1 yr gamr]); 
    x = x + aj*dx; 
    gama = gama + aj*dgama; 
    zeta = zeta + aj*delz;
    y = Ane*x - bne; 
    ym = y.*gama; 
    sumy = sum(ym);
    iter = iter + 1;
end
xqr = x;
end

3.4 示例

考虑以下约束最小化问题：
Minimize (f(x) = 2x_1 + x_2^2)
Subject to (h(x) = x_1^2 + x_2^2 - 8 = 0)，(0 \leq x_1 \leq 4)，(1 \leq x_2 \leq 5)

定义目标函数和约束函数

function fv = fun(x)
fv = [2*x(1) + x(2) ^2;       % objective function
      x(1)^2 + x(2)^2 - 8;  % equality constraint: h(x) = 0
      x(1);             % inequality constraint: g1(x) >= 0
      -x(1) + 4;            % inequality constraint: g2(x) >= 0
      x(2) - 1;             % inequality constraint: g3(x) >= 0
      -x(2) + 5];           % inequality constraint: g4(x) >= 0
end

定义梯度函数

function dfv = dfun(x)
df = [2 2*x(2)];   % gradient of objective function
dh = [2*x(1) 2*x(2)];  % gradient of equality constraint (h(x) = 0)
dg = [1 0;-1 0;0 1;0 -1]; % gradient of inequality constraint (gi(x) >= 0)
dfv = [df' dh' dg'];
end

执行 SQP 算法

x0 = [2 3]'; 
lam0 = 1; 
mu0 = 3; 
crit = 1e-6;
[xopt,fopt,iter] = sqpopt(@fun,@dfun,x0,lam0,mu0,crit)

执行结果如下：

xopt =
    2.6458
    1.0000
fopt =
    6.2915
iter =
     8

3.5 SQP 算法流程图

graph TD;
    A[初始化x0, μ0, k = 0, Z0 = In] --> B[计算gk, Aek, Aik, vk, wk];
    B --> C[求解二次规划问题得到δx, 计算λk+1, ˆμk+1, αk];
    C --> D[更新δx = αkδx, xk+1 = xk + δx];
    D --> E{||δx|| ≤ ε?};
    E -- 是 --> F[停止计算];
    E -- 否 --> G[计算γk, θ, ηk, Zk+1];
    G --> H[k = k + 1];
    H --> B;

4. 直接搜索方法

4.1 循环坐标法

4.1.1 算法原理

在循环坐标搜索方法中，沿着每个坐标方向进行搜索以找到最小值。设 (u_k) 是沿着坐标方向 (k) 的单位向量。确定使 (f(\alpha) = f(x + \alpha u_k)) 最小的 (\alpha_k) 值，并在沿着方向 (k) 的搜索结束时移动到新点 (x + \alpha_k u_k)。将所有可能方向的搜索作为一个阶段。当梯度似乎为零时，搜索最小值的过程终止。

4.1.2 算法步骤

选择一个起始点 (x_1) 并计算 (f_1 = f(x_1))。设置 (k = 1)，(x = x_1) 和 (f = f_1)。
确定使 (f(\alpha) = f(x + \alpha u_k)) 最小的 (\alpha_k)。
更新当前点 (x = x + \alpha_k u_k) 并计算 (f = f(x))。
如果 (k < n)，设置 (k = k + 1) 并返回步骤 2。
指定方向 (d = x - x_1) 并确定使 (f(x + \alpha_d)) 最小的 (\alpha_d)。
移动到新点 (x = x + \alpha_d d) 并计算 (f = f(x))。
如果 (|d| > \varepsilon) 或 (|f - f_1| > \varepsilon)，设置 (x_1 = x)，(f_1 = f) 并返回步骤 2。

4.1.3 MATLAB 实现

function [xopt,fopt,iter] = cycopt(fcyc,x0,crit)
stsize = 0.1; 
x = x0; 
xk = x; 
nv = length(x); 
fk = fcyc(xk); 
fold = fk; 
iter = 0;
while (1)
    iter = iter + 1;
    for j = 1:nv
        h(j) = 0;
    end
    for k = 1:nv
        for j = 1:nv
            d(j) = 0;
        end
        d(k) = 1; 
        aut = 0;
        [stfit,fvfit] = quadfit(fcyc,aut,fk,stsize,xk,d,crit);
        xk(k) = xk(k) + stfit; 
        h(k) = stfit; 
        fk = fvfit;
    end
    aut = 0; 
    [stfit,fvfit] = quadfit(fcyc,aut,fk,stsize,xk,h,crit);
    for j = 1:nv
        xk(j) = xk(j) + stfit * h(j);
    end
    if (abs(fk-fold) <= crit)
        break;
    end
    fold = fk;
end
xopt = xk; 
fopt = fcyc(xk);
end

4.1.4 示例

考虑 Powell 函数：
(f(x) = (x_1 + 10x_2)^2 + 5(x_3 - x_4)^2 + (x_2 - 2x_3)^4 + 10(x_1 - x_4)^4)

x0 = [-3 -1 0 1]; 
crit = 1e-4;
fc = @(x) (x(1)+10*x(2))^2+5*(x(3)-x(4))^2+(x(2)-2*x(3))^4+10*(x(1)-x(4))^4;
[xopt,fopt,iter] = cycopt(fc,x0,crit)

执行结果如下：

xopt =
    0.1042   -0.0097    0.0543    0.0566
fopt =
   3.1864e-04
iter =
    14

4.2 Hooke - Jeeves 模式搜索法

4.2.1 算法原理

在 Hooke - Jeeves 模式搜索方法中，选择一个初始步长 (\alpha) 并从给定的起始点开始搜索。设 (u_k) 是沿着坐标方向 (x_k) 的单位向量。评估 (x + \alpha u_k) 处的函数值，如果函数值减小，则更新当前点 (x) 为 (x + \alpha u_k)。如果函数值不减小，则计算 (f(x - \alpha u_k))。如果该函数值减小，则将 (x) 更新为 (x - \alpha u_k)。

4.2.2 算法步骤

选择一个起始点 (x_1)，一个初始步长 (\alpha)，一个步长缩减参数 (c) 和一个加速参数 (\beta)。
设置 (x_Q = x_1) 并在点 (x_Q) 周围进行搜索。如果函数值减小，则将生成的新点定义为新的起始点 (x_Q) 并将旧的起始点重命名为 (x_Q’)。如果搜索不成功，则将步长缩减为 (\alpha = c\alpha)。
找到点 (x_P = x_Q + \beta(x_Q - x_Q’)) 并在点 (x_P) 周围进行探索。
如果函数值减小，则将生成的新点定义为新的起始点 (x_Q)，将旧的起始点重命名为 (x_Q’) 并返回步骤 3。
如果搜索不成功，则在点 (x_Q) 周围进行探索。将步长缩减为 (\alpha = c\alpha) 直到收敛。将新点作为新的起始点并返回步骤 3。

4.2.3 MATLAB 实现

function [xopt,fopt,iter] = hjopt(fhj,x0,crit)
inist = 0.5; 
fr = 0.125; 
x = x0; 
xb = x; 
n = length(x0); 
stsize = inist;
nc = 0; % number of contraction steps
nb = 0; % number of base changes
np = 0; % number of pattern moves
fb = fhj(xb); 
iter = 1;
while (1)
    fk = fb; 
    [fk, x] = search(fhj, n, stsize, fk, x);
    if (fb-fk > crit)
        while (2)
            icv = 0;
            for j = 1:n
                cp = x(j); 
                x(j) = 2*cp - xb(j); 
                xb(j) = cp;
            end
            fb = fk; 
            fk = fhj(x); 
            [fk, x] = search(fhj, n, stsize, fk, x);
            if (fb-fk <= crit)
                x = xb;
                if (nb > 1) 
                    if abs(fk - fold) < crit
                        icv = 1; 
                        break;
                    end
                end
                nb = nb + 1; 
                break;
            end
            np = np + 1;
        end
        if (icv == 1)
            break;
        end
    else
        fold = fk;
        if (stsize < crit)
            break;
        end
        stsize = fr*stsize; 
        nc = nc + 1;
    end
    iter = iter + 1;
end
xopt = x; 
fopt = fk;
end

function [fk, xk] = search(fhj, n, stsize, fk, x)
for k = 1:n
    cpt = x(k); 
    x(k) = cpt + stsize;
    f = fhj(x);
    if (f < fk)
        fk = f;
    else
        x(k) = cpt - stsize; 
        f = fhj(x);
        if (f < fk )
            fk = f;
        else
            x(k) = cpt;
        end
    end
end
xk = x;
end

4.2.4 示例

考虑 Rosenbrock 函数：
(f(x) = 100(x_1^2 - x_2)^2 + (1 - x_1)^2)

f = @(x) 100*(x(1)^2 - x(2))^2 + (1 - x(1))^ 2;
x0 = [-1 1]; 
crit = 1e-6;
[xopt,fopt,iter] = hjopt(f,x0,crit)

执行结果如下：

xopt =
    1.0078    1.0156
fopt =
   6.0860e-05
iter =
    10

4.3 Hooke - Jeeves 模式搜索法流程图

graph TD;
    A[选择起始点x1, 初始步长α, 步长缩减参数c, 加速参数β] --> B[设置xQ = x1, 搜索xQ];
    B --> C{函数值减小?};
    C -- 是 --> D[定义新点为xQ, 旧点为xQ'];
    C -- 否 --> E[α = cα];
    D --> F[找到xP = xQ + β(xQ - xQ'), 探索xP];
    F --> G{函数值减小?};
    G -- 是 --> D;
    G -- 否 --> H[探索xQ, α = cα直到收敛];
    H --> F;

5. 总结

本文介绍了几种在 MATLAB 中实现的优化算法，包括广义简约梯度法（GRG）、序列二次规划法（SQP）、循环坐标法和 Hooke - Jeeves 模式搜索法。通过具体的代码示例和操作步骤，展示了如何使用这些算法解决有约束和无约束的优化问题。每种算法都有其适用场景和优缺点，在实际应用中可以根据具体问题选择合适的算法。

6. 算法对比分析

6.1 算法特点对比

算法名称	特点	适用场景	优点	缺点
广义简约梯度法（GRG）	通过定义梯度函数，使用 `grgopt` 函数求解优化问题	有约束的优化问题	实现相对简单，能有效处理一定规模的约束优化问题	对于复杂的约束条件和目标函数，收敛速度可能较慢
序列二次规划法（SQP）	将有约束的优化问题转化为一系列二次规划问题求解	有约束的非线性优化问题	收敛速度较快，能处理复杂的约束条件	计算量较大，需要计算目标函数和约束函数的梯度
循环坐标法	沿着每个坐标方向进行搜索以找到最小值	无约束或简单约束的优化问题	算法简单，易于实现	收敛速度慢，对于高维问题效率较低
Hooke - Jeeves 模式搜索法	通过选择初始步长，在起始点周围进行搜索和模式移动	无约束或简单约束的优化问题	不需要计算梯度，对目标函数的要求较低	收敛速度不稳定，可能需要较多的迭代次数

6.2 示例问题对比

以之前的示例问题为例，对比不同算法的性能：
- Powell 函数 ：使用循环坐标法求解，迭代 14 次得到结果 (f_{opt}=3.1864\times10^{-4})。
- Rosenbrock 函数 ：使用 Hooke - Jeeves 模式搜索法求解，迭代 10 次得到结果 (f_{opt}=6.0860\times10^{-5})。
- 约束最小化问题 ：使用序列二次规划法求解，迭代 8 次得到结果 (f_{opt}=6.2915)。

从这些示例可以看出，不同算法在不同问题上的表现有所差异。对于复杂的约束问题，SQP 算法可能更具优势；而对于简单的无约束问题，循环坐标法和 Hooke - Jeeves 模式搜索法可能更易于实现。

7. 实际应用建议

7.1 问题分析

在选择优化算法之前，需要对问题进行全面的分析，包括目标函数的性质（如线性、非线性、凸性等）、约束条件的类型（等式约束、不等式约束）以及问题的维度等。

7.2 算法选择

线性问题 ：如果目标函数和约束条件都是线性的，可以优先考虑使用线性规划算法。
非线性问题 ：
- 对于有约束的非线性问题，序列二次规划法（SQP）通常是一个不错的选择。
- 对于无约束或简单约束的非线性问题，可以尝试使用循环坐标法或 Hooke - Jeeves 模式搜索法。
- 如果问题的规模较小且目标函数的梯度计算较为复杂，可以考虑使用广义简约梯度法（GRG）。

7.3 参数调整

在使用优化算法时，需要合理调整算法的参数，如起始点、步长、收敛准则等。不同的参数设置可能会影响算法的收敛速度和结果的准确性。可以通过多次试验来找到合适的参数值。

8. 结论

本文详细介绍了几种常见的优化算法在 MATLAB 中的实现，包括广义简约梯度法（GRG）、序列二次规划法（SQP）、循环坐标法和 Hooke - Jeeves 模式搜索法。通过具体的代码示例和操作步骤，展示了如何使用这些算法解决不同类型的优化问题。在实际应用中，需要根据问题的特点选择合适的算法，并合理调整参数，以提高算法的效率和结果的准确性。同时，还可以结合其他优化技术和算法，进一步优化问题的求解过程。