集总式新安江代码(二、三水源)python

首先设置模型参数

import numpy as np
import pandas as pd
import math
df = pd.read_excel(r"C:\Users\Admin\Desktop\pet.xlsx")
df = df.fillna(0)
print(df)
WM,WUM,WLM,WDM,b,K,C,FC,F,T,CG,L,CR = 120,15,85,20,0.3,0.95,0.14,15,2856,24,0.9,5,0.5    # F单位为km^2 T单位为h
Sm,KI,KG,EX,FR1,CI,Q = 20,0.2,0.5,1.5,0.4175,0.9,5.6
WMM = (1+b)*WM
U = F / (3.6 * T)
SMM = Sm * (1+EX)

我的数据是excel数据其中有流域面平均雨量和蒸发资料

首先是设置一个for循坏我一共5000多行降雨数据,设置一个for循坏

for i in range(0,5479):
# 计算土湿
    if i == 0:
        df.loc[i, "WU"] = 0
        df.loc[i,"WL"] = 2.2
        df.loc[i, "WD"] = 20
        df.loc[i, "EP"] = K * df.loc[i, "E0"]
        if df.loc[i, "WU"] + df.loc[i, "P"] >= df.loc[i, "EP"]:
            df.loc[i, "EU"] = df.loc[i, "EP"]
            df.loc[i, "EL"] = 0
            df.loc[i, "ED"] = 0
        if df.loc[i, "WU"] + df.loc[i, "P"] < df.loc[i, "EP"] and df.loc[i, "WL"] >= C * WLM:
            df.loc[i, "EU"] = df.loc[i, "WU"] + df.loc[i, "P"]
            df.loc[i, "EL"] = (df.loc[i, "EP"] - df.loc[i, "EU"]) * df.loc[i, "WL"] / WLM
            df.loc[i, "ED"] = 0
        if df.loc[i, "WU"] + df.loc[i, "P"] < df.loc[i, "EP"] and df.loc[i, "WL"] < C * WLM and df.loc[i, "WL"] >= C * (df.loc[i, "EP"] - df.loc[i, "EU"]):
            df.loc[i, "EU"] = df.loc[i, "WU"] + df.loc[i, "P"]
            df.loc[i, "EL"] = C * (df.loc[i, "EP"] - df.loc[i, "EU"])
            df.loc[i, "ED"] = 0
        if df.loc[i, "WU"] + df.loc[i, "P"] < df.loc[i, "EP"] and df.loc[i, "WL"] < C * (
                df.loc[i, "EP"] - df.loc[i, "EU"]):
            df.loc[i, "EU"] = df.loc[i, "WU"] + df.loc[i, "P"]
            df.loc[i, "EL"] = df.loc[i, "WL"]
            df.loc[i, "ED"] = C * (df.loc[i, "EP"] - df.loc[i, "EU"]) - df.loc[i, "EL"]
        df.loc[i, "E"] = df.loc[i, "EU"] + df.loc[i, "EL"] + df.loc[i, "ED"]
        df.loc[i, "PE"] = df.loc[i, "P"] - df.loc[i, "E"]
        df.loc[i, "W"] = df.loc[i, "WU"] + df.loc[i, "WL"] + df.loc[i, "WD"]
    else:
        dw = df.loc[i-1,"PE"] - df.loc[i-1,"R"] #这边需要连接下面R,否则会报错
        df.loc[i,"WU"] = df.loc[i-1,"WU"] + dw
        df.loc[i,"WL"] = df.loc[i-1,"WL"]
        df.loc[i,"WD"] = df.loc[i-1,"WD"]
        if df.loc[i,"WU"] < 0 :
            df.loc[i,"WL"] = df.loc[i-1,"WL"] + dw + df.loc[i-1,"WU"]
            df.loc[i,"WD"] = df.loc[i-1,"WD"]
            df.loc[i,"WU"] = 0
            if df.loc[i,"WL"] < 0:
                df.loc[i,"WD"] = df.loc[i-1,"WD"] + dw + df.loc[i-1,"WU"] + df.loc[i-1,"WL"]
                df.loc[i,"WU"] = 0
                df.loc[i,"WL"] = 0
                if df.loc[i,"WD"] < 0:
                    df.loc[i,"WU"] = 0
                    df.loc[i,"WD"] = 0
                    df.loc[i,"WL"] = 0
        if df.loc[i,"WU"] > WUM:
            df.loc[i,"WU"] = WUM
            df.loc[i,"WL"] = df.loc[i-1,"WL"] + dw - WUM + df.loc[i-1,"WU"]
            df.loc[i,"WD"] = df.loc[i-1,"WD"]
            if df.loc[i,"WL"] > WLM:
                df.loc[i,"WU"] = WUM
                df.loc[i,"WL"] = WLM
                df.loc[i,"WD"] = df.loc[i-1,"WD"] + dw - WLM - WUM + df.loc[i-1,"WU"] + df.loc[i-1,"WL"]
        df.loc[i,"EP"] = K * df.loc[i,"E0"]
        if df.loc[i,"WU"] + df.loc[i,"P"] >= df.loc[i,"EP"]:
            df.loc[i,"EU"] = df.loc[i,"EP"]
            df.loc[i,"EL"] = 0
            df.loc[i,"ED"] = 0
        if df.loc[i,"WU"] + df.loc[i,"P"] < df.loc[i,"EP"] and df.loc[i,"WL"] >= C * WLM:
            df.loc[i,"EU"] = df.loc[i,"WU"] + df.loc[i,"P"]
            df.loc[i,"EL"] = (df.loc[i,"EP"] - df.loc[i,"EU"]) * df.loc[i,"WL"] / WLM
            df.loc[i,"ED"] = 0
        if df.loc[i,"WU"] + df.loc[i,"P"] < df.loc[i,"EP"] and df.loc[i,"WL"] < C * WLM and df.loc[i,"WL"] >= C * (df.loc[i,"EP"]-df.loc[i,"EU"]):
            df.loc[i,"EU"] = df.loc[i,"WU"] + df.loc[i,"P"]
            df.loc[i,"EL"] = C * (df.loc[i,"EP"] - df.loc[i,"EU"])
            df.loc[i,"ED"] = 0
        if df.loc[i,"WU"] + df.loc[i,"P"] < df.loc[i,"EP"] and df.loc[i,"WL"] < C * (df.loc[i,"EP"]-df.loc[i,"EU"]):
            df.loc[i,"EU"] = df.loc[i,"WU"] + df.loc[i,"P"]
            df.loc[i,"EL"] = df.loc[i,"WL"]
            df.loc[i,"ED"] = C * (df.loc[i,"EP"] - df.loc[i,"EU"]) - df.loc[i,"EL"]
        df.loc[i,"E"] = df.loc[i,"EU"] + df.loc[i,"EL"] + df.loc[i,"ED"]
        df.loc[i,"PE"] = df.loc[i,"P"] - df.loc[i,"E"]
        df.loc[i,"W"] = df.loc[i,"WU"] + df.loc[i,"WL"] + df.loc[i,"WD"]

再根据自己的要求计算R

# 计算R
    a = WMM * (1- math.pow((1-df.loc[i,"W"]/WM),1/(1+b)))
    if i ==0:
        if df.loc[i,"PE"] > 0:
            if a + df.loc[i,"PE"] <= WMM:
                df.loc[i,"R"] = df.loc[i,"PE"] + df.loc[i,"W"] - WM + WM * math.pow(1-((df.loc[i,"PE"]+a)/WMM),b+1)  # df["R"][i] = df["PE"][i] + df["W"][i] - WM + WM * math.pow(1 - ((df["PE"][i] + a) / WMM), b + 1)
            if a + df.loc[i,"PE"] > WMM:
                df.loc[i,"R"] = df.loc[i,"PE"] - (WM - df.loc[i,"W"])
        if df.loc[i,"R"] < 0:
            df.loc[i,"R"] = 0
    else:
        if df.loc[i,"PE"] > 0:
            if a + df.loc[i, "PE"] <= WMM:
                 df.loc[i, "R"] = df.loc[i, "PE"] + df.loc[i, "W"] - WM + WM * math.pow(1 - ((df.loc[i, "PE"] + a) / WMM),b + 1)  # df["R"][i] = df["PE"][i] + df["W"][i] - WM + WM * math.pow(1 - ((df["PE"][i] + a) / WMM), b + 1)
            if a + df.loc[i, "PE"] > WMM:
                df.loc[i, "R"] = df.loc[i, "PE"] - (WM - df.loc[i, "W"])
        if df.loc[i, "R"] < 0:
            df.loc[i, "R"] = 0

在计算产流面积

# 计算流域产流面积
    if i == 0:
        if df.loc[i,"R"] > 0 :
            df.loc[i,"FR"] = df.loc[i,"R"] / df.loc[i,"PE"]
            if df.loc[i,"FR"] >1:
                df.loc[i,"FR"] =1

        else:
            df.loc[i,"FR"] = FR1
    else:
        if df.loc[i, "R"] > 0:
            df.loc[i, "FR"] = df.loc[i, "R"] / df.loc[i, "PE"]
            if df.loc[i, "FR"] > 1:
                df.loc[i, "FR"] = 1
        else:
            df.loc[i, "FR"] = df.loc[i-1,"FR"]

这边是二水源代码

# 划分二水源
#     if i == 0:
#         if df.loc[i, "R"] > 0:
#             df.loc[i, "FR"] = df.loc[i, "R"] / df.loc[i, "PE"]
#             if df.loc[i,"FR"] > 1:
#                 df.loc[i,"FR"] = 1
#             if df.loc[i, "FR"] > 0:
#                 if df.loc[i, "PE"] > FC:
#                     df.loc[i, "RS"] = (df.loc[i, "PE"] - FC) * df.loc[i, "FR"]
#                     df.loc[i, "RG"] = df.loc[i, "R"] - df.loc[i, "RS"]
#                 if df.loc[i, "PE"] <= FC:
#                     df.loc[i, "RS"] = 0
#                     df.loc[i, "RG"] = df.loc[i, "R"]
#             if df.loc[i, "FR"] < 0:
#                 df.loc[i, "FR"] = 0
#                 df.loc[i, "RS"] = 0
#                 df.loc[i, "RG"] = 0
#         if df.loc[i, "R"] == 0:
#             df.loc[i, "RS"] = 0
#             df.loc[i, "RG"] = 0
#     else:
#         if df.loc[i,"R"] > 0 :
#             df.loc[i,"FR"] = df.loc[i,"R"] / df.loc[i,"PE"]
#             if df.loc[i,"FR"] > 0:
#                 if df.loc[i,"PE"] > FC:
#                     df.loc[i,"RS"] = (df.loc[i,"PE"] - FC) * df.loc[i,"FR"]
#                     df.loc[i,"RG"] = df.loc[i,"R"] - df.loc[i,"RS"]
#                 if df.loc[i,"PE"] <= FC:
#                     df.loc[i,"RS"] = 0
#                     df.loc[i,"RG"] = df.loc[i,"R"]
#             if df.loc[i,"FR"] < 0:
#                 df.loc[i,"FR"] = 0
#                 df.loc[i,"RS"] = 0
#                 df.loc[i,"RG"] = 0
#         if df.loc[i,"R"] == 0:
#             df.loc[i,"RS"] = 0
#             df.loc[i,"RG"] = 0

三水源

# 划分三水源
    if i == 0:     # 第一时段计算
        df.loc[i, "S1"] = 3.4525
        if df.loc[i,"PE"] > 0:
            AU = SMM * (1 - math.pow((1-(((df.loc[i,"S1"]*FR1)/df.loc[i,"FR"])/Sm)),1/(1+EX)))
            if df.loc[i,"PE"] + AU < SMM:
                df.loc[i,"RS"] = df.loc[i,"FR"] * (df.loc[i,"PE"] +(df.loc[i,"S1"] * FR1)/df.loc[i,"FR"] - Sm +
                                                 Sm *math.pow((1-(df.loc[i,"PE"] +AU)/SMM),1+EX))
            if df.loc[i,"PE"] + AU >= SMM:
                df.loc[i,"RS"] = df.loc[i,"FR"] * (df.loc[i,"PE"] + (df.loc[i,"S1"] * FR1)/df.loc[i,"FR"] - Sm)
            S = (df.loc[i,"S1"] * FR1)/df.loc[i,"FR"] + (df.loc[i,"R"] - df.loc[i,"RS"])/df.loc[i,"FR"]
            df.loc[i,"RI"] = KI * S *df.loc[i,"FR"]
            df.loc[i,"RG"] = KG * S *df.loc[i,"FR"]
            df.loc[i+1,"S1"] = S*(1-KI-KG)
        else:
            S = (df.loc[i, "S1"] * FR1) / df.loc[i, "FR"]
            df.loc[i + 1, "S1"] = S * (1 - KG - KI)
            df.loc[i, "RS"] = 0
            df.loc[i, "RG"] = KI*S * df.loc[i,"FR"]
            df.loc[i, "RI"] = KG * S *df.loc[i,"FR"]
    else: #其余时段计算
        if df.loc[i, "PE"] > 0:
            AU = SMM * (1 - math.pow((1-(((df.loc[i,"S1"]*df.loc[i-1,"FR"])/df.loc[i,"FR"])/Sm)),1/(1+EX)))
            if df.loc[i,"PE"] + AU < SMM:
                df.loc[i,"RS"] = df.loc[i,"FR"] * (df.loc[i,"PE"] +(df.loc[i,"S1"] * df.loc[i-1,"FR"])/df.loc[i,"FR"] - Sm +
                                                 Sm *math.pow((1-(df.loc[i,"PE"] +AU)/SMM),1+EX))
            if df.loc[i,"PE"] + AU >= SMM:
                df.loc[i,"RS"] = df.loc[i,"FR"] * (df.loc[i,"PE"] + (df.loc[i,"S1"] * df.loc[i-1,"FR"])/df.loc[i,"FR"] - Sm)
            S = (df.loc[i,"S1"] * df.loc[i-1,"FR"])/df.loc[i,"FR"] + (df.loc[i,"R"] - df.loc[i,"RS"])/df.loc[i,"FR"]
            df.loc[i,"RI"] = KI * S *df.loc[i,"FR"]
            df.loc[i,"RG"] = KG * S *df.loc[i,"FR"]
            df.loc[i+1,"S1"] = S*(1-KI-KG)
        else:
            S = (df.loc[i, "S1"] * df.loc[i-1,"FR"]) / df.loc[i, "FR"]
            df.loc[i + 1, "S1"] =S * (1 - KG - KI)
            df.loc[i, "RS"] = 0
            df.loc[i, "RG"] = KG *S *df.loc[i,"FR"]
            df.loc[i, "RI"] = KI * S *df.loc[i,"FR"]

二水源汇流

# # 计算二水源汇流
#     if i == 0:
#         df.loc[i, "Qs"] = df.loc[i, "RS"] * U
#         df.loc[i,"Qg"] = 5.6
#         df.loc[i, "QT"] = df.loc[i, "Qs"] + df.loc[i, "Qg"]
#         df.loc[i,"Qt"] = 5.6
#     else:
#         df.loc[i,"Qs"] = df.loc[i,"RS"] * U
#         df.loc[i,"Qg"] = CG * df.loc[i-1, "Qg"] + (1-CG) * df.loc[i,"RG"] * U
#         df.loc[i,"QT"] = df.loc[i,"Qs"] + df.loc[i,"Qg"]
#         df.loc[i,"Qt"] = CR * df.loc[i-1,"Qt"] + (1 - CR) * df.loc[i-L,"QT"]

三水源汇流

# 计算三水源汇流
    if i == 0:
        df.loc[i,"QS"] = df.loc[i,"RS"] * U
        df.loc[i,"QI"] = 1/3 * Q
        df.loc[i,"QG"] = 1/3 * Q
        df.loc[i, "QT"] = df.loc[i, "QS"] + df.loc[i, "QI"] + df.loc[i, "QG"]
    if i >0:
        df.loc[i, "QS"] = df.loc[i, "RS"] * U
        df.loc[i,"QI"] = CI * df.loc[i-1,"QI"] + (1-CI)*df.loc[i,"RI"]*U
        df.loc[i,"QG"] = CG* df.loc[i-1,"QG"]+(1-CG)*df.loc[i,"RG"] *U
        df.loc[i,"QT"] = df.loc[i,"QS"]+df.loc[i,"QI"]+df.loc[i,"QG"]
    if i >=0 and i <= L:
        df.loc[i,"Qt"] = Q
    if i >= L+1:
        df.loc[i,"Qt"] = df.loc[i-1,"Qt"]*CR+(1-CR) * df.loc[i-L,"QT"]

最后写入excel表,保留4位小数

ds = pd.DataFrame(df)
ds.to_excel(r"C:\Users\Admin\Desktop\新安江1.xlsx",float_format='%.4f')

最后三水源结果

写到这感觉自己写的太麻烦了,希望以后能够改进吧

这边蒸散发还是借鉴了一个哥们的代码,链接在这边:(26条消息) Python数据分析实例,利用Pandas建立流域三层蒸发和蓄满产流模型_三层蒸发模型代码_Joyonlyonly的博客-CSDN博客

  • 13
    点赞
  • 70
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 12
    评论
新安江模型是一种常用的水资源管理模型,用于预测和分析水文和水资源系统。在Python中,可以使用以下代码实现一个简单的新安江模型程序: ```python # 导入所需的库 import numpy as np # 定义新安江模型的参数 initial_storage = 100 # 初始水库蓄水量(单位:百万立方米) inflow = np.array([50, 60, 70, 80, 90, 100, 110, 120]) # 每个时期的入库流量(单位:百万立方米) demand = np.array([20, 25, 30, 35, 40, 45, 50, 55]) # 每个时期的需水量(单位:百万立方米) # 创建一个空数组来存储每个时期的蓄水量和开闸流量 storage = np.zeros(len(inflow)) # 蓄水量 release = np.zeros(len(inflow)) # 开闸流量 # 计算每个时期的蓄水量和开闸流量 for i in range(len(inflow)): if i == 0: storage[i] = initial_storage + inflow[i] - demand[i] else: storage[i] = storage[i-1] + inflow[i] - demand[i] if storage[i] < 0: release[i] = demand[i] - storage[i] storage[i] = 0 # 打印每个时期的蓄水量和开闸流量 for i in range(len(inflow)): print("第", i+1, "时期的蓄水量为:", storage[i], "百万立方米") print("第", i+1, "时期的开闸流量为:", release[i], "百万立方米") ``` 运行该代码,可以得到每个时期的蓄水量和开闸流量。这个简单的新安江模型程序可以作为进一步研究和优化水资源管理的基础。注意,该模型是根据给定的入库流量和需水量进行模拟计算,并没有考虑其他影响因素,如降雨量、蒸发量等。在实际应用中,通常需要添加更多因素和参数来提高模型的准确性和可靠性。
评论 12
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Peanut-uu

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值