基于改进鲸鱼算法的门式起重机主梁可靠度优化设计

门式起重机主梁可靠度优化设计，这个名字看着长，实际做起来最头疼的不是优化算法本身，而是"可靠度约束"怎么算、怎么嵌进优化循环。主梁要轻，就要把截面尺寸压下来；但载荷、材料强度、弹性模量全是随机量，截面尺寸一旦过小，失效概率就上去了。用改进鲸鱼算法PWSDWOA来搜最优截面，等于在"轻量化"和"可靠"之间跑优化，目标明确，但路不好走。

这篇复现文章，我尽量把整个链条写清楚：从主梁可靠度问题的数学模型，到原始WOA为什么需要改进，再到PWSDWOA的混沌初始化、自适应权重、差分变异三个策略，最后给出Matlab代码骨架和一组能直接参考的算例参数。适合正在做智能算法应用、结构可靠性优化方向的同学，也适合需要把材料力学问题转成优化模型的工程师。

1. 门式起重机主梁可靠度优化设计：问题本质与数学模型

1.1 主梁设计变量与目标函数

我先说清楚优化问题的"输入"是什么。门式起重机主梁普遍采用焊接工字形或箱形截面，我这里复现时以工字形焊接截面为算例。设计变量取四个基本尺寸：

• h：腹板高度
• tw：腹板厚度
• bf：翼缘宽度
• tf：翼缘厚度

目标函数是主梁单位长度的截面面积：

[
A(h, tw, bf, tf) = h \cdot tw + 2 \cdot bf \cdot tf
]

如果材料密度取 (\rho)，那么每延米重量就是 (\rho \cdot A)。优化目标取A最小，等价于主梁最轻。这个目标函数本身很简单，难点在于约束条件不是常规的"应力不超过许用值"，而是"强度可靠度达标"和"刚度可靠度达标"，所以整个问题叠加了两层计算：一层是优化搜索，一层是每次搜索都要执行的可靠度分析。

1.2 随机变量与功能函数

工程中的载荷、材料属性都有离散性。我按常见处理方式，把吊重轮压Q、材料屈服强度fy、弹性模量E都看作随机变量，并且假设服从正态分布。算例里我用的统计参数是：

• Q：均值按额定起重量对应的跨中轮压，取313.6 kN，变异系数0.1
• fy：均值按Q345钢取345 MPa，变异系数0.07
• E：均值取206 GPa，变异系数0.05

主梁跨中受集中轮压时，最大弯矩为 (M = QL / 4)，其中L是跨度。截面模量 (W = I / (h/2 + tf))，(I) 是主梁截面惯性矩。于是强度功能函数可以写成：

[
g_1 = fy - \frac{M}{W}
]

刚度功能函数以跨中挠度为衡量指标：

[
g_2 = [\delta] - \frac{QL^3}{48EI}
]

其中允许挠度 ([\delta]) 我按常规取 (L/700)。(g_1)、(g_2) 大于0说明安全，小于0说明失效。由于Q、fy、E是随机变量，(g_1)、(g_2) 自然也是随机变量，所以光看均值不够，要看它们落在失效区域的概率。

1.3 可靠指标与失效概率

可靠指标 (\beta) 和失效概率 (P_f) 之间的关系是：

[
P_f = \Phi(-\beta)
]

如果功能函数g服从正态分布，(\beta = \mu_g / \sigma_g)。但实际问题里g往往不是严格正态，所以工程上常用一次二阶矩方法，在"验算点"附近线性化处理。我用的是HL-RF迭代格式，也叫改进一次二阶矩法。大致思路：

1. 从随机变量均值点出发；
2. 计算当前点的功能函数值 (g(x)) 和梯度 (\nabla g(x))；
3. 根据梯度方向修正验算点位置；
4. 计算对应的 (\beta)，重复迭代直到 (\beta) 变化量小于收敛阈值。

这个计算过程在Matlab里写成循环非常快，但要注意一点：梯度不能全用数值差分，否则会慢到让人失去耐心。具体做法后面专门讲。

1.4 优化模型与约束条件

把可靠度要求写成约束，完整的优化模型长这样：

[
\min A(h, tw, bf, tf)
]

