数学建模C题Python实现:从线性规划到粒子群与聚类预测
2026/9/13 5:57:05 网站建设 项目流程

简介:针对2024年全国大学生数学建模竞赛C题,这份资源围绕农作物种植策略提供了一整套可落地的分析方案,涵盖建模思路、程序源码和优秀论文,适合参赛学生、赛事指导老师以及想入门数学建模的读者使用。代码按赛题问题拆分:前两问均有独立求解脚本,第三问综合运用线性回归、非线性回归、Pearson相关性检验和KMeans聚类算法,对预期销售量、种植成本和销售单价进行统计建模,并对地块—大棚做分类处理,还附带整理好的地块信息、作物信息、相关性信息及2023年汇总补充数据,基本还原了从数据处理到结果输出的完整链路。资源共28个文件,包括12个Excel数据文件、8个Python程序脚本、2个PDF赛题与论文、2个Markdown说明文档及4个代码备份文件,压缩包整体仅2.21MB,轻量易用,方便按需查阅。当前已有133人学习,适合需要复现国赛获奖思路、参考建模代码或快速了解农作物种植优化问题的读者。

1. 从赛题到代码:一份完整 C 题方案的解构

2024 年国赛 C 题把「农作物种植策略」做成了一个相当实际的优化问题:村里有平旱地、梯田、山坡地、水浇地,还有大棚,十几种作物要在 2024—2030 年间逐年排产,产量、售价、成本还会波动。这份资源里真正值钱的不是那篇论文,而是六个按问题顺序编号的 Python 求解文件:从 Q1_1.py、Q1_2.py 的线性规划起步,到 Q2.py 的随机优化,再到问题三的 pearson.py、线性/非线性回归脚本和 KMeans 聚类,完整覆盖了确定性优化、不确定性建模、相关性检验和聚类四类数学建模模型。它适合两类人:正在备战数学建模竞赛、想找一套可复现代码路径的选手;以及用 Python 处理生产调度、想把优化工具链串起来落地的开发者。下文按赛题推进顺序逐层拆这套代码。

2. 地块与作物的约束系统:Q1_1.py 和 Q1_2.py 的线性规划实现

2.1 从赛题数据到决策变量

问题一的核心假设是未来产量、售价、成本相对稳定,因此方案可以在同一套约束下逐年复用。数据基础是资源包里的「地块信息.xlsx」和「作物信息.xlsx」,前者记录每块地的类型与面积,后者记录每种作物的亩产量、种植成本、销售单价。建模的第一步是把自然语言描述转成线性规划的标准形式。

决策变量定义为:

# x[i * n_crop + j] 表示第 i 块地上种植第 j 种作物的面积 # 索引展开方式:地块优先,同一地块内的作物连续排列 n_land = len(land) n_crop = len(crop) x = [0] * (n_land * n_crop)

变量采用「地块优先 + 块内作物连续」的扁平排列,是因为scipy.optimize.linprog只接受一维决策变量,后续构造约束矩阵时需要按i * n_crop + j这个索引规则定位。如果换成二维变量,矩阵构造会直观一些,但还原成线性规划标准形式时需要多一层 reshape,反而容易出错。

2.2 目标函数与约束矩阵的构造

目标函数是总利润最大化。单块地上某作物的利润按「面积 × 亩产量 × 销售单价 − 面积 × 种植成本」计算,因此目标系数向量里的每一项都是单位面积净利润。linprog只做最小化,所以把利润向量整体取负。

