24年全国数模B题第三问优秀论文（B196）复现---python实现---附代码

本文链接：https://blog.csdn.net/ABCabc0286/article/details/144545100

（1）求解思路：

（1）决策树构建：确定了32决策节点，8个零件，3个半成品，1个成品的各种策略

反馈循环机制剖析：存在五个反馈循环机制，包括三个半成品不合格品拆解为待检测零件的循环、一个成品不合格品拆解为待检测半成品的循环以及一个客户退回产品重新归类为成品次品的循环，这些循环显著增加了计算复杂度。

（2）确定成本利润的计算公式，然后模拟一下。

（3）基于蚁群算法进行优化（引入精英）：

编码决策：将 32 个决策节点编码为 16 位二进制的染色体形式，每位对应特定阶段决策。

初始化种群：随机生成一定规模个体组成初始种群，每个个体代表一种生产决策组合。

适应度计算：针对种群中每个个体，算出利润，作为个体适应度值。

选择操作：采用轮盘赌选择方式，依据个体适应度计算被选中概率，（适应度高，利润大，选中的概率也就越大）按概率选择父代个体。

（引入精英策略进一步优化，每一次迭代完之后强制保存适应度最高的前10%的个体）

交叉操作：对父代个体进行单点交叉（随机更换决策），随机确定交叉点产生子代个体，模拟基因重组，保留优良基因特征并增加种群多样性。

变异操作：按设定变异概率对子代个体基因位随机变异（也是更换决策选择），探索更多决策空间，避免过早陷入局部最优。

迭代优化：重复上述 3 - 6 步，多代迭代进化，使种群向利润更高方向发展，最终选最后一代种群中适应度最高个体对应的决策组合为最优生产决策方案。

（2）决策树绘制：

from graphviz import Digraph

# 创建一个有向图
dot = Digraph(comment='Production Decision Flowchart', encoding='utf-8')
dot.node_attr['fontname'] = 'SimHei'  # 设置节点字体为SimHei
dot.edge_attr['fontname'] = 'SimHei'  # 设置边的字体为SimHei

# 定义零件节点
parts = ['零件{}'.format(i + 1) for i in range(8)]
semi_products = ['半成品1', '半成品2', '半成品3']
finished_product = '成品'
user = '用户'

# 增加表示拆解后的零件节点
disassembled_parts = ['拆解零件']

# 添加零件节点
for part in parts:
    dot.node(part, shape='oval')

# 添加半成品节点
for semi in semi_products:
    dot.node(semi, shape='box')

# 添加成品和用户节点
dot.node(finished_product, shape='box')
dot.node(user, shape='oval')

# 添加拆解后的零件节点
for dis_part in disassembled_parts:
    dot.node(dis_part, shape='oval')

# 添加决策节点和箭头（边）
# 零件到半成品的决策
for i in range(len(parts)):
    # 检测 -> 合格路径，设置黑色箭头和黑色文字
    dot.edge(parts[i], semi_products[0 if i < 3 else (1 if i < 6 else 2)], label='检测 -> 合格', color='black', fontcolor='black')
    # 不检测路径，设置红色箭头和红色文字
    dot.edge(parts[i], semi_products[0 if i < 3 else (1 if i < 6 else 2)], label='不检测', color='red', fontcolor='red')

# 半成品之间的决策，按照顺序流转，前一个半成品指向后一个半成品
for i in range(len(semi_products) - 1):
    dot.edge(semi_products[i], semi_products[i + 1], label='检测 -> 合格', color='black', fontcolor='black')
    dot.edge(semi_products[i], semi_products[i + 1], label='不检测', color='red', fontcolor='red')
# 最后一个半成品指向成品
dot.edge(semi_products[-1], finished_product, label='检测 -> 合格', color='black', fontcolor='black')
dot.edge(semi_products[-1], finished_product, label='不检测', color='red', fontcolor='red')

# 成品到用户的决策
dot.edge(finished_product, user, label='合格', color='black', fontcolor='black')

# 处理成品不合格情况，增加是否拆解判断
dot.edge(finished_product, disassembled_parts[0], label='不合格 -> 拆解 -> 检测', color='red', fontcolor='red')
dot.edge(finished_product, disassembled_parts[0], label='不合格 -> 直接报废', color='gray', fontcolor='gray')

# 处理半成品不合格情况，增加是否拆解判断
for i, semi in enumerate(semi_products):
    dot.edge(semi, disassembled_parts[0], label='不合格 -> 拆解 -> 检测', color='red', fontcolor='red')
    # 不合格但不拆解时直接指向成品
    dot.edge(semi, finished_product, label='不合格 -> 继续流转（不拆解）', color='orange', fontcolor='orange')

# 用户反馈不合格产品返回后的处理
dot.edge(user, finished_product, label='不合格（退换产品） -> 检测', color='red', fontcolor='red')

# 渲染并保存图像到images文件夹
dot.render('images/production_decision_flowchart', view=True)

（3）决策数直接求解：

import graphviz
from itertools import product
import matplotlib.pyplot as plt
from matplotlib import rcParams
import os

# 设置中文字体
rcParams['font.sans-serif'] = ['SimHei']  # 用于显示中文标签
rcParams['axes.unicode_minus'] = False  # 用于正常显示负号