约束条件：

• (\beta_{strength} \ge \beta_{target})
• (\beta_{stiffness} \ge \beta_{target})
• 边界约束：(h \in [h_{min}, h_{max}])，(tw \in [tw_{min}, tw_{max}])，(bf \in [bf_{min}, bf_{max}])，(tf \in [tf_{min}, tf_{max}])
• 局部稳定约束，比如腹板高厚比和翼缘宽厚比不能超限

目标可靠指标 (\beta_{target}) 我取3.2，对应失效概率约 (6.9 \times 10^{-4})，对起重机械来说属于比较常规的要求。如果做更严格的设计，可以取3.7或者更高。

在实际用智能算法求解时，这种带可靠度约束的问题很难直接处理，我用的是外点罚函数法：把约束不满足的量乘以一个较大系数加到目标函数里。这样算法的搜索方向能自然倾向于满足可靠度约束的区域。

2. PWSDWOA改良思路：从原始WOA到混沌差分混合

2.1 原始WOA的三类捕食行为

鲸鱼算法WOA是Mirjalili在2016年提出的群智能算法，模仿座头鲸用气泡网捕食的行为。它的核心公式分三部分：

第一是包围猎物。当算法认为当前最优解是猎物位置时，其他个体向最优个体靠近：

[
D = |C \cdot X^*(t) - X(t)|
]

[
X(t+1) = X^*(t) - A \cdot D
]

其中A和C是系数，A随迭代次数从2线性降到0。第二是气泡网攻击，使用螺旋更新：

[
X(t+1) = D' \cdot e^{bl} \cdot \cos(2\pi l) + X^*(t)
]

第三是随机搜索。当|A|大于等于1时，个体不再跟随最优解，而是随机选一个同伴作为参考。这三类行为交替出现，理论上具备全局探索和局部开发的能力。

2.2 为什么原始WOA在可靠度优化里不够用

直接拿原始WOA跑主梁可靠度优化，实测下来有几个明显短板。

第一，初始种群用rand生成伪随机数，碰到底层约束强的工程问题时，大量初始个体落在不可行域，算法前几十代基本都在"往可行域里爬"，浪费迭代次数。第二，WOA后期A系数趋近0，种群会快速聚向当前最优个体。一旦这个最优解卡在局部，比如某组截面只满足了强度却没满足刚度，整个种群就被"带偏"了，很难跳出来。第三，可靠度函数内部嵌套迭代，目标函数整体是近似非光滑的，WOA在这种地形上容易反复震荡，收敛不稳定。

2.3 三个改进方向

PWSDWOA本质上是在原始WOA基础上做三处手术：

第一，初始化阶段用分段线性混沌映射PWLCM生成初始种群，替代伪随机数。混沌序列在[0,1]区间上遍历性更好，初始个体能更均匀地铺满解空间，从而减少前期"爬行"耗时。

第二，位置更新时引入自适应惯性权重w。迭代初期w设大一点，保留较强的全局探索能力；迭代后期w变小，让算法在局部做精细化搜索。这个思路和PSO的惯性权重很类似，但放到WOA里效果依然明显。

第三，引入差分进化DE的变异与选择机制。对部分个体执行 (V = X_{r1} + F \cdot (X_{r2} - X_{r3})) 的差分变异，再用贪心策略保留更优个体。这样每当种群趋于同质化时，就会有新个体注入，避免早熟。

所以在这套复现里，PWSDWOA的含义可以理解为：PWS（分段线性混沌映射）+ DWOA（融合差分进化的鲸鱼算法）。它是把WOA在"初始化、位置更新、种群多样性"三个环节同时做了改进。

3. PWSDWOA算法实现细节与Matlab代码骨架

3.1 PWLCM混沌初始化

分段线性混沌映射的迭代公式是：

[
x_{n+1} =
\begin{cases}
x_n / p, & 0 < x_n < p \
(1-x_n) / (1-p), & p \le x_n < 1
\end{cases}
]

其中p一般取0.4到0.5。不同p对序列分布有一定影响，实测取p=0.4时均匀性较好。生成N组、每组dim维的混沌序列后，再把值映射到设计变量的上下界：