import pandas as pd from scipy.optimize import linprog def build_q1_model(land, crop): n_land, n_crop = len(land), len(crop) c = [] for _, ld in land.iterrows(): for _, cr in crop.iterrows(): profit = cr["亩产量"] * cr["销售单价"] - cr["种植成本"] c.append(-profit) # 最小化负利润 = 最大化利润 A_ub, b_ub = [], [] # 约束:每块地上所有作物的种植面积之和 <= 该地块面积 for i in range(n_land): row = [0] * (n_land * n_crop) for j in range(n_crop): row[i * n_crop + j] = 1 A_ub.append(row) b_ub.append(land.loc[i, "面积"]) # 变量下界为 0,上界为该地块面积,保证单个变量不会超过地块总面积 bounds = [(0, land.loc[i, "面积"]) for i in range(n_land) for _ in range(n_crop)] res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method="highs") return res

这里的method="highs"是 SciPy 1.7 之后默认启用的 HiGHS 求解器,适合中等规模的线性规划;如果直接用旧代码里的method="simplex",在数千个变量的场景下迭代步数和数值误差都会明显变大。约束矩阵 A_ub 的每一行只在一块地对应的 n_crop 个位置上填 1,本质是块内加和约束,这在实际求解中比单独对每块地每种作物分别设限更接近赛题的「可以少种但不可以多种」逻辑。

2.3 Q1_2.py 在 Q1_1 之上叠加的时间维约束

Q1_1.py 求解的是稳定假设下的一期最优布局,而 Q1_2.py 解决的是问题一的第二问:方案要逐年滚动到 2030 年,同一块地上不能连续多年种植同一种作物。此时决策变量需要扩成三维,x[t][i][j]表示第 t 年第 i 块地种第 j 种作物的面积,并额外追加重茬约束。

T = 7 # 2024 到 2030 def add_rotation_constraint(A_ub, b_ub, n_land, n_crop, T, forbid_pairs): for t in range(1, T): # 从第二年开始检查上一年的种植情况 for i in range(n_land): for j in forbid_pairs.get(i, []): row = [0] * (T * n_land * n_crop) # 上一年的面积 + 这一年的面积 <= 上一年面积的阈值 row[(t-1) * n_land * n_crop + i * n_crop + j] = 1 row[t * n_land * n_crop + i * n_crop + j] = 0.2 A_ub.append(row) b_ub.append(0) return A_ub, b_ub

这段约束的含义是:如果第 t−1 年某地块种了易重茬障碍的作物,第 t 年的同作物面积不能超过上一年的 0.2 倍,相当于允许少量留茬,但不允许连续性大规模种植。forbid_pairs需要根据作物信息里的「重茬敏感」标记手动配置,不同作物对连作的容忍度差异很大,这属于赛题数据里不会直接给出、但实际求解时必须补全的工程细节。

Q1_1 和 Q1_2 的关系可以简单理解成:前者是「一年定终身」,后者是「七年滚动排产」。从代码量看,Q1_2 只多了时间索引和轮作约束,但求解规模扩大了近 7 倍,这也是我建议用highs而不是老求解器的直接原因。

3. 不确定性建模:Q2.py 与改进粒子群优化算法(IPSO)

3.1 为什么线性规划在问题二失效

问题二把产量、售价、成本同时变成随机变量。如果继续用线性规划,就得把每个随机参数替换成期望值,但期望值解在真实波动下往往违背面积约束或产生明显亏损。原因是线性规划天然假设参数确定,一旦目标函数里的系数带随机性,最优解对参数扰动非常敏感,尤其是产量和价格呈负相关时(丰年产量高但价格低),直接用均值建模会系统性高估收益。

资源包里的论文标题给出了作者团队的解法:基于粒子群优化算法的华北农作物种植策略模型。代码里对应的是 Q2.py,它引入三个情景——丰年、平年、歉年,每个情景有独立的产量折减系数和价格波动系数,目标函数改为多年期望利润。

