241 lines
13 KiB
Python
241 lines
13 KiB
Python
import pulp
|
||
import pandas as pd
|
||
from tabulate import tabulate
|
||
import os
|
||
|
||
|
||
def calculate_nutrients_from_formula(formula_proportions, ingredients):
|
||
"""
|
||
计算市面配方营养总量，作为需求校准
|
||
:param formula_proportions: 配方比例字典，如 {'玉米': 0.59, '豆粕': 0.15}
|
||
:param ingredients: 原料数据字典，含营养值/价格
|
||
:return: 营养总量字典，如 {'蛋白': 17.76, '能量': 11.56}
|
||
"""
|
||
nutrients = list(set(key for ing in ingredients.values() for key in ing if key not in ['价格']))
|
||
total_nutrients = {}
|
||
for nut in nutrients:
|
||
total = sum(formula_proportions.get(i, 0) * ingredients.get(i, {}).get(nut, 0) for i in formula_proportions)
|
||
total_nutrients[nut] = round(total, 2)
|
||
return total_nutrients
|
||
|
||
|
||
def calculate_cost_from_formula(formula_proportions, ingredients):
|
||
"""
|
||
计算市面配方成本，作为参考
|
||
:param formula_proportions: 配方比例字典，如 {'玉米': 0.59}
|
||
:param ingredients: 原料数据字典，含价格
|
||
:return: 成本字典，含总价格和一斤价格
|
||
"""
|
||
total_price = sum(
|
||
formula_proportions.get(i, 0) * ingredients.get(i, {}).get('价格', 0) for i in formula_proportions)
|
||
return {
|
||
'总价格': round(total_price, 2),
|
||
'一斤饲料价格': round(total_price / 2000, 2) # 吨=2000斤
|
||
}
|
||
|
||
|
||
def optimize_feed(requirements, ingredients, result_fields, optimization_type='min'):
|
||
"""
|
||
优化饲料配方，计算最低/最高目标值的原料比例
|
||
:param requirements: 猪营养需求字典，如 {'蛋白_下限': 17, '能量_下限': 12.5, '纤维_上限': 6}
|
||
:param ingredients: 原料数据字典，如 {'构树叶_发酵': {'蛋白': 15, '能量': 9, '纤维': 15, '价格': 600}}
|
||
:param result_fields: 结果字段列表，如 ['比例', '成本']
|
||
:param optimization_type: 'min'（最小化，如成本）或 'max'（最大化，如消化率）
|
||
:return: 字典，含结果字段和值，如 {'构树叶_发酵': {'比例': 14.29, '成本': 85.74}, '总价格': 1989.12}
|
||
"""
|
||
sense = pulp.LpMinimize if optimization_type == 'min' else pulp.LpMaximize
|
||
model = pulp.LpProblem("Feed_Optimization", sense)
|
||
|
||
# 决策变量：每种原料比例（0-1，单位：比例）
|
||
x = pulp.LpVariable.dicts("比例", ingredients, lowBound=0, upBound=1)
|
||
|
||
# 目标函数：最小化总价格（纯价格，无额外权重）
|
||
model += pulp.lpSum([x[i] * ingredients[i]['价格'] for i in ingredients]), "总价格"
|
||
|
||
# 约束1：总比例=100%
|
||
model += pulp.lpSum([x[i] for i in ingredients]) == 1, "总比例"
|
||
|
||
# 约束2：营养需求（下限/上限）
|
||
for req, value in requirements.items():
|
||
if '下限' in req:
|
||
nutrient = req.replace('_下限', '')
|
||
# 营养总和≥下限（如蛋白≥17%）
|
||
model += pulp.lpSum([x[i] * ingredients[i][nutrient] for i in ingredients]) >= value, req
|
||
elif '上限' in req:
|
||
nutrient = req.replace('_上限', '')
|
||
# 营养总和≤上限（如纤维≤6%）
|
||
model += pulp.lpSum([x[i] * ingredients[i][nutrient] for i in ingredients]) <= value, req
|
||
|
||
# 求解：用CBC极限精度，深度优先，限时10秒
|
||
try:
|
||
solver = pulp.PULP_CBC_CMD(msg=True,
|
||
options=['primalT', '1e-12', 'dualT', '1e-12', 'maxIt', '10000000', 'presolve', 'on',
|
||
'strategy', '2', 'randomC', '123', 'sec', '10'])
|
||
model.solve(solver)
|
||
except AttributeError:
|
||
print("警告：CBC不可用，使用默认求解器，精度可能稍低")
|
||
model.solve()
|
||
|
||
# 检查是否找到最优解
|
||
if pulp.LpStatus[model.status] != 'Optimal':
|
||
print("约束值（调试用）：")
|
||
for name, constraint in model.constraints.items():
|
||
print(f"{name}: {constraint.value()}")
|
||
return {"错误": "无解！检查数据（纤维/单宁超？蛋白/氨基酸不足？）"}
|
||
|
||
# 构建结果
|
||
result = {}
|
||
total_price = pulp.value(model.objective)
|
||
result['总价格'] = round(total_price, 2)
|
||
result['一斤饲料价格'] = round(total_price / 2000, 2) # 吨=2000斤
|
||
|
||
# 每种原料的比例（%）和成本（元/吨）
|
||
for i in ingredients:
|
||
result[i] = {field: round(pulp.value(x[i]) * 100, 2) if field == '比例' else round(
|
||
pulp.value(x[i]) * ingredients[i]['价格'], 2) for field in result_fields}
|
||
|
||
# 计算营养总量
|
||
nutrients = list(set(key for ing in ingredients.values() for key in ing if key not in ['价格']))
|
||
result['营养总量'] = {}
|
||
for nut in nutrients:
|
||
total = sum(pulp.value(x[i]) * ingredients[i].get(nut, 0) for i in ingredients)
|
||
result['营养总量'][nut] = round(total, 2)
|
||
|
||
# 检查超标（宁可无解也不超）
|
||
for req, value in requirements.items():
|
||
if '上限' in req:
|
||
nut = req.replace('_上限', '')
|
||
if result['营养总量'][nut] > value + 0.01: # 容忍小误差
|
||
return {"错误": f"{nut}超上限: {result['营养总量'][nut]} > {value} (可能浮点误差)"}
|
||
|
||
# 存CSV到results/
|
||
os.makedirs('results', exist_ok=True)
|
||
df_formula = pd.DataFrame(
|
||
[(k, v['比例'], v['成本']) for k, v in result.items() if k not in ['总价格', '一斤饲料价格', '营养总量']],
|
||
columns=['原料', '比例', '成本'])
|
||
df_formula.to_csv('results/优化配方结果.csv', index=False)
|
||
|
||
# 营养偏差表
|
||
nutrient_table = []
|
||
for nut, value in result['营养总量'].items():
|
||
req_key = f"{nut}_下限" if f"{nut}_下限" in requirements else f"{nut}_上限"
|
||
req_type = '下限' if f"{nut}_下限" in requirements else '上限'
|
||
req_value = requirements.get(req_key, None)
|
||
deviation = value - req_value if req_key.endswith('_下限') else req_value - value
|
||
deviation_pct = (deviation / req_value * 100) if req_value != 0 else 0
|
||
status = '达标' if (req_key.endswith('_下限') and value >= req_value) or (
|
||
req_key.endswith('_上限') and value <= req_value) else '不达标'
|
||
nutrient_table.append([nut, value, req_value, req_type, deviation, deviation_pct, status])
|
||
df_nutrients = pd.DataFrame(nutrient_table,
|
||
columns=['营养', '实际值', '需求值', '需求类型', '偏差', '偏差比例(%)', '状态'])
|
||
df_nutrients.to_csv('results/营养总量与需求偏差.csv', index=False)
|
||
|
||
return result
|
||
|
||
|
||
if __name__ == "__main__":
|
||
# 原料数据（单位：%或MJ/kg，价格元/吨，2025年9月农业农村部/饲料市场信息网）
|
||
ingredients = {
|
||
'构树叶_未加工': {'蛋白': 13, '能量': 8.5, '纤维': 25, '赖氨酸': 0.4, '蛋氨酸': 0.08, '钙': 1.8, '磷': 0.15,
|
||
'单宁': 2.5, '价格': 500},
|
||
'构树叶_发酵': {'蛋白': 15, '能量': 9, '纤维': 15, '赖氨酸': 0.5, '蛋氨酸': 0.1, '钙': 1.5, '磷': 0.2,
|
||
'单宁': 1.5, '价格': 600},
|
||
'玉米胚芽_压榨': {'蛋白': 11, '能量': 13, '纤维': 7, '赖氨酸': 0.4, '蛋氨酸': 0.2, '钙': 0.1, '磷': 0.3,
|
||
'单宁': 0, '价格': 1700},
|
||
'麦麸': {'蛋白': 15, '能量': 10, '纤维': 12, '赖氨酸': 0.4, '蛋氨酸': 0.2, '钙': 0.1, '磷': 1.0, '单宁': 0,
|
||
'价格': 1300},
|
||
'米糠': {'蛋白': 13, '能量': 11, '纤维': 10, '赖氨酸': 0.4, '蛋氨酸': 0.2, '钙': 0.1, '磷': 1.5, '单宁': 0,
|
||
'价格': 1100},
|
||
'牧草_苜蓿干': {'蛋白': 17, '能量': 8.5, '纤维': 20, '赖氨酸': 0.6, '蛋氨酸': 0.2, '钙': 1.2, '磷': 0.3,
|
||
'单宁': 0.5, '价格': 800},
|
||
'玉米': {'蛋白': 9, '能量': 13.5, '纤维': 2.5, '赖氨酸': 0.3, '蛋氨酸': 0.17, '钙': 0.02, '磷': 0.3, '单宁': 0,
|
||
'价格': 2367},
|
||
'豆粕': {'蛋白': 46, '能量': 10.5, '纤维': 4, '赖氨酸': 2.8, '蛋氨酸': 0.6, '钙': 0.3, '磷': 0.7, '单宁': 0,
|
||
'价格': 3300},
|
||
'植物油_豆油': {'蛋白': 0, '能量': 36, '纤维': 0, '赖氨酸': 0, '蛋氨酸': 0, '钙': 0, '磷': 0, '单宁': 0,
|
||
'价格': 7500},
|
||
'赖氨酸_L盐': {'蛋白': 78, '能量': 0, '纤维': 0, '赖氨酸': 78, '蛋氨酸': 0, '钙': 0, '磷': 0, '单宁': 0,
|
||
'价格': 15000},
|
||
'蛋氨酸_DL': {'蛋白': 99, '能量': 0, '纤维': 0, '赖氨酸': 0, '蛋氨酸': 99, '钙': 0, '磷': 0, '单宁': 0,
|
||
'价格': 25000},
|
||
'预混料_维生素微量元素': {'蛋白': 0, '能量': 0, '纤维': 0, '赖氨酸': 0, '蛋氨酸': 0, '钙': 10, '磷': 5,
|
||
'单宁': 0, '价格': 3000},
|
||
'鱼粉': {'蛋白': 60, '能量': 12, '纤维': 1, '赖氨酸': 3.5, '蛋氨酸': 1.5, '钙': 5, '磷': 3, '单宁': 0,
|
||
'价格': 5000},
|
||
'骨粉': {'蛋白': 0, '能量': 0, '纤维': 0, '赖氨酸': 0, '蛋氨酸': 0, '钙': 30, '磷': 15, '单宁': 0,
|
||
'价格': 1500},
|
||
'食盐': {'蛋白': 0, '能量': 0, '纤维': 0, '赖氨酸': 0, '蛋氨酸': 0, '钙': 0, '磷': 0, '单宁': 0,
|
||
'价格': 500}
|
||
}
|
||
|
||
# 育肥猪需求（30-115kg，NY/T 65-2004+丹麦标准，单位：%或MJ/kg）
|
||
requirements = {
|
||
'蛋白_下限': 17,
|
||
'能量_下限': 12.5,
|
||
'纤维_上限': 6,
|
||
'赖氨酸_下限': 0.8,
|
||
'蛋氨酸_下限': 0.25,
|
||
'钙_下限': 0.5,
|
||
'磷_下限': 0.4,
|
||
'单宁_上限': 0.5
|
||
}
|
||
|
||
# 优化配方
|
||
result = optimize_feed(requirements, ingredients, ['比例', '成本'], 'min')
|
||
|
||
# 美化输出：优化配方
|
||
print("\n=== 优化配方结果 ===")
|
||
table_data = [(k, v['比例'], v['成本']) for k, v in result.items() if
|
||
k not in ['总价格', '一斤饲料价格', '营养总量']]
|
||
headers = ['原料', '比例 (%)', '成本 (元/吨)']
|
||
print(tabulate(table_data, headers=headers, tablefmt='grid', floatfmt='.2f'))
|
||
print(f"总价格: {result['总价格']:.2f} 元/吨")
|
||
print(f"一斤饲料价格: {result['一斤饲料价格']:.2f} 元")
|
||
|
||
# 营养总量与需求偏差
|
||
print("\n=== 营养总量与需求偏差 ===")
|
||
nutrient_table = []
|
||
for nut, value in result['营养总量'].items():
|
||
req_key = f"{nut}_下限" if f"{nut}_下限" in requirements else f"{nut}_上限"
|
||
req_type = '下限' if f"{nut}_下限" in requirements else '上限'
|
||
req_value = requirements.get(req_key, None)
|
||
deviation = value - req_value if req_key.endswith('_下限') else req_value - value
|
||
deviation_pct = (deviation / req_value * 100) if req_value != 0 else 0
|
||
status = '达标' if (req_key.endswith('_下限') and value >= req_value) or (
|
||
req_key.endswith('_上限') and value <= req_value) else '不达标'
|
||
nutrient_table.append([nut, value, req_value, req_type, deviation, deviation_pct, status])
|
||
print(tabulate(nutrient_table, headers=['营养', '实际值', '需求值', '需求类型', '偏差', '偏差比例(%)', '状态'],
|
||
tablefmt='grid', floatfmt='.2f'))
|
||
|
||
# 市面配方（2025年养猪网，21-35kg调整至30-115kg）
|
||
market_formula = {
|
||
'玉米': 0.59,
|
||
'麦麸': 0.13,
|
||
'豆粕': 0.15,
|
||
'鱼粉': 0.06,
|
||
'骨粉': 0.015,
|
||
'食盐': 0.005
|
||
}
|
||
|
||
# 市面配方营养和成本
|
||
market_nutrients = calculate_nutrients_from_formula(market_formula, ingredients)
|
||
market_cost = calculate_cost_from_formula(market_formula, ingredients)
|
||
|
||
print("\n=== 市面配方营养总量 ===")
|
||
print(tabulate([(k, v) for k, v in market_nutrients.items()], headers=['营养', '值'], tablefmt='grid',
|
||
floatfmt='.2f'))
|
||
pd.DataFrame([(k, v) for k, v in market_nutrients.items()], columns=['营养', '值']).to_csv(
|
||
'results/市面配方营养总量.csv', index=False)
|
||
|
||
print("\n=== 市面配方成本与偏差 ===")
|
||
cost_table = [
|
||
['总价格 (元/吨)', market_cost['总价格'], result['总价格'], market_cost['总价格'] - result['总价格']],
|
||
['一斤饲料价格 (元)', market_cost['一斤饲料价格'], result['一斤饲料价格'],
|
||
market_cost['一斤饲料价格'] - result['一斤饲料价格']]
|
||
]
|
||
print(tabulate(cost_table, headers=['项', '市面配方', '优化配方', '偏差'], tablefmt='grid', floatfmt='.2f'))
|
||
pd.DataFrame(cost_table, columns=['项', '市面配方', '优化配方', '偏差']).to_csv('results/市面配方成本与偏差.csv',
|
||
index=False)
|
||
|
||
print("\n所有结果已存至 results/ 文件夹")
|