[
x = lb + (ub - lb) \cdot seq
]

Matlab里的初始化代码可以这样写：

matlab复制function X = pwlcm_init(N, dim, lb, ub, p)
    % PWLCM混沌序列初始化种群
    if nargin < 5
        p = 0.4;
    end
    seq = zeros(N, dim);
    for d = 1:dim
        x = rand; % 混沌初值，避免取0或p附近
        for i = 1:N
            if x < p
                x = x / p;
            else
                x = (1 - x) / (1 - p);
            end
            seq(i, d) = x;
        end
    end
    % 映射到变量边界
    X = lb + (ub - lb) .* seq;
end

这里有几个细节容易踩坑：初值x不能正好落在0、p或1上，否则序列会退化；每维最好独立生成一条混沌序列，不要用同一序列转置填充，否则变量之间会引入虚假相关性。我一开始图省事用了一维序列扩维，后来发现初始种群在二维剖面上有明显对角聚集，换成逐维生成后分布正常了。

3.2 自适应惯性权重

我把自适应权重加在包围猎物的更新公式里，位置更新改成：

[
X(t+1) = w \cdot X^*(t) - A \cdot D
]

权重w按迭代次数非线性递减：

[
w = w_{min} + (w_{max} - w_{min}) \cdot e^{-\lambda (t/T)^2}
]

我常用的参数是 (w_{max}=0.9)，(w_{min}=0.4)，(\lambda=4)。这个权重的http作用在于：前期w大，个体受最优解"牵引"的力度强，能维持较大步长去探索陌生区域；后期w小，个体更依赖当前A、C带来的精细调节，收敛更平稳。

实现时就是在每次迭代开始前先更新w，然后在包围猎物的公式里把原来的 (X^) 替换成 (w \cdot X^)。注意螺旋更新和随机搜索这两支我没有加权重，不然会过度收缩搜索范围。

3.3 差分变异与选择

差分变异主要借鉴DE/rand/1策略：

[
V_i = X_{r1} + F \cdot (X_{r2} - X_{r3})
]

其中r1、r2、r3是从种群中随机挑选的三个不同个体，F是缩放因子，一般取0.5。并不是每个个体都做变异，我设置了一个概率CR，只有 rand < CR 才执行。变异后还要对越界的变量重新拉回边界内。

选择环节用贪心策略：如果变异后的V_i适应度比当前个体X_i好，就替换X_i；否则保留原个体。这里的"适应度"就是目标函数加罚函数后的值。

matlab复制% 差分变异与选择
function X_new = de_perturb(X, i, F, CR, lb, ub, obj_func)
    N = size(X, 1);
    if rand < CR
        idx = randperm(N, 3);
        % 保证不是当前个体
        while any(idx == i)
            idx = randperm(N, 3);
        end
        V = X(idx(1), :) + F * (X(idx(2), :) - X(idx(3), :));
        V = max(min(V, ub), lb);
        if obj_func(V) < obj_func(X(i, :))
            X_new = V;
        else
            X_new = X(i, :);
        end
    else
        X_new = X(i, :);
    end
end

注意randperm选三个索引时，如果种群很小，要避免r1、r2、r3和当前个体重合，否则变异方向会退化。实际算例中N取50，碰撞概率很低，但代码里还是加上while判断更稳妥。

3.4 PWSDWOA主循环

整个PWSDWOA的主循环不复杂，核心骨架如下：