# 定义各零配件、半成品、成品相关参数
part_info = {
    1: {'次品率': 0.1, '购买单价': 2, '检测成本': 1},
    2: {'次品率': 0.1, '购买单价': 8, '检测成本': 1},
    3: {'次品率': 0.1, '购买单价': 12, '检测成本': 2},
    4: {'次品率': 0.1, '购买单价': 2, '检测成本': 1},
    5: {'次品率': 0.1, '购买单价': 8, '检测成本': 1},
    6: {'次品率': 0.1, '购买单价': 12, '检测成本': 2},
    7: {'次品率': 0.1, '购买单价': 8, '检测成本': 1},
    8: {'次品率': 0.1, '购买单价': 12, '检测成本': 2}
}
semi_info = {
    1: {'次品率': 0.1, '装配成本': 8, '检测成本': 4, '拆解费用': 6},
    2: {'次品率': 0.1, '装配成本': 8, '检测成本': 4, '拆解费用': 6},
    3: {'次品率': 0.1, '装配成本': 8, '检测成本': 4, '拆解费用': 6}
}
finished_product_info = {'次品率': 0.1, '装配成本': 8, '检测成本': 6, '拆解费用': 10, '市场售价': 200, '调换损失': 40}

# 假设每种零配件购买数量为15000
part_num = 1666
strategies = []
results = []

# 生成所有可能的决策策略组合（每个零件、半成品是否检测，半成品是否拆解，成品是否检测、是否拆解的所有情况）
for part_decisions in product([0, 1], repeat=8):  # 零件是否检测决策，0表示不检测，1表示检测
    for semi_detect_decisions in product([0, 1], repeat=3):  # 半成品是否检测决策
        for semi_disassemble_decisions in product([0, 1], repeat=3):  # 半成品是否拆解决策
            for finished_product_detect_decision in [0, 1]:  # 成品是否检测决策
                for finished_product_disassemble_decision in [0, 1]:  # 成品是否拆解决策
                    strategies.append((part_decisions, semi_detect_decisions, semi_disassemble_decisions,
                                       finished_product_detect_decision, finished_product_disassemble_decision))

# 输出策略的总数
print(len(strategies))

# 计算每种策略下的成本、收益和利润，并找出最优策略以及利润最大的前十种策略
profits = []
for strategy in strategies:
    part_detections, semi_detect_decisions, semi_disassemble_decisions, finished_product_detect_decision, \
    finished_product_disassemble_decision = strategy
    # 零件购买成本
    part_purchase_cost = sum(part_num * part_info[i + 1]['购买单价'] for i in range(8))
    # 零件检测成本
    part_detection_cost = sum(part_num * part_info[i + 1]['检测成本'] * part_detections[i] for i in range(8))

    # 计算各半成品有效零件数量（考虑零件检测情况）
    semi_available_parts = []
    for semi_idx in range(3):
        available_parts = []
        start_idx = semi_idx * 3
        for i in range(3):
            actual_idx = start_idx + i
            if actual_idx < 8:  # 增加索引范围判断，防止越界
                if part_detections[actual_idx] == 0:
                    available_parts.append(part_num)
                else:
                    available_parts.append(part_num * (1 - part_info[actual_idx + 1]['次品率']))
            else:
                break
        if available_parts:  # 确保列表不为空再取最小值
            semi_available_parts.append(min(available_parts))

    # 半成品装配成本
    semi_assembly_cost = sum(semi_info[j + 1]['装配成本'] * semi_available_parts[j] if j < len(semi_available_parts) else 0 for j in range(3))
    # 半成品检测成本
    semi_detection_cost = sum(semi_info[j + 1]['检测成本'] * semi_detect_decisions[j] * semi_available_parts[j] if j < len(semi_available_parts) else 0
                              for j in range(3))
    # 半成品拆解成本（考虑是否拆解决策及相关条件）
    semi_disassembly_cost = sum(semi_info[j + 1]['拆解费用'] * semi_disassemble_decisions[j] * semi_info[j + 1]['次品率'] *
                                part_detections[(j * 3)] * semi_available_parts[j] if j < len(semi_available_parts) else 0
                                for j in range(3))

    # 成品装配成本（考虑半成品情况）
    finished_product_assembly_cost = 0
    if semi_available_parts:
        available_assembly_rates = [
            x * (1 - semi_info[i + 1]['次品率'] * semi_detect_decisions[i]) if i < len(semi_available_parts) else 1
            for i, x in enumerate(semi_available_parts)
        ]
        finished_product_assembly_cost = finished_product_info['装配成本'] * min(available_assembly_rates)

    # 成品检测成本（考虑是否检测决策及半成品情况）
    finished_product_detection_cost = 0
    if semi_available_parts and finished_product_detect_decision:
        available_detection_rates = [
            x * (1 - semi_info[i + 1]['次品率'] * semi_detect_decisions[i]) if i < len(semi_available_parts) else 1
            for i, x in enumerate(semi_available_parts)
        ]
        finished_product_detection_cost = finished_product_info['检测成本'] * min(available_detection_rates)

    # 成品拆解成本（考虑是否拆解决策、次品率及相关条件）
    finished_product_disassembly_cost = 0
    if semi_available_parts and finished_product_disassemble_decision:
        available_disassembly_rates = [
            x * (1 - semi_info[i + 1]['次品率'] * semi_detect_decisions[i]) if i < len(semi_available_parts) else 1
            for i, x in enumerate(semi_available_parts)
        ]
        finished_product_disassembly_cost = finished_product_info['拆解费用'] * finished_product_info['次品率'] * part_detections[0] * min(
            available_disassembly_rates)

    # 用户调换成本（考虑成品次品率及相关情况）
    user_exchange_cost = 0
    if semi_available_parts:
        available_exchange_rates = [
            x * (1 - semi_info[i + 1]['次品率'] * semi_detect_decisions[i]) if i < len(semi_available_parts) else 1
            for i, x in enumerate(semi_available_parts)
        ]
        user_exchange_cost = finished_product_info['调换损失'] * finished_product_info['次品率'] * min(available_exchange_rates)

    # 收益计算（考虑成品情况）
    revenue = 0
    if semi_available_parts:
        available_revenue_rates = [
            x * (1 - semi_info[i + 1]['次品率'] * semi_detect_decisions[i]) if i < len(semi_available_parts) else 1
            for i, x in enumerate(semi_available_parts)
        ]
        revenue = finished_product_info['市场售价'] * (1 - finished_product_info['次品率']) * min(available_revenue_rates)

    # 总成本
    total_cost = part_purchase_cost + part_detection_cost + semi_assembly_cost + semi_detection_cost + \
                 semi_disassembly_cost + finished_product_assembly_cost + finished_product_detection_cost + \
                 finished_product_disassembly_cost + user_exchange_cost
    # 利润
    profit = revenue - total_cost
    profits.append(profit)
    results.append((strategy, profit))

