Water Level Management of Great Lakes
Greetings, readers. Today, I would like to delve into my recent research project on the water level management of the Great Lakes. This study focuses on developing a model to enhance the efficiency of dam utilization, thereby aligning with the varied interests of different stakeholders.
Problem 1: Stakeholder Identification and Nash Equilibrium
To commence, this study approached the issue from the perspective of game theory, aiming to identify the key stakeholders involved in water level management, such as safety, shipping, environmentalists, fishery, and tourism sectors. By seeking equilibrium points centered around the Nash equilibrium, we determined the relative weights of these stakeholders across different months. Employing a multi-objective programming model, we ultimately derived the optimal water levels for the five lakes over twelve months.
Problem 2: Network Flow Model and PID Parameter Adjustment
Next, the research abstracted the Great Lakes system into a network flow diagram and established a dynamic algorithm for adjusting PID (Proportional-Integral-Derivative) parameters in conjunction with optimization techniques. We modeled the lake water as an equivalent cone to facilitate volume calculations. By fine-tuning the gain coefficients for proportion, integration, and differentiation, the model aims to stabilize water levels and inflows while controlling outflows to minimize deviations from the optimal water levels. The findings indicate that this control algorithm enables rapid and stable attainment of the target water levels.
Problem 3: Sensitivity and Efficiency Analysis
The study further evaluated the sensitivity of the water level control algorithm and assessed its efficiency in aligning with stakeholders' interests. When applied to the 2017 water level and flow data, the algorithm demonstrated an improved capacity to maintain water levels closer to the optimal levels compared to actual 2017 data. Overall, the algorithm exhibits high flexibility and better satisfies the stakeholders' interests.
Problem 4: Sensitivity Analysis under External Factors
Additionally, the research conducted sensitivity analyses of the models under various conditions such as precipitation and ice blockage. The models were then tested under the comprehensive influence of climate and natural factors. Using the Monte Carlo algorithm, the Sobol index was calculated under these conditions. The results showed that the model effectively mitigates significant instabilities under different external factors, maintaining water levels within a reasonable and optimal range with high sensitivity.
Problem 5: Lake Ontario Specific Model
Lastly, based on the water level requirements of the stakeholders specified in the appendix, the study recalibrated the weights and utilized the initial model to determine the unique optimal water level for Lake Ontario. By analyzing recent precipitation, temperature, and runoff data from institutional databases, the study explored potential water level changes in Lake Ontario. It was concluded that the control algorithm proposed in this paper can effectively regulate the water level at Moses Saunders Dam.
In conclusion, this comprehensive study offers a robust framework for managing the water levels of the Great Lakes, addressing the multifaceted interests of various stakeholders. If you have any questions or insights, please feel free to share them in the comments section.
纳什均衡和多目标规划
首先从博弈论的角度出发，识别利益相关者，并尝试找到围绕纳什均衡的平衡点。其次，我们确定了五个不同的利益相关者（安全、航运、环保人士、渔民、旅游业），并建立了不同月份利益相关者的权重。最后，在应用多目标规划模型后，得到了五大湖在十二个月中的最佳水位。
import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm
# 利益相关者权重矩阵 (5个利益相关者 x 12个月)
weights = np.array([
[0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.8, 0.7, 0.6, 0.5], # 安全
[0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4], # 航运
[0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3], # 环保人士
[0.5, 0.6, 0.7, 0.8, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2], # 渔民
[0.6, 0.7, 0.8, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1] # 旅游业
])
# 动态权重调整函数
def dynamic_weights(month, historical_data):
# 根据历史数据调整权重
adjusted_weights = weights.copy()
for i in range(5):
adjusted_weights[i, month] *= (1 + norm.pdf(historical_data[month], loc=50, scale=10))
return adjusted_weights
# 目标函数：最大化所有利益相关者的利益
def objective_function(water_levels, weights):
total_benefit = 0
for i in range(5): # 5个利益相关者
for j in range(12): # 12个月
total_benefit += weights[i, j] * water_levels[j]
return -total_benefit # 最小化负的总利益
# 约束条件：水位范围在0到100之间
def constraint(water_levels):
return water_levels - 0 # 水位 >= 0
def constraint2(water_levels):
return 100 - water_levels # 水位 <= 100
# 约束条件：水位变化率
def constraint3(water_levels):
return np.diff(water_levels) - 10 # 水位变化率 <= 10
def constraint4(water_levels):
return 10 - np.diff(water_levels) # 水位变化率 >= -10
# 初始猜测值
initial_guess = np.full(12, 50) # 初始水位为50
# 历史数据（假设）
historical_data = np.random.normal(50, 10, 12)
# 动态权重调整
adjusted_weights = dynamic_weights(np.arange(12), historical_data)
# 约束条件
constraints = [
{'type': 'ineq', 'fun': constraint},
{'type': 'ineq', 'fun': constraint2},
{'type': 'ineq', 'fun': constraint3},
{'type': 'ineq', 'fun': constraint4}
]
# 优化
result = minimize(objective_function, initial_guess, args=(adjusted_weights,), method='SLSQP', constraints=constraints)
# 输出结果
print("最佳水位 (12个月):", result.x)
print("总利益:", -result.fun)网络流图和PID动态调整
将五大湖系统抽象为网络流图模型，并结合优化建立了PID参数的动态调整算法。其次，将湖水建模为一个等效圆锥体，以计算其体积。通过调整比例、积分和微分的增益系数，保持水位和流入量在稳定水平，控制流出量，以最小化水位与最佳水位之间的偏差。最终结果表明，在该控制算法下，湖水位可以快速且稳定地达到目标水位。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.stats import norm
# 湖泊参数
r = 1000 # 底面半径 (m)
h_target = 50 # 目标水位 (m)
h_initial = 40 # 初始水位 (m)
flow_in = 100 # 流入量 (m^3/s)
flow_out_initial = 50 # 初始流出量 (m^3/s)
# 初始PID控制参数
Kp = 1.0 # 比例增益
Ki = 0.1 # 积分增益
Kd = 0.01 # 微分增益
# 时间步长
dt = 1 # 时间步长 (s)
# 初始化变量
h = h_initial
flow_out = flow_out_initial
integral_error = 0
previous_error = 0
# 记录水位和流出量
h_history = []
flow_out_history = []
pid_history = []
# 模拟时间
T = 1000 # 总时间 (s)
# 自适应PID控制参数调整函数
def adaptive_pid(error_history):
# 使用最小二乘法调整PID参数
def objective_function(params):
Kp, Ki, Kd = params
total_error = 0
for i in range(len(error_history)):
error = error_history[i]
integral_error = np.sum(error_history[:i+1]) * dt
derivative_error = (error - previous_error) / dt
control_output = Kp * error + Ki * integral_error + Kd * derivative_error
total_error += abs(control_output)
return total_error
initial_guess = [Kp, Ki, Kd]
result = minimize(objective_function, initial_guess, method='Nelder-Mead')
return result.x
# 模拟过程
error_history = []
for t in range(T):
# 计算当前水位
volume = (1/3) * np.pi * r**2 * h
# 计算误差
error = h_target - h
error_history.append(error)
# 自适应调整PID参数
if t % 100 == 0:
Kp, Ki, Kd = adaptive_pid(error_history)
pid_history.append((Kp, Ki, Kd))
# 计算PID控制输出
integral_error += error * dt
derivative_error = (error - previous_error) / dt
flow_out = Kp * error + Ki * integral_error + Kd * derivative_error
# 更新水位
dh = (flow_in - flow_out) * dt
h += dh
# 记录数据
h_history.append(h)
flow_out_history.append(flow_out)
# 更新前一个误差
previous_error = error
# 绘制结果
plt.figure(figsize=(12, 12))
plt.subplot(3, 1, 1)
plt.plot(h_history, label='Water Level')
plt.axhline(y=h_target, color='r', linestyle='--', label='Target Water Level')
plt.xlabel('Time (s)')
plt.ylabel('Water Level (m)')
plt.legend()
plt.subplot(3, 1, 2)
plt.plot(flow_out_history, label='Flow Out')
plt.xlabel('Time (s)')
plt.ylabel('Flow Out (m^3/s)')
plt.legend()
plt.subplot(3, 1, 3)
plt.plot([p[0] for p in pid_history], label='Kp')
plt.plot([p[1] for p in pid_history], label='Ki')
plt.plot([p[2] for p in pid_history], label='Kd')
plt.xlabel('Time (s)')
plt.ylabel('PID Parameters')
plt.legend()
plt.tight_layout()
plt.show()蒙特卡罗参数调整
分析水位控制算法的敏感性，并评估利益的效率。在将该算法应用于2017年的水位和流量数据后，结果显示该算法能够使水位更接近最佳水位。总体而言，该算法具有高灵活性，能够更好地满足利益相关者的利益，优于2017年的实际数据。
# 敏感性分析
def sensitivity_analysis(test_data, Kp_range, Ki_range, Kd_range):
results = []
for Kp in Kp_range:
for Ki in Ki_range:
for Kd in Kd_range:
h = test_data['water_level'].iloc[0]
flow_out = test_data['flow_out'].iloc[0]
integral_error = 0
previous_error = 0
h_history = []
for t in range(len(test_data)):
error = h_target - h
integral_error += error * dt
derivative_error = (error - previous_error) / dt
flow_out = Kp * error + Ki * integral_error + Kd * derivative_error
dh = (test_data['flow_in'].iloc[t] - flow_out) * dt
h += dh
h_history.append(h)
previous_error = error
mse = np.mean((np.array(h_history) - h_target) ** 2)
mae = np.mean(np.abs(np.array(h_history) - h_target))
results.append((Kp, Ki, Kd, mse, mae))
return results
# 参数范围
Kp_range = np.linspace(0.5, 1.5, 5)
Ki_range = np.linspace(0.05, 0.15, 5)
Kd_range = np.linspace(0.005, 0.015, 5)
# 敏感性分析结果
sensitivity_results = sensitivity_analysis(test_data, Kp_range, Ki_range, Kd_range)
# 结果分析
best_params = min(sensitivity_results, key=lambda x: x[3]) # 选择MSE最小的参数组合
print("Best Parameters (Kp, Ki, Kd):", best_params[:3])
print("Best MSE:", best_params[3])
print("Best MAE:", best_params[4])Sobel指数计算
对降水、冰封等条件下的模型进行了敏感性分析。然后在气候和自然因素的综合影响下进行测试，并使用蒙特卡洛算法计算上述条件下的Sobol指数。结果表明，该模型能在不同外部因素下避免显著的不稳定因素，并在高敏感度下维持水位在合理和最佳范围内。
# Sobol指数分析
def model(X):
results = []
for x in X:
Kp, Ki, Kd, precipitation, ice_cover = x
h = test_data['water_level'].iloc[0]
flow_out = test_data['flow_out'].iloc[0]
integral_error = 0
previous_error = 0
h_history = []
for t in range(len(test_data)):
error = h_target - h
integral_error += error * dt
derivative_error = (error - previous_error) / dt
flow_out = Kp * error + Ki * integral_error + Kd * derivative_error
dh = (test_data['flow_in'].iloc[t] - flow_out) * dt
h += dh
h_history.append(h)
previous_error = error
mse = np.mean((np.array(h_history) - h_target) ** 2)
results.append(mse)
return np.array(results)安大略湖特定分析
根据附录中现有利益相关者的水位需求，重新设定权重，并使用第一个问题模型获得安大略湖的唯一最佳水位。随后，根据最近几周从机构数据库下载的降水、温度和径流数据，分析了安大略湖可能的水位变化。在摩西·桑德斯大坝使用本文设计的控制算法能够有效调节安大略湖的水位。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import minimize
# 模拟数据（假设）
stakeholder_data = pd.read_csv('stakeholder_water_level_requirements.csv') # 利益相关者水位需求数据
climate_data = pd.read_csv('recent_climate_data.csv') # 最近几周的降水、温度和径流数据
# 湖泊参数
r = 1000 # 底面半径 (m)
h_initial = 40 # 初始水位 (m)
flow_in_initial = 100 # 初始流入量 (m^3/s)
flow_out_initial = 50 # 初始流出量 (m^3/s)
# 初始PID控制参数
Kp = 1.0 # 比例增益
Ki = 0.1 # 积分增益
Kd = 0.01 # 微分增益
# 时间步长
dt = 1 # 时间步长 (s)
# 初始化变量
h = h_initial
flow_in = flow_in_initial
flow_out = flow_out_initial
integral_error = 0
previous_error = 0
# 记录水位和流出量
h_history = []
flow_out_history = []
pid_history = []
# 自适应PID控制参数调整函数
def adaptive_pid(error_history):
def objective_function(params):
Kp, Ki, Kd = params
total_error = 0
for i in range(len(error_history)):
error = error_history[i]
integral_error = np.sum(error_history[:i+1]) * dt
derivative_error = (error - previous_error) / dt
control_output = Kp * error + Ki * integral_error + Kd * derivative_error
total_error += abs(control_output)
return total_error
initial_guess = [Kp, Ki, Kd]
result = minimize(objective_function, initial_guess, method='Nelder-Mead')
return result.x
# 多目标优化函数
def multi_objective_optimization(stakeholder_data):
# 设定权重
weights = stakeholder_data['weight'].values
target_levels = stakeholder_data['target_level'].values
def objective_function(h):
total_error = 0
for i in range(len(stakeholder_data)):
error = target_levels[i] - h
total_error += weights[i] * abs(error)
return total_error
initial_guess = h_initial
result = minimize(objective_function, initial_guess, method='Nelder-Mead')
return result.x
# 计算唯一最佳水位
h_target = multi_objective_optimization(stakeholder_data)
print("Optimal Water Level (h_target):", h_target)
# 模拟过程
error_history = []
for t in range(len(climate_data)):
# 计算当前水位
volume = (1/3) * np.pi * r**2 * h
# 计算误差
error = h_target - h
error_history.append(error)
# 自适应调整PID参数
if t % 100 == 0:
Kp, Ki, Kd = adaptive_pid(error_history)
pid_history.append((Kp, Ki, Kd))
# 计算PID控制输出
integral_error += error * dt
derivative_error = (error - previous_error) / dt
flow_out = Kp * error + Ki * integral_error + Kd * derivative_error
# 更新水位
dh = (climate_data['flow_in'].iloc[t] - flow_out) * dt
h += dh
# 记录数据
h_history.append(h)
flow_out_history.append(flow_out)
# 更新前一个误差
previous_error = error
# 绘制结果
plt.figure(figsize=(12, 12))
plt.subplot(3, 1, 1)
plt.plot(h_history, label='Water Level')
plt.axhline(y=h_target, color='r', linestyle='--', label='Optimal Water Level')
plt.xlabel('Time (s)')
plt.ylabel('Water Level (m)')
plt.legend()
plt.subplot(3, 1, 2)
plt.plot(flow_out_history, label='Flow Out')
plt.xlabel('Time (s)')
plt.ylabel('Flow Out (m^3/s)')
plt.legend()
plt.subplot(3, 1, 3)
plt.plot([p[0] for p in pid_history], label='Kp')
plt.plot([p[1] for p in pid_history], label='Ki')
plt.plot([p[2] for p in pid_history], label='Kd')
plt.xlabel('Time (s)')
plt.ylabel('PID Parameters')
plt.legend()
plt.tight_layout()
plt.show()
# 分析最近几周的降水、温度和径流数据对水位的影响
def analyze_climate_impact(climate_data, h_history):
# 计算水位变化
h_change = np.diff(h_history)
# 计算降水、温度和径流对水位变化的影响
precipitation_impact = np.corrcoef(climate_data['precipitation'][1:], h_change)[0, 1]
temperature_impact = np.corrcoef(climate_data['temperature'][1:], h_change)[0, 1]
runoff_impact = np.corrcoef(climate_data['runoff'][1:], h_change)[0, 1]
return precipitation_impact, temperature_impact, runoff_impact
# 结果分析
precipitation_impact, temperature_impact, runoff_impact = analyze_climate_impact(climate_data, h_history)
print("Precipitation Impact:", precipitation_impact)
print("Temperature Impact:", temperature_impact)
print("Runoff Impact:", runoff_impact)
