2024年第五届“华数杯”全国大学生数学建模竞赛 B题详细思路代码

没有更新完之前，专栏价格为59，更新完毕之后恢复到99. 专栏内包含2024年所有数学建模比赛思路和代码，有些重要比赛着重更新（华数杯、国赛、美赛），小比赛可能会有chatgpt4更新，只需订阅一次。有些文章没有完整代码，请到专栏内查找最新代码和思路。如果比赛结束后没有更新代码（可能会有事情来不及更新）赛后我会统一退款。

2024年第五届“华数杯”专栏地址：https://blog.csdn.net/m0_52343631/category_12482955.html?spm=1001.2014.3001.5482

目录

2024年第五届“华数杯”专栏地址：https://blog.csdn.net/m0_52343631/category_12482955.html?spm=1001.2014.3001.5482

摘　要

一、问题重述

二、模型的假设

三、问题的分析

五、模型建立与求解

六、模型总结

6.2 模型评价

参考文献

附　录

---

2023年优秀论文

摘　要

本文结合 Kubelka-Munk 光学模型 以及 Kubelka-Munk 单常数理论 建立出合适的

关于不透明制品配色的优化模型并使用计算机方法进行配色以提高实际生产效率．首

先使用 一元线性回归分析拟合参数 ，接着利用 SLSQP 和 随机化 方法求解 单目标优化

问题 得到与目标样本色差最小的配色方案，同时在考虑节约成本和资源以及批量配色

的要求下进一步改进模型求解 多目标优化问题 结合 遗传算法 和 TOPSIS 评价方法得到

最佳的多个配色方案，为实际生产提供重要的理论参考．

针对问题一

首先根据附件 2 ，绘制出红、黄、蓝三种着色剂在不同波长下的 𝐾 / 𝑆

与浓度关系图原始数据散点图后，观察分析图像，发现图中的绝大部分点均位于一条

直线附近，符合 线性回归 的特点．接着调用 Python 里 Scipy 库中的 linregress 线性回

归方法 求解拟合本题数据．建立 一元单变量线性回归模型 如： 𝑦 = 𝑎𝑥 +

𝑏

，其中 𝑦 代表

𝐾 / 𝑆 ,

𝑥

代表单位换算后的浓度如

0 . 5% 也就是 0.005.

最终答案详见

表 2

( 第

6 页 ).

针对问题二 首先依据问题一求出的在不同的波长

𝜆

下三种着色剂吸收系数 𝐾

与散热系数 𝑆 的比值 𝐾 / 𝑆 与浓度的函数关系式，根据已知条件和假设以及

K − M 光

学模型 和 K −

M

单常数理论 先求解配色而成的混色物的 𝑅 值；接着利用

辛普森积分

方法 按照 𝑑 ( 𝜆 )

=20

𝑛𝑚

为间隔将光谱波长

𝜆 从 400 𝑛𝑚 到

700 𝑛𝑚 进行积分求和 从而计算

得到目标样本的三刺激值和配方而成的混色物的三刺激值；最后可根据 CIELAB 色彩

空间的总色差计算方法建立出计算材料与目标样的色差的关系式．设定判定方案相同

条件，建立以 三种着色剂浓度为自变量 的 以色差最小为优化目标 以 浓度范围 、 色差要

求 和 方案不同 为约束条件的 单目标非线性优化模型 ，对每一个目标样本使用 SLSQP

方法和 随机化算法 从初始随机点不断进行迭代计算配出与目标样本的色差最为接近的

10 个不同配方 . 最终计算得到配方如图 4- 图 8( 第 10-11 页 ) ．

针对问题三

本题实质上是要求在色差接近的情况下，保证成本尽可能少．首先

在问题二优化模型的基础上建立以 三种着色剂浓度为自变量 的 以色差最小和成本尽可

能小为优化目标 的以 浓度范围 、 色差要求 和 方案不同 为约束条件 多目标非线性优化模

型 ，首先利用 遗传算法 直接求解该模型，接着利用 动态线性加权 的方法 , 首先去除色

