利用pyswmm+pymoo完成基于swmm模型与NSGA-Ⅱ算法的城市系统排水多目标优化（2024.3.9更新）

已于 2024-03-09 11:20:54 修改
阅读量3k 收藏 41
点赞数 24
文章标签： 算法 机器学习 人工智能 python
于 2023-04-28 22:09:51 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/mjzjz_c/article/details/130435178
版权

目录

前言

比较闲暇的时间，为了保持点编程的工作量，帮忙解决了一个小小的编程问题。由于本人并非该专业的学生，对于专业问题就门外汉了，本文主要以解决编程问题的思维来记录这个过程。

原问题（门外汉的表述）

背景

为了更好的解决城市洪涝问题，在某小区内铺设海绵城市设施（LID设施）并检验其效果，根据遗传算法的筛选，找到在不同降雨条件下小区中布置排水设施的最优方案。

优化目标

  1. 排水设施（LID）成本最优：即费用最少
  2. LID效果最佳：即径流削减率最大
  3. LID景观效果最好：即美观程度最佳

约束条件

另行讨论。

算法思路

a) 寻找目标函数
b) 初始化
c) 进行适应度概率计算
d) 优胜劣汰（保留优项）
e) 重组、交叉、变异
f) 化为表现型，适应度概率计算
g) 迭代至满足要求
h) 按照适应度高低输出结果
在这里插入图片描述

参考案例

a) 将子汇水区看作是基因，用1和0来表示子汇水区上面是否有排水设施（LID设施）。比如设定1表示该子汇水区存在LID设施（LID类型可能有好几种），0则代表没有设施。
b) 利用swmm软件计算初设种群的目标函数作没有优化时的对照。
c) 交叉、重组、变异之后得到子代。
交叉：例如设定交叉概率为0.5，当P>0.5时，子汇水区1与子汇水区2的LID布设方式进行交换。P<0.5时，保留原状。
变异：例如设定变异概率为0.5，当P<0.5时，对子汇水区1中的LID类型做变换。当0.8>P>0.5时，子汇水区上的LID设施不做改变。当P>0.8时，子汇水区1插入（或增设）一种LID设施。
d) 父代与子代合并，得到2n个个体。
e) 2n个个体进行非支配排序。
f) 计算拥挤度，选择拥挤度较大的目标进行探索。
g) 根据计算出的排序以及拥挤度筛选新的父代。
h) 直至达到设定的循环代数。找到pareto曲线，这条曲线上面的点都是最优解。

在这里插入图片描述

预期结果

最后要求得到pareto前沿曲线，即第一序列曲线（最优前沿曲线），并找几个方案进行比较，可以找最极端的方案，分析最特殊的点（例如两头的点）如：当完全不考虑成本时，或当追求成本最小时，可能lid表现并不很好，但还是有一定效果。

问题提炼（编程角度）

三个min/max优化函数

LID成本
在这里插入图片描述
设置LID前后效果对比
在这里插入图片描述
在这里插入图片描述

三个约束（实际还添加了优化函数2的非负约束）

在这里插入图片描述

解决问题

因为流程大致是修改输入inp文件、通过inp文件模拟swmm得到输出out文件、根据out文件判断子代是否保留等、再修改输入inp文件，不断重复。swmm是个开源程序，手动实现迭代进化肯定不现实，一开始把目标放在了matswmm上。然而所需的几个关键指标不能通过matswmm直接进行提取和修改，非常繁琐，则考虑通过pyswmm实现。此外由于该流程如何实现NSGA-Ⅱ不是重点，则直接使用pymoo库实现算法。

pyswmm

GitHub地址在这：pymoo官方仓库
此外还有官方的简单教程：官方文档
在本问题中需要提取并修改LID的面积，仿真，得到子汇水区的径流量，得到节点的超载时间
这里先给出test文件代码，里面是进行pyswmm调用的测试，部分地方给出了注释。

"""
@File:test.py
@Author:MJZJZ
@Version:1.0
@Date:2023/04/16
@Time:10:55
"""
import numpy as np
# pyswmm需要的module
from pyswmm import Simulation, Nodes, Subcatchments, LidGroups

x = [2.66, 0.49, 0.68, 0.23, 0.06, 0.04, 0.44, 0.04, 0.03, 0.07, 0.72, 4.81, 0.47, 0.34,
     3.29, 0.01, 0.09, 0.21, 0.49, 2.94, 0.26, 3.01, 0.06, 0.05, 2.34, 0.20, 0.02, 0.67,
     0.13, 0.19, 0.00, 0.02, 0.17, 5.07, 0.28, 2.62, 0.09, 4.41, 2.68, 0.20, 6.28, 2.75,
     0.10, 10.27, 5.00, 1.57, 0.78, 0.02, 0.50, 0.07, 0.92, 0.03, 6.03, 0.40, 0.20, 0.80,
     0.01, 1.11, 0.02, 0.00, 0.07]  # 测试