# 找出利润最大的前十种策略的索引
top_10_indices = sorted(range(len(profits)), key=lambda k: profits[k], reverse=True)[:10]

# 输出利润最大的前十种策略（输出具体策略内容）
print("利润最大的前十种策略如下：")
for index in top_10_indices:
    part_detections, semi_detect_decisions, semi_disassemble_decisions, finished_product_detect_decision, \
    finished_product_disassemble_decision = strategies[index]
    print(f"策略 {index}: 零件检测决策 {part_detections}, 半成品检测决策 {semi_detect_decisions}, 半成品拆解决策 {semi_disassemble_decisions}, "
          f"成品检测决策 {finished_product_detect_decision}, 成品拆解决策 {finished_product_disassemble_decision}")

# 绘制柱状图展示利润最大的前十种策略的利润情况并保存到images文件夹，进行图像美化
top_10_profits = [profits[i] for i in top_10_indices]
plt.figure(figsize=(10, 6))  # 设置图像大小，宽度10，高度6
# 设置柱状图颜色、边框颜色、边框宽度等，使外观更美观
bars = plt.bar(range(10), top_10_profits, color='skyblue', edgecolor='black', linewidth=1)
for bar in bars:
    bar.set_alpha(0.8)  # 设置柱子透明度，增加层次感

plt.xlabel('策略序号（前十）', fontsize=12)  # 设置x轴标签及字体大小
plt.ylabel('利润', fontsize=12)  # 设置y轴标签及字体大小
plt.title('利润最大的前十种策略利润情况', fontsize=14)  # 设置标题及字体大小
plt.xticks(range(10), [f"策略 {i}" for i in top_10_indices], rotation=45)  # 设置x轴刻度标签及旋转角度
plt.tick_params(axis='both', which='major', labelsize=10)  # 设置坐标轴刻度字体大小

# 添加数据标签到柱子上，显示具体利润数值
for i, v in enumerate(top_10_profits):
    plt.text(i, v + 0.05 * max(top_10_profits), str(round(v, 2)), ha='center', va='bottom', fontsize=10)

plt.grid(axis='y', linestyle='--', alpha=0.7)  # 添加y轴网格线，设置样式和透明度
# 添加水平参考线（均值线等，这里以平均利润为例），可以帮助对比数据
mean_profit = sum(top_10_profits) / len(top_10_profits)
plt.axhline(y=mean_profit, color='r', linestyle='--', label='平均利润')
plt.legend()  # 显示图例

# 保存柱状图到images文件夹，需确保文件夹已存在
save_path = os.path.join('images', 'top_10_profits_bar_chart.png')
plt.savefig(save_path)

# 绘制决策树（简单示例，仅展示结构，可根据实际进一步完善节点内容等），进行图像美化
dot = graphviz.Digraph(comment='Decision Tree', format='png')
dot.attr('node', shape='box', style='filled', fillcolor='lightgray', fontname='SimHei')  # 设置节点通用样式，如形状、填充色、字体

# 添加根节点（这里简单示例，可细化）
dot.node('Root', '生产决策流程', color='blue', fontcolor='white')

# 添加零件相关节点及边（检测决策）
for i in range(8):
    part_node_name = f'Part_{i + 1}'
    part_detect_name = f'{part_node_name}_Detect'
    part_no_detect_name = f'{part_node_name}_No_Detect'
    dot.node(part_node_name, f'零件{i + 1}', color='green', fontcolor='white')
    dot.edge('Root', part_node_name, color='black')
    dot.node(part_detect_name, f'零件{i + 1}检测', color='lightgreen', fontcolor='black')
    dot.edge(part_node_name, part_detect_name, color='black')
    dot.node(part_no_detect_name, f'零件{i + 1}不检测', color='darkgreen', fontcolor='black')
    dot.edge(part_node_name, part_no_detect_name, color='black')

    # 基于零件检测情况添加指向对应半成品检测节点的边（这里简单示意，实际逻辑更复杂，需结合具体流向）
    for j in range(3):
        semi_node_name = f'Semi_{j + 1}'
        semi_detect_name = f'{semi_node_name}_Detect'
        dot.edge(part_detect_name, semi_detect_name, color='gray')