差和成本的量纲将二者 归一化 ，接着设置权重 𝛼 改进优化模型，按照一定间隔从 0-1

遍历权重，对每一个目标样本使用 SLSQP 方法和 随机化算法 从初始随机点不断进行

迭代计算该权重下的最优的 10 个不同配方．为选取最佳权重 𝛼 和相对应的最优方案，

利用 TOPSIS 评价法对所得到的方案进行评分 , 依据评分选择得分高的方案．由于遗

传算法运行效率远优于第二种方法，且解集沿帕累托前沿均匀分布，避免了对方案评

分的主观性，最终选择 遗传算法 计算答案如图 9- 图 12( 第 12-13 页 ).

针对问题四

本题增加了使用配色的着色剂越少越好的条件，设所用的配色剂种

类数为 𝛽 ，由题目条件可建立以 三种着色剂浓度为自变量 的 以色差最小、成本尽可能

小以及种类尽可能少为优化目标 的 多目标非线性优化模型 ，同样利用 遗传算法 从初始

种群不断进化，交叉，变异求解得到最优解．计算答案如图 15- 图 19( 第 16-17 页 )

关键词： K-M 光学模型 回归分析 单 / 多目标优化 遗传算法

𝑇𝑂𝑃𝑆𝐼𝑆 评价

一、问题重述

1.1 问题背景

日常生活中色彩斑斓的不透明有色制品深受大众喜爱，它们由多种着色剂染色而

成．在实际生产中，许多相关配色行业还是以传统的人工配色为主，随着社会和经济

的发展，有色制品产品品质的要求日益提高．然而人工配色存在着主观性强、效率低

下等局限性，依靠人工配色已经远远不能满足生产需要．

随着计算机技术的日益成熟和计算能力的提高，依据相关光学和色彩理论的基础，

人们开始尝试利用计算机来处理色彩和图形设计相关的任务．相较于传统的人工配色，

采用计算机方法配色不仅能节省大量的人力，物力，财力而且也节约资源，对社会可持

续发展有着深远影响．在此大背景下，基于计算机配色理论以及相关的光学模型知识

研究如何通过计算机方法来实现不透明制品的配色已经成为一个重要的研究课题，对

社会发展具有重要意义．

1.2 问题提出

对于不透明材料而言，大部分光线会被其表面吸收或反射，其中吸收系数 𝐾 与散

热系数 𝑆 的比值与反射率 𝑅 之间存在一定的数学关系，为简化计算，最终配色后的颜

色的反射率以 20 𝑛𝑚 为间隔的光谱数据表示．针对某一不透明制品，已知红、黄、蓝 3

种着色剂在不同浓度下不同波长的 𝐾 /

𝑆

值以及基地材料在不同波长下的 𝐾 / 𝑆 值（见

附件 2 ）．色差计算中出现的三刺激值

XYZ

计算方法如下：

二、模型的假设

1 、假设着色剂和基底材料的吸收系数 𝐾 与散热系数 𝑆 的比值在同一波长下保持

稳定，不受其他材料影响．

2 、假设材料均匀，忽略材料边缘的影响．

3 、假设入射光为光谱均匀的漫射光．

4 、在 CIELAB 均匀色彩空间下计算色差．

5 、假设着色剂对物体均匀着色．

6 、假设色料的散射系数 𝑆 相对于基材的散热系数 𝑆 0 可忽略不记．

7 、忽略同色异谱现象．

三、问题的分析

3.1 问题一的分析

本题是一个关于拟合和回归分析的问题，需要在不同的波长下得到三种着色剂吸

收系数 𝐾 与散热系数 𝑆 的比值 𝐾 / 𝑆 与浓度的函数关系式 . 首先可以先根据题目所给数

据绘出原始数据散点图，然后分析图像得出数据可能符合的回归模型，最后再利用最

小二乘法的思想和 python 中自带的相关回归分析函数来求解本题．

3.2 问题二的分析

本题是建立在问题一求出不同的波长下得到三种着色剂吸收系数 𝐾 与散热系数

