45、MATLAB 优化算法：从 Rosenbrock 到遗传算法的实现与应用

perl8

于 2025-11-22 12:57:20 发布

阅读量13

点赞数

CC 4.0 BY-SA版权

分类专栏： MATLAB赋能化工计算文章标签： MATLAB 优化算法 Rosenbrock方法

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

MATLAB赋能化工计算专栏收录该内容

48 篇文章 ¥499.90

订阅专栏¥69.90

会员秒杀 ¥9.9 重磅福利

超级会员免费看

MATLAB 优化算法：从 Rosenbrock 到遗传算法的实现与应用

1. Rosenbrock 方法

Rosenbrock 方法在每一步沿着 n 个正交方向进行搜索，然后确定新的正交方向。以下是该方法的步骤：
1. 设置方向向量为单位向量，选择起始点 $x_0$、初始步长 $\alpha_i$、步长扩展参数 $\beta$（$\beta>1$）和步长缩减参数 $c$（$c<1$），并初始化相关变量。
2. 根据函数值比较更新步长和点的位置。若新点函数值更小，则扩展步长；否则，缩减步长。
3. 构建新的线性无关方向，计算其模长。若满足停止准则，则停止计算；否则，使用 Gram - Schmidt 过程形成新的正交方向。
4. 进入下一阶段计算，以新的坐标方向作为起始向量。

以下是实现 Rosenbrock 方法的 MATLAB 函数 rbopt ：

function [xopt,fopt,istag] = rbopt(fros,x0,crit)
% rbopt.m: minimization by the Rosenbrock’s method
% Inputs:
%   fros: objective functions
%   x0: starting point
%   crit: stopping criterion
% Outputs:
%   xopt: optimal point
%   fopt: objective function value at x=xopt
%   istag: number of stages
% Example:
% x0 = [-3 -1 0 1]; crit = 1e-4;
% f = @(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,istag] = rbopt(f,x0,crit)     
n = length(x0);  % number of variables
mi = 7; stsize = 1; alpha = 3; beta = 0.5; xt = x0; 
% Initial vectors along coordinates
for j = 1:n
    for k = 1:n, v(j,k) = 0; end
    v(j,j) = 1;
end
ft = fros(xt); istag = 0;

while (1)
    istag = istag + 1;
    % Initialization
    for j = 1:n, stp(j) = stsize; stj(j) = 0; sfv(j,1) = 0; sfv(j,2) = 0; end
    iter = 0; icont = 0; idn = 0;
    while (2)
        iter = iter + 1;
        if iter > n, iter = 1; end
        for j = 1:n, x(j) = xt(j) + stp(iter)*v(j,iter); end
        f = fros(x);
        if f < ft
            stj(iter) = stj(iter) + stp(iter); % successful
            for j = 1:n, xt(j) = x(j); end
            ft = f; idn = 0; stp(iter) = alpha*stp(iter);
            if sfv(iter,1) == 0, sfv(iter,1) = 1; icont = icont + 1; end
        else % failure
            stp(iter) = -beta*stp(iter);
            if sfv(iter,2) == 0, sfv(iter,2) = 1; icont = icont + 1; end
            idn = idn + 1;
        end
        if icont == 2*n, break; end
        if idn > 50, mi = mi + 1; return; end
    end
    if istag == 1, fold = ft;
    else
        if abs(ft - fold) < crit, mi = mi + 1; display('Convergence achieved.'); break; end
    end
    fold = ft;
    % Construct new vectors
    for k = 1:n 
        for j = 1:n, u(k,j) = 0; end
    end
    for i = 1:n
        for j = 1:n
            for k = 1:n, u(k,i) = u(k,i) + stj(j)*v(k,j); end
        end
    end
    % Orthogonalization by Gram-Schmidt procedure
    for i = 1:n
        if i > 1
            for j = 1:n, v(j,i) = 0; end
            for k = 1:i-1
                c = 0;
                for j = 1:n, c = c + u(j,i)*v(j,k); end
                for j = 1:n, v(j,i) = v(j,i) + c*v(j,k); end
            end
            for j = 1:n, u(j,i) = u(j,i) - v(j,i); end
        end
        c = 0;
        for j = 1:n, c = c + u(j,i)*u(j,i); end
        c = sqrt(c);
        % Take new step length as the length of the first vector
        if i == 1
            stsize = c;
            if stsize < crit, mi = mi + 1; break; end
            if c < crit % Orthogonality is lost: reset vectors.
                for j = 1:n
                    for k = 1:n, v(j,k) = 0; end
                    v(j,j) = 1;  break;
                end
            end
            for j = 1:n, v(j,i) = u(j,i)/c; end
        end
    end
