49、基于遗传算法的优化方法及MATLAB实现

算法笑匠

于 2025-11-21 11:04:43 发布

阅读量9

点赞数

CC 4.0 BY-SA版权

分类专栏： MATLAB工程应用精解文章标签：遗传算法 MATLAB 优化方法

本文链接：https://blog.youkuaiyun.com/1a2s3d4f5g/article/details/155214405

MATLAB工程应用精解专栏收录该内容

53 篇文章 ¥499.90

订阅专栏¥69.90

会员秒杀 ¥9.9 重磅福利

超级会员免费看

基于遗传算法的优化方法及MATLAB实现

1. 引言

在优化问题中，传统方法往往在处理复杂的目标函数和约束条件时存在一定的局限性。而遗传算法（Genetic Algorithm，GA）作为一种基于生物进化原理的优化方法，为解决这些复杂问题提供了新的思路。它适用于目标函数和/或约束函数高度非线性、不可微或不连续的优化问题，并且变量可以是连续的、离散的或两者的组合。

2. 遗传算法基础

2.1 遗传算法原理

遗传算法基于生物进化和适者生存的原则。与传统方法从单个点开始并在每次迭代中生成单个点不同，遗传算法从一组随机点（称为种群）开始，并在算法的每次迭代（代）中生成一组随机点（称为个体）。遗传算法使用三种算子来改进种群并收敛到最终的解种群，分别是选择、交叉和变异：
- 选择算子 ：在种群中识别父代。
- 交叉算子 ：应用于随机选择并匹配的两个父代，以产生子代。
- 变异算子 ：应用于个体，以在种群中产生多样性。

这些操作在每一代中执行，当种群无法再改进或达到最大迭代次数时，算法停止。由于遗传算法是一种随机优化方法，每次运行可能会产生略有不同的解。

2.2 MATLAB中的遗传算法函数

MATLAB提供了两个基于遗传算法的函数： ga 和 gamultiobj 。
- ga 函数 ：用于解决单目标优化问题，形式如下：

minimize
x
f(x)
subject to:
lb ≤ x ≤ ub
Ceq(x) = 0    (nonlinear equality constraints)
C(x) ≤ 0   (nonlinear inequality constraints)
Aeqx = beq (linear equality constraints)
Ax ≤ b   (linear inequality constraints)

基本命令为：

[x, f] = ga(@UserFunction, nvars, A, b, Aeq, beq, lb, ub, @NonLinConstr, options)

其中：
- UserFunction ：计算目标函数的函数名。
- nvars ：变量的数量。
- A 和 b ：线性不等式约束的系数。
- Aeq 和 beq ：线性等式约束的系数。
- lb 和 ub ：变量的下界和上界向量。
- NonLinConstr ：定义非线性不等式和等式约束的函数。
- options ：设置 gaoptimset 中描述的参数。

gamultiobj 函数 ：用于解决多目标优化问题，形式如下：

minimize
x
{f1, f2, …, fm}
subject to:
lb ≤ x ≤ ub
Aeqx = beq      (linear equality constraints)
Ax ≤ b           (linear inequality constraints)

MATLAB函数为：

[x, f] = gamultiobj(@UserFunction, nvars, A, b, Aeq, beq, lb, ub, options)

其中：
- x ：获得的Pareto解。
- f ：对应的目标函数值。
- 其他参数与 ga 函数类似。

2.3 约束处理

在 ga 和 gamultiobj 函数中，变量可以是连续的或二进制的（0或1）。但在当前版本的MATLAB中，当 'PopulationType' 选项设置为 'bitstring' （二进制）时，线性和非线性约束可能不满足。为了解决这个问题，可以将约束以不可行性约束的形式添加到目标函数中。即当约束被违反时，将一个大的惩罚项添加到目标函数中，否则惩罚项设置为零。惩罚方法的表达式为：

minimize
x
f + P[max(0, (Ax - b)) +  max(0, C(x)) + (Aeqx - beq)2 + Ceq
2 ]

其中 P 是一个大的惩罚参数。

3. 示例分析

3.1 示例1：背包问题（单目标，二进制变量）

考虑一个背包问题，通过最小化以下目标函数来确定七个二进制变量：

minimize f(x) = -(12x1 + 12x2 + 9x3 + 15x4 + 90x5 + 26x6 + 112x7) + P max[0, (3x1 + 4x2 + 3x3 + 3x4 + 15x5 + 13x6 + 16x7 - 35)]

其中 P 设置为100。MATLAB代码如下：