𝑆 的比值 𝐾 / 𝑆 与浓度的函数关系式后，根据已知条件和假设以及 K − M 光学模型 和

K − M 单常数理论 先求解配色而成的混色物的 𝑅 值后，接着计算出目标样本的三刺激

值和配方而成的混色物的三刺激值，最后可根据

CIELAB

色彩空间的总色差计算方法

建立出计算材料与目标样的色差的关系式．于是便建立以 三种着色剂浓度为自变量 的

以色差最小为优化目标 的 单目标非线性优化模型 ，可以调用 Python 中 scipy 库中的相

关函数对每一个目标样本进行求解得到合适的 10 个不同配方．

3.3 问题三的分析

本题给出了基地材料重量和色母粒价格，要求在问题 2 的基础上，考虑成本控制

和批量配色改进优化模型解决问题．实质上是要求在色差接近的情况下，保证成本尽

可能少．于是便建立以 三种着色剂浓度为自变量 的 以色差最小和成本尽可能小为优化

目标 的 多目标非线性优化模型 ，可采用 遗传算法 从初始种群不断进行进化，交叉，变

异求解此模型也可以利用 动态线性加权 的方法 , 首先去除色差和成本的量纲将二者 归

一化 ，接着设置权重 𝛼 改进优化模型将其转化为 单目标非线性优化模型 ，按照一定间

隔从 0-1 遍历权重，同时调用 Python 中 scipy 库中的相关函数对每一个目标样本进行

求解得到在该权重下的合适的 10 个不同配方．最后可利用 TOPSIS 法对所得到的方

案进行评分 , 最后选择得分高的方案．

3.4 问题四的分析

本题在问题三基础上增加了使用配色的着色剂越少越好的条件，设所用的配色剂

种类数为 𝛽 ，由题目条件可建立以 三种着色剂浓度为自变量 的 以色差最小、成本尽可

五、模型建立与求解

5.1 问题一的模型建立与求解

5.1.1 问题一模型建立

剂在波长 以下以求解红色着色剂在波长 𝜆 为 400 𝑛𝑚 时 𝐾 / 𝑆 与浓度的函数关系，黄色着色

𝜆 为 420 𝑛𝑚 时 𝐾 / 𝑆 与浓度的函数关系，以及蓝色着色剂在波长 𝜆 为 440 𝑛𝑚

时 𝐾 / 𝑆 与浓度的函数关系为例进行模型建立与求解说明，其余的函数关系式求解过程

均类似．

首先根据附件 2 绘出三种情况下的原始数据散点图：

问题二要求在满足 (7) 式约束条件下，运用优化模型配出与目标样本的色差最为接

近的 10 个不同配方，求解与第一个目标样本的色差最为接近的 10 个不同配方步骤如

下：

Step 1:

定义 配色方案相同的条件

• 若每个颜色的相对浓度差都在 5% 以内，则认为是同一种方案．

• 若全部颜色的相对浓度差的和在 10% 以内，则认为是同一种方案．

• 若每个颜色的浓度差不超过 0.0001 ，则认为是同一种方案 .

为了求解得到 10 个不同的配方，需要将每一次得到的方案与前面得到的符合题意的

答案进行比较，方案不能相同．

Step 2:

进行第一次最优解搜索．选择利用 python 中自带的 SLSQP 方法进行

最优化求解此问题．开始的第一次最优解搜索从一个预先计算的粗略值 [²] 开始，利用

SLSQP 方法不断进行迭代最后得出第一次的最优解．

Step 3:

为保证找到不同的方案，后续的最优解方案搜索利用 随机化思想 ，每一

次从 随机点 开始搜索 , 利用 SLSQP 方法不断进行 迭代计算 得出最优解，每一次得到的

方案均与前面已经得到的方案进行比较，如果满足配色方案相同的条件，则放弃此次

方案，开始新的搜索；否则，将此方案加入已选方案中，更新约束条件，继续搜索下

一个最优解，直至找出 10 个不同的方案．算法简要流程如图 3 ：

六、模型总结