# 添加半成品相关节点及边（检测、拆解决策）
for j in range(3):
    semi_node_name = f'Semi_{j + 1}'
    semi_detect_name = f'{semi_node_name}_Detect'
    semi_no_detect_name = f'{semi_node_name}_No_Detect'
    dot.node(semi_node_name, f'半成品{j + 1}', color='orange', fontcolor='white')
    dot.edge('Root', semi_node_name, color='black')
    dot.node(semi_detect_name, f'半成品{j + 1}检测', color='yellow', fontcolor='black')
    dot.edge(semi_node_name, semi_detect_name, color='black')
    dot.node(semi_no_detect_name, f'半成品{j + 1}不检测', color='darkorange', fontcolor='black')
    dot.edge(semi_node_name, semi_no_detect_name, color='black')

    semi_disassemble_name = f'{semi_node_name}_Disassemble'
    semi_no_disassemble_name = f'{semi_node_name}_No_Disassemble'
    dot.node(semi_disassemble_name, f'半成品{j + 1}拆解', color='brown', fontcolor='black')
    dot.edge(semi_detect_name, semi_disassemble_name, color='black')
    dot.node(semi_no_disassemble_name, f'半成品{j + 1}不拆解', color='pink', fontcolor='black')
    dot.edge(semi_detect_name, semi_no_disassemble_name, color='black')

    # 基于半成品情况添加指向成品检测节点的边（这里简单示意，实际逻辑更复杂，需结合具体流向）
    finished_product_detect_name = 'Finished_Product_Detect'
    dot.edge(semi_no_disassemble_name, finished_product_detect_name, color='gray')

# 添加成品相关节点及边（检测、拆解决策）
finished_product_detect_name = 'Finished_Product_Detect'
finished_product_no_detect_name = 'Finished_Product_No_Detect'
dot.node(finished_product_detect_name, '成品检测', color='cyan', fontcolor='white')
dot.edge('Root', finished_product_detect_name, color='black')
dot.node(finished_product_no_detect_name, '成品不检测', color='magenta', fontcolor='white')
dot.edge(finished_product_detect_name, finished_product_no_detect_name, color='black')

finished_product_disassemble_name = 'Finished_Product_Disassemble'
finished_product_no_disassemble_name = 'Finished_Product_No_Disassemble'
dot.node(finished_product_disassemble_name, '成品拆解', color='lime', fontcolor='white')
dot.edge(finished_product_detect_name, finished_product_disassemble_name, color='black')
dot.node(finished_product_no_disassemble_name, '成品不拆解', color='gold', fontcolor='white')
dot.edge(finished_product_detect_name, finished_product_no_disassemble_name, color='black')

# 保存决策树图像到images文件夹（需提前创建该文件夹）
dot.render('images/decision_tree', view=True)

（4）遗传算法直接求解：

import random
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
# 设置中文字体为黑体（需确保系统中有黑体字体文件，可根据实际情况调整字体）
mpl.rcParams['font.sans-serif'] = ['SimHei']
# 解决负号显示问题
mpl.rcParams['axes.unicode_minus'] = False
# 零配件、半成品、成品相关参数定义
part_info = {
    1: {'次品率': 0.1, '购买单价': 2, '检测成本': 1},
    2: {'次品率': 0.1, '购买单价': 8, '检测成本': 1},
    3: {'次品率': 0.1, '购买单价': 12, '检测成本': 2},
    4: {'次品率': 0.1, '购买单价': 2, '检测成本': 1},
    5: {'次品率': 0.1, '购买单价': 8, '检测成本': 1},
    6: {'次品率': 0.1, '购买单价': 12, '检测成本': 2},
    7: {'次品率': 0.1, '购买单价': 8, '检测成本': 1},
    8: {'次品率': 0.1, '购买单价': 12, '检测成本': 2}
}
semi_info = {
    1: {'次品率': 0.1, '装配成本': 8, '检测成本': 4, '拆解费用': 6},
    2: {'次品率': 0.1, '装配成本': 8, '检测成本': 4, '拆解费用': 6},
    3: {'次品率': 0.1, '装配成本': 8, '检测成本': 4, '拆解费用': 6}
}
finished_product_info = {'次品率': 0.1, '装配成本': 8, '检测成本': 6, '拆解费用': 10, '市场售价': 200, '调换损失': 40}
# 假设每种零配件购买数量
part_num = 15000
# 种群规模
pop_size = 100
# 染色体长度（对应决策数量）
chromosome_length = 24
# 最大迭代次数
max_generation = 500
# 变异概率
mutation_prob = 0.01


# 初始化种群，生成pop_size个长度为chromosome_length的二进制染色体个体
def initialize_population():
    population = []
    for _ in range(pop_size):
        individual = [random.randint(0, 1) for _ in range(chromosome_length)]
        population.append(individual)
    return population