function Example13_20
    nvars = 7;
    options = gaoptimset;
    options = gaoptimset(options, 'PopulationType', 'bitstring');
    options = gaoptimset(options, 'PopulationSize', 100);
    options = gaoptimset(options, 'PlotFcns', {@gaplotbestf, @gaplotbestindiv});
    [x, f] = ga(@KnapsackObj, nvars, [], [], [], [], [], [], [], options);
    disp([repmat('x(', 7,1) num2str((1:7)') repmat(') = ', 7,1) int2str(x')])
    disp(['Total value of placed objects = $' num2str(-f, '%6.2f')])
end

function f = KnapsackObj(x)
    P = 100;
    c = 3*x(1)+4*x(2)+3*x(3)+3*x(4)+15*x(5)+13*x(6)+16*x(7)-35;
    f = -(12*x(1)+12*x(2)+9*x(3)+15*x(4)+90*x(5)+26*x(6)+112*x(7))+P*max(0,c);
end

运行该程序，得到的结果与之前的示例一致。

3.2 示例2：两杆桁架问题（单目标，连续变量）

考虑两杆桁架问题，通过 ga 函数求解，同时对变量进行缩放以提高数值精度。MATLAB代码如下：

function Example13_21
    global sigma
    nvars = 3;
    lb = [100, 0, 0];
    ub = [300, 20, 20];
    sigma = 1e5;
    options = gaoptimset;
    options = gaoptimset(options, 'PopulationSize', 60);
    options = gaoptimset(options, 'Generations', 200);
    options = gaoptimset(options, 'Display', 'final');
    options = gaoptimset(options, 'PlotFcns', {@gaplotbestindiv, @gaplotscorediversity});
    [x, V] = ga(@TwoBarTrussObj, nvars, [], [], [], [], lb, ub, @TwoBarTrussNonLinCon, options);
    disp(['y = ' num2str(x(1)) ' cm x1 = ' num2str(x(2)) ' cm^2 x2= ' num2str(x(3)) ' cm^2 Volume = ' num2str(V*1e6) ' cm^3'])
end

function f = TwoBarTrussObj(x)
    x(1) = x(1)/1e2;
    x(2:3) = x(2:3)/1e4;
    % Values scaled
    f = x(2)*sqrt(16+x(1)^2)+x(3)*sqrt(1+x(1)^2);
end

function [C, Ceq] = TwoBarTrussNonLinCon(x)
    global sigma
    x(1) = x(1)/1e2;
    x(2:3) = x(2:3)/1e4;
    % Values scaled
    C(1) = 20*sqrt(16+x(1)^2)-sigma*x(1)*x(2);
    C(2) = 80*sqrt(1+x(1)^2)-sigma*x(1)*x(3);
    Ceq = [];
end

运行该程序，结果与之前的示例相当吻合。由于遗传算法的随机性，每次运行 ga 可能会产生略有不同的结果，因此通常建议多次运行并选择目标函数值最小的解。

3.3 示例3：两杆桁架问题（多目标，连续变量）

将两杆桁架问题转换为多目标形式，使用 gamultiobj 函数求解。目标函数如下：

minimize {f1, f2}
subject to:
0.05 ≤ (x1, x2) ≤ 0.25
1 ≤ y ≤ 3
80√(1 + y²)/(σyx2) - 1 ≤ 0
f2 = 20√(16 + y²)/(σyx1) - 1
where f1 = x1√(16 + y²) + x2√(1 + y²)

MATLAB代码如下：

function Example13_22
    global P sigma
    P = 100;
    sigma = 1e5;
    nvars = 3;
    lb = [0.05, 0.05, 1];
    ub = [0.25, 0.25, 3];
    options = gaoptimset;
    options = gaoptimset(options, 'Vectorized', 'on');
    [x, V] = gamultiobj(@TwoBarTrussTwoObj, nvars, [], [], [], [], lb, ub, options);
    [f , C] = TwoBarTrussTwoObj(x);
    disp(C)
    plot(V(:,1),V(:,2), 'k*')
    xlabel('f_1')
    ylabel('f_2')
    grid on
end

function [f, C] = TwoBarTrussTwoObj(x)
    global P sigma
    f1(:,1) = x(:,1).*sqrt(16+x(:,3).^2)+x(:,2).*sqrt(1+x(:,3).^2);
    f1(:,2) = 20*sqrt(16+x(:,3).^2)./(sigma*x(:,1).*x(:,3))-1;
    C (:,1) = 80*sqrt(1+x(:,3).^2)./(sigma*x(:,2).*x(:,3))-1;
    f(:,1) = f1(:,1)+P*max(0,C(:,1));
    f(:,2) = f1(:,2)+P*max(0,C(:,1));