end
xopt = xt; fopt = ft;
end

示例：求函数 $f(x) = x_1^2 + 2x_2^2 + 2x_1x_2$ 的最小值，使用 $x_0 = [0.5, 1]$ 和 $crit = 1\times10^{-6}$。

f = @(x) x(1)^2 + 2*x(2)^2 + 2*x(1)*x(2); x0 = [0.5 1]; crit = 1e-6;
[xopt,fopt,istag] = rbopt(f,x0,crit)

2. Nelder 和 Mead 的单纯形法

在 n 维空间中，n + 1 个点构成一个单纯形。该方法通过反射、扩展或收缩操作来更新单纯形的顶点，以寻找函数的最小值。以下是具体步骤：
1. 选择起始点 $x_1$ 并计算函数值 $f_1$，设置 $k = 0$。
2. 生成单纯形并计算各顶点的函数值，$k = k + 1$。
3. 确定最高、第二高和最小函数值及其对应的坐标。若满足收敛条件，则停止计算。
4. 计算不包含最高点的 n 个点的平均值 $x_m$ 和反射点 $x_r$。根据反射点的函数值决定后续操作。
5. 若反射点函数值小于最小函数值，则进行扩展操作；否则，根据不同情况进行收缩操作。
6. 若满足收敛条件，则停止计算；否则，更新单纯形并继续迭代。

以下是实现 Nelder 和 Mead 单纯形法的 MATLAB 函数 nmopt ：

function [xopt,fopt,iter] = nmopt(fun,x0,crit)
% nmopt.m: minimization by the Nelder and Mead's simplex method 
% Inputs:
%   fun: objective functions
%   x0: starting point
%   crit: stopping criterion
% Outputs:
%   xopt: optimal point
%   fopt: objective function value at x=xopt
%   iter: number of iterations
% Example:
% x0 = [-3 -1 0 1]; crit = 1e-4;
% f = @(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] = nmopt(f,x0,crit)
rf = 1; % reflection parameter
ef = 2; % expansion parameter
cf = 0.5; % contraction parameter
sf = 0.5; % scale parameter
n = length(x0); beta = 1; n1 = n+1; iter = 0; 
x = x0; fv(1) = fun(x);
for j = 1:n, pt(1,j) = x(j); end
% Create simplex and evaluate function values
for k = 2:n1
    u = zeros(1,n); u(1,k-1) = 1;
    pt(k,:) = pt(1,:) + beta*u; x = pt(k,:); fv(k) = fun(x);