# 计算个体适应度（利润）
def calculate_fitness(individual):
    # 零件购买成本
    part_purchase_cost = sum(part_num * part_info[i + 1]['购买单价'] for i in range(8))
    # 零件检测成本
    part_detection_cost = sum(part_num * part_info[i + 1]['检测成本'] * individual[i] for i in range(8))

    # 计算各半成品有效零件数量
    semi_available_parts = []
    for semi_idx in range(3):
        available_parts = []
        start_idx = semi_idx * 3
        for i in range(3):
            actual_idx = start_idx + i
            if actual_idx < 8:  # 增加判断，防止索引超出part_info的键范围
                if individual[actual_idx] == 0:
                    available_parts.append(part_num)
                else:
                    available_parts.append(part_num * (1 - part_info[actual_idx + 1]['次品率']))
            else:
                break
        if available_parts:
            semi_available_parts.append(min(available_parts))

    # 半成品装配成本
    semi_assembly_cost = sum(semi_info[j + 1]['装配成本'] * semi_available_parts[j] if j < len(semi_available_parts) else 0 for j in range(3))
    # 半成品检测成本
    semi_detection_cost = sum(semi_info[j + 1]['检测成本'] * individual[8 + j] * semi_available_parts[j] if j < len(semi_available_parts) else 0
                              for j in range(3))
    # 半成品拆解成本
    semi_disassembly_cost = sum(semi_info[j + 1]['拆解费用'] * semi_info[j + 1]['次品率'] * individual[8 + j] *
                                individual[start_idx] * semi_available_parts[j] if j < len(semi_available_parts) else 0
                                for j in range(3))

    # 成品装配成本
    finished_product_assembly_cost = 0
    if semi_available_parts:
        available_assembly_rates = [
            x * (1 - semi_info[i + 1]['次品率'] * individual[8 + i]) if i < len(semi_available_parts) else 1
            for i, x in enumerate(semi_available_parts)
        ]
        finished_product_assembly_cost = finished_product_info['装配成本'] * min(available_assembly_rates)

    # 成品检测成本
    finished_product_detection_cost = 0
    if semi_available_parts and individual[17]:
        available_detection_rates = [
            x * (1 - semi_info[i + 1]['次品率'] * individual[8 + i]) if i < len(semi_available_parts) else 1
            for i, x in enumerate(semi_available_parts)
        ]
        finished_product_detection_cost = finished_product_info['检测成本'] * min(available_detection_rates)

    # 成品拆解成本
    finished_product_disassembly_cost = 0
    if semi_available_parts and individual[17]:
        available_disassembly_rates = [
            x * (1 - semi_info[i + 1]['次品率'] * individual[8 + i]) if i < len(semi_available_parts) else 1
            for i, x in enumerate(semi_available_parts)
        ]
        finished_product_disassembly_cost = finished_product_info['拆解费用'] * finished_product_info['次品率'] * individual[0] * min(
            available_disassembly_rates)

    # 用户调换成本
    user_exchange_cost = 0
    if semi_available_parts:
        available_exchange_rates = [
            x * (1 - semi_info[i + 1]['次品率'] * individual[8 + i]) if i < len(semi_available_parts) else 1
            for i, x in enumerate(semi_available_parts)
        ]
        user_exchange_cost = finished_product_info['调换损失'] * finished_product_info['次品率'] * min(available_exchange_rates)

    # 收益计算
    revenue = 0
    if semi_available_parts:
        available_revenue_rates = [
            x * (1 - semi_info[i + 1]['次品率'] * individual[8 + i]) if i < len(semi_available_parts) else 1
            for i, x in enumerate(semi_available_parts)
        ]
        revenue = finished_product_info['市场售价'] * (1 - finished_product_info['次品率']) * min(available_revenue_rates)

    # 总成本
    total_cost = part_purchase_cost + part_detection_cost + semi_assembly_cost + semi_detection_cost + \
                 semi_disassembly_cost + finished_product_assembly_cost + finished_product_detection_cost + \
                 finished_product_disassembly_cost + user_exchange_cost
    # 利润（适应度）
    fitness = revenue - total_cost
    return fitness


# 选择操作（轮盘赌选择）
def selection(population):
    fitness_values = [calculate_fitness(ind) for ind in population]
    total_fitness = sum(fitness_values)
    selection_probs = [f / total_fitness for f in fitness_values]
    selected_indices = np.random.choice(len(population), size=pop_size, p=selection_probs, replace=True)
    selected_population = [population[i] for i in selected_indices]
    return selected_population


# 交叉操作（单点交叉）
def crossover(parent1, parent2):
    crossover_point = random.randint(1, chromosome_length - 1)
    child1 = parent1[:crossover_point] + parent2[crossover_point:]
    child2 = parent2[:crossover_point] + parent1[crossover_point:]
    return child1, child2


# 变异操作
def mutation(individual):
    for i in range(len(individual)):
        if random.random() < mutation_prob:
            individual[i] = 1 - individual[i]
    return individual


def genetic_algorithm():
    all_fitnesses = []  # 用于存储每次运行遗传算法得到的每代适应度列表
    num_runs = 10  # 多次运行遗传算法的次数，可根据实际调整
    for _ in range(num_runs):
        population = initialize_population()
        best_fitness_list = []
        for _ in range(max_generation):
            new_population = []
            selected_population = selection(population)
            for i in range(0, pop_size, 2):
                parent1 = selected_population[i]
                parent2 = selected_population[i + 1]
                child1, child2 = crossover(parent1, parent2)
                child1 = mutation(child1)
                child2 = mutation(child2)
                new_population.extend([child1, child2])
            population = new_population
            # 找到当前代的最优个体及适应度并记录
            best_individual = max(population, key=calculate_fitness)
            best_fitness = calculate_fitness(best_individual)
            best_fitness_list.append(best_fitness)
        all_fitnesses.append(best_fitness_list)
    return all_fitnesses