import numpy as np from sko.PSO import PSO def expected_profit(x, land, crop, years=7): X = x.reshape(years, len(land), len(crop)) scenarios = {"丰年": (1.10, 0.95), # 产量上浮 10%,价格回落 5% "平年": (1.00, 1.00), "歉年": (0.85, 1.08)} # 产量下降 15%,价格上浮 8% prob = {"丰年": 0.3, "平年": 0.5, "歉年": 0.2} total = 0.0 for t in range(years): for name, (y_f, p_f) in scenarios.items(): profit_t = 0.0 for i in range(len(land)): for j in range(len(crop)): profit_t += ( X[t, i, j] * crop.loc[j, "亩产量"] * y_f * crop.loc[j, "销售单价"] * p_f - X[t, i, j] * crop.loc[j, "种植成本"] ) total += prob[name] * profit_t return -total + penalty(X, land) # 负号让 PSO 从最大化转最小化

expected_profit的目标值是对七个年份、三个情景的加权求和。penalty(X, land)是面积越界和轮作约束的罚函数实现,我通常把越界面积平方乘一个大系数加回去,而不是直接截断变量,因为截断会破坏 PSO 搜索的连续性。

3.2 粒子群参数怎么设置

pso = PSO(func=expected_profit, n_dim=years * len(land) * len(crop), pop=40, max_iter=200, lb=0, ub=land["面积"].max(), w=0.8, c1=1.4, c2=1.4) pso.run() best_x = pso.gbest_x

这里的pop=40是粒子数,max_iter=200是最大迭代轮数,lb=0ub=land["面积"].max()给出变量上下界。w=0.8是惯性权重,控制粒子保持原有速度的程度;c1=1.4是自我认知系数,c2=1.4是社会认知系数,两者相等会让粒子在局部搜索和全局收敛之间保持平衡。如果赛题时间紧张,可以把max_iter降到 100,但代价是结果不稳定。对比线性规划的单次确定解,PSO 要跑 200 × 40 次目标函数评估,Q2.py 运行十几分钟是正常的。

与 Lingo 或 Matlab 里的ga工具箱相比,sko.PSO不用写完整的约束矩阵,只需要把硬约束塞进目标函数。代价是每次运行结果有随机性,我建议运行三次取最优,同时记录每次的期望利润浮动范围,作为模型鲁棒性的参考指标。

4. 销售量预测链路:pearson.py、Linear_Regression.py 与 nonlinear_regression.py 的组合

4.1 先检验相关性,再决定用哪种模型

问题三需要预测每种作物的预期销售量,但不同作物的销售量随时间呈现出完全不同的形态:有的近似直线上升,有的先增后平,有的周期性波动。pearson.py的作用是在做回归之前先回答一个问题:这组年份和销售量之间到底存不存在线性相关?如果皮尔逊系数绝对值低于 0.8 或 p 值大于 0.05,直接套线性回归就是无效建模,必须进入非线性分支。

from scipy.stats import pearsonr def check_linearity(data, crop_name, year_col="年份", sales_col="预期销售量"): sub = data[data["作物名称"] == crop_name] r, p = pearsonr(sub[year_col], sub[sales_col]) # r 在 [-1, 1],p 小于 0.05 说明相关性显著 return r, p for crop_name in data["作物名称"].unique(): r, p = check_linearity(data, crop_name) branch = "linear" if abs(r) > 0.8 and p < 0.05 else "nonlinear" print(f"{crop_name}: r={r:.3f}, p={p:.3f}, -> {branch}")

赛题场景里样本量通常很小,一个作物往往只有几年数据,p 值的参考意义要打折。我的处理方式是同时看 r 和残差图:如果 r 大但数据点明显弯成弧形,仍然归入非线性分支;如果 r 中等但 p 值极小,说明线性趋势存在但噪声占比高,应该线性与非线性各跑一版对比。

4.2 线性回归用 statsmodels 而不是 sklearn

Linear_Regression.py做的是普通最小二乘回归,但有一个容易被忽略的细节:sklearn 的LinearRegression不输出系数显著性、R² 和 F 统计量,用来预测没问题,用来解释模型就缺了诊断信息。比赛论文里需要写「模型显著性强」,所以代码里更适合用 statsmodels。

import statsmodels.api as sm def fit_ols(x, y): X = sm.add_constant(x) # 加截距项 model = sm.OLS(y, X).fit() return model.params, model.rsquared, model.pvalues

add_constant是必要步骤,缺失时模型会被强行过原点,这对销售量数据是错误假设。model.rsquared就是拟合优度,论文里报这个值;model.pvalues用于判断年份变量是否显著,通常小于 0.05 才会保留。

4.3 非线性回归的曲线族选择

nonlinear_regression.py对应三种最常用的曲线形态:

from scipy.optimize import curve_fit def quadratic(x, a, b, c): return a * x**2 + b * x + c def exp_growth(x, a, b, c): return a * np.exp(b * (x - x[0])) + c def log_growth(x, a, b, c): return a * np.log(b * (x - x[0]) + 1) + c

quadratic适合先升后降或持续加速的作物;exp_growth适合处于快速扩张期的品种;log_growth适合市场趋于饱和、增长逐步放缓的品种。curve_fit返回的popt是最优参数数组,pcov是参数协方差矩阵,对角线开方就是参数标准误。如果拟合后某个参数的标准误比参数本身还大,说明曲线选型有问题,应该换一个函数再拟合,而不是直接采用结果。

线性回归与非线性回归之间没有绝对优劣,关键在前期判断。pearson.py 负责把「能不能用线性」这个问题量化,后两个脚本负责在各自分支里产出预测值,这条链路先检验再建模的顺序,比拿到数据直接 fit 更经得起复现和答辩。

5. 地块与大棚的差异化处理:Q3_cluster.py 的 KMeans 与 Q3_ABCD_EF.py 的手动二分类

5.1 为什么问题三要先聚类再预测

问题三的数据里,每种作物的预期销售量、销售单价、种植成本各有各的波动节奏,把所有作物放进同一个回归模型会造成严重的异方差问题,比如大宗粮食作物数量级是蔬菜的几十倍,一起拟合时蔬菜的预测误差会被粮食作物的残差稀释掉。Q3_cluster.py的思路是先按销售行为把作物分成三类,让同类作物拥有相近的参数尺度,再用类别数据单独做回归。

因为资源附带了「相关性信息_类别1.xlsx」「相关性信息_类别2.xlsx」「相关性信息_类别3.xlsx」三个输出文件,说明 Q3_cluster.py 实际使用的就是三分类。聚类特征选取的是预期销售量、种植成本、销售单价三个维度的均值向量。

from sklearn.cluster import KMeans from sklearn.preprocessing import StandardScaler def cluster_crops(df, feature_cols, n_clusters=3): X = df[feature_cols].values X = StandardScaler().fit_transform(X) # 消除量纲差异 km = KMeans(n_clusters=n_clusters, random_state=42, n_init=10) df["类别"] = km.fit_predict(X) return df

StandardScaler这一步非常关键。预期销售量的方差远大于种植成本,如果不做标准化,聚类结果基本只被销售量主导,类别 2 和类别 3 的区别会变成单纯的「量级差异」而不是「行为差异」。n_init=10是让 KMeans 用 10 个不同的随机中心点各跑一次,取最优结果,避免初始化导致的局部最优。类别数 3 不是随便定的,可以用肘部法则验证:画类别数从 2 到 6 的 SSE 曲线,选拐点。

5.2 类别特征如何对应到回归

聚类完成后,输出文件「相关性信息_类别x.xlsx」里每一类包含一组作物。接下来的回归不是对所有作物混在一起做,而是每个类别单独执行第 4 章的 pearson + 回归流程。类别之间的差异往往体现在曲线形态上:类别 1 可能是量价齐升的经济作物,类别 2 是价格平稳的粮食作物,类别 3 是价格随季节波动大的果蔬类。这样分类后,回归得到的 R² 会比整体拟合高出不少,论文里说明「先聚类再分类预测」也就是在讲这一步。

5.3 露天地块与大棚的拆分

Q3_ABCD_EF.py是一个手动二分类脚本,逻辑比 KMeans 更直接:文件名里的 A/B/C/D 对应赛题中的四类露天地块,E/F 对应大棚地块。大棚作物的生长环境可控,销售量与温度的耦合程度低于露天作物,混合建模会引入难以解释的干扰项。常见做法是先用地块类型字段把数据拆成「露天」和「大棚」两份,再分别走预测流程。

def split_land_category(land): outdoor = ["平旱地", "梯田", "山坡地", "水浇地"] land["地块大类"] = land["地块类型"].apply( lambda t: "露天" if t in outdoor else "大棚") return land[land["地块大类"] == "露天"], land[land["地块大类"] == "大棚"]

这里apply配合 lambda 做映射,比逐行for循环快一个数量级。拆分之后两份数据各自做一次聚类和回归,输出对应为「result3(地块-大棚).xlsx」和「result3(聚类).xlsx」。两个结果文件不是替代关系:前者回答地块层面的分配约束,后者回答作物层面的销售量预测,最终方案要把两者合并:先用聚类结果算出每类作物的预测销售量,再用地块-大棚分组把销售量落实到具体地块面积上。

6. 环境复现、运行顺序与结果交叉验证

6.1 依赖库安装与版本坑

资源标注要求 Python 3.11,核心依赖是 pandas、numpy、sko、matplotlib。其中sko是 scikit-opt 的缩写,PyPI 包名是scikit-opt,直接pip install sko会装错包。

pip install pandas numpy scipy scikit-learn scikit-opt matplotlib openpyxl

openpyxl是 pandas 读写 .xlsx 的底层引擎,不装会在pd.read_excel时报 ImportError。scikit-opt 目前对 numpy 2.x 的兼容性有历史问题,建议先安装 numpy==1.26.4 再装 scikit-opt,避免出现module 'numpy' has no attribute 'float'这类报错。statsmodels 只在回归脚本里用到,如果只跑 Q3 的线性回归分支,可以用 sklearn 替代,但建议还是装上。

6.2 脚本执行顺序

顺序脚本输入输出
1Q1_1.py地块信息.xlsx、作物信息.xlsxresult1_1.xlsx
2Q1_2.py地块信息.xlsx、作物信息.xlsxresult1_2.xlsx
3Q2.py地块信息.xlsx、作物信息.xlsx、汇总补充信息.xlsxresult2.xlsx
4Q3_cluster.py相关性信息.xlsx相关性信息_类别1/2/3.xlsx
5pearson.py相关性信息_类别1/2/3.xlsx相关系数与分支标记
6Linear_Regression.py + nonlinear_regression.py类别数据预期销售量预测值
7Q3_ABCD_EF.py地块信息.xlsx、作物信息.xlsx、预测值result3(地块-大棚).xlsx、result3(聚类).xlsx

6.3 一个值得养成的验证技巧

def verify_area_constraint(solution, land, eps=1e-6): # solution 是 predict 或 优化结果,行是地块,列是作物 total_planted = solution.sum(axis=1) # 每块地的总种植面积 exceed = total_planted - land["面积"] return exceed[exceed > eps]

这套资源的问题在于脚本多、中间文件多,最容易出的错不是模型错误,而是某个中间文件被上一步覆盖,导致后续脚本读到错位数据。我的习惯是跑完 Q3_cluster.py 和 Q3_ABCD_EF.py 之后,第一时间用上面的函数检查「result3(聚类).xlsx」里每块地的种植面积总和是否仍然满足约束。如果超过,问题大概率出在聚类后的作物归属发生了错位;如果没超过,但方案利润异常低,再回去查回归脚本里是不是把某个类别的预测销售量串到了另一个类别。复现这套代码时把握好「输入文件路径正确」「输出文件不互相覆盖」「约束校验通过」这三条,整套流程基本不会出问题。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询