end

运行该程序，得到Pareto最优解，并绘制Pareto前沿图。需要验证约束条件是否满足，在本示例中， C < 0 。

3.4 示例4：两杆桁架问题（单目标，连续和离散变量）

假设两杆桁架问题具有混合连续 - 离散变量，变量 x1 和 x2 为连续变量，变量 y 为离散变量，取值为 {1, 1.3, 1.6, 1.9, 2.2, 2.5, 2.8} 。目标函数为：

minimize 
f = x1√(16 + y²) + x2√(1 + y²)
subject to:
0 ≤ (x1, x2) ≤ 0.25
80√(1 + y²)/(σyx2) - 1 ≤ 0
20√(16 + y²)/(σyx1) - 1 ≤ 0

为了处理连续变量，将其通过二进制变量表示。对于连续变量 x ，定义在闭区间 [lb, ub] 内，可以表示为：

x = lb + (ub - lb)/(2^N - 1) * D(s)
where D(s) = ∑(i = 0 to N - 1) 2^i * si, si = 0 or 1

其中 N 是所需的最小二进制变量数量，通过以下公式计算：

N = 1 + floor(d + log10(ub - lb) / log10(2))

d 是所需的有效小数位数。

MATLAB代码如下：

function Example13_23
    nvars = 47;
    PopulationSize = 100;
    Generations = 100;
    options = gaoptimset;
    options = gaoptimset(options, 'PopulationType', 'bitstring');
    options = gaoptimset(options, 'PopulationSize', PopulationSize);
    options = gaoptimset(options, 'Generations', Generations);
    options = gaoptimset(options, 'PlotFcns', {@gaplotbestindiv, @gaplotscorediversity});
    [x, fval] = ga(@TwoBarTrussTwoObjGa, nvars, [], [], [], [], [], [], [], options);
    finalx1 = GetRealFromBinary(x, 1, 22, 0, 0.25);
    finalx2 = GetRealFromBinary(x, 23, 44, 0, 0.25);
    finaly = GetDiscreteFromBinary(x, 45, 47);
    disp(['x1 = ' num2str(finalx1*1e4) ' cm^2 x2 = ' num2str(finalx2*1e4) ' cm^2 y = ' num2str(finaly*100) ' cm'])
    disp(['V = ' num2str(fval*1e6)])
end

function f = TwoBarTrussTwoObjGa(x)
    P = 100;
    sigma = 1e5;
    var(1) = GetRealFromBinary(x, 1, 22, 0, 0.25);
    var(2) = GetRealFromBinary(x, 23, 44, 0, 0.25);
    var(3) = GetDiscreteFromBinary(x, 45, 47);
    f1 = var(1)*sqrt(16+var(3)^2)+var(2)*sqrt(1+var(3)^2);
    C1 = 80*sqrt(1+var(3)^2)/(sigma*var(2)*var(3))-1;
    C2 = 20*sqrt(16+var(3)^2)/(sigma*var(1)*var(3))-1;
    f = f1+P*(max(0,C1)+max(0,C2));
end

function x = GetRealFromBinary(y, startbit, endbit, lowerlimit, upperlimit)
    temp = y(startbit:endbit);
    ii = length(temp);
    newint = 0;
    for i = 1:ii
        newint = newint + temp(i)*2^(i - 1);
    end
    maxbitvalue = 2^(length(temp)) - 1;
    x = lowerlimit + (upperlimit - lowerlimit)/maxbitvalue * newint;
end

function x = GetDiscreteFromBinary(y, startbit, endbit)
    choice = [1, 1.3, 1.6, 1.9, 2.2, 2.5, 2.8];
    temp = y(startbit:endbit);
    ii = length(temp);
    newint = 0;
    for i = 1:ii
        newint = newint + temp(i)*2^(i - 1);
    end
    if newint == 0
        x = choice(1);
    else
        x = choice(newint);
    end
end

运行该程序，得到优化结果。由于变量 y 被限制为离散值，该解可能比之前的示例稍差。

4. MATLAB优化工具箱函数总结

