import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize
import seaborn as sns
# 设置全局样式(使用 seaborn 主题)
sns.set_theme(style="whitegrid")
# ======================
# 1. 数据模拟与预处理
# ======================
def simulate_leakage_data(days=365):
"""模拟校园供水系统漏水量数据"""
np.random.seed(2025) # 固定随机种子保证可复现
# 基础漏水量(正态分布)
base_leakage = np.random.normal(loc=100, scale=20, size=days)
# 添加季节性波动(夏季漏水量更高)
seasonal_effect = 1 + 0.3 * np.sin(2 * np.pi * np.arange(days) / 365 + np.pi / 2)
leakage = base_leakage * seasonal_effect
# 添加随机事件(暴雨、管道老化等)
event_days = np.random.choice(range(days), size=int(days * 0.1), replace=False)
leakage[event_days] *= np.random.uniform(1.5, 3.0, size=len(event_days))
# 创建日期索引
dates = pd.date_range(start='2024-01-01', periods=days)
return pd.Series(leakage, index=dates, name='漏水量')
# 生成模拟数据
leakage_data = simulate_leakage_data()
print("模拟漏水量数据统计摘要:")
print(leakage_data.describe())
# ======================
# 2. 维修决策优化模型
# ======================
def repair_optimization_model(leakage, water_price, repair_cost):
"""
维修决策优化模型(0-1规划)
:param leakage: 每日漏水量序列
:param water_price: 水价(元/吨)
:param repair_cost: 单次维修成本(元)
:return: 最优维修决策, 维修次数, 总损失金额
"""
n_days = len(leakage)
# 目标函数:最小化总损失金额
def objective(x):
# x: 决策变量 (0/1 表示每天是否维修)
repair_expense = repair_cost * x
leakage_loss = water_price * leakage * (1 - x)
return np.sum(repair_expense + leakage_loss)
# 约束条件:无特殊约束
constraints = []
# 变量边界 (0 <= x_i <= 1)
bounds = [(0, 1) for _ in range(n_days)]
# 初始解(全部不维修)
x0 = np.zeros(n_days)
# 求解优化问题
res = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=constraints)
# 提取最优解(四舍五入为0或1)
optimal_decision = np.round(res.x).astype(int)
# 计算维修次数和总损失金额
repair_count = np.sum(optimal_decision)
total_loss = objective(optimal_decision)
return optimal_decision, repair_count, total_loss
# ======================
# 3. 三维曲面分析
# ======================
def plot_3d_surface_analysis(leakage):
"""绘制维修次数和损失金额关于水价和维修成本的三维曲面图"""
# 设置参数范围
water_prices = np.linspace(5.5, 6.5, 20) # 水价范围:5.5-6.5元/吨
repair_costs = np.linspace(500, 1500, 20) # 维修成本范围:500-1500元/次
# 创建网格存储结果
repair_counts = np.zeros((len(water_prices), len(repair_costs)))
total_losses = np.zeros((len(water_prices), len(repair_costs)))
# 计算所有组合的结果
for i, wp in enumerate(water_prices):
for j, rc in enumerate(repair_costs):
_, repair_count, total_loss = repair_optimization_model(leakage.values, wp, rc)
repair_counts[i, j] = repair_count
total_losses[i, j] = total_loss
# 创建三维网格
WP, RC = np.meshgrid(water_prices, repair_costs)
# 创建三维画布
fig = plt.figure(figsize=(18, 12))
# 维修次数曲面图
ax1 = fig.add_subplot(211, projection='3d')
surf1 = ax1.plot_surface(WP, RC, repair_counts.T, cmap='viridis',
edgecolor='none', alpha=0.8, rstride=1, cstride=1)
ax1.set_xlabel('水价 (元/吨)', fontsize=12, labelpad=10)
ax1.set_ylabel('维修成本 (元/次)', fontsize=12, labelpad=10)
ax1.set_zlabel('维修次数', fontsize=12, labelpad=10)
ax1.set_title('维修次数与水价和维修成本的关系', fontsize=16, pad=15)
ax1.view_init(elev=30, azim=45)
fig.colorbar(surf1, ax=ax1, shrink=0.6, aspect=10, pad=0.1, label='维修次数')
# 总损失金额曲面图
ax2 = fig.add_subplot(212, projection='3d')
surf2 = ax2.plot_surface(WP, RC, total_losses.T, cmap='plasma',
edgecolor='none', alpha=0.8, rstride=1, cstride=1)
ax2.set_xlabel('水价 (元/吨)', fontsize=12, labelpad=10)
ax2.set_ylabel('维修成本 (元/次)', fontsize=12, labelpad=10)
ax2.set_zlabel('总损失金额 (元)', fontsize=12, labelpad=10)
ax2.set_title('总损失金额与水价和维修成本的关系', fontsize=16, pad=15)
ax2.view_init(elev=30, azim=45)
fig.colorbar(surf2, ax=ax2, shrink=0.6, aspect=10, pad=0.1, label='总损失金额 (元)')
plt.tight_layout()
plt.savefig('repair_3d_analysis.png', dpi=300, bbox_inches='tight')
plt.show()
return WP, RC, repair_counts, total_losses
# ======================
# 4. 维修计划可视化
# ======================
def visualize_repair_schedule(leakage, optimal_decision, water_price, repair_cost):
"""可视化维修安排计划"""
# 创建时间序列
dates = leakage.index
n_days = len(dates)
# 计算每日损失成本
daily_loss = np.where(optimal_decision == 1, repair_cost, water_price * leakage)
# 创建子图
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 10), sharex=True)
# 漏水量和维修决策
ax1.bar(dates, leakage, color='skyblue', alpha=0.7, label='日漏水量')
ax1.plot(dates, leakage.rolling(7).mean(), 'r-', linewidth=2, label='7日移动平均')
repair_dates = dates[optimal_decision == 1]
ax1.scatter(repair_dates, leakage[optimal_decision == 1],
color='red', s=80, zorder=5, label='维修日')
ax1.set_ylabel('漏水量 (吨)', fontsize=12)
ax1.set_title(f'漏水量监测与维修决策 (水价={water_price}元/吨, 维修成本={repair_cost}元/次)', fontsize=14)
ax1.legend(loc='upper left')
ax1.grid(True, linestyle='--', alpha=0.7)
# 维修状态
ax2.fill_between(dates, 0, 1, where=optimal_decision == 1,
color='red', alpha=0.3, label='维修日')
ax2.fill_between(dates, 0, 1, where=optimal_decision == 0,
color='blue', alpha=0.3, label='非维修日')
ax2.set_xlabel('日期', fontsize=12)
ax2.set_ylabel('维修状态', fontsize=12)
ax2.set_yticks([])
ax2.legend(loc='upper left')
ax2.grid(True, linestyle='--', alpha=0.7)
# 添加月度分隔线
months = pd.date_range(start=dates[0], end=dates[-1], freq='MS')
for month in months:
ax1.axvline(month, color='gray', linestyle='--', alpha=0.5)
ax2.axvline(month, color='gray', linestyle='--', alpha=0.5)
plt.tight_layout()
plt.savefig('repair_schedule.png', dpi=300, bbox_inches='tight')
plt.show()
# 创建损失金额数据框
loss_df = pd.DataFrame({
'日期': dates,
'漏水量': leakage.values,
'维修决策': optimal_decision,
'日损失金额': daily_loss
})
# 按月汇总
monthly_summary = loss_df.resample('M', on='日期').agg({
'漏水量': 'sum',
'维修决策': 'sum',
'日损失金额': 'sum'
}).rename(columns={
'漏水量': '月漏水量',
'维修决策': '月维修次数',
'日损失金额': '月损失金额'
})
print("\n月度维修与损失摘要:")
print(monthly_summary)
return loss_df, monthly_summary
# ======================
# 5. 执行主程序
# ======================
def main():
"""主执行函数"""
print("=" * 70)
print("校园供水系统漏损维修决策优化系统")
print("=" * 70)
# 1. 模拟漏水量数据
print("\n步骤1: 模拟漏水量数据...")
leakage_data = simulate_leakage_data()
# 2. 三维曲面分析
print("\n步骤2: 进行水价与维修成本的三维曲面分析...")
WP, RC, repair_counts, total_losses = plot_3d_surface_analysis(leakage_data)
# 3. 特定参数分析(文档中的案例)
water_price_case = 6.0
repair_cost_case = 1000
print(f"\n步骤3: 分析特定案例 (水价={water_price_case}元/吨, 维修成本={repair_cost_case}元/次)...")
optimal_decision, repair_count, total_loss = repair_optimization_model(
leakage_data.values, water_price_case, repair_cost_case
)
print(f"维修次数: {repair_count}次")
print(f"总损失金额: {total_loss:.2f}元")
# 4. 可视化维修安排
print("\n步骤4: 可视化维修安排...")
loss_df, monthly_summary = visualize_repair_schedule(
leakage_data, optimal_decision, water_price_case, repair_cost_case
)
# 5. 结果总结
print("\n" + "=" * 70)
print("优化结果总结")
print("=" * 70)
print(f"1. 全年总漏水量: {leakage_data.sum():.2f} 吨")
print(f"2. 最优维修次数: {repair_count} 次")
print(f"3. 总损失金额: {total_loss:.2f} 元")
print(f"4. 平均每次维修节省损失: {leakage_data.sum() * water_price_case / repair_count:.2f} 元/次")
# 保存结果到CSV
loss_df.to_csv('daily_repair_loss.csv', index=False)
monthly_summary.to_csv('monthly_repair_summary.csv')
print("\n结果已保存至 daily_repair_loss.csv 和 monthly_repair_summary.csv")
if __name__ == "__main__":
main()
输出可视化