遗传算法应用

最新推荐文章于 2024-05-04 10:52:19 发布

若枫林涧

最新推荐文章于 2024-05-04 10:52:19 发布

阅读量442

点赞数 1

文章标签：算法机器学习 python

本文链接：https://blog.csdn.net/qq_51765896/article/details/132876372

版权

遗传算法

遗传算法简介

遗传算法是一种模拟自然选择过程的搜索算法。它是启发式算法的一种，常用于解决优化和搜索问题。

遗传算法的工作原理

初始化: 创建一个随机的解集，称为种群。
选择: 评估种群中每个解的适应度，并根据其适应度选择解。
交叉: 随机选择两个解，并组合它们以产生一个新的解。
变异: 以小的概率随机修改新的解。
评估: 评估新解的适应度。
替换: 用新解替换老的解。
终止: 如果满足某些条件（例如达到最大迭代次数或解不再变化），则停止算法。

遗传算法的应用

遗传算法在各种领域都有应用，如：

函数优化
机器学习
路径规划
游戏策略
…

下面我们来看一个非线性规划问题：

$\begin{align*} \text{max} \quad & f(x) = 2x_1 + 3x_1^2 + 3x_2 + x_2^2 + x_3 \\ \text{s.t.} \quad & \begin{cases} x_1 + 2x_1^2 + x_2 + 2x_2^2 + x_3 \leq 10 \\ x_1 + x_1^2 + x_2 + x_2^2 - x_3 \leq 50 \\ 2x_1 + x_1^2 + 2x_2 + x_3 \leq 40 \\ x_1^2 + x_3 = 2 \\ x_1 + 2x_2 \geq 1 \\ x_1 \geq 0 \\ x_2, x_3 \text{ 无约束} \end{cases} \end{align*}$
首先使用matlab的工具箱进行求解：

clear,clc
close all
% 定义目标函数  
% 目标函数为 -(2*x(1) + 3*x(1)^2 + 3*x(2) + x(2)^2 + x(3))  
% x 是一个向量，包含三个元素 x(1), x(2), x(3)  
fun = @(x) -(2*x(1) + 3*x(1)^2 + 3*x(2) + x(2)^2 + x(3));  
  
% 线性不等式约束  
% Aineq 是一个矩阵，定义了约束条件的系数  
% bineq 是一个向量，定义了约束条件的右侧值  
% 例如，第一个约束条件为 1*x(1) + 1*x(2) + 1*x(3) <= 10  
Aineq = [1 1 1;   
         1 1 -1;   
         2 2 1;   
         1 2 0];  
bineq = [10; 50; 40; 1];  
  
% 边界条件  
% lb 是一个向量，定义了每个变量的下界  
% ub 是一个向量，定义了每个变量的上界  
lb = [0, -inf, -inf];  
ub = [inf, inf, inf];  
  
% 使用fmincon解决问题  
% x0 是初始解，这里选择 [0, 0, 0]  
% options 定义了优化选项，这里选择显示迭代过程  
x0 = [0, 0, 0];  
options = optimoptions('fmincon','Display','iter');  
[x,fval] = fmincon(fun, x0, Aineq, bineq, [], [], lb, ub, @nonlin_constraints, options);  
  
% 显示结果  
disp('Solution using MATLAB:');  
disp(x);  
disp('Objective function value:');  
disp(-fval);  
  
% 非线性约束条件  
% 该函数定义了非线性约束条件，返回一个向量 c 和一个标量 ceq  
% c 是一个向量，定义了不等式约束条件  
% ceq 是一个标量，定义了等式约束条件  
function [c, ceq] = nonlin_constraints(x)  
    c(1) = x(1) + 2*x(1)^2 + x(2) + 2*x(2)^2 + x(3) - 10;  
    c(2) = x(1) + x(1)^2 + x(2) + x(2)^2 - x(3) - 50;  
    c(3) = 2*x(1) + x(1)^2 + 2*x(2) + x(3) - 40;  
    ceq(1) = x(1)^2 + x(3) - 2;  
end

求得结果：
Solution using MATLAB:
2.3333 -0.6667 -3.4444
Objective function value:
16.0000