MATLAB函数	描述
`bintprog`	二进制整数规划求解器
`fgoalattain`	多目标目标达成求解器
`fminbnd`	固定区间上单变量函数的最小值
`fmincon`	约束非线性多变量函数的最小值
`fminimax`	极小极大求解器
`fminsearch`	无约束多变量函数的最小值
`fminunc`	无约束多变量函数的最小值
`fseminf`	半无限约束多变量非线性函数的最小值
`linprog`	线性规划求解器
`lsqcurvefit`	非线性最小二乘曲线拟合求解器
`lsqnonlin`	非线性最小二乘求解器
`quadprog`	二次规划求解器
`ga`	实现遗传算法，求解单目标函数的最小值
`gamultiobj`	实现遗传算法，求解多目标函数的最小值

5. 遗传算法优化流程

graph TD;
    A[初始化种群] --> B[评估适应度];
    B --> C{是否满足停止条件};
    C -- 是 --> D[输出最优解];
    C -- 否 --> E[选择操作];
    E --> F[交叉操作];
    F --> G[变异操作];
    G --> B;

通过以上内容，我们可以看到遗传算法在解决复杂优化问题方面的强大能力，以及MATLAB提供的相关函数为我们实现这些算法提供了便利。在实际应用中，可以根据具体问题选择合适的优化方法和函数。

6. 优化问题练习案例

6.1 两杆机构平衡问题

考虑一个两杆机构，受到外部力作用，通过最小化势能函数来获得机构的平衡状态。势能函数为：

minimize PE = -P1*L1*sin(α) - P2*L2*sin(β) + 2*P3*[L1*(1 - cos(α)) + L2*(1 - cos(β))]

操作步骤如下：
1. 创建势能函数的等高线和曲面图。
2. 使用 fminunc 函数，以初始点 (α, β) = (1, 1) 进行求解，得到 α = 0.4636 rad ， β = 0.1651 rad ， PE = -0.3789 。

6.2 水道设计优化问题

对于水道的设计，目标是在固定横截面积的情况下，最大化流量。流量与湿周的倒数成正比，湿周 p = c + 2*h/sin(θ) ，横截面积 A = c*h + h^2*cot(θ) 。要最小化的函数为：

minimize 1/p = (A/(h - h*cot(θ) + 2*h/sin(θ)))^(-1)

操作步骤如下：
1. 创建目标函数的等高线和曲面图。
2. 使用 fminsearch 函数，以初始点 (h, θ) = (1, 1) 进行求解，得到 h = 2.3172 ， θ = 60° 。

6.3 两弹簧系统平衡问题

两弹簧系统在受力后达到平衡，势能函数为：

PE = 0.5*k1*x1^2 + 0.5*k2*(x2 - x1)^2 - F*x2

操作步骤如下：
1. 创建势能函数的等高线和曲面图。
2. 使用 fminsearch 函数，以初始点 (x1, x2) = (1, 1) 进行求解，得到 x1 = 4.1289 ， x2 = 0 ， PE = -15.2802 。

6.4 加工操作优化问题

加工操作中，平均总生产时间 T = tm + tm*tc/Tl + taux ，其中 tm 为切削时间， Tl 为刀具寿命， tc 为换刀时间， taux 为辅助时间。切削时间 tm = π*D*L/(1000*V*f) ，刀具寿命受切削速度 V 、进给率 f 和切削深度 d 影响。
操作步骤如下：
1. 创建系统的等高线和曲面图。
2. 使用 fmincon 函数，以初始点 (V, f) = (10, 1) 进行求解，得到 T = 3.942 ， V = 100.46 ， f = 2.000 。

6.5 弹簧 - 质量系统平衡问题

对于弹簧 - 质量系统，通过最小化势能函数 PE = 0.5*X'*K*X - X'*P 来求解平衡位置，其中 K 为刚度矩阵， P 为外力向量。
操作步骤如下：
使用 fminunc 函数，以初始点 (x1, x2) = (1, 1) 进行求解，得到 x1 = 4.7059 ， x2 = 10.5882 ， PE = 529.41 。

6.6 弹簧 - 重量系统平衡问题

弹簧 - 重量系统中，通过最小化势能函数来求解平衡位置。势能函数为：

PE = 0.5*∑(i = 1 to 6) Ki*(ΔLi)^2 + ∑(j = 1 to 5) Wj*Yj+1

操作步骤如下：
使用 fminunc 函数，以初始点 (X2, X3, X4, X5, X6, Y2, Y3, Y4, Y5, Y6) = (7.5, 15, 22.5, 30, 37.5, 0, 0, 0, 0, 0) 进行求解。

6.7 管状柱屈曲问题

管状柱的屈曲载荷 F = π*a*E*R^b*t^(4 - b)/(4*L^2) ，通过曲线拟合实验数据来确定未知常数 a 和 b 。
操作步骤如下：
使用 lsqcurvefit 函数，以初始点 (a, b) = (1, 1) 进行求解，得到 a = -3.3182 ， b = 3.1588 。