# 读取inp文件作为仿真对象
with Simulation(r'C:\Users\13618\Desktop\inp\P=10.inp') as sim:
    LidGroups(sim)['H00000'][0].unit_area = x[0]  # 设定inp文件里的LID面积
    LidGroups(sim)['H00000'][1].unit_area = x[1]  # 该文件里不同subcatchments（子汇水区）的LID项数不一样 手动添加所有项
    LidGroups(sim)['H00001'][0].unit_area = x[2]  # 对应inp文件的[LID_USAGE]
    LidGroups(sim)['H00001'][1].unit_area = x[3]
    LidGroups(sim)['H00002'][0].unit_area = x[4]
    LidGroups(sim)['H00004'][0].unit_area = x[5]
    LidGroups(sim)['H00005'][0].unit_area = x[6]
    LidGroups(sim)['H00006'][0].unit_area = x[7]
    LidGroups(sim)['H00007'][0].unit_area = x[8]
    LidGroups(sim)['H00008'][0].unit_area = x[9]
    LidGroups(sim)['H00009'][0].unit_area = x[10]
    LidGroups(sim)['H00009'][1].unit_area = x[11]
    LidGroups(sim)['H00010'][0].unit_area = x[12]
    LidGroups(sim)['H00011'][0].unit_area = x[13]
    LidGroups(sim)['H00012'][0].unit_area = x[14]
    LidGroups(sim)['H00013'][0].unit_area = x[15]
    LidGroups(sim)['H00014'][0].unit_area = x[16]
    LidGroups(sim)['H00015'][0].unit_area = x[17]
    LidGroups(sim)['H00015'][1].unit_area = x[18]
    LidGroups(sim)['H00016'][0].unit_area = x[19]
    LidGroups(sim)['H00017'][0].unit_area = x[20]
    LidGroups(sim)['H00018'][0].unit_area = x[21]
    LidGroups(sim)['H00019'][0].unit_area = x[22]
    LidGroups(sim)['H00020'][0].unit_area = x[23]
    LidGroups(sim)['H00021'][0].unit_area = x[24]
    LidGroups(sim)['H00021'][1].unit_area = x[25]
    LidGroups(sim)['H00022'][0].unit_area = x[26]
    LidGroups(sim)['H00023'][0].unit_area = x[27]
    LidGroups(sim)['H00024'][0].unit_area = x[28]
    LidGroups(sim)['H00025'][0].unit_area = x[29]
    LidGroups(sim)['H00026'][0].unit_area = x[30]
    LidGroups(sim)['H00027'][0].unit_area = x[31]
    LidGroups(sim)['H00029'][0].unit_area = x[32]
    LidGroups(sim)['H00029'][1].unit_area = x[33]
    LidGroups(sim)['H00031'][0].unit_area = x[34]
    LidGroups(sim)['H00032'][0].unit_area = x[35]
    LidGroups(sim)['H00032'][1].unit_area = x[36]
    LidGroups(sim)['H00034'][0].unit_area = x[37]
    LidGroups(sim)['H00035'][0].unit_area = x[38]
    LidGroups(sim)['H00036'][0].unit_area = x[39]
    LidGroups(sim)['H00036'][1].unit_area = x[40]
    LidGroups(sim)['H00037'][0].unit_area = x[41]
    LidGroups(sim)['H00038'][0].unit_area = x[42]
    LidGroups(sim)['H00038'][1].unit_area = x[43]
    LidGroups(sim)['H00039'][0].unit_area = x[44]
    LidGroups(sim)['H00041'][0].unit_area = x[45]
    LidGroups(sim)['H00044'][0].unit_area = x[46]
    LidGroups(sim)['H00046'][0].unit_area = x[47]
    LidGroups(sim)['H00047'][0].unit_area = x[48]
    LidGroups(sim)['H00047'][1].unit_area = x[49]
    LidGroups(sim)['H00048'][0].unit_area = x[50]
    LidGroups(sim)['H00048'][1].unit_area = x[51]
    LidGroups(sim)['H00049'][0].unit_area = x[52]
    LidGroups(sim)['H00051'][0].unit_area = x[53]
    LidGroups(sim)['H00052'][0].unit_area = x[54]
    LidGroups(sim)['H00053'][0].unit_area = x[55]
    LidGroups(sim)['H00053'][1].unit_area = x[56]
    LidGroups(sim)['H00054'][0].unit_area = x[57]
    LidGroups(sim)['H00054'][1].unit_area = x[58]
    LidGroups(sim)['H00055'][0].unit_area = x[59]
    LidGroups(sim)['H00056'][0].unit_area = x[60]
    # 开始仿真
    for step in enumerate(sim):
        pass  # 若需要每步的仿真都修改 则在这里添加代码
    # 本问题只需完整仿真完 直接pass即可
    # 仿真结束
    lid_area = np.array(np.zeros(61))  # 便于后续使用所有子汇水区LID面积 放在数组里
    lid_area[0] = LidGroups(sim)['H00000'][0].unit_area  # 仿真后提取各子汇水区LID面积（实际则为上面设定的x，这是为了方便测试）
    lid_area[1] = LidGroups(sim)['H00000'][1].unit_area
    lid_area[2] = LidGroups(sim)['H00001'][0].unit_area
    lid_area[3] = LidGroups(sim)['H00001'][1].unit_area
    lid_area[4] = LidGroups(sim)['H00002'][0].unit_area
    lid_area[5] = LidGroups(sim)['H00004'][0].unit_area
    lid_area[6] = LidGroups(sim)['H00005'][0].unit_area
    lid_area[7] = LidGroups(sim)['H00006'][0].unit_area
    lid_area[8] = LidGroups(sim)['H00007'][0].unit_area
    lid_area[9] = LidGroups(sim)['H00008'][0].unit_area
    lid_area[10] = LidGroups(sim)['H00009'][0].unit_area
    lid_area[11] = LidGroups(sim)['H00009'][1].unit_area
    lid_area[12] = LidGroups(sim)['H00010'][0].unit_area
    lid_area[13] = LidGroups(sim)['H00011'][0].unit_area
    lid_area[14] = LidGroups(sim)['H00012'][0].unit_area
    lid_area[15] = LidGroups(sim)['H00013'][0].unit_area
    lid_area[16] = LidGroups(sim)['H00014'][0].unit_area
    lid_area[17] = LidGroups(sim)['H00015'][0].unit_area
    lid_area[18] = LidGroups(sim)['H00015'][1].unit_area
    lid_area[19] = LidGroups(sim)['H00016'][0].unit_area
    lid_area[20] = LidGroups(sim)['H00017'][0].unit_area
    lid_area[21] = LidGroups(sim)['H00018'][0].unit_area
    lid_area[22] = LidGroups(sim)['H00019'][0].unit_area
    lid_area[23] = LidGroups(sim)['H00020'][0].unit_area
    lid_area[24] = LidGroups(sim)['H00021'][0].unit_area
    lid_area[25] = LidGroups(sim)['H00021'][1].unit_area
    lid_area[26] = LidGroups(sim)['H00022'][0].unit_area
    lid_area[27] = LidGroups(sim)['H00023'][0].unit_area
    lid_area[28] = LidGroups(sim)['H00024'][0].unit_area
    lid_area[29] = LidGroups(sim)['H00025'][0].unit_area
    lid_area[30] = LidGroups(sim)['H00026'][0].unit_area
    lid_area[31] = LidGroups(sim)['H00027'][0].unit_area
    lid_area[32] = LidGroups(sim)['H00029'][0].unit_area
    lid_area[33] = LidGroups(sim)['H00029'][1].unit_area
    lid_area[34] = LidGroups(sim)['H00031'][0].unit_area
    lid_area[35] = LidGroups(sim)['H00032'][0].unit_area
    lid_area[36] = LidGroups(sim)['H00032'][1].unit_area
    lid_area[37] = LidGroups(sim)['H00034'][0].unit_area
    lid_area[38] = LidGroups(sim)['H00035'][0].unit_area
    lid_area[39] = LidGroups(sim)['H00036'][0].unit_area
    lid_area[40] = LidGroups(sim)['H00036'][1].unit_area
    lid_area[41] = LidGroups(sim)['H00037'][0].unit_area
    lid_area[42] = LidGroups(sim)['H00038'][0].unit_area
    lid_area[43] = LidGroups(sim)['H00038'][1].unit_area
    lid_area[44] = LidGroups(sim)['H00039'][0].unit_area
    lid_area[45] = LidGroups(sim)['H00041'][0].unit_area
    lid_area[46] = LidGroups(sim)['H00044'][0].unit_area
    lid_area[47] = LidGroups(sim)['H00046'][0].unit_area
    lid_area[48] = LidGroups(sim)['H00047'][0].unit_area
    lid_area[49] = LidGroups(sim)['H00047'][1].unit_area
    lid_area[50] = LidGroups(sim)['H00048'][0].unit_area
    lid_area[51] = LidGroups(sim)['H00048'][1].unit_area
    lid_area[52] = LidGroups(sim)['H00049'][0].unit_area
    lid_area[53] = LidGroups(sim)['H00051'][0].unit_area
    lid_area[54] = LidGroups(sim)['H00052'][0].unit_area
    lid_area[55] = LidGroups(sim)['H00053'][0].unit_area
    lid_area[56] = LidGroups(sim)['H00053'][1].unit_area
    lid_area[57] = LidGroups(sim)['H00054'][0].unit_area
    lid_area[58] = LidGroups(sim)['H00054'][1].unit_area
    lid_area[59] = LidGroups(sim)['H00055'][0].unit_area
    lid_area[60] = LidGroups(sim)['H00056'][0].unit_area
    # 下面是优化函数1的计算 需要各类型LID的面积分别求和
    a_lswd = lid_area[1] + lid_area[3] + lid_area[5] + lid_area[6] + lid_area[7] \
             + lid_area[10] + lid_area[12] + lid_area[13] + lid_area[14] + lid_area[17] \
             + lid_area[19] + lid_area[20] + lid_area[21] + lid_area[24] + lid_area[26] \
             + lid_area[27] + lid_area[28] + lid_area[32] + lid_area[34] + lid_area[39] \
             + lid_area[42] + lid_area[46] + lid_area[47] + lid_area[48] + lid_area[50] \
             + lid_area[52] + lid_area[53] + lid_area[54] + lid_area[55] + lid_area[57] \
             + lid_area[59] + lid_area[60]

    b_xcsld = lid_area[0] + lid_area[2] + lid_area[11] + lid_area[18] + lid_area[25] \
              + lid_area[33] + lid_area[35] + lid_area[37] + lid_area[38] + lid_area[40] \
              + lid_area[41] + lid_area[43] + lid_area[44] + lid_area[45]

    d_tslm = lid_area[4] + lid_area[8] + lid_area[9] + lid_area[15] + lid_area[16] \
             + lid_area[22] + lid_area[23] + lid_area[29] + lid_area[30] + lid_area[31] \
             + lid_area[36] + lid_area[49] + lid_area[51] + lid_area[56] + lid_area[58]

    # max R 径流量
    # 目前的径流量runoff只在subcatchment的statistic函数里保存 需要单独提取出来 单位不一致 需要单独换算
    # 详见https://github.com/OpenWaterAnalytics/pyswmm/discussions/416 开发团队已表示pyswmm2将会认真考虑该问题
    subcatchment_runoff = np.array(np.zeros(57))
    for i in range(10):
        subcatchment_runoff[i] = Subcatchments(sim)["H0000" + str(i)].statistics['runoff'] / Subcatchments(sim)[
            "H0000" + str(i)].area / 10
    for i in range(10, 57):
        subcatchment_runoff[i] = Subcatchments(sim)["H000" + str(i)].statistics['runoff'] / Subcatchments(sim)[
            "H000" + str(i)].area / 10
    R = 0
    # global Pre_runoff
    # for i in range(57):
    #     R = R + (subcatchment_runoff[i] - Pre_runoff[i]) / Pre_runoff[i]
    # Pre_runoff = subcatchment_runoff
    for i in range(57):
        R = R + subcatchment_runoff[i]

        print(R)

    # min T 节点超载时间
    # 节点超载时间surcharge_duration同样在node的statistic函数里保存 这些都在官网的文档里有描述函数用法
    node_surcharged_duration = np.array(np.zeros(57))
    for i in range(1, 10):
        node_surcharged_duration[i - 1] = Nodes(sim)[str("J00" + str(i))].statistics['surcharge_duration']
    for i in range(10, 58):
        node_surcharged_duration[i - 1] = Nodes(sim)["J0" + str(i)].statistics['surcharge_duration']
    T = np.sum(node_surcharged_duration) / 57
    
    # 计算三个优化函数的测试
    f1 = 221.64 * a_lswd + 110.78 * b_xcsld + 303.13 * d_tslm
    f2 = R / 4037.39 - 1
    f3 = T
    # 一些测试
    print("\narea1:{} ".format(lid_area[0]))
    print("area2:{} ".format(lid_area[1]))
    print("area3:{}\n ".format(lid_area[2]))

    print(f1)
    print(f2)
    print(f3)
    # 更新 养成好习惯 关闭仿真
    sim.close()