我们再调用python的scipy库进行求解：

from scipy.optimize import minimize, NonlinearConstraint
import numpy as np

# 定义目标函数
def objective(x):
    return -(2*x[0] + 3*x[0]**2 + 3*x[1] + x[1]**2 + x[2])

# 定义非线性约束
def constraint_ineq(x):
    return [
        x[0] + 2*x[0]**2 + x[1] + 2*x[1]**2 + x[2] - 10,
        x[0] + x[0]**2 + x[1] + x[1]**2 - x[2] - 50,
        2*x[0] + x[0]**2 + 2*x[1] + x[2] - 40,
        x[0] + 2*x[1] - 1
    ]

def constraint_eq(x):
    return x[0]**2 + x[2] - 2

nonlinear_constraints = NonlinearConstraint(constraint_ineq, -np.inf, 0)
equality_constraint = NonlinearConstraint(constraint_eq, 0, 0)

# 边界条件
bounds = [(0, np.inf), (None, None), (None, None)]

# 初始猜测
x0 = [0, 0, 0]

# 使用scipy的minimize函数求解
solution_python_final = minimize(objective, x0, constraints=[nonlinear_constraints, equality_constraint], bounds=bounds, method='SLSQP')

print(f"x1 = {solution_python_final.x[0]:.2f}")
print(f"x2 = {solution_python_final.x[1]:.2f}")
print(f"x3 = {solution_python_final.x[2]:.2f}")
print(f"目标函数的最大值为：f(x) = {-solution_python_final.fun:.2f}")

所得结果：
x1 = 2.33
x2 = -0.67
x3 = -3.44
目标函数的最大值为：f(x) = 16.00

matlab和python调库得到的结果一致，接下来，我们使用遗传算法求解这个问题，调用matlab的遗传算法工具箱ga：

clc, clear
a=[-1 -2 0;-1 0 0];b=[-1;0];
[x,y]=ga(@ycfun1,3,a,b,[],[],[],[],@ycfun2);
x, y=-y

function y=ycfun1(x)  %x为行向量
c1=[2 3 1]; c2=[3 1 0];
y= c1* x' + c2* x'.^2; y=-y;
end

function [f,g]=ycfun2(x)
f=[x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10
   x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50
   2*x(1)+x(1)^2+2*x(2)+x(3)-40];
g=x(1)^2+x(3)-2;
end

所得结果：

x =[1.6987 1.0807 -0.8855]

y =15.5787

可见与优化工具箱求解出来的结果相近。

下面考虑一个NP-hard问题——旅行商问题

TSP的整数线性规划模型

旅行商问题（Traveling Salesman Problem, TSP）是经典的组合优化问题之一。给定一组城市和每对城市之间的距离，TSP的目标是找到一个最短的可能访问所有城市并返回起始城市的路径。

给定：

$n$ ：城市的数量。
$C$ ：城市的集合， $\{1,2,\ldots,n\}$ 。
$d_{ij}$ ：城市 $i$ 和城市 $j$ 之间的距离。
$x_{ij}$ ：决策变量，如果路径从城市 $i$ 到城市 $j$ ，则 $x_{ij} = 1$ ，否则 $x_{ij} = 0$ 。

目标函数：

$\min Z = \sum\limits _{i=1}^{n} \sum \limits _{j=1, j\neq i}^{n} d_{ij} x_{ij}$

约束条件：

每个城市都被访问一次（每个城市只有一个进入的边和一个离开的边）：

$\sum \limits _{i=1, i\neq j}^{n} x_{ij} = 1, \quad \forall j \in C$
$\sum \limits _{j=1, j\neq i}^{n} x_{ij} = 1, \quad \forall i \in C$

子旅程约束（为了消除不允许的子旅程）：

引入一个新的连续变量 $u_i$ 来表示城市 $i$ 在旅程中的顺序。然后添加以下约束：

$u_i - u_j + nx_{ij} \leq n-1, \quad \forall 1 \leq i \neq j \leq n$
$\leq u_i \leq n, \quad \forall 2 \leq i \leq n$

整数约束：