end
while (1)
    nloop = 0;
    while (2)
        nloop = nloop + 1;
        if (nloop == 2), break; end
        iter = iter + 1;
        % Find highest, lease, and second highest values
        [fh, indh] = max(fv); [fl, indl] = min(fv); fs = fv(indl); inds = indl;
        for k = 1:n1
            if (k ~= indh)
                if (fs < fv(k)), fs = fv(k); inds = k; end
            end
        end
        % Find the average of n points xm
        for k = 1:n
            xm(k) = 0;
            for j = 1:n1
                if (j ~= indh), xm(k) = xm(k) + pt(j,k); end
            end
            xm(k) = xm(k)/n;
        end
        % Reflection procedure
        xr = xm + rf*(xm - pt(indh,:)); fr = fun(xr);
        if (fr >= fl & fr <= fs) % accept reflection
            fv(indh) = fr; pt(indh,:) = xr; break;
        end
        % Expansion procedure
        if (fr < fl )
            xe = xr + ef*(xr - xm); fe = fun(xe);
            if (fe < fl) % accept expansion
                fv(indh) = fe; pt(indh,:) = xe; break;
            else % accept reflection
                fv(indh) = fr; pt(indh,:) = xr; break;
            end
        end
        % Contraction procedure
        if (fr > fh)
            xc = xm + cf*(pt(indh,:) - xm); [fc] = fun(xc);
            if (fc <= fh) % accept contraction
                fv(indh) = fc; pt(indh,:) = xc; break;
            end
        elseif (fr > fs & fr  <= fh)
            xc = xm + cf*(xr - xm); fc = fun(xc);
            if (fc <= fr) % accept contraction
                fv(indh) = fc; pt(indh,:) = xc; break;
            end
        end
        % Scaling
        for k = 1:n1
            if (k ~= indl)
                for j = 1:n, x(j) = sf *pt(k,j) + (1 - sf) *pt(indl,j); 
                pt(k,j) = x(j); end
                fv(k) = fun(x);
            end
        end
    end
    sigma = std(fv); avg = mean(fv);
    if (sigma <= crit)
        inc = inc + 1;
        if (inc == 2), break; end
    else
        inc = 0;
    end
end
xopt = pt(indl,:); fopt = fl;
end

示例：求 Powell 函数的最小值，使用 $x_0 = [-3, -1, 0, 1]$ 和 $crit = 1\times10^{-6}$。

x0 = [-3 -1 0 1]; crit = 1e-6;
f = @(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] = nmopt(f,x0,crit)

3. 模拟退火（SA）方法

模拟退火方法借鉴了金属退火过程中原子从高能态向低能态转变的原理。在优化问题中，通过一定的概率接受函数值增加的点，从而避免陷入局部最优。以下是该方法的步骤：
1. 选择起始点 $x$ 并计算函数值 $f$，设置 $x_{min} = x$ 和 $f_{min} = f$。
2. 选择步长向量 $s$ 和步长缩减参数 $r_s$。
3. 定义接受率向量 $a$，初始值均为 1。
4. 选择起始温度 $T$ 和温度缩减因子 $r_T$，设置 $k = 1$。
5. 设置温度为 $r_TT$ 和步长为 $r_ss_T$，在每个温度下进行 $N_T$ 次迭代，每次迭代包含 $N_C$ 个循环。
6. 生成随机数 $r$，计算新点 $x_s$。若新点超出边界，则调整其坐标。
7. 计算新点的函数值 $f_s$。若 $f_s \leq f$，则接受新点；否则，以一定概率接受新点。
8. 根据接受情况更新接受率向量 $a$。
9. 更新 $k$ 并返回步骤 5。
10. 根据接受率调整步长。

以下是实现模拟退火方法的 MATLAB 函数 saopt ：

function [xopt,fopt,iter] = saopt(fun,T,xl,xu,rp,rs,crit)
% saopt.m: minimization by simulated annealing (SA) method
% Inputs:
%   fun: objective function
%   T: initial temperature
%   xl,xu: lower and upper bounds on x
%   rp: temperature reduction factor
%   rs: step reduction parameter
%   crit: stopping criterion
% Outputs:
%   xopt: optimal point
%   fopt: objective function value at x=xopt
%   iter: number of cycle iterations
% Example:
%   f = @(x) -cos(5*sqrt(sum(x-5).^2)) + 0.1*sum(x-5).^2;
%   T = 100; xl = -10*[1 1]; xu = 10*[1 1]; rp = 0.8; rs = 0.9; crit = 1e-8;
%   [xopt,fopt,iter] = saopt(f,T,xl,xu,rp,rs,crit)
% Initialization
sf = 2; stp = 1; np = 10; nc = 20; nt = 2e4; ir = 16; rcrit = 1e-10; n = 
length(xu);
% Feasible starting point
for j = 1:n, x(j) = xl(j) + rand*(xu(j) - xl(j)); xs(j) = x(j); xmin(j) = 
x(j); end
f = fun(x); fmin = f; fold = f; citer = 0;
% Set step sizes, step factors and acceptance ratios
for j = 1:n, stsize(j) = stp; ar(j) = 1; end
while (1)
    for piter = 1:np % temperature step loop
        for ic = 1:nc % % search cycles: search along coordinate 