pymoo

GitHub地址：pymoo官方仓库
官网：官方文档
用途可以说非常多了，本人基本是对着官方文档，自定义ElementwiseProblem，调用了NSGA-Ⅱ算法，将输出的Design space values、Objective spaces values、Constraint values分别保存下来。测试用例基本是按官网给的，这里就不贴了。

完整代码供参考

注：完整代码的目标函数进行过修改，原目标函数的实现可以看上面pyswmm部分，修改后的具体目标函数就不放了，感谢认真读者的提醒。另，pyswmm库也进行过改动，作者水平也有所提高，重新写一遍代码可读性和简洁性肯定更高一些，这里保留原代码仅为了留有纪念，些许笨拙地方还请海涵。

"""
@File:swmm_nsga2.py
@Author:MJZJZ
@Version:1.0
@Date:2023/04/03
@Time:14:31
"""
import numpy as np

from pyswmm import Simulation, Nodes, Subcatchments, LidGroups
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.core.problem import ElementwiseProblem
from pymoo.optimize import minimize
from pymoo.visualization.scatter import Scatter


# Pre_runoff = [28.84793721, 31.06321707, 11.96966298, 29.81131282, 22.76600575, 20.5654663,
#               32.03460237, 11.54445679, 11.39371698, 29.31124099, 26.40671222, 23.23278734,
#               29.09242243, 14.4631284, 13.2157434, 29.24380799, 26.48086273, 19.48246667,
#               29.14060795, 11.70136631, 12.75527874, 29.25169541, 26.61230667, 27.28638322,
#               26.29139334, 12.81874815, 13.18556032, 18.39743068, 26.2929191, 29.45805644,
#               29.632263, 24.50601238, 22.87947793, 26.18072453, 30.5640098, 28.75785111,
#               29.25407989, 28.78694888, 29.23303827, 30.4157768, 29.4229525, 29.88916333,
#               30.6586755, 29.72818431, 29.14815945, 32.80225754, 29.03260243, 15.24535308,
#               22.10563715, 28.67348727, 29.36437599, 31.8935914, 29.07932792, 23.60885133,
#               23.93630402, 29.05591441, 29.04791563]