matlab复制function [best_x, best_f] = pwsdwoa(obj_func, dim, lb, ub, N, MaxIter)
    % 参数设置
    p = 0.4;
    F = 0.5;
    CR = 0.5;
    w_max = 0.9;
    w_min = 0.4;
    lambda = 4;

    % PWLCM混沌初始化
    X = pwlcm_init(N, dim, lb, ub, p);
    f = zeros(N, 1);
    for i = 1:N
        f(i) = obj_func(X(i, :));
    end
    [best_f, idx_best] = min(f);
    best_x = X(idx_best, :);

    % 主迭代
    for t = 1:MaxIter
        a = 2 - 2 * t / MaxIter;
        w = w_min + (w_max - w_min) * exp(-lambda * (t / MaxIter)^2);

        for i = 1:N
            r = rand;
            A = 2 * a * r - a;
            C = 2 * r;
            p_rand = rand;

            if p_rand < 0.5
                if abs(A) < 1
                    % 包围猎物
                    D = abs(C .* best_x - X(i, :));
                    X_new = w .* best_x - A .* D;
                else
                    % 随机搜索
                    rand_idx = randi(N);
                    D = abs(C .* X(rand_idx, :) - X(i, :));
                    X_new = X(rand_idx, :) - A .* D;
                end
            else
                % 螺旋更新
                l = (rand - 0.5) * 2;
                D_best = abs(best_x - X(i, :));
                X_new = D_best .* exp(3 * l) .* cos(2 * pi * l) + best_x;
            end

            % 边界约束
            X_new = max(min(X_new, ub), lb);

            % 差分变异与贪心选择
            X_new = de_perturb([X; X_new], i, F, CR, lb, ub, obj_func);
            X(i, :) = X_new;

            % 立即更新该个体的适应度
            f(i) = obj_func(X(i, :));
        end

        % 更新全局最优
        [f_min, idx_min] = min(f);
        if f_min < best_f
            best_f = f_min;
            best_x = X(idx_min, :);
        end
    end
end

注意de_perturb函数里我把X_new临时拼到了X矩阵末尾，再从中提取变异参考个体，其实更好的做法是把当前个体跳过，直接在原种群选r1、r2、r3。上面代码重点是表达思路，直接照搬的话建议把随机索引冲突判断写完善。

4. 可靠度计算与适应度函数怎么搭

4.1 截面几何参数计算

主梁的截面面积、惯性矩、截面模量最好单独写成函数，这既是工程计算的公摊模块，也方便后续把工字形改成箱形截面时只替换这个文件。

matlab复制function [A, I, W] = cross_section_geom(h, tw, bf, tf)
    % 工字形截面几何参数
    A = h * tw + 2 * bf * tf;
    % 惯性矩：腹板矩形 + 两个翼缘矩形绕自身形心及平移项
    I_w = tw * h^3 / 12;
    I_f = 2 * (bf * tf^3 / 12 + bf * tf * (h/2 + tf/2)^2);
    I = I_w + I_f;
    % 截面模量，按最外缘纤维
    W = I / (h/2 + tf);
end

这里惯性矩的翼缘部分我保留了翼缘自身的惯性矩项 (bf \cdot tf^3 / 12)。虽然多数情况下这一项比平移项小两个数量级，但严谨起见保留没有坏处。注意如果主梁是箱形截面，这个函数要换成两根腹板和上下翼缘的组合，逻辑类似。

4.2 HL-RF法求可靠指标

可靠度计算是整个优化里最耗时的环节。我采用HL-RF迭代求可靠指标，但为了避免数值梯度拖慢速度，我直接推导了功能函数对随机变量的解析梯度。

以强度功能函数为例：

[
g_1 = fy - \frac{Q \cdot L}{4W}
]

其中只有fy和Q是随机变量，W是设计变量的确定性函数。梯度为：

[
\frac{\partial g_1}{\partial fy} = 1
]

[
\frac{\partial g_1}{\partial Q} = -\frac{L}{4W}
]

刚度功能函数：

[
g_2 = [\delta] - \frac{Q \cdot L^3}{48 E I}
]

梯度为：

[
\frac{\partial g_2}{\partial Q} = -\frac{L^3}{48 E I}
]

[
\frac{\partial g_2}{\partial E} = \frac{Q \cdot L^3}{48 E^2 I}
]

有了解析梯度，HL-RF迭代的代码可以写得很干净：