direction
            for k = 1:n
                xs(k) = x(k) + (2*rand - 1)*stsize(k);
                if (xs(k) < xl(k)) | (xs(k) > xu(k)), xs(k) = xl(k) + 
rand*(xu(k) - xl(k));
                end
                fs = fun(xs);
                    if fs <= f % point is accepted: update xmin and fmin
                    x(k) = xs(k); f = fs;
                    if fs < fmin, xmin = xs; fmin = fs; end
                    else
                        p = exp((f - fs)/T);
                        if rand < p, x(k) = xs(k); f = fs;
                        else % point rejected
                            xs(k) = x(k); ar(k) = ar(k) - 1/nc;
                        end
                    end
            end
        end
        % Adjust step so that about half the points are accepted.
        for j = 1:n
            if ar(j) > 0.6, stsize(j) = stsize(j)*(1 + sf*(ar(j) - 0.6)/0.4);
            elseif ar(j) < 0.4, stsize(j) = stsize(j)/(1 + sf*(0.4 
- ar(j))/0.4);
            end
            if stsize(j) > xu(j) - xl(j), stsize(j) = xu(j) - xl(j); end
            ar(j) = 1;
        end
    end
    stsize = stp*(xu - xl); fcrit = crit + rcrit*abs(fmin);
    if (fmin <= fold) & (fold - fmin < fcrit)
        citer = citer + 1;
        if citer >= ir, break; end
    else, citer = 0;
    end
    % Reduce temperature
    T = rp*T; stp = rs*stp; x = xmin; f = fmin; fold = f;
end
xopt = xmin; fopt = fmin; iter = citer;
end

示例：求波纹弹簧函数 $f(x) = -\cos(5\sqrt{(x_1 - 5)^2 + (x_2 - 5)^2}) + 0.1((x_1 - 5)^2 + (x_2 - 5)^2)$ 的最小值，设置 $T = 100$，$r_p = 0.8$，$r_s = 0.9$，$crit = 1\times10^{-8}$。

function fv = fcy(x)
% Corrugated spring function
C = 0; for j = 1:2, C = C + (x(j) - 5)^2; end
fv = -cos(5*sqrt(C)) + 0.1*C;
end

T = 100; xl = -10*[1 1]; xu = 10*[1 1]; rp = 0.8; rs = 0.9; crit = 1e-8;
[xopt,fopt,iter] = saopt(@fcy,T,xl,xu,rp,rs,crit)

各方法对比表格

方法	优点	缺点	适用场景
Rosenbrock 方法	利用正交方向搜索，能有效避免局部最优	计算复杂度较高	适用于多维优化问题
Nelder 和 Mead 的单纯形法	简单易懂，不需要计算梯度	收敛速度较慢	适用于函数梯度难以计算的情况
模拟退火方法	能跳出局部最优，全局搜索能力强	计算时间长	适用于复杂的全局优化问题

模拟退火方法流程图

graph TD;
    A[选择起始点 x 并计算 f] --> B[选择步长向量 s 和步长缩减参数 rs];
    B --> C[定义接受率向量 a];
    C --> D[选择起始温度 T 和温度缩减因子 rT];
    D --> E[设置温度和步长];
    E --> F[生成随机数 r 计算新点 xs];
    F --> G{xs 是否超出边界};
    G -- 是 --> H[调整 xs 坐标];
    G -- 否 --> I[计算 fs];
    H --> I;
    I --> J{fs <= f};
    J -- 是 --> K[更新 x 和 fmin];
    J -- 否 --> L[计算接受概率 p];
    L --> M{rand < p};
    M -- 是 --> K;
    M -- 否 --> N[更新接受率向量 a];
    K --> O{是否完成 NC 个循环};
    N --> O;
    O -- 否 --> F;
    O -- 是 --> P[调整步长];
    P --> Q{是否完成 NT 次迭代};
    Q -- 否 --> E;
    Q -- 是 --> R[降低温度];
    R --> S{是否满足收敛条件};
    S -- 否 --> E;
    S -- 是 --> T[输出结果];

