import numpy as np
import matplotlib.pyplot as plt
def interpolate_storage_to_level(storage, curve):
"""根据库容插值计算水位"""
levels = np.array(list(curve.keys()))
storages = np.array(list(curve.values()))
if storage <= min(storages):
return min(levels)
elif storage >= max(storages):
return max(levels)
else:
# 线性插值计算水位
idx = np.searchsorted(storages, storage)
x0, x1 = storages[idx - 1], storages[idx]
y0, y1 = levels[idx - 1], levels[idx]
return y0 + (storage - x0) * (y1 - y0) / (x1 - x0)
def optimize_reservoir_operation(
num_stages, dt, inflow_rates, safe_discharge_limit,
water_level_curve, initial_storage, final_storage
):
"""水库优化调度主函数"""
# 离散化库容状态
storage_levels = np.array(sorted(water_level_curve.values()))
storage_levels_count = len(storage_levels)
# 检查初始和最终库容是否在范围内
if initial_storage not in storage_levels:
raise ValueError(f"初始库容 {initial_storage} 不在离散化状态列表中")
if final_storage not in storage_levels:
raise ValueError(f"最终库容 {final_storage} 不在离散化状态列表中")
# 初始化动态规划数组
dp_cost = np.full((num_stages + 1, storage_levels_count), np.inf)
dp_policy = np.full((num_stages + 1, storage_levels_count), -1)
# 设置边界条件
final_storage_idx = np.argmin(np.abs(storage_levels - final_storage))
dp_cost[num_stages, final_storage_idx] = 0
# 逆序动态规划求解
for t in range(num_stages - 1, -1, -1):
for i in range(storage_levels_count):
current_storage = storage_levels[i]
# 计算可能的下泄流量范围
max_discharge = min(safe_discharge_limit, inflow_rates[t] * 2)
discharge_range = np.arange(0, max_discharge + 1, 10)
for discharge in discharge_range:
# 计算下一状态的库容
next_storage = current_storage + (inflow_rates[t] - discharge) * dt
# 查找最接近的离散化库容状态
next_idx = np.argmin(np.abs(storage_levels - next_storage))
next_actual_storage = storage_levels[next_idx]
# 计算误差
storage_error = np.abs(next_actual_storage - next_storage)
# 如果误差在可接受范围内,更新DP表
if storage_error < 1e6: # 1e6 m³的容差
cost = discharge ** 2 + dp_cost[t + 1, next_idx]
if cost < dp_cost[t, i]:
dp_cost[t, i] = cost
dp_policy[t, i] = discharge
# 回溯最优策略
optimal_discharges = []
optimal_storages = [initial_storage]
current_storage_idx = np.argmin(np.abs(storage_levels - initial_storage))
# 修正:num_stages应该等于时段数,而不是时刻数
for t in range(num_stages):
discharge = dp_policy[t, current_storage_idx]
optimal_discharges.append(discharge)
# 计算下一时刻的库容
next_storage = optimal_storages[-1] + (inflow_rates[t] - discharge) * dt
optimal_storages.append(next_storage)
next_storage_idx = np.argmin(np.abs(storage_levels - next_storage))
current_storage_idx = next_storage_idx
# 计算对应的水位
water_levels = [interpolate_storage_to_level(s, water_level_curve) for s in optimal_storages]
return {
'optimal_discharges': optimal_discharges,
'optimal_storages': optimal_storages[1:],
'water_levels': water_levels[1:],
'total_cost': dp_cost[0, np.argmin(np.abs(storage_levels - initial_storage))]
}
def plot_results(inflow_rates, optimal_discharges, optimal_storages, water_levels, dt):
"""可视化优化结果"""
# 修正:时段数应该等于len(inflow_rates)
num_stages = len(inflow_rates)
time_hours = np.arange(num_stages + 1) * dt / 3600 # 转换为小时
plt.figure(figsize=(15, 10))
# 绘制入库流量和下泄流量
plt.subplot(3, 1, 1)
plt.bar(time_hours[:-1], inflow_rates, width=dt / 3600 / 1.5, label='入库流量', color='blue')
plt.bar(time_hours[:-1], optimal_discharges, width=dt / 3600 / 3, label='下泄流量', color='orange')
plt.axhline(y=2010, color='r', linestyle='--', label='安全泄量限制')
plt.xlabel('时间 (小时)')
plt.ylabel('流量 (m³/s)')
plt.title('水库调度流量过程')
plt.legend()
plt.grid(True)
# 绘制库容变化
plt.subplot(3, 1, 2)
plt.plot(time_hours, optimal_storages, 'o-', label='库容变化')
plt.xlabel('时间 (小时)')
plt.ylabel('库容 (m³)')
plt.title('水库库容变化过程')
plt.legend()
plt.grid(True)
# 绘制水位变化
plt.subplot(3, 1, 3)
plt.plot(time_hours, water_levels, 'o-', label='水位变化', color='green')
plt.xlabel('时间 (小时)')
plt.ylabel('水位 (m)')
plt.title('水库水位变化过程')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
def main():
"""主函数,整合所有功能"""
# 水库参数设置
# 修正:时段数应为4,5个时刻对应4个时段
num_stages = 4
dt = 10800 # 时段长度(秒)
# 修正:将时段长度转换为小时,以便正确计算m³/s·h
dt_hours = dt / 3600
# 5个时刻对应的4个时段平均入库流量(m³/s)
inflow_rates = np.array([5403, 8125, 6575, 4188])
safe_discharge_limit = 2010 # 下游安全泄量限制(m³/s)
# 水位-库容曲线 (水位(m) -> 库容(m³))
water_level_curve = {
110.00: 60500,
110.80: 68900,
110.90: 69950,
111.00: 71000,
111.50: 76500,
111.80: 79800,
112.00: 82000,
112.60: 88600,
112.70: 89700,
113.00: 93000,
113.70: 102100,
113.80: 103400,
114.00: 106000,
114.25: 109370
}
# 初始和最终条件
initial_storage = 60500 # 初始库容(m³)
final_storage = 109370 # 最终库容(m³)
try:
# 执行优化计算
results = optimize_reservoir_operation(
num_stages, dt, inflow_rates, safe_discharge_limit,
water_level_curve, initial_storage, final_storage
)
# 输出结果
print("最优下泄流量计划 (m³/s):", results['optimal_discharges'])
print("各时段末库容 (m³):", results['optimal_storages'])
print("各时段末水位 (m):", results['water_levels'])
print("总目标函数值 (下泄流量平方和):", results['total_cost'])
# 可视化结果
plot_results(
inflow_rates,
results['optimal_discharges'],
[initial_storage] + results['optimal_storages'],
[interpolate_storage_to_level(initial_storage, water_level_curve)] + results['water_levels'],
dt
)
except Exception as e:
print(f"计算过程中发生错误: {e}")
if __name__ == "__main__":
main()