6.2 模型评价

6.2.1 模型优点

增加了方案相同的判定条件从而有助于增加方案间差异性

如下：

• 若每个颜色的相对浓度差都在 5% 以内，则认为是同一种方案．

• 若全部颜色的相对浓度差的和在 10% 以内，则认为是同一种方案．

• 若每个颜色的浓度差不超过 0.0001 ，则认为是同一种方案 .

为了求解得到 10 个不同的配方，需要将每一次得到的方案与前面得到的符合题意的

答案进行比较得到不同方案，有助于避免实际生产过程中由于配色方案的过于接近产

生的可能不良影响．

计算可信度大

优化模型迭代计算中计算大量散点得到最终配方．以下用图 23 说

明 问题二 算法迭代计算得到目标样本一和目标样本二的十个配方的过程：

由图 24 和图 25 两图直观分析，利用二次函数模型在图像上拟合效果更好．但考

虑到红、黄、蓝三种着色剂在 𝜆 从 400 𝑛𝑚 - 700 𝑛𝑚 不同波长下 𝐾 / 𝑆 与浓度的关系整体

上均保持线性关系，于是在关于红色着色剂在 700 𝑛𝑚𝐾 / 𝑆 与浓度的关系上，本文仍然

采用线性回归模型进行拟合和计算．

可推广至极小区域里最优化搜索

由对图 23 的分析可知，所有符合目标样本的

10 个不同配方的点均在各自同一个小区域里，可以分析得知，该区域为最佳配色方案

区域．后续可以根据实际情况进行优化，在满足实际要求的情况下在由这 10 个不同配

方的点构成的最佳配色方案区域 , 运用 3 维点云处理的聚类分析方法 进行进一步搜索

最终得到最优的配色方案．

搜索方法改进

搜索最优解过程中可以使用图论思想，将搜索过程中满足条件的

优解视作一个点，点权为色差．若两个解相似，则在图上连边．为得到十个不同的方

案，在图上 求解大小为十的最小点权和点独立集 即可解决问题．但实际上求解最小点

权和点独立集难度很大而且点的规模过大，故本文未采取这种方法．

参考文献

[1] 姜鹏飞 . 计算机配色理论及算法的研究 [D/OL] . 硕士 , 中原工学院 , 2016.

[2] 周春霞，唐正宁 . 包装印刷专色油墨的计算机配色理论研究 . 包装工程 , (5):121–123,

2006.

[3] 王林吉 . 基于 CIELAB 均匀颜色空间和聚类算法的混纺测色研究 [D] . 硕士 , 浙江理

工大学 , 2010.

附　录