4. 遗传算法（GA）

遗传算法模拟自然进化和选择过程，以寻找最优解。它是一种全局优化方法，利用随机信息进行搜索。以下是基本遗传算法的步骤：
1. 设置初始种群 $P(0)$，评估初始解的适应度，设置 $t = 0$。
2. 使用遗传算子生成子代，添加随机生成的新种群，并再次评估种群适应度。
3. 确定哪些成员将进入下一代（竞争选择），选择种群 $P(t + 1)$。
4. 如果未达到收敛条件，更新 $t$ 并返回步骤 2。

在遗传算法中，问题通常表述为最大化函数 $f(x)$，同时满足 $l_i \leq x_i \leq u_i$（$i = 1, \ldots, n$）。每个变量被表示为 $m$ 位二进制数，通过将变量的可行区间划分为 $2^m - 1$ 个区间来实现。变量 $x_i$ 的区间步长 $s_i$ 为：
[s_i = \frac{u_i - l_i}{2^m - 1}]

初始种群的每个成员是一个长度为 $n * m$ 位的字符串。在评估阶段，提取每个成员的变量值，计算函数值 $f_1, f_2, \ldots$，这些值被称为适应度值。繁殖过程包括创建交配池、交叉和变异操作。缩放适应度值的总和 $S$ 计算如下：
[S = \sum_{j = 1}^{z} \hat{f}_j, \quad \hat{f}_j = Q + R \frac{f_j - f_l}{f_h - f_l}, \quad Q = 0.1, \quad R = 0.9]

使用轮盘赌选择来复制最适应的成员进行繁殖。交叉操作在交配父池上进行以产生子代。变异操作随机改变预期适应度。再次评估种群，存储最高适应度值和对应的变量值作为 $f_{max}$ 和 $x_{max}$，完成一代的计算。如果达到预定的代数，则终止计算过程。

以下是实现遗传算法的 MATLAB 函数 gaopt ：

function [xopt,fopt,iter] = gaopt(fun,xl,xu,nb,ps,ng,mp)
% gaopt.m: maximization by the genetic algorithm 
% Inputs:
%   fun: objective functions
%   xl,xu: lower and upper bounds on x
%   nb: number of binary digits
%   ps: population size
%   ng: number of generations
%   mp: mutation probability
% Outputs:
%   xopt: optimal point
%   fopt: objective function value at x=xopt
%   iter: number of function evaluations
% Example:
%   xl=-6*[1 1 1 1]; xu=6*[1 1 1 1]; nb=8; ps=40; ng=50; mp=0.02;
%   f = @(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] = gaopt(f,xl,xu,nb,ps,ng,mp)
n = length(xu);
Pn = round(rand(ps,n*nb)); % initial population
for kg = 1:ng
    ft = fitness(fun,ps,nb,n,xl,xu,Pn); [fmax,kmax] = max(ft);
    if (kg == 1)
        fbest = fmax;
        for k = 1:n
            ab = num2str(Pn(kmax,nb*(k-1)+1:nb*k));
            adec = bin2dec(ab); % onvert binary string to decimal integer
            xbest(k) = xl(k) + (xu(k)-xl(k))/(2^nb - 1)*adec;
        end
    else
        if (fmax > fbest)
            fbest = fmax;
            for k = 1:n
                ab = num2str(Pn(kmax,nb*(k-1)+1:nb*k)); adec = bin2dec(ab);
                xbest(k) = xl(k) + (xu(k)-xl(k))/(2^nb - 1)*adec;
            end
        end
    end
    [frk,ft] = frank(ps, ft);  % scaling
    for k = 1: ps % shuffling: roulette wheel
        ik = frk(k); ptemp =  ps; ktemp = k; am(ik) = 2*(ptemp + 1 - ktemp)/(ptemp*(ptemp + 1));
    end
    for k = 2: ps, ptemp = am(k - 1); am(k) = ptemp + am(k); end
    %
    for k = 1:ps % selection
        kstr1 = roul(ps,am); kstr2 = roul(ps,am); kchd = cros(Pn,nb,n,kstr1, kstr2);
        kchd = mutat(kchd,mp,nb,n);
        for j = 1: nb*n, Pn(k, j) = kchd(j); end
    end