def plot_fitness_change(all_fitnesses):
    generations = range(1, max_generation + 1)
    fig, ax = plt.subplots()
    # 添加背景方格
    ax.grid(which='both', color='gray', linestyle=':', linewidth=0.5)

    # 计算平均适应度，确保维度匹配
    average_fitness_list = np.mean(all_fitnesses, axis=0)

    # 分段设置折线图颜色，这里简单示例分为三段，可根据实际情况调整逻辑
    segment_size = max_generation // 3
    colors = ['blue', 'green', 'red']
    for i in range(len(colors)):
        start = i * segment_size
        end = (i + 1) * segment_size if i < len(colors) - 1 else max_generation
        ax.plot(generations[start:end], average_fitness_list[start:end], color=colors[i], label=f'阶段{i + 1}')

    ax.set_xlabel('迭代次数')
    ax.set_ylabel('最优个体平均适应度')
    ax.set_title('迭代过程中最优个体平均适应度变化情况')
    # 添加颜色条，这里以平均适应度列表的值作为颜色映射依据，需根据实际合理调整
    norm = plt.Normalize(min(average_fitness_list), max(average_fitness_list))
    sm = plt.cm.ScalarMappable(cmap='viridis', norm=norm)
    sm.set_array([])
    cbar = fig.colorbar(sm, ax=ax)
    cbar.set_label('平均适应度数值')
    plt.legend()
    plt.savefig('images/最优个体平均适应度变化情况.png')
    plt.show()


def plot_3d_dummy(all_fitnesses):
    import numpy as np
    from mpl_toolkits.mplot3d import Axes3D

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    # 种群规模范围，这里只有一个固定的种群规模值，为了和网格数据生成逻辑匹配，仍使用数组形式
    pop_sizes = np.array([pop_size])
    # 迭代次数范围
    generations = np.arange(1, max_generation + 1)
    # 生成网格数据
    pop_sizes_grid, generations_grid = np.meshgrid(pop_sizes, generations)

    # 确保all_fitnesses是二维数组形式（如果不是，进行转换）
    all_fitnesses = np.array(all_fitnesses)
    num_runs, num_generations = all_fitnesses.shape
    # 将二维的all_fitnesses列表转换为适合绘制散点图的一维数组形式
    # 这里调整为确保数据正确平铺展开
    profits = all_fitnesses.flatten()

    # 调整x、y坐标数据形状使其满足广播机制要求
    # 重复pop_sizes_grid对应每个迭代次数 * 运行次数，得到形状为(num_runs * max_generation,)的数组
    x_data = np.repeat(pop_sizes_grid.flatten(), num_runs)
    # 对应每个种群规模重复迭代次数范围 * 运行次数，得到形状为(num_runs * max_generation,)的数组
    y_data = np.tile(generations, (num_runs, 1)).flatten()

    # 使用散点图绘制，颜色根据利润（适应度）值映射
    sc = ax.scatter(x_data, y_data, profits, c=profits, cmap='viridis', s=50)

    # 明确标注x、y、z轴含义
    ax.set_xlabel('初始种群个数')
    ax.set_ylabel('迭代次数')
    ax.set_zlabel('利润（适应度）')
    ax.set_title('种群个数、迭代次数与利润（适应度）关系三维示例图像')

    # 添加颜色条说明散点颜色对应的数值含义
    fig.colorbar(sc, ax=ax, label='利润（适应度）数值')

    plt.savefig('images/三维示例图像.png')
    plt.show()


# 运行遗传算法并获取所有适应度数据
all_fitnesses = genetic_algorithm()

# 绘制迭代过程中最优个体变化情况图像
plot_fitness_change(all_fitnesses)

# 绘制三维示例图像
plot_3d_dummy(all_fitnesses)

（5）遗传算法+精英策略：

import random
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import rcParams
import os

# 设置中文字体，解决中文显示问题
rcParams['font.sans-serif'] = ['SimHei']
rcParams['axes.unicode_minus'] = False

# 零配件、半成品、成品相关参数定义
part_info = {
    1: {'次品率': 0.1, '购买单价': 2, '检测成本': 1},
    2: {'次品率': 0.1, '购买单价': 8, '检测成本': 1},
    3: {'次品率': 0.1, '购买单价': 12, '检测成本': 2},
    4: {'次品率': 0.1, '购买单价': 2, '检测成本': 1},
    5: {'次品率': 0.1, '购买单价': 8, '检测成本': 1},
    6: {'次品率': 0.1, '购买单价': 12, '检测成本': 2},
    7: {'次品率': 0.1, '购买单价': 8, '检测成本': 1},
    8: {'次品率': 0.1, '购买单价': 12, '检测成本': 2}
}
semi_info = {
    1: {'次品率': 0.1, '装配成本': 8, '检测成本': 4, '拆解费用': 6},
    2: {'次品率': 0.1, '装配成本': 8, '检测成本': 4, '拆解费用': 6},
    3: {'次品率': 0.1, '装配成本': 8, '检测成本': 4, '拆解费用': 6}
}
finished_product_info = {'次品率': 0.1, '装配成本': 8, '检测成本': 6, '拆解费用': 10, '市场售价': 200, '调换损失': 40}
# 假设每种零配件购买数量
part_num = 15000
# 最大迭代次数
max_generation = 30
# 变异概率
mutation_prob = 0.01
# 精英个体数量占种群规模的比例，可根据实际情况调整
elite_ratio = 0.1

# 初始化种群，生成指定个数的长度为chromosome_length的二进制染色体个体
def initialize_population(pop_size, chromosome_length):
    population = []
    for _ in range(pop_size):
        individual = [random.randint(0, 1) for _ in range(chromosome_length)]
        population.append(individual)
    return population