from collections import defaultdict
import pandas as pd
import numpy as np
from scipy import stats
import scipy.integrate
import scipy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mplfonts.bin.cli import init
init()
from mplfonts import use_font
use_font("Noto Serif CJK SC")
colors = ["红", "黄", "蓝"]
df = pd.read_excel("附件2.xlsx", skiprows=1)
df = df.rename(columns={"Unnamed: 0": "颜色"})
r_idx = df[df["颜色"] == "红"].index[0]
y_idx = df[df["颜色"] == "黄"].index[0]
b_idx = df[df["颜色"] == "蓝"].index[0]
default_idx = df[df["颜色"] == "基材"].index[0]
r_df = df[r_idx:y_idx].copy()
y_df = df[y_idx:b_idx].copy()
b_df = df[b_idx:].copy()
default_df = df[default_idx:r_idx].copy()
r_df["颜色"] = "红"
y_df["颜色"] = "黄"
b_df["颜色"] = "蓝"
df_all = pd.concat([r_df, y_df, b_df], ignore_index=True)
dfs = {"红": r_df, "黄": y_df, "蓝": b_df}
params = defaultdict(dict)
param_A_mat = np.array([])
param_B_mat = np.array([])
# AX + B
def solve(color: str, i: int):
slope, intercept, r_value, p_value, std_err =
stats.linregress(
dfs[color]["浓度（%）"] / 100, dfs[color][f"{i}nm"]
) # /100 是将浓度百分比转化为小数
params[color][i] = [slope, intercept]
def calc_Lab(XYZ1: np.array):
X0 = 94.83Y0 = 100.00
Z0 = 107.38
XYZ0 = np.array([X0, Y0, Z0])
X1 = XYZ1[0]
Y1 = XYZ1[1]
Z1 = XYZ1[2]
magic = 0.008856
if np.min(XYZ1 / XYZ0) > magic:
L = 116 * (Y1 / Y0) ** (1 / 3) - 16
a = 500 * ((X1 / X0) ** (1 / 3) - (Y1 / Y0) ** (1 / 3))
b = 200 * ((Y1 / Y0) ** (1 / 3) - (Z1 / Z0) ** (1 / 3))
else:
L = 903.3 * (Y1 / Y0)
a = 3893.5 * (X1 / X0 - Y1 / Y0)
b = 1557.4 * (Y1 / Y0 - Z1 / Z0)
# print(L,a,b,X1,Y1,Z1,"ASKD")
return np.array([L, a, b])
def calc_delta(XYZ1: np.array, XYZ2: np.array):
"""计算色差"""
Lab1 = calc_Lab(XYZ1)
Lab2 = calc_Lab(XYZ2)
return np.sum((Lab1 - Lab2) ** 2) ** 0.5
R_examples = defaultdict(list)
weights = [[], [], []]
df = pd.read_excel("附件3.xlsx", skiprows=1)
for id, dat in df.iterrows():
idx = int(dat["样本"])
for i in range(400, 720, 20):
R_examples[idx].append(dat[f"{i}nm"])
R_examples[idx] = np.array(R_examples[idx])
df = pd.read_excel("附件1.xlsx", skiprows=2)
for id, dat in df.iterrows():
weights[0].append(dat["S(λ)x ̅(λ)"])
weights[1].append(dat["S(λ)y ̅(λ)"])
weights[2].append(dat["S(λ)z ̅(λ)"])
weights[0] = np.array(weights[0])
weights[1] = np.array(weights[1])
weights[2] = np.array(weights[2])
weights = np.array(weights).T
# weights 是一个 (720-400)/20 行 3 列的矩阵
def get_D(idx: int):
R = R_examples[idx]
K = 2 * R**2 / (R**2 - 1)
return np.diag(K)
def ks_to_R(ks: np.array) -> np.array:
return 1 + ks - np.sqrt(ks * 2 + ks**2)
def R_to_ks(R: np.array) -> np.array:
return (1 - R) * (1 - R) / (2 * R)
def get_default_ks():
return np.array([default_df[f"{lmd}nm"][0] for lmd in
range(400, 720, 20)])
def get_special_ks(idx: int):
return R_to_ks(R_examples[idx])
def get_init_c(idx: int):
TE = weights.T
D = get_D(idx)
phi = param_B_mat
TED = TE @ D
ks_default = get_default_ks()
ks_special = get_special_ks(idx)
A = TED @ phi
A = np.linalg.inv(A)
C = A @ TED
C = C @ (ks_special - ks_default)
C = np.maximum(C, 0)
return C
xxs = np.arange(400,720,20)
def calc_XYZ(R: np.array) -> np.array:
"""返回 array[X,Y,Z]"""
k = 0.1
y_mat = R * weights.T
return k * np.array(
[
scipy.integrate.simpson(y=y_mat[i], x=xxs, dx=20)
for i in range(3)
]
)
def calc_R(C: np.array):
ks_d = np.array([default_df[f"{lmd}nm"][0] for lmd in
range(400, 720, 20)]).T
ks_all = np.sum(C * (param_A_mat * C + param_B_mat),
axis=1) + ks_d
return ks_to_R(ks_all)
XYZ_examples = {idx: calc_XYZ(R_examples[idx]) for idx in
R_examples.keys()}
def calc_loss(idx: int, R: np.array) -> float:
return calc_delta(XYZ_examples[idx], calc_XYZ(R))
def calc_loss_by_C(C: np.array, idx: int) -> float:
res = calc_loss(idx, calc_R(C))
return res
def solve_first(idx: int, C1: float, C2: float, C3: float):
C = np.array([C1, C2, C3])
return calc_loss_by_C(C, idx)
def constraints_0(params, idx):
"""限制色差 < 1"""
ret = 1 - calc_loss_by_C(params, idx)
return [ret, params[0], params[1], params[2]]
def is_same(C1: np.array, C2: np.array) -> bool:
"""如果相似则返回 True,否则返回 False"""
return (
np.max(np.abs((C1 - C2)) / C1) < 0.05
or np.max(np.abs(C1 - C2)) < 0.0001
or np.sum(np.abs((C1 - C2)) / C1) < 0.10
)
now_sols = []
def constraints_1(params):
"""限制与之前的方案都不相同"""
for sol in now_sols:
if is_same(params, sol):
return -1
return 1
def draw_graph(idx):
# 生成示例数据
l = 10
x = np.linspace(0, 0.10, l)
y = np.linspace(0, 0.10, l)
z = np.linspace(0, 0.10, l)
X, Y, Z = np.meshgrid(x, y, z)
w = np.zeros((l, l, l))
for i in range(l):
for j in range(l):
for k in range(l):
w[i][j][k] = calc_loss_by_C(
np.array([X[i][j][k], Y[i][j][k], Z[i][j][k]]),
idx
)
# Creating figure
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
img = ax.scatter3D(X, Y, Z, c=w, cmap=plt.hot())
ax.set_xlabel("C0")
ax.set_ylabel("C1")
ax.set_zlabel("C2")
fig.colorbar(img)
plt.title(f"目标{idx}")
plt.savefig(f"p3_{idx}.png")
ANS = defaultdict(list)
def solve_second():
global now_sols
failed = 0
fail_log = []
for idx in range(1, 11):
now_sols = []
initial_params = np.array(get_init_c(idx))
sol_num = 0
iter_num = 0
def find_sol(begin) -> bool:
print(f"Searching solution {iter_num} times for
{idx}...")
result = scipy.optimize.minimize(
calc_loss_by_C,
begin,
args=(idx,),
bounds=((0, 1), (0, 1), (0, 1)),
constraints=[
{"type": "ineq", "fun": constraints_0,
"args": (idx,)},
{"type": "ineq", "fun": constraints_1},
],
method="SLSQP",
)
if result.success:
ANS[idx].append((result.fun, result.x))
now_sols.append(result.x)
return result.success
if find_sol(initial_params):
sol_num += 1
while sol_num < 20 and iter_num < 200:
ret = find_sol(np.random.rand(len(initial_params))) #
[0,1) 的随机数入口
if ret:
sol_num += 1
iter_num += 1
if sol_num >= 10:
print(f"Finished {idx}. All success.\n\n")
else:
print(f"Failed to finish {idx}.")
failed += 1
fail_log.append((idx, sol_num))
print("All finished.")
print(f"Failed {failed} times. Fail log:")
for idx, num in fail_log:
print(f"Failed to finish {idx}. Found {num}
solutions.")
print("\n\n")
output = ""
for key in ANS.keys():
ANS[key].sort(key=lambda x: x[0])
output += f"目标{key}配方:\n"
for i in range(min(10, len(ANS[key]))):
output += f"{i + 1}: Delta =
{ANS[key][i][0].round(4)}, C = {[x.round(4) for
x in ANS[key][i][1]]}\n"
output += "\n"
with open("result2.txt", "w+", encoding="utf-8") as f:
f.write(output)
if __name__ == "__main__":
for i in range(400, 720, 20):
solve("红", i)
solve("黄", i)
solve("蓝", i)
param_A_mat = np.array(
[[params[c][i][0] for i in range(400, 720, 20)] for c
in colors]
).T
param_B_mat = np.array(
[[params[c][i][1] for i in range(400, 720, 20)] for c
in colors]
).T
# print(XYZ_examples)
solve_second()
for i in range(1,11):
draw_graph(i)
✡✝ ✆