matlab复制function beta = hlrf_beta(mu, sigma, x0, g, grad_g, max_iter, tol)
    % mu, sigma: 随机变量均值与标准差
    % x0: 迭代初值，一般取mu
    % g: 功能函数句柄
    % grad_g: 梯度函数句柄，返回列向量
    x = x0;
    for k = 1:max_iter
        g_val = g(x);
        grad = grad_g(x);
        % 计算alpha方向
        alpha = -sigma .* grad / norm(sigma .* grad);
        % 计算当前可靠指标
        beta = (g_val + grad' * (mu - x)) / norm(sigma .* grad);
        % 更新验算点
        x_new = mu + beta * alpha .* sigma;
        if norm(x_new - x) < tol
            x = x_new;
            break;
        end
        x = x_new;
    end
end

注意这个迭代式的符号处理要小心。功能函数失效域定义为g<0，计算时要把梯度方向、alpha方向统一。我实际调试时遇到过beta为正但验算点反复振荡的情况，后来发现是alpha符号反了。如果你从其他资料里抄了公式，建议先拿一个一维问题验证收敛性。

4.3 适应度函数与罚函数

适应度函数要把目标值和约束罚项整合起来：

matlab复制function fitness = obj_fun(x)
    h = x(1); tw = x(2); bf = x(3); tf = x(4);

    % 几何参数
    [A, I, W] = cross_section_geom(h, tw, bf, tf);

    % 随机变量统计参数
    mu_Q = 313.6e3;  % N
    sigma_Q = mu_Q * 0.1;
    mu_fy = 345e6;   % Pa
    sigma_fy = mu_fy * 0.07;
    mu_E = 206e9;    % Pa
    sigma_E = mu_E * 0.05;

    L = 30;  % 跨度 m

    % 强度可靠指标
    g1 = @(r) r(2) - (r(1) * L) / (4 * W);  % r = [Q, fy]
    grad_g1 = @(r) [ -L / (4 * W); 1 ];
    beta_strength = hlrf_beta([mu_Q; mu_fy], [sigma_Q; sigma_fy], [mu_Q; mu_fy], g1, grad_g1, 100, 1e-6);

    % 刚度可靠指标
    g2 = @(r) (L / 700) - (r(1) * L^3) / (48 * r(2) * I); % r = [Q, E]
    grad_g2 = @(r) [ -L^3 / (48 * r(2) * I); (r(1) * L^3) / (48 * r(2)^2 * I) ];
    beta_stiffness = hlrf_beta([mu_Q; mu_E], [sigma_Q; sigma_E], [mu_Q; mu_E], g2, grad_g2, 100, 1e-6);

    beta_target = 3.2;

    % 罚函数
    penalty = 1e4;
    g_penalty = max(0, beta_target - beta_strength) + max(0, beta_target - beta_stiffness);

    fitness = A + penalty * g_penalty;
end

这里有几个单位细节必须统一。长度用m，力用N，应力用Pa，这样I的单位是m^4，W的单位是m^3。如果长度用mm，力用N，算出来的应力单位是MPa，虽然数值上也没错，但和材料强度单位保持一致才不会出低级错误。我建议全部用国际单位制，后面画图时如果需要mm再转。

罚系数penalty的取值直接影响优化结果。太小时算法会允许部分解不满足可靠度；太大时目标函数中截面积的影响被压缩，搜索会偏向"只要可靠度满足就行"，导致收敛慢。我试了几组，penalty取1e4到1e5之间比较合适。还有一种做法是随迭代次数逐步增大penalty，前期允许算法在不可行域探索，后期强制收紧约束，这也能用。

4.4 主程序整体调用结构

主程序main.m结构不复杂，关键是输出内容要完整，方便后处理：

matlab复制% 主脚本：PWSDWOA优化门式起重机主梁截面
clc; clear; close all;

dim = 4;
lb = [1.0; 0.006; 0.3; 0.008];   % h, tw, bf, tf 单位m
ub = [2.0; 0.016; 0.6; 0.024];

N = 50;
MaxIter = 200;

[best_x, best_f] = pwsdwoa(@obj_fun, dim, lb, ub, N, MaxIter);

fprintf('最优截面：h=%.3f m, tw=%.3f m, bf=%.3f m, tf=%.3f m\n', ...
    best_x(1), best_x(2), best_x(3), best_x(4));
fprintf('最小截面面积：%.4f m^2\n', best_f);

% 最终验算可靠度
[A_opt, I_opt, W_opt] = cross_section_geom(best_x(1), best_x(2), best_x(3), best_x(4));
% 再调用一次可靠度计算，输出beta

跑之前建议先在命令行直接调用obj_fun一个随机点，确保能正常返回值。接着再跑10代左右，观察适应度是否下降，确认通道畅通后再跑完整迭代。不要一上来就500代，出了问题排查太慢。

5. 常见问题、调参经验与结果验证

5.1 可靠度迭代不收敛

HL-RF迭代遇到非线性较强的功能函数时，偶尔会振荡甚至发散。我遇到最多的情况是随机变量变异系数取得过大，比如载荷变异系数超过0.2，验算点迭代就容易在几个点之间来回跳。解决思路有两个：

一个是给迭代步长加阻尼，把验算点更新公式改成：

[
x_{new} = x + \eta \cdot (x_{HLRF} - x)
]

(\eta)取0.5左右，能显著抑制振荡。代价是收敛变慢，但可靠度计算本身很快，可以接受。

另一个是用带均值修正的HL-RF变体，或者干脆改用蒙特卡洛抽样做校验。不过蒙特卡洛在优化循环里太慢，一般只用于最后验证结果。

5.2 PWSDWOA参数设置建议

种群N和迭代次数MaxIter不必一味求大。这个算例里设计变量只有4个，N取40到60足够；迭代次数150到300可以看到明显的收敛趋势。如果变量增加到8到12个，比如加入主梁跨度和腹板间距，N可以加到80到100，MaxIter相应增加到500。

差分变异参数F和CR对结果有直接影响。F太小，变异步长短，跳不出局部；F太大，变异个体随机性过强，优秀基因容易被冲散。我测试下来F=0.5、CR=0.5是一个比较稳妥的组合。如果发现收敛曲线波动剧烈，把CR降到0.3；如果发现算法停滞，把F提到0.8。

混沌映射参数p建议固定为0.4，不用来回调。如果你把p改成0.9，序列分布会偏向[1-p, 1]区间，初始种群在边界附近堆积，效果反而不如rand。

5.3 怎么判断结果合理

优化跑完不要急着采信，至少做三件事验证。

第一，用蒙特卡洛模拟对最优解做一次可靠度复核。在最优截面下，抽样Q、fy、E各10万次，统计g1和g2小于0的比例，再换算成失效概率和beta。如果MCS算出的beta和目标beta相差超过0.1，说明HL-RF计算或罚函数实现可能有问题。

第二，比较最优设计变量和边界约束的关系。如果某个变量卡在边界上，比如tf取到了上界，要思考是不是设计范围给小了，或者这个变量对可靠度的贡献必须靠增大截面来满足。边界解不是错误，但需要解释。

第三，画收敛曲线。我用原始WOA和PWSDWOA各跑一遍对比，PWSDWOA通常能在前30代就下降到接近最优的水平，而原始WOA到100代还在震荡。这类对比图在写研究报告或论文时也是很好的素材。

5.4 耗时优化技巧

可靠度优化最头疼的是计算时间。每次适应度调用都要做两组HL-RF迭代，每组迭代还要算功能函数和梯度，计算量比普通无约束优化大得多。两个技巧非常管用。

第一个是缓存优化后的结果。如果种群中有重复或接近重复的设计变量组合，直接返回已算好的适应度，避免重复计算可靠度。可以用一个简单的容器保存最近几百组解及其适应度，每组解先查表再计算。

第二个是用parfor并行计算种群个体的适应度。种群中每个个体的可靠度计算互相独立，天然适合并行。但注意Matlab的parfor要求内部函数都能正常序列化，obj_fun里不要依赖全局变量，所有参数都通过函数入参传递。跑多核的时候，我的经验是从4并行核起步，提升非常明显。

另外，HL-RF的容差tol也不要设得过小。优化过程中取1e-4完全够用，只有最终验证时才取1e-6。这样整个优化过程能快将近一倍。

最后再说一点调试心得：不要一上来就追求复杂的箱形截面、多随机变量、多目标函数。先拿工字形简化模型、两个功能函数跑通闭环，确认算法和可靠度求解模块都正常，再逐步增加复杂度。实际工程里，主梁可能还要处理局部屈曲、疲劳、动态挠度等约束，每加一个约束都要重新评估罚函数和变量边界，这比算法本身更考验耐心。