# 用ElementwiseProblem自定义优化问题
class MyProblem(ElementwiseProblem):
    #  初始化
    def __init__(self):
        super().__init__(n_var=61,  # 自变量个数
                         n_obj=3,  # 优化目标函数个数
                         n_ieq_constr=4,  # 约束个数
                         xl=0.0,  # 自变量约束 这里上下限一样就直接用一个数来表示 否则可以用数组的形式
                         xu=300.0)

    #  重写问题 这里直接在_evaluate下写pyswmm的部分
    def _evaluate(self, x, out, *args, **kwargs):
        with Simulation(r'C:\Users\13618\Desktop\inp\P=1.inp') as sim:
            LidGroups(sim)['H00000'][0].unit_area = x[0]  # 设定inp文件里的LID面积
            LidGroups(sim)['H00000'][1].unit_area = x[1]  # 该文件里不同subcatchments（子汇水区）的LID项数不一样 手动添加所有项
            LidGroups(sim)['H00001'][0].unit_area = x[2]  # 对应inp文件的[LID_USAGE]
            LidGroups(sim)['H00001'][1].unit_area = x[3]
            LidGroups(sim)['H00002'][0].unit_area = x[4]
            LidGroups(sim)['H00004'][0].unit_area = x[5]
            LidGroups(sim)['H00005'][0].unit_area = x[6]
            LidGroups(sim)['H00006'][0].unit_area = x[7]
            LidGroups(sim)['H00007'][0].unit_area = x[8]
            LidGroups(sim)['H00008'][0].unit_area = x[9]
            LidGroups(sim)['H00009'][0].unit_area = x[10]
            LidGroups(sim)['H00009'][1].unit_area = x[11]
            LidGroups(sim)['H00010'][0].unit_area = x[12]
            LidGroups(sim)['H00011'][0].unit_area = x[13]
            LidGroups(sim)['H00012'][0].unit_area = x[14]
            LidGroups(sim)['H00013'][0].unit_area = x[15]
            LidGroups(sim)['H00014'][0].unit_area = x[16]
            LidGroups(sim)['H00015'][0].unit_area = x[17]
            LidGroups(sim)['H00015'][1].unit_area = x[18]
            LidGroups(sim)['H00016'][0].unit_area = x[19]
            LidGroups(sim)['H00017'][0].unit_area = x[20]
            LidGroups(sim)['H00018'][0].unit_area = x[21]
            LidGroups(sim)['H00019'][0].unit_area = x[22]
            LidGroups(sim)['H00020'][0].unit_area = x[23]
            LidGroups(sim)['H00021'][0].unit_area = x[24]
            LidGroups(sim)['H00021'][1].unit_area = x[25]
            LidGroups(sim)['H00022'][0].unit_area = x[26]
            LidGroups(sim)['H00023'][0].unit_area = x[27]
            LidGroups(sim)['H00024'][0].unit_area = x[28]
            LidGroups(sim)['H00025'][0].unit_area = x[29]
            LidGroups(sim)['H00026'][0].unit_area = x[30]
            LidGroups(sim)['H00027'][0].unit_area = x[31]
            LidGroups(sim)['H00029'][0].unit_area = x[32]
            LidGroups(sim)['H00029'][1].unit_area = x[33]
            LidGroups(sim)['H00031'][0].unit_area = x[34]
            LidGroups(sim)['H00032'][0].unit_area = x[35]
            LidGroups(sim)['H00032'][1].unit_area = x[36]
            LidGroups(sim)['H00034'][0].unit_area = x[37]
            LidGroups(sim)['H00035'][0].unit_area = x[38]
            LidGroups(sim)['H00036'][0].unit_area = x[39]
            LidGroups(sim)['H00036'][1].unit_area = x[40]
            LidGroups(sim)['H00037'][0].unit_area = x[41]
            LidGroups(sim)['H00038'][0].unit_area = x[42]
            LidGroups(sim)['H00038'][1].unit_area = x[43]
            LidGroups(sim)['H00039'][0].unit_area = x[44]
            LidGroups(sim)['H00041'][0].unit_area = x[45]
            LidGroups(sim)['H00044'][0].unit_area = x[46]
            LidGroups(sim)['H00046'][0].unit_area = x[47]
            LidGroups(sim)['H00047'][0].unit_area = x[48]
            LidGroups(sim)['H00047'][1].unit_area = x[49]
            LidGroups(sim)['H00048'][0].unit_area = x[50]
            LidGroups(sim)['H00048'][1].unit_area = x[51]
            LidGroups(sim)['H00049'][0].unit_area = x[52]
            LidGroups(sim)['H00051'][0].unit_area = x[53]
            LidGroups(sim)['H00052'][0].unit_area = x[54]
            LidGroups(sim)['H00053'][0].unit_area = x[55]
            LidGroups(sim)['H00053'][1].unit_area = x[56]
            LidGroups(sim)['H00054'][0].unit_area = x[57]
            LidGroups(sim)['H00054'][1].unit_area = x[58]
            LidGroups(sim)['H00055'][0].unit_area = x[59]
            LidGroups(sim)['H00056'][0].unit_area = x[60]
            # 开始仿真
            for step in enumerate(sim):
                pass  # 若需要每步的仿真都修改 则在这里添加代码
            # 本问题只需完整仿真完 直接pass即可
            # 仿真结束
            lid_area = np.array(np.zeros(61))  # 便于后续使用所有子汇水区LID面积 放在数组里
            lid_area[0] = LidGroups(sim)['H00000'][0].unit_area  # 仿真后提取各子汇水区LID面积
            lid_area[1] = LidGroups(sim)['H00000'][1].unit_area
            lid_area[2] = LidGroups(sim)['H00001'][0].unit_area
            lid_area[3] = LidGroups(sim)['H00001'][1].unit_area
            lid_area[4] = LidGroups(sim)['H00002'][0].unit_area
            lid_area[5] = LidGroups(sim)['H00004'][0].unit_area
            lid_area[6] = LidGroups(sim)['H00005'][0].unit_area
            lid_area[7] = LidGroups(sim)['H00006'][0].unit_area
            lid_area[8] = LidGroups(sim)['H00007'][0].unit_area
            lid_area[9] = LidGroups(sim)['H00008'][0].unit_area
            lid_area[10] = LidGroups(sim)['H00009'][0].unit_area
            lid_area[11] = LidGroups(sim)['H00009'][1].unit_area
            lid_area[12] = LidGroups(sim)['H00010'][0].unit_area
            lid_area[13] = LidGroups(sim)['H00011'][0].unit_area
            lid_area[14] = LidGroups(sim)['H00012'][0].unit_area
            lid_area[15] = LidGroups(sim)['H00013'][0].unit_area
            lid_area[16] = LidGroups(sim)['H00014'][0].unit_area
            lid_area[17] = LidGroups(sim)['H00015'][0].unit_area
            lid_area[18] = LidGroups(sim)['H00015'][1].unit_area
            lid_area[19] = LidGroups(sim)['H00016'][0].unit_area
            lid_area[20] = LidGroups(sim)['H00017'][0].unit_area
            lid_area[21] = LidGroups(sim)['H00018'][0].unit_area
            lid_area[22] = LidGroups(sim)['H00019'][0].unit_area
            lid_area[23] = LidGroups(sim)['H00020'][0].unit_area
            lid_area[24] = LidGroups(sim)['H00021'][0].unit_area
            lid_area[25] = LidGroups(sim)['H00021'][1].unit_area
            lid_area[26] = LidGroups(sim)['H00022'][0].unit_area
            lid_area[27] = LidGroups(sim)['H00023'][0].unit_area
            lid_area[28] = LidGroups(sim)['H00024'][0].unit_area
            lid_area[29] = LidGroups(sim)['H00025'][0].unit_area
            lid_area[30] = LidGroups(sim)['H00026'][0].unit_area
            lid_area[31] = LidGroups(sim)['H00027'][0].unit_area
            lid_area[32] = LidGroups(sim)['H00029'][0].unit_area
            lid_area[33] = LidGroups(sim)['H00029'][1].unit_area
            lid_area[34] = LidGroups(sim)['H00031'][0].unit_area
            lid_area[35] = LidGroups(sim)['H00032'][0].unit_area
            lid_area[36] = LidGroups(sim)['H00032'][1].unit_area
            lid_area[37] = LidGroups(sim)['H00034'][0].unit_area
            lid_area[38] = LidGroups(sim)['H00035'][0].unit_area
            lid_area[39] = LidGroups(sim)['H00036'][0].unit_area
            lid_area[40] = LidGroups(sim)['H00036'][1].unit_area
            lid_area[41] = LidGroups(sim)['H00037'][0].unit_area
            lid_area[42] = LidGroups(sim)['H00038'][0].unit_area
            lid_area[43] = LidGroups(sim)['H00038'][1].unit_area
            lid_area[44] = LidGroups(sim)['H00039'][0].unit_area
            lid_area[45] = LidGroups(sim)['H00041'][0].unit_area
            lid_area[46] = LidGroups(sim)['H00044'][0].unit_area
            lid_area[47] = LidGroups(sim)['H00046'][0].unit_area
            lid_area[48] = LidGroups(sim)['H00047'][0].unit_area
            lid_area[49] = LidGroups(sim)['H00047'][1].unit_area
            lid_area[50] = LidGroups(sim)['H00048'][0].unit_area
            lid_area[51] = LidGroups(sim)['H00048'][1].unit_area
            lid_area[52] = LidGroups(sim)['H00049'][0].unit_area
            lid_area[53] = LidGroups(sim)['H00051'][0].unit_area
            lid_area[54] = LidGroups(sim)['H00052'][0].unit_area
            lid_area[55] = LidGroups(sim)['H00053'][0].unit_area
            lid_area[56] = LidGroups(sim)['H00053'][1].unit_area
            lid_area[57] = LidGroups(sim)['H00054'][0].unit_area
            lid_area[58] = LidGroups(sim)['H00054'][1].unit_area
            lid_area[59] = LidGroups(sim)['H00055'][0].unit_area
            lid_area[60] = LidGroups(sim)['H00056'][0].unit_area
            # 下面是优化函数1的计算 需要各类型LID的面积分别求和
            a_lswd = lid_area[1] + lid_area[3] + lid_area[5] + lid_area[6] + lid_area[7] \
                     + lid_area[10] + lid_area[12] + lid_area[13] + lid_area[14] + lid_area[17] \
                     + lid_area[19] + lid_area[20] + lid_area[21] + lid_area[24] + lid_area[26] \
                     + lid_area[27] + lid_area[28] + lid_area[32] + lid_area[34] + lid_area[39] \
                     + lid_area[42] + lid_area[46] + lid_area[47] + lid_area[48] + lid_area[50] \
                     + lid_area[52] + lid_area[53] + lid_area[54] + lid_area[55] + lid_area[57] \
                     + lid_area[59] + lid_area[60]

            b_xcsld = lid_area[0] + lid_area[2] + lid_area[11] + lid_area[18] + lid_area[25] \
                      + lid_area[33] + lid_area[35] + lid_area[37] + lid_area[38] + lid_area[40] \
                      + lid_area[41] + lid_area[43] + lid_area[44] + lid_area[45]

            d_tslm = lid_area[4] + lid_area[8] + lid_area[9] + lid_area[15] + lid_area[16] \
                     + lid_area[22] + lid_area[23] + lid_area[29] + lid_area[30] + lid_area[31] \
                     + lid_area[36] + lid_area[49] + lid_area[51] + lid_area[56] + lid_area[58]

            # max R 径流量
            # 目前的径流量runoff只在subcatchment的statistic函数里保存 需要单独提取出来 单位不一致 需要单独换算
            # 详见https://github.com/OpenWaterAnalytics/pyswmm/discussions/416 开发团队已表示pyswmm2将会认真考虑该问题
            subcatchment_runoff = np.array(np.zeros(57))
            for i in range(10):
                subcatchment_runoff[i] = Subcatchments(sim)["H0000" + str(i)].statistics['runoff'] / Subcatchments(sim)[
                    "H0000" + str(i)].area / 10
            for i in range(10, 57):
                subcatchment_runoff[i] = Subcatchments(sim)["H000" + str(i)].statistics['runoff'] / Subcatchments(sim)[
                    "H000" + str(i)].area / 10
            R = 0
            # global Pre_runoff
            # for i in range(57):
            #     R = R + (subcatchment_runoff[i] - Pre_runoff[i]) / Pre_runoff[i]
            # Pre_runoff = subcatchment_runoff
            for i in range(57):
                R = R + subcatchment_runoff[i]
            R = R / 1640 - 1

            # min T 节点超载时间
            # 节点超载时间surcharge_duration同样在node的statistic函数里保存 这些都在官网的文档里有描述函数用法
            node_surcharged_duration = np.array(np.zeros(57))
            for i in range(1, 10):
                node_surcharged_duration[i - 1] = Nodes(sim)[str("J00" + str(i))].statistics['surcharge_duration']
            for i in range(10, 58):
                node_surcharged_duration[i - 1] = Nodes(sim)["J0" + str(i)].statistics['surcharge_duration']
            T = np.sum(node_surcharged_duration) / 30

            # 注意优化目标函数与约束都只能是min的形式 这里上面的R已经全部添加了负号进行计算
            f1 = 221.64 * a_lswd + 110.78 * b_xcsld + 303.13 * d_tslm
            f2 = R
            f3 = T
            # 约束 同样是min的形式
            g1 = a_lswd - 5500
            g2 = b_xcsld - 14600
            g3 = d_tslm - 3800
            g4 = R
            # 将优化目标函数与约束分别赋值给F、G
            out["F"] = [f1, f2, f3]
            out["G"] = [g1, g2, g3, g4]
            # 更新 养成好习惯 关闭仿真
    		sim.close()
			#如需单步内仿真多个inp，则需多次使用sim和sim.close()