$x_{ij} \in \{0,1\}, \quad \forall i,j \in C, i \neq j$
$u_i \text{ 是整数}, \quad \forall i \in C$

接下来我们使用遗传算法来求解这个问题：

# 可视化版本
import numpy as np
import time
import matplotlib.pyplot as plt
import imageio
import io


# 城市坐标数据
cities = [
    (565, 575), (25, 185), (345, 750), (945, 685), (845, 655),
    (880, 660), (25, 230), (525, 1000), (580, 1175), (650, 1130),
    (1605, 620), (1220, 580), (1465, 200), (1530, 5), (845, 680),
    (725, 370), (145, 665), (415, 635), (510, 875), (560, 365),
    (300, 465), (520, 585), (480, 415), (835, 625), (975, 580),
    (1215, 245), (1320, 315), (1250, 400), (660, 180), (410, 250),
    (420, 555), (575, 665), (1150, 1160), (700, 580), (685, 595),
    (685, 610), (770, 610), (795, 645), (720, 635), (760, 650),
    (475, 960), (95, 260), (875, 920), (700, 500), (555, 815),
    (830, 485), (1170, 65), (830, 610), (605, 625), (595, 360),
    (1340, 725), (1740, 245)
]


# 计算所有城市之间的距离矩阵
def distance_matrix(cities):
    cities_array = np.array(cities)
    diff = cities_array[:, np.newaxis, :] - cities_array[np.newaxis, :, :]
    dist_matrix = np.linalg.norm(diff, axis=-1)
    return dist_matrix

# 评估路线的适应度（总距离）
def fitness(route, dist_matrix):
    shifted_route = np.roll(route, -1)
    return -np.sum(dist_matrix[route, shifted_route])

# 基于适应度概率选择两个父代
def select(population, dist_matrix):
    fit_values = np.array([fitness(route, dist_matrix) for route in population])
    fit_values -= np.max(fit_values)  # 为了数值稳定性
    probs = np.exp(fit_values) / np.sum(np.exp(fit_values))
    selected_indices = np.random.choice(len(population), size=2, p=probs)
    return [population[index] for index in selected_indices]

# 交叉：结合两个父代以产生一个孩子
def crossover(parent1, parent2):
    size = len(parent1)
    start, end = sorted(np.random.choice(np.arange(size), size=2, replace=False))
    child = [-1] * size
    for i in range(start, end + 1):
        child[i] = parent1[i]
    pointer = 0
    for i in range(size):
        if child[i] == -1:
            while parent2[pointer] in child:
                pointer += 1
            child[i] = parent2[pointer]
    return np.array(child)

# 变异：随机交换路线中的两个城市
def mutate(route):
    i, j = np.random.choice(np.arange(len(route)), size=2, replace=False)
    new_route = route.copy()
    new_route[i], new_route[j] = new_route[j], new_route[i]
    return new_route

# 绘制路径
def plot_route(cities, route, generation, frames):
    route_cities = [cities[i] for i in route]
    route_cities.append(route_cities[0])  # 添加第一个城市以闭合路径
    xs, ys = zip(*route_cities)
    
    plt.clf()  # 清除先前的图形
    plt.scatter(xs, ys, c='red', marker='o')  # 绘制城市的位置
    plt.plot(xs, ys, c='blue')  # 绘制路径
    plt.title(f"Generation: {generation+1}")
    plt.pause(0.1)  # 暂停以创建动画效果
    
    # 将当前帧保存到frames列表
    fig = plt.gcf()
    buf = io.BytesIO()
    fig.savefig(buf, format='png')
    buf.seek(0)
    frames.append(imageio.imread(buf))

# 遗传算法主函数
def genetic_algorithm(cities, pop_size, generations):
    num_cities = len(cities)
    dist_matrix = distance_matrix(cities)
    
    # 初始化种群
    population = [np.random.permutation(num_cities) for _ in range(pop_size)]
    
    plt.ion()  # 打开交互模式
    frames = []  # 保存每一帧的图像
    for generation in range(generations):
        new_population = []
        for i in range(pop_size):
            parents = select(population, dist_matrix)
            child = crossover(parents[0], parents[1])
            child = mutate(child)
            new_population.append(child)
        population = new_population
        
        # 绘制当前最佳路径
        best_route = max(population, key=lambda x: fitness(x, dist_matrix))
        plot_route(cities, best_route, generation, frames)
    
    plt.ioff()  # 关闭交互模式

    # 使用imageio将frames保存为GIF
    imageio.mimsave('route_evolution.gif', frames, duration=0.1)
    
    return best_route, -fitness(best_route, dist_matrix)