6.8 药物成分含量问题

药物中某成分的含量随时间变化，假设非线性模型为 Y = a + (0.51 - a)*exp(-b*(t - 8)) 。
操作步骤如下：
使用 lsqcurvefit 函数，以初始点 (a, b) = (1, 1) 进行求解，得到 a = 0.3918 ， b = 0.1394 。

6.9 车辆运动问题

车辆的加速度 a(t) = C*t^2 + D*t + 2*a0 ，位置 x(t) = A*t^4 + B*t^3 + a0*t^2 。
操作步骤如下：
使用 lsqnonlin 函数，以初始点 (A, B) = (0, 0) 进行求解，得到车辆速度关于时间的方程 v(t) = -0.032693*t^3 + 0.2515*t^2 + 1.22*t ，误差为 5.6533 。

6.10 圆孔定位问题

对于圆孔，通过 rn = sqrt((xc - xn)^2 + (yc - yn)^2) 来确定圆心位置 (xc, yc) 和半径 r 。
操作步骤如下：
使用模拟数据，通过MATLAB表达式 r = ro + 0.02*randn(N, 1) 生成数据，其中 ro = 0.075 ， N = 40 ， xc = 1.25 ， yc = 1.5 ，求解得到圆心和半径。

6.11 最小球体问题

找到包含四个给定点的最小球体的圆心和半径。球体半径公式为 R^2 = (xc - x)^2 + (yc - y)^2 + (zc - z)^2 。
操作步骤如下：
使用 fmincon 函数，以初始点 (xc, yc, zc, R) = (4, 4, 4, 4) 进行求解，得到 (xc, yc, zc, R) = (-0.5, -0.5, 2.5, 4.5552) 。

6.12 三杆桁架问题

三杆桁架的目标函数为其加载节点的垂直挠度：

minimize f = P*h/(E*(1/(x1 + 12*x2)))

约束条件包括应力约束和设计变量的上下界。
操作步骤如下：
使用 fmincon 函数，以初始点 (x1, x2) = (0, 0) 进行求解，得到 x1 = 8.3177 ， x2 = 26.924 ， f = 0.021554 。

6.13 电阻桥网络问题

电阻桥网络的总功率 P = ∑(n = 1 to 5) In^2*Rn ，约束条件为电流范围 1 ≤ In ≤ 2 。
操作步骤如下：
使用 fmincon 函数，以初始点 (R1, R2, R3, R4, R5) = (1, 1, 1, 1, 1) 进行求解，得到 [R1, R2, R3, R4, R5] = [3, 1, 1, 2, 0.5] 欧姆，总功率耗散为 12 W 。

6.14 生产运输优化问题

公司有 m 个生产设施和 n 个仓库，目标是最小化运输成本：

minimize ∑(i = 1 to m) ∑(j = 1 to n) (cij*xij + dij*xij^2)

约束条件包括生产能力和仓库需求。
操作步骤如下：
使用 quadprog 函数，以给定的初始点进行求解，得到最优的生产和运输数量。

6.15 两杆桁架多目标问题

两杆桁架的多目标问题，目标函数为重量和应力相关的函数，约束条件包括应力范围和设计变量的上下界。
操作步骤如下：
使用 fseminf 函数，以初始点 (x1, x2) = (0.1, 0.1) 进行求解，得到 x1 = 0.8088 ， x2 = 0.45058 ， w = 12.6431 。

7. 总结

在实际应用中，需要根据问题的特点，如目标函数的性质（线性、非线性、单目标、多目标等）、约束条件的类型（等式约束、不等式约束等）以及变量的类型（连续、离散、二进制等），选择合适的优化方法和函数。同时，由于遗传算法等方法具有随机性，多次运行并选择最优结果是一个不错的策略。

graph LR;
    A[问题定义] --> B{问题类型};
    B -- 单目标无约束 --> C[fminsearch/fminunc];
    B -- 单目标约束 --> D[fmincon/ga];
    B -- 多目标 --> E[fgoalattain/gamultiobj];
    B -- 线性规划 --> F[linprog];
    B -- 二次规划 --> G[quadprog];
    B -- 曲线拟合 --> H[lsqcurvefit/lsqnonlin];

这个流程图展示了根据问题类型选择合适MATLAB优化函数的过程。希望这些内容能帮助大家在遇到优化问题时，能够更加从容地选择和应用合适的方法。