用Python模拟湖羊养殖场：从数学建模到生产计划优化实战（附完整代码）

用Python模拟湖羊养殖场从数学建模到生产计划优化实战附完整代码湖羊养殖作为现代农业的重要组成部分其生产效率直接关系到养殖场的经济效益。传统的人工管理方式往往难以应对复杂的生产调度问题而数学建模与编程技术的结合为这一领域带来了全新的解决方案。本文将带领读者从零开始构建一个完整的湖羊养殖场模拟系统通过Python实现从基础数据建模到高级生产优化的全流程。1. 环境准备与基础数据建模在开始编码之前我们需要明确湖羊养殖的基本参数和规则。根据行业标准湖羊的生命周期可分为以下几个关键阶段# 湖羊生命周期各阶段时长天 LIFE_STAGES { mating: 20, # 交配期 pregnancy: 149, # 孕期 lactation: 40, # 哺乳期 resting: 20 # 休整期 } # 羊栏容量限制 PEN_CAPACITY { empty: 14, # 空怀休整期 ram: 4, # 非交配期种公羊 mating: 14, # 交配期母羊 pregnant: 8, # 怀孕期 lactating: 6, # 哺乳期(母羊羔羊) fattening: 14 # 育肥期羔羊 }建立基础母羊类时我们需要跟踪每只羊的当前状态和状态持续时间class Ewe: def __init__(self, ewe_id): self.id ewe_id self.state resting # 初始状态为休整期 self.state_days 0 # 当前状态已持续天数 self.pregnant False # 是否怀孕 self.lambs 0 # 产羔数量 def update(self): 更新母羊状态 self.state_days 1 if self.state mating and self.state_days LIFE_STAGES[mating]: self._transition_from_mating() elif self.state pregnancy and self.state_days LIFE_STAGES[pregnancy]: self._transition_from_pregnancy() # 其他状态转换逻辑... def _transition_from_mating(self): 从交配期转换 self.state pregnancy if random.random() 0.85 else resting self.state_days 0 self.pregnant self.state pregnancy2. 养殖场核心模拟引擎养殖场模拟引擎需要处理羊群的日常状态更新、羊栏分配以及生产记录等功能。我们采用面向对象的方式构建这个系统class SheepFarm: def __init__(self, num_ewes, num_rams, total_pens112): self.ewes [Ewe(i) for i in range(num_ewes)] self.rams num_rams self.pens { mating: [], # 交配羊栏 pregnant: [], # 怀孕羊栏 lactating: [], # 哺乳羊栏 fattening: [], # 育肥羊栏 resting: [] # 休整羊栏 } self.day 0 self.total_pens total_pens self.used_pens 0 self.records [] def daily_update(self): 每日更新养殖场状态 self.day 1 self._update_sheep_states() self._allocate_pens() self._record_daily_stats() def _update_sheep_states(self): 更新所有羊的状态 for ewe in self.ewes: ewe.update() def _allocate_pens(self): 分配羊栏 # 重置每日羊栏使用 self.pens {k: [] for k in self.pens} self.used_pens 0 # 按状态分组母羊 grouped_ewes defaultdict(list) for ewe in self.ewes: grouped_ewes[ewe.state].append(ewe) # 分配交配羊栏 mating_groups self._group_ewes( grouped_ewes.get(mating, []), PEN_CAPACITY[mating] ) self.pens[mating] mating_groups self.used_pens len(mating_groups) # 分配其他状态羊栏... def _group_ewes(self, ewes, capacity): 将母羊分组以适应羊栏容量 return [ewes[i:icapacity] for i in range(0, len(ewes), capacity)]3. 蒙特卡洛模拟与不确定性处理现实养殖中存在诸多不确定性因素我们需要通过蒙特卡洛模拟来评估不同生产计划的效果。以下是处理随机性的关键函数def simulate_pregnancy(ewe): 模拟受孕过程考虑85%的成功率 success random.random() 0.85 if success: # 孕期在147-150天之间波动 actual_pregnancy random.randint(147, 150) return actual_pregnancy return None def simulate_lambing(): 模拟产羔过程 lamb_count 0 rand_val random.random() # 产羔数量概率分布 if rand_val 0.7: # 70%概率产2只 lamb_count 2 elif rand_val 0.85: # 15%概率产1只 lamb_count 1 else: # 15%概率产3只或更多 lamb_count 3 int(random.expovariate(0.5)) # 模拟羔羊死亡率(3%) survived sum(1 for _ in range(lamb_count) if random.random() 0.03) return survived def adjust_fattening_period(lactation_days): 根据哺乳期调整育肥期 base_fattening 210 # 基础育肥期 deviation lactation_days - 40 # 与基准40天的偏差 return base_fattening - deviation * 24. 生产计划优化算法基于约束的优化算法可以帮助我们找到最佳的生产计划。以下是遗传算法的实现示例def genetic_algorithm_optimization(farm, generations100, pop_size50): 遗传算法优化生产计划 # 初始化种群 population [] for _ in range(pop_size): plan { batch_interval: random.randint(15, 30), ewes_per_batch: random.randint(30, 60), lactation_days: random.randint(35, 45) } population.append(plan) # 进化循环 for gen in range(generations): # 评估适应度 fitness [] for plan in population: score evaluate_plan(farm, plan) fitness.append(score) # 选择 selected tournament_selection(population, fitness) # 交叉和变异 new_population [] while len(new_population) pop_size: parent1, parent2 random.sample(selected, 2) child crossover(parent1, parent2) child mutate(child) new_population.append(child) population new_population # 返回最佳个体 best_idx np.argmax([evaluate_plan(farm, p) for p in population]) return population[best_idx] def evaluate_plan(farm, plan): 评估生产计划的适应度 # 模拟该计划下的养殖场运行 temp_farm copy.deepcopy(farm) total_days 365 * 2 # 模拟两年 output 0 for _ in range(total_days): temp_farm.daily_update() # 记录出栏数量 output temp_farm.records[-1][lambs_sold] # 惩罚羊栏不足 if temp_farm.used_pens temp_farm.total_pens: output - (temp_farm.used_pens - temp_farm.total_pens) * 3 return output5. 结果可视化与分析数据可视化是理解模拟结果的关键。我们使用Matplotlib创建多种图表来分析养殖场运营情况def visualize_farm_performance(farm, days365): 可视化养殖场性能指标 # 准备数据 pen_usage [r[used_pens] for r in farm.records[:days]] lamb_production [r[lambs_born] for r in farm.records[:days]] sales [r[lambs_sold] for r in farm.records[:days]] # 创建图表 plt.figure(figsize(15, 10)) # 羊栏使用情况 plt.subplot(3, 1, 1) plt.plot(pen_usage, label羊栏使用数) plt.axhline(yfarm.total_pens, colorr, linestyle--, label羊栏总数) plt.title(每日羊栏使用情况) plt.legend() # 羔羊生产情况 plt.subplot(3, 1, 2) plt.bar(range(days), lamb_production, label每日产羔数) plt.title(每日羔羊生产情况) # 出栏情况 plt.subplot(3, 1, 3) plt.plot(np.cumsum(sales), label累计出栏数) plt.title(累计出栏羊只数量) plt.legend() plt.tight_layout() plt.show() def analyze_optimization_results(results): 分析优化结果 df pd.DataFrame(results) # 计算关键指标 metrics { 年均出栏数: df[annual_output].mean(), 羊栏利用率: df[pen_utilization].mean(), 最大连续缺口: df[max_shortage].max(), 平均每日损失: df[daily_loss].mean() } # 创建热力图展示参数组合效果 pivot_table df.pivot_table( valuesannual_output, indexbatch_interval, columnsewes_per_batch, aggfuncmean ) plt.figure(figsize(10, 8)) sns.heatmap(pivot_table, annotTrue, fmt.0f, cmapYlGnBu) plt.title(不同参数组合下的年均出栏数) plt.xlabel(每批次基础母羊数) plt.ylabel(批次间隔(天)) plt.show() return metrics6. 完整系统集成与实战案例将上述模块整合为一个完整的养殖场管理系统我们可以进行实际案例分析class FarmManagementSystem: 湖羊养殖场综合管理系统 def __init__(self, config_fileNone): self.farm None self.current_plan None self.history [] if config_file: self.load_config(config_file) def initialize_farm(self, num_ewes, num_rams, total_pens): 初始化养殖场 self.farm SheepFarm(num_ewes, num_rams, total_pens) return self.farm def run_simulation(self, days365, planNone): 运行模拟 if not self.farm: raise ValueError(养殖场未初始化) if plan: self.current_plan plan self._apply_plan(plan) for _ in range(days): self.farm.daily_update() self.history.append({ days: days, plan: copy.deepcopy(self.current_plan), results: copy.deepcopy(self.farm.records) }) return self._summarize_results() def _apply_plan(self, plan): 应用生产计划 # 实现计划应用到养殖场的逻辑 pass def _summarize_results(self): 汇总模拟结果 if not self.history: return None last_run self.history[-1][results] df pd.DataFrame(last_run) summary { total_lambs_born: df[lambs_born].sum(), total_lambs_sold: df[lambs_sold].sum(), avg_pen_usage: df[used_pens].mean(), max_pen_usage: df[used_pens].max(), pen_shortage_days: sum(df[used_pens] self.farm.total_pens), annual_output: df[lambs_sold].sum() * (365 / len(df)) } return summary # 实战案例演示 if __name__ __main__: # 初始化管理系统 farm_system FarmManagementSystem() # 创建养殖场实例(112个标准羊栏) farm_system.initialize_farm(num_ewes500, num_rams10, total_pens112) # 运行基础模拟(1年) print(运行基础模拟...) base_results farm_system.run_simulation(days365) print(f基础方案年化出栏数: {base_results[annual_output]:.0f}只) # 优化生产计划 print(\n优化生产计划...) best_plan genetic_algorithm_optimization(farm_system.farm) print(f最佳生产计划参数: {best_plan}) # 应用优化后的计划 farm_system.initialize_farm(num_ewes500, num_rams10, total_pens112) optimized_results farm_system.run_simulation(days365, planbest_plan) print(f优化后年化出栏数: {optimized_results[annual_output]:.0f}只) # 结果可视化 visualize_farm_performance(farm_system.farm)7. 高级主题与扩展方向在基础系统之上我们可以进一步探索以下高级主题多目标优化同时考虑出栏数量、羊栏利用率、饲料成本等多个目标def multi_objective_optimization(farm, generations100): 多目标优化 # 使用NSGA-II算法 creator.create(FitnessMulti, base.Fitness, weights(1.0, -1.0)) creator.create(Individual, list, fitnesscreator.FitnessMulti) toolbox base.Toolbox() toolbox.register(attr_interval, random.randint, 15, 30) toolbox.register(attr_ewes, random.randint, 30, 60) toolbox.register(individual, tools.initCycle, creator.Individual, (toolbox.attr_interval, toolbox.attr_ewes), n1) toolbox.register(population, tools.initRepeat, list, toolbox.individual) def evaluate(individual): plan {batch_interval: individual[0], ewes_per_batch: individual[1]} output, loss evaluate_plan_multi(farm, plan) return output, loss toolbox.register(evaluate, evaluate) toolbox.register(mate, tools.cxBlend, alpha0.5) toolbox.register(mutate, tools.mutGaussian, mu0, sigma1, indpb0.2) toolbox.register(select, tools.selNSGA2) pop toolbox.population(n50) algorithms.eaMuPlusLambda(pop, toolbox, mu50, lambda_100, cxpb0.7, mutpb0.3, ngengenerations, statsNone, halloffameNone, verboseFalse) return pop机器学习预测使用历史数据训练预测模型优化生产决策def train_production_predictor(history_data): 训练生产预测模型 # 准备数据 X [] y [] for record in history_data: X.append([ record[plan][batch_interval], record[plan][ewes_per_batch], record[plan][lactation_days], record[results][avg_pen_usage] ]) y.append(record[results][annual_output]) X np.array(X) y np.array(y) # 训练随机森林回归模型 model RandomForestRegressor(n_estimators100, random_state42) model.fit(X, y) return model def predict_optimal_parameters(model, farm): 使用模型预测最优参数 # 生成候选参数 candidates [] for interval in range(15, 31): for ewes in range(30, 61): candidates.append([interval, ewes, 40, 0]) # 0为占位符 candidates np.array(candidates) # 预测各候选的表现 predictions model.predict(candidates) best_idx np.argmax(predictions) return { batch_interval: candidates[best_idx, 0], ewes_per_batch: candidates[best_idx, 1] }分布式模拟使用多进程加速蒙特卡洛模拟def parallel_monte_carlo(farm, plan, simulations1000, days365): 并行蒙特卡洛模拟 def worker(_): temp_farm copy.deepcopy(farm) temp_farm.apply_plan(plan) for _ in range(days): temp_farm.daily_update() return temp_farm.records[-1][total_lambs_sold] with Pool(processescpu_count()) as pool: results pool.map(worker, range(simulations)) return { mean_output: np.mean(results), std_output: np.std(results), min_output: np.min(results), max_output: np.max(results), distribution: results }