# 开始计时
start_time = time.time()

best_route, best_distance = genetic_algorithm(cities, 50, 500)

# 停止计时
end_time = time.time()
elapsed_time = end_time - start_time

print(f"Best Route: {best_route}")
print(f"Best Distance: {best_distance}")
print(f"Operation Hours: {elapsed_time:.2f} s")

运行结果：
Best Route: [17 2 16 20 41 6 1 29 22 19 49 28 12 13 51 10 11 27 26 25 46 15 45 24
3 50 32 42 14 23 5 4 47 37 36 39 38 43 33 34 35 48 44 18 40 7 8 9 31 0 21 30]
Best Distance: 9275.342356906101
Operation Hours: 228.96 s

在这里插入图片描述

接下来，我们调用python的deap库解决：

import numpy as np
import random
from deap import algorithms, base, creator, tools

# 城市坐标数据
cities = [
    (565, 575), (25, 185), (345, 750), (945, 685), (845, 655),
    (880, 660), (25, 230), (525, 1000), (580, 1175), (650, 1130),
    (1605, 620), (1220, 580), (1465, 200), (1530, 5), (845, 680),
    (725, 370), (145, 665), (415, 635), (510, 875), (560, 365),
    (300, 465), (520, 585), (480, 415), (835, 625), (975, 580),
    (1215, 245), (1320, 315), (1250, 400), (660, 180), (410, 250),
    (420, 555), (575, 665), (1150, 1160), (700, 580), (685, 595),
    (685, 610), (770, 610), (795, 645), (720, 635), (760, 650),
    (475, 960), (95, 260), (875, 920), (700, 500), (555, 815),
    (830, 485), (1170, 65), (830, 610), (605, 625), (595, 360),
    (1340, 725), (1740, 245)
]

# 计算城市间的距离矩阵
dist_matrix = [[np.linalg.norm(np.array(city1) - np.array(city2))
                for city2 in cities] for city1 in cities]

# 计算路径总距离
def evaluate(individual):
    return sum([dist_matrix[individual[i-1]][individual[i]] for i in range(len(individual))]),


# 设置优化目标：最小化路径距离
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", list, fitness=creator.FitnessMin)

toolbox = base.Toolbox()
toolbox.register("indices", random.sample, range(len(cities)), len(cities))
toolbox.register("individual", tools.initIterate,
                 creator.Individual, toolbox.indices)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

toolbox.register("mate", tools.cxOrdered)
toolbox.register("mutate", tools.mutShuffleIndexes, indpb=0.05)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("evaluate", evaluate)

population = toolbox.population(n=300)

# 算法参数
ngen = 1000
crossover_prob = 0.75
mutation_prob = 0.2

# 执行遗传算法
algorithms.eaSimple(population, toolbox, crossover_prob, mutation_prob, ngen)

# 输出结果
top_individual = tools.selBest(population, k=1)[0]
print('Optimal Route:', top_individual)
print('Optimal Distance:', evaluate(top_individual)[0])

运行结果：
Optimal Route: [0, 43, 45, 36, 33, 34, 35, 38, 39, 37, 47, 23, 4, 14, 5, 3, 24, 11, 27, 26, 25, 46, 12, 13,
51, 10, 50, 32, 42, 9, 8, 7, 40, 18, 44, 2, 16, 20, 41, 6, 1, 29, 28, 15, 49, 19, 22, 30, 17, 21, 31, 48]
Optimal Distance: 7797.455992319521

可见调库所得结果要好得多，并且运行时间短了数倍。

原因在于我们自己写的遗传算法没有仔细调参数，而库函数的参数已经是优化的比较好的了。当然自己写的算法也有好处，那就是自由度高，比如我想加上可视化操作就比较容易。