end
xopt = xbest; fopt = fbest; iter = ng*ps;
end

function ft = fitness(fun,ps,nb,n,xl,xu,Pn)  % function evaluation 
for k = 1:ps
    for j = 1:n
        a = num2str(Pn(k,nb*(j-1) + 1:nb*j)); adec = bin2dec(a);
        x(j) = xl(j) + (xu(j)-xl(j))/(2^nb - 1)*adec;
    end
    f = fun(x); ft(k) = f;
end
end

function [frk,ft] = frank(ps, ft) 
for k = 1:ps, frk(k) = k; end
for k = 1:ps-1
    temp = ft(k); ktemp = k;
    for j = k+1:ps
        if (ft(j) > temp), temp = ft(j); ktemp = j; end
    end
    jtemp = frk(ktemp); frk(ktemp) = frk(k); frk(k) = jtemp;
    ftemp = ft(ktemp); ft(ktemp) = ft(k); ft(k) = ftemp;
end
end

function kstr = roul(ps,am)  % shuffling operation
temp = rand;
for k = 1:ps
    while (1)
        if (temp > am(k)), break; end
        kstr = k; return;
    end
end
end

function kchd = cros(Pn,nb,n,kstr1, kstr2)  % crossover operation
for j = 1:n
    jp = nb*(j - 1); nc = floor(rand*nb + 0.5);
    if (nc == 0)
        for k = 1:nb, kchd(k+jp) = Pn(kstr2, k+jp); end
    elseif (nc == nb)
        for k = 1:nb, kchd(k+jp) = Pn(kstr1, k+jp); end
    else
        for k = 1: nc, kchd(k+jp) = Pn(kstr1, k+jp); end
        for k = nc+1:nb, kchd(k+jp) = Pn(kstr2, k+jp); end
    end
end
end

function [kchd] = mutat(kchd,mp,nb,n)  % mutation operation
for k = 1:nb*n
    if rand < mp
        kchd(k) = ~kchd(k);
    end
end
end

示例：求函数的最大值，使用 $x_l = -6 [1 1 1 1]$，$x_u = 6 [1 1 1 1]$，$nb = 8$，$ps = 40$，$ng = 50$，$mp = 0.02$。

xl=-6*[1 1 1 1]; xu=6*[1 1 1 1]; nb=8; ps=40; ng=50; mp=0.02;
f = @(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] = gaopt(f,xl,xu,nb,ps,ng,mp)

遗传算法流程图

graph TD;
    A[初始化种群 P(0)] --> B[评估初始解适应度];
    B --> C{t = 0};
    C --> D[使用遗传算子生成子代];
    D --> E[添加随机种群并评估适应度];
    E --> F[竞争选择，确定 P(t + 1)];
    F --> G{是否收敛};
    G -- 否 --> H[t = t + 1];
    H --> D;
    G -- 是 --> I[输出结果];

各优化算法的性能对比

算法	收敛速度	全局搜索能力	计算复杂度	对初始点的依赖
Rosenbrock 方法	中等	一般	较高	有一定依赖
Nelder 和 Mead 的单纯形法	较慢	一般	较低	有一定依赖
模拟退火方法	较慢	强	较高	依赖较小
遗传算法	中等	强	较高	依赖较小

总结

本文介绍了四种优化算法：Rosenbrock 方法、Nelder 和 Mead 的单纯形法、模拟退火方法和遗传算法。每种算法都有其独特的优点和适用场景。Rosenbrock 方法利用正交方向搜索，适用于多维优化问题；Nelder 和 Mead 的单纯形法简单易懂，适合函数梯度难以计算的情况；模拟退火方法能跳出局部最优，适用于复杂的全局优化问题；遗传算法模拟自然进化过程，具有较强的全局搜索能力。通过 MATLAB 代码实现了这些算法，并给出了具体的示例。在实际应用中，应根据问题的特点选择合适的优化算法。