# 计算个体适应度（利润）
def calculate_fitness(individual):
    try:
        # 零件购买成本
        part_purchase_cost = sum(part_num * part_info[i + 1]['购买单价'] for i in range(8))
        # 零件检测成本
        part_detection_cost = sum(part_num * part_info[i + 1]['检测成本'] * individual[i] for i in range(8))

        # 计算各半成品有效零件数量
        semi_available_parts = []
        for semi_idx in range(3):
            available_parts = []
            start_idx = semi_idx * 3
            for i in range(3):
                actual_idx = start_idx + i
                if actual_idx < 8:  # 增加判断，防止索引超出part_info的键范围
                    if individual[actual_idx] == 0:
                        available_parts.append(part_num)
                    else:
                        available_parts.append(part_num * (1 - part_info[actual_idx + 1]['次品率']))
                else:
                    break
            if available_parts:
                semi_available_parts.append(min(available_parts))

        # 半成品装配成本
        semi_assembly_cost = sum(semi_info[j + 1]['装配成本'] * semi_available_parts[j] if j < len(semi_available_parts) else 0 for j in range(3))
        # 半成品检测成本
        semi_detection_cost = sum(semi_info[j + 1]['检测成本'] * individual[8 + j] * semi_available_parts[j] if j < len(semi_available_parts) else 0
                                  for j in range(3))
        # 半成品拆解成本
        semi_disassembly_cost = sum(semi_info[j + 1]['拆解费用'] * individual[11 + j] * semi_info[j + 1]['次品率'] *
                                    individual[start_idx] * semi_available_parts[j] if j < len(semi_available_parts) else 0
                                    for j in range(3))

        # 成品装配成本
        finished_product_assembly_cost = 0
        if semi_available_parts:
            available_assembly_rates = [
                x * (1 - semi_info[i + 1]['次品率'] * individual[8 + i]) if i < len(semi_available_parts) else 1
                for i, x in enumerate(semi_available_parts)
            ]
            finished_product_assembly_cost = finished_product_info['装配成本'] * min(available_assembly_rates)

        # 成品检测成本
        finished_product_detection_cost = 0
        if semi_available_parts and individual[14]:
            available_detection_rates = [
                x * (1 - semi_info[i + 1]['次品率'] * individual[8 + i]) if i < len(semi_available_parts) else 1
                for i, x in enumerate(semi_available_parts)
            ]
            finished_product_detection_cost = finished_product_info['检测成本'] * min(available_detection_rates)

        # 成品拆解成本
        finished_product_disassembly_cost = 0
        if semi_available_parts and individual[15]:
            available_disassembly_rates = [
                x * (1 - semi_info[i + 1]['次品率'] * individual[8 + i]) if i < len(semi_available_parts) else 1
                for i, x in enumerate(semi_available_parts)
            ]
            finished_product_disassembly_cost = finished_product_info['拆解费用'] * finished_product_info['次品率'] * individual[0] * min(
                available_disassembly_rates)

        # 用户调换成本
        user_exchange_cost = 0
        if semi_available_parts:
            available_exchange_rates = [
                x * (1 - semi_info[i + 1]['次品率'] * individual[8 + i]) if i < len(semi_available_parts) else 1
                for i, x in enumerate(semi_available_parts)
            ]
            user_exchange_cost = finished_product_info['调换损失'] * finished_product_info['次品率'] * min(available_exchange_rates)

        # 收益计算
        revenue = 0
        if semi_available_parts:
            available_revenue_rates = [
                x * (1 - semi_info[i + 1]['次品率'] * individual[8 + i]) if i < len(semi_available_parts) else 1
                for i, x in enumerate(semi_available_parts)
            ]
            revenue = finished_product_info['市场售价'] * (1 - finished_product_info['次品率']) * min(available_revenue_rates)

        # 总成本
        total_cost = part_purchase_cost + part_detection_cost + semi_assembly_cost + semi_detection_cost + \
                     semi_disassembly_cost + finished_product_assembly_cost + finished_product_detection_cost + \
                     finished_product_disassembly_cost + user_exchange_cost
        # 利润（适应度）
        fitness = revenue - total_cost
    except ZeroDivisionError:
        fitness = 0
    except:
        fitness = 0
    return fitness


# 选择操作（轮盘赌选择）
def selection(population):
    fitness_values = [calculate_fitness(ind) for ind in population]
    if len(fitness_values) == 0:
        return population
    elif all(f == 0 for f in fitness_values):
        selection_probs = [1 / len(fitness_values)] * len(fitness_values)
    else:
        total_fitness = sum(fitness_values)
        if total_fitness == 0:
            selection_probs = [1 / len(fitness_values)] * len(fitness_values)
        else:
            selection_probs = [f / total_fitness for f in fitness_values]
            # 归一化处理，确保概率之和为1
            prob_sum = sum(selection_probs)
            if prob_sum!= 0:
                selection_probs = [p / prob_sum for p in selection_probs]
    selected_indices = np.random.choice(len(population), size=len(population), p=selection_probs, replace=True)
    selected_population = [population[i] for i in selected_indices]
    return selected_population


# 交叉操作（单点交叉）
def crossover(parent1, parent2):
    crossover_point = random.randint(1, chromosome_length - 1)
    child1 = parent1[:crossover_point] + parent2[crossover_point:]
    child2 = parent2[:crossover_point] + parent1[crossover_point:]
    return child1, child2


# 变异操作
def mutation(individual):
    for i in range(len(individual)):
        if random.random() < mutation_prob:
            individual[i] = 1 - individual[i]
    return individual


# 遗传算法主函数，添加精英策略
def genetic_algorithm(pop_sizes):
    all_results = []
    for pop_size in pop_sizes:
        # 对于每个种群规模，进行多次遗传算法运行并记录结果
        results_for_size = []
        for _ in range(5):  # 这里可以调整针对每个种群规模运行的次数，取平均等操作更稳定
            population = initialize_population(pop_size, chromosome_length)
            for generation in range(max_generation):
                # 选择精英个体
                fitness_values = [calculate_fitness(ind) for ind in population]
                elite_size = int(pop_size * elite_ratio)
                elite_indices = np.argsort(fitness_values)[-elite_size:]
                elites = [population[i] for i in elite_indices]

                # 剩余种群进行选择、交叉、变异操作
                remaining_population = [population[i] for i in range(len(population)) if i not in elite_indices]
                selected_population = selection(remaining_population)
                new_population = []
                for i in range(0, len(selected_population), 2):
                    parent1 = selected_population[i]
                    if i + 1 < len(selected_population):
                        parent2 = selected_population[i + 1]
                        child1, child2 = crossover(parent1, parent2)
                        child1 = mutation(child1)
                        child2 = mutation(child2)
                        new_population.extend([child1, child2])
                    else:
                        # 如果只剩下一个个体，无法进行交叉操作，直接将其加入新种群（也可以考虑其他处理方式，比如复制一份加入等）
                        new_population.append(parent1)

                # 合并精英个体和新生成个体形成新一代种群
                if len(new_population) < pop_size:
                    # 如果新生成个体数量不足种群规模，补充精英个体至种群规模
                    new_population.extend(elites[:pop_size - len(new_population)])
                else:
                    # 如果新生成个体数量足够，随机替换部分个体为精英个体
                    replace_indices = random.sample(range(pop_size), elite_size)
                    for i, idx in enumerate(replace_indices):
                        new_population[idx] = elites[i]
                population = new_population

                # 添加判断，如果种群为空，赋予一个默认个体（这里简单示例为全0个体，可根据实际调整）
                if not population:
                    default_individual = [0] * chromosome_length
                    population = [default_individual]

                best_individual = max(population, key=calculate_fitness)
                best_fitness = calculate_fitness(best_individual)
                results_for_size.append(best_fitness)
        all_results.append(results_for_size)
    return all_results


# 将二进制个体转换为中文描述的策略，添加类型判断确保参数符合要求
def binary_to_chinese_strategy(individual):
    if isinstance(individual, list):
        part_decision_str = "零件检测决策：" + ", ".join(
            ["检测" if bit == 1 else "不检测" for bit in individual[:8]])
        semi_detect_decision_str = "半成品检测决策：" + ", ".join(
            ["检测" if bit == 1 else "不检测" for bit in individual[8:11]])
        semi_disassemble_decision_str = "半成品拆解决策：" + ", ".join(
            ["拆解" if bit == 1 else "不拆解" for bit in individual[11:14]])
        finished_product_detect_decision_str = "成品检测决策：" + ("检测" if individual[14] == 1 else "不检测")
        finished_product_disassemble_decision_str = "成品拆解决策：" + ("拆解" if individual[15] == 1 else "不拆解")
        return part_decision_str + "\n" + semi_detect_decision_str + "\n" + semi_disassemble_decision_str + "\n" + \
               finished_product_detect_decision_str + "\n" + finished_product_disassemble_decision_str
    else:
        return "传入参数类型不符合要求，无法生成中文描述策略"


# 设置种群规模的范围，可根据实际情况调整范围和步长
pop_sizes = range(0, 101, 10)
chromosome_length = 16 #  二进制长度
results = genetic_algorithm(pop_sizes)

# 绘制三维图像展示种群个数、迭代次数和利润（适应度）关系，并保存到images文件夹
os.makedirs('images', exist_ok=True)  # 确保images文件夹存在
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
colors = plt.cm.jet(np.linspace(0, 1, len(pop_sizes)))  # 生成与种群规模数量对应的颜色序列
for i, pop_size in enumerate(pop_sizes):
    mean_results = np.mean(results[i])  # 计算每个种群规模下多次运行的平均利润
    ax.scatter([pop_size] * max_generation, list(range(1, max_generation + 1)), [mean_results] * max_generation, c=[colors[i]])

ax.set_xlabel('种群个数')
ax.set_ylabel('迭代次数')
ax.set_zlabel('利润（适应度）')
ax.set_title('种群个数、迭代次数与利润关系')
plt.savefig('images/population_iteration_fitness_3d.png')
plt.show()

# 运行遗传算法并输出结果（示例取最后一个种群规模下的结果展示）
best_individual, best_fitness = genetic_algorithm([pop_sizes[-1]])[0][-1], np.mean(
    genetic_algorithm([pop_sizes[-1]])[0])  # 取最后一个种群规模下最后一次运行结果作为示例展示，也可根据需求调整
# 添加类型判断，确保传入binary_to_chinese_strategy函数的参数符合要求
if isinstance(best_individual, list):
    print("最优决策组合（染色体表示）:", best_individual)
    print("最大利润（适应度）:", best_fitness)
    print("最优策略详情（中文描述）：\n", binary_to_chinese_strategy(best_individual))
else:
    print("获取到的最优个体类型不符合要求，无法输出中文描述策略。最优决策组合（染色体表示）:", best_individual)
    print("最大利润（适应度）:", best_fitness)