# 照搬pymoo即可
problem = MyProblem()

algorithm = NSGA2(pop_size=100)  # 种群大小100

res = minimize(problem,
               algorithm,
               ("n_gen", 500),  # 迭代500代
               verbose=True,  # true则在终端显示每代部分参数
               # save_history=True,
               seed=1)  # 这里并未改动具体的变异率等，详细参数可以转到NSGA2的类型声明查看
# val = [e.opt.get("F")[0] for e in res.history]
# import matplotlib.pyplot as plt
# plt.plot(np.arange(len(val)), val)
# plt.show()
print(res.exec_time)  # 打印耗时
np.savetxt(r'C:\Users\13618\Desktop\design.txt', res.X)  # 保存最后pop_size大小的结果 分别是design spaces values 即自变量的值
np.savetxt(r'C:\Users\13618\Desktop\objective.txt', res.F)  # objective spaces values 即优化目标函数的值
np.savetxt(r'C:\Users\13618\Desktop\constraint.txt', res.G)  # constraint values 即约束的值
# 用三维散点图展示结果
plot = Scatter()
plot.add(res.F, edgecolor="red", facecolor="none")
plot.show()

总结

对编程能力提升倒是不大，主要还是解决问题的思维锻炼，很多问题并没有现成的途径直接解决，去官网仔细阅读函数说明和examples是非常有帮助的。另外pyswmm对比matswmm真是好用极了，确实还得找仍在更新维护的工具才顺手，原开发团队很快也会发布pyswmm2，应该上手更方便些。
如果这篇文章能帮助你，感谢你的点赞收藏支持，欢迎该方面问题的任何交流。

mjzjz_c
关注
  • 24
    点赞
  • 踩
  • 41
    收藏
    觉得还不错? 一键收藏
  • 16
    评论
swmm源代码及代码路线图和中文使用手册
07-09
swmm管网的源代码，包括核心算法及界面；附代码路线图文档，供代码解读这查看；中文使用手册，详细描述SWMM操作使用方法，供初学者使用。
SWMM5中文操作手册
08-12
SWMM5中文操作手册，从简单入手，一步步的教使用SWMM5
Python库 | pyswmm-0.6.2.tar.gz
03-07
python库。 资源全名：pyswmm-0.6.2.tar.gz
探索强大的Python接口：PySWMM
最新发布
gitblog_00055的博客
05-26 268
探索强大的Python接口：PySWMM 项目地址:https://gitcode.com/pyswmm/pyswmm 项目介绍 PySWMM是面向Python的EPA Stormwater Management Model（SWMM5）的软件包，提供了一个用于创建、操纵和研究复杂网络结构、动态特性和功能的强大工具。这个开源项目使得SWMM模型的操作变得更加便捷，为工程师、模型构建者和研究人员打开...
python预测控制_【模型工具】耦合python和 SWMM的城市排水系统模型预测算法
weixin_39819671的博客
12-11 748
原标题：【模型工具】耦合python和 SWMM的城市排水系统模型预测算法原文信息题目：Leveraging open source software and parallel computing for model predictive control of urban drainage systems using EPA-SWMM5作者：Jeffrey M. Sadlera, Jonathan...
用python解析SWMM的.OUT输出文件
MissWhite的博客
04-24 2164
SWMM中.OUT输出文件为二进制文件，用记事本打开会呈现如下图所示的乱码状态。但是.OUT文件中有很多重要结果信息（例如：节点全模拟时段水深、侧向来水、出口总流量等数据），用官方提供的pyswmm库对.OUT文件进行解析。将二进制数据输出为Excel。文件主要包括三种类型数据，节点类型数据、管道类型数据、汇水区类型数据。节点类型包括检查井、排口、调蓄池、分流器；管道类型包括管道、泵、孔口、堰、出口；汇水区包括汇水分区。
关于读取swmm的out文件的另一条路
weixin_44276743的博客
06-15 718
关于SWMM的out文件的处理方式？
多目标遗传算法NSGA-II原理详解及算法实现
热门推荐
不曾走远的博客
03-11 4万+
非支配排序遗传算法NSGA (Non-dominated Sorting Genetic Algorithms)是一种基于Pareto最优概念的遗传算法，它是众多的多目标优化遗传算法中体现Goldberg思想最直接的方法。该算法就是在基本遗传算法的基础上，对选择再生方法进行改进:将每个个体按照它们的支配与非支配关系进行分层，再做选择操作，从而使得该算法在多目标优化方面得到非常满意的结果。
SWMM对象模型
nekewang的博客
09-26 6905
图1-1  典型的城市排水系统图（来自EPA） 上图描绘了典型的城市排水系统所包含的要素。SWMM将该系统进行概化为几个主要环境体之间的一系列水流和物质流。这些环境体包括如下： 大气环境体产生降水并将污染物沉积到陆地表面环境体 陆地表面环境体以雨或雪的形式接收来自大气环境体的降水。它再以三种形式输出：1）蒸发回到大气环境体；2）渗透到下一层地表环境体；3）流动过程中表面径流和污染物流（传...
swmm5中文操作手册
08-03
学习过程中顺便翻译了SWMM的用户手册，限于翻译水平，译文中不足之处在所难免，不足之处敬请见谅，希望对大家有用，谢谢！
JAVA 对SWMM模型DLL调用
05-23
java调用SWMM模型的DLL 实现通过INP（输入文件）生成 rpt (报告文件)out（输出文件），并读取输出文件里面的数据。
swmm模型软件
01-22
（Environmental Protection Agency，环境保护署） SWMM（storm water management model，暴雨洪水管理模型）是一个动态的降水-径流模拟模型，主要用于模拟城市某一单一降水事件或长期的水量和水质模拟。其径流模块部分综合处理各子流域所发生的降水，径流和污染负荷。其汇流模块部分则通过管网、渠道、蓄水和处理设施、水泵、调节闸等进行水量传输。该模型可以跟踪模拟不同时间步长任意时刻每个子流域所产生径流的水质和水量，以及每个管道和河道中水的流量、水深及水质等情况。 SWMM模型的核心水力模块主要包括径流模块(runoff block)、输送模块（transport block）、扩展传送模块（extend transport block）、调蓄/处理模块（storage/treatment block）和受纳水体模块 receiving water block等。 此外，SWMM还包括多个服务模块，主要有降雨模块（rain block），统计模块（statistics block）、绘图模块（graph block），联合模块（combine block）、运行模块（temp block）和执行模块（executive block）等。
EPA SWMM 暴雨管理模型
10-29
EPA（Environmental Protection Agency，环境保护署） SWMM（storm water management model，暴雨洪水管理模型）是一个开源的降水-径流模拟模型，主要用于模拟城市某一单一降水事件或长期的水量和水质模拟。其径流模块部分综合处理各子流域所发生的降水，径流和污染负荷。其汇流模块部分则通过管网、渠道、蓄水和处理设施、水泵、调节闸等进行水量传输。该模型可以跟踪模拟不同时间步长任意时刻每个子流域所产生径流的水质和水量，以及每个管道和河道中水的流量、水深及水质等情况。 SWMM自1971年开发以来，已经经历过多次升级。在世界范围内广泛应用于城市地区的暴雨洪水、合流式下水道、排污管道以及其它排水系统的规划、分析和设计，在其它非城市区域也有广泛的应用。当前最新版本5.0是在以前版本基础上进行了全新升级的结果，可以在Windows操作系统下运行SWMM5提供了一个宽松的综合性环境，可以对研究区输入的数据进行编辑、模拟水文、水力和水质情况，并可以用多种形式对结果进行显示，包括对排水区域和系统输水路线进行彩色编码，提供结果的时间序列曲线和图表、坡面图以及统计频率的分析结果。 最新的版本开发由国家环境保护署国家风险管理研究中心实验室下属的供水和水资源研究中心资助，同时也得到了来自CDM咨询公司的协助。
提取cad管线信息到swmm模型工具
02-17
该cad插件可以提取dwg图纸中的管网信息，并直接转成swmm模型的inp文件格式，提取的内容包括管线的拓扑关系及管径、管长、标高、坐标；会生成inp文件，可直接拖入swmm中显示。包含[OUTFALLS]，[JUNCTIONS]，[CONDUITS]，[XSECTIONS]，[COORDINATES]，[VERTICES]部分。
基于NSGA-Ⅱ算法的雨水管网优化模型研究_毕业论文.pdf
10-29
基于NSGA-Ⅱ算法的雨水管网优化模型研究是雨水管网优化设计的一种新方法，该方法通过使用快速非支配排序遗传算法（NSGA-Ⅱ）和SWMM软件建立雨水管网模型，实现雨水管网优化设计。该方法可以解决雨水管网设计中存在的...
基于SWMM的城市雨洪管理模型构建及应用研究.ppt
05-28
基于SWMM的城市雨洪管理模型构建及应用研究.ppt基于SWMM的城市雨洪管理模型构建及应用研究.ppt基于SWMM的城市雨洪管理模型构建及应用研究.ppt基于SWMM的城市雨洪管理模型构建及应用研究.ppt基于SWMM的城市雨洪管理...
元胞自动机模型代码matlab-SWMM-2DCA:此代号命名为雨水管理模型二维元胞自动机(SWMM-2DCA)，用于模拟城市排水系统与地表径
06-17
元细胞自动机模型代码matlab SWMM-2DCA 此代号命名为雨水管理模型二维元胞自动机 (SWMM-2DCA)，用于模拟城市排水系统与地表径流之间的洪水特性和流动相互作用。 代码是我的 MATLAB。 有关科学和技术信息，请参阅
Python库 | SWMM5-5.1.12-cp37-cp37m-win_amd64.whl
02-28
python库。 资源全名：SWMM5-5.1.12-cp37-cp37m-win_amd64.whl
基于swmm的城市雨洪管理模型构建及应用研究.pdf
10-16
《基于SWMM的城市雨洪管理模型构建及应用研究》是一篇关于城市雨洪管理模型构建和应用的研究论文。文章主要通过利用Storm Water Management Model（SWMM）进行城市雨洪管理模型的构建和应用研究。 在论文中，研究者首先对SWMM进行了简要介绍，SWMM是一种广泛应用于城市雨水系统规划和设计的模型工具，可以模拟城市中雨水径流的产生和传输。 接着，研究者基于实际城市的数据和特征，构建了城市雨洪管理模型。该模型包括了城市的地理信息、建筑物分布、道路网络等，同时考虑了降雨过程和排水系统的运行规则。通过将这些元素结合起来，模型可以模拟出城市内雨水径流的路径和演变过程。 在模型的应用研究中，研究者利用了实际工程案例进行了仿真实验，并对模型的预测效果进行了验证。通过对比模型预测结果与实际观测数据，研究者发现模型的预测准确性较高，能够较好地模拟城市雨洪管理过程。 最后，研究者还对模型的一些应用进行了探讨，如可以用于城市雨洪管理规划、城市水资源利用优化等。他们提出了一些建议和对策，以改善城市雨洪管理模型的应用。 总之，《基于SWMM的城市雨洪管理模型构建及应用研究》这篇论文通过利用SWMM模型进行城市雨洪管理的研究，并深入探讨了模型的构建和应用。研究结果对于城市雨洪管理工作具有一定的指导意义，同时也为城市规划和设计提供了可靠的工具和参考。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 16
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

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

余额充值