import numpy as np from app.infra.epanet import run_project from app.native.wndb.core.database import ChangeSet from app.native.wndb.model.demands import get_demand, set_demand from app.native.wndb.model.options import get_option from app.native.wndb.model.patterns import get_pattern, set_pattern from app.native.wndb.model.pumps import get_pump from app.native.wndb.model.reservoirs import get_reservoir from app.native.wndb.model.status import get_status, set_status from app.native.wndb.model.tanks import get_tank, set_tank from app.native.wndb.model.times import get_time, set_time # from get_real_status import * from datetime import datetime from math import modf import json import pytz import time from uuid import UUID import logging from app.infra.db.postgresql.scada import ( ScadaElementMappings, load_realtime_element_mappings, ) from app.domain.time import parse_beijing_time, parse_clock_duration_seconds from app.native.wndb.core.connection import project_transaction from app.native.wndb.core.database import refresh_materialized_views_after_commit from app.infra.db.timescaledb.internal_queries import ( InternalQueries as TimescaleInternalQueries, ) from app.services.scheme_management import create_analysis_run, update_analysis_run from app.infra.db.timescaledb.internal_queries import ( InternalStorage as TimescaleInternalStorage, ) logging.basicConfig( level=logging.INFO, format="%(asctime)s - %(levelname)s - %(message)s" ) RESERVOIR_BASIC_HEIGHT = 250.35 def _primary_demand(demand_set: dict) -> dict: """Return sequence-zero demand, creating it when a junction has none.""" demands = demand_set.setdefault("demands", []) if not demands: demands.append({"demand": 0.0, "pattern": None, "category": None}) return demands[0] def _primary_demand_pattern(demand_set: dict) -> str: """Return the first configured pattern in deterministic sequence order.""" pattern = next( ( demand.get("pattern") for demand in demand_set.get("demands", []) if demand.get("pattern") ), None, ) if pattern is None: raise ValueError( f"Junction {demand_set.get('junction')!r} has no demand pattern" ) return str(pattern) def query_corresponding_element_id_and_query_id(name: str) -> ScadaElementMappings: """Return an immutable project-local SCADA-to-model mapping snapshot.""" return load_realtime_element_mappings(name) def get_pattern_index(cur_datetime: str, pattern_time_step: float) -> int: """ 根据给定的日期时间字符串，计算并返回对应的模式索引。 :param cur_datetime: str, 当前的日期时间字符串，格式为“YYYY-MM-DD HH:MM:SS”。 :return: int, 基于预定义的时间步长 PATTERN_TIME_STEP。 """ str_format = "%Y-%m-%d %H:%M:%S" dt = datetime.strptime(cur_datetime, str_format) hr = dt.hour mnt = dt.minute i = int((hr * 60 + mnt) / pattern_time_step) return i def get_pattern_index_str(current_time: str, pattern_time_step: float) -> str: """ 根据当前时间获取时间步长的模式索引，并将其格式化为“HH:MM:00”字符串。 :param current_time: str, 当前时间，格式为"YYYY-MM-DD HH:MM:SS" :return: str， 以“HH:MM:00”格式返回 """ i = get_pattern_index(current_time, pattern_time_step) [minN, hrN] = modf(i * pattern_time_step / 60) minN_str = str(int(minN * 60)) minN_str = minN_str.zfill(2) hrN_str = str(int(hrN)) hrN_str = hrN_str.zfill(2) str_i = "{}:{}:00".format(hrN_str, minN_str) return str_i def from_seconds_to_clock(secs: int) -> str: """ 从秒格式化为“HH:MM:00”字符串 :param secs: int，秒 :return: str， 以“HH:MM:00”格式返回 """ hrs = int(secs / 3600) minutes = int((secs - hrs * 3600) / 60) seconds = secs - hrs * 3600 - minutes * 60 hrs_str = str(hrs).zfill(2) minutes_str = str(minutes).zfill(2) seconds_str = str(seconds).zfill(2) str_clock = "{}:{}:{}".format(hrs_str, minutes_str, seconds_str) return str_clock def convert_time_format(original_time: str) -> str: """ 格式转换，将带时区的 ISO 8601 / RFC3339 时间转为北京时间的“YYYY-MM-DD HH:MM:SS” :param original_time: str，带显式时区的时间 :return: str，“2024-04-13 08:00:00”格式的时间 """ normalized_time = parse_beijing_time( original_time, field_name="modify_pattern_start_time" ) return normalized_time.replace(microsecond=0).strftime("%Y-%m-%d %H:%M:%S") def _apply_valve_control( project_name: str, valve_control: dict[str, dict] ) -> None: """Apply explicit valve status, setting, and opening controls.""" for valve_name, control in valve_control.items(): valve_status = get_status(project_name, valve_name) if "status" in control: valve_status["status"] = control["status"] if "setting" in control: valve_status["setting"] = control["setting"] if "k" in control: valve_k = control["k"] if valve_k == 0: valve_status["status"] = "CLOSED" else: valve_status["setting"] = 0.1036 * pow(valve_k, -3.105) cs = ChangeSet() cs.append(valve_status) set_status(project_name, cs) # 2025/01/11 def run_simulation( name: str, simulation_type: str, modify_pattern_start_time: str, modify_total_duration: int = 0, modify_reservoir_head_pattern: dict[str, list] = None, modify_tank_initial_level: dict[str, float] = None, modify_junction_base_demand: dict[str, float] = None, modify_junction_damand_pattern: dict[str, list] = None, modify_fixed_pump_pattern: dict[str, list] = None, modify_variable_pump_pattern: dict[str, list] = None, modify_valve_opening: dict[str, float] = None, scheme_type: str = None, scheme_name: str = None, result_db_name: str = None, valve_control: dict[str, dict] = None, scheme_username: str = "system", scheme_detail: dict | None = None, scada_mappings: ScadaElementMappings | None = None, ) -> UUID | None: """ 传入需要修改的参数，改变数据库中对应位置的值，然后计算，返回结果 :param name: 模型名称，数据库中对应的名字 :param simulation_type: 模拟的类型，realtime为实时模拟，修改原数据库；extended为多步长模拟，需要复制数据库 :param modify_pattern_start_time: 模拟开始时间,格式为'2024-11-25T09:00:00+08:00' :param modify_total_duration: 模拟总历时，秒 :param modify_reservoir_head_pattern: dict中包含多个水库模式，str为水库head_pattern的id，list为修改后的head_pattern :param modify_tank_initial_level: dict中包含多个水塔，str为水塔的id，float为修改后的initial_level :param modify_junction_base_demand: dict中包含多个节点，str为节点的id，float为修改后的base_demand :param modify_junction_damand_pattern: dict中包含多个节点模式，str为节点demand_pattern的id，list为修改后的demand_pattern :param modify_fixed_pump_pattern: dict中包含多个水泵模式，str为工频水泵的id，list为修改后的pattern :param modify_variable_pump_pattern: dict中包含多个水泵模式，str为变频水泵的id，list为修改后的pattern :param modify_valve_opening: dict中包含多个阀门开启度，str为阀门的id，float为修改后的阀门开启度 :param valve_control: dict中可分别指定阀门的status、setting和k；存在时优先于modify_valve_opening :param scheme_type: 模拟方案类型 :param scheme_name：模拟方案名称 :return: """ # 记录开始时间 time_cost_start = time.perf_counter() print("name", name) print("simulation_type", simulation_type) print("modify_pattern_start_time", modify_pattern_start_time) print("modify_total_duration", modify_total_duration) print("modify_reservoir_head_pattern", modify_reservoir_head_pattern) print("modify_tank_initial_level", modify_tank_initial_level) print("modify_junction_base_demand", modify_junction_base_demand) print( "{} -- Hydraulic simulation started.".format( datetime.now(pytz.timezone("Asia/Shanghai")).strftime("%Y-%m-%d %H:%M:%S") ) ) # 判断是实时模拟还是多步长模拟 # if simulation_type.upper() == 'REALTIME': # 实时模拟(修改原数据库) # name_c = name # elif simulation_type.upper() == 'EXTENDED': # 扩展模拟(复制数据库) # name_c = '_'.join([name, 'c']) # if have_project(name_c): # delete_project(name_c) # copy_project(name, name_c) # 备份项目 # else: # raise Exception('Incorrect simulation type, choose in (realtime, extended)') name_c = name with project_transaction(name_c): dic_time = get_time(name_c) print(dic_time) # 获取水力模拟步长，如’0:15:00‘ hydraulic_timestep = dic_time["HYDRAULIC TIMESTEP"] # 转换为分钟浮点数，兼容 EPANET 的 H:MM 和 H:MM:SS 写法 pattern_time_step = ( parse_clock_duration_seconds( hydraulic_timestep, field_name="HYDRAULIC TIMESTEP", ) / 60 ) project_scada = scada_mappings or load_realtime_element_mappings(name_c) # 对输入的时间参数进行处理 pattern_start_time = convert_time_format(modify_pattern_start_time) # 获取模拟开始时间是对应pattern的第几个数 modify_index = get_pattern_index(pattern_start_time, pattern_time_step) # 遍历水泵的pattern_id，并根据输入的pump_pattern修改pattern的值 # for pump_pattern_id in pump_pattern_ids: # # 检查pump_pattern中pump_pattern_id对应的第一个频率值是否为有效数字（非空、非NaN）。如果该值有效，则继续执行代码块。 # if not np.isnan(modify_pump_pattern[pump_pattern_id][0]): # # 取出数据库中的pattern # pump_pattern = get_pattern(name_c, get_pump(name_c, pump_pattern_id)['pattern']) # # 替换数据库中的pattern为modify_pump_pattern # pump_pattern['factors'][modify_index: modify_index + len(modify_pump_pattern[pump_pattern_id])] \ # = modify_pump_pattern[pump_pattern_id] # cs = ChangeSet() # cs.append(pump_pattern) # set_pattern(name_c, cs) # 修改模拟开始的时间 str_pattern_start = get_pattern_index_str( convert_time_format(modify_pattern_start_time), pattern_time_step ) dic_time = get_time(name_c) dic_time["PATTERN START"] = str_pattern_start dic_time["DURATION"] = from_seconds_to_clock(modify_total_duration) if simulation_type.upper() == "REALTIME": dic_time["DURATION"] = 0 cs = ChangeSet() cs.operations.append(dic_time) set_time(name_c, cs) # 根据SCADA实时数据进行修改，如果没有对应的SCADA数据，如未来的时间点，则不改变pg数据库的数据 if project_scada.reservoirs: # reservoirs_id = {'ZBBDJSCP000002': '2497', 'R00003': '2571'} # 1.获取reservoir的SCADA数据,形式如{'2497': '3.1231', '2571': '2.7387'} reservoir_SCADA_data_dict = TimescaleInternalQueries.query_scada_by_ids_time( device_ids=list(project_scada.reservoirs.values()), query_time=modify_pattern_start_time, db_name=name, ) # 2.构建出新字典，形式如{'ZBBDJSCP000002': '3.1231', 'R00003': '2.7387'} reservoir_dict = { key: reservoir_SCADA_data_dict[value] for key, value in project_scada.reservoirs.items() } # 3.修改reservoir液位模式 for reservoir_name, value in reservoir_dict.items(): if value and float(value) != 0: # 先根据reservoir获取对应的pattern，再对pattern进行修改 reservoir_pattern = get_pattern( name_c, get_reservoir(name_c, reservoir_name)["pattern"] ) reservoir_pattern["factors"][modify_index] = ( float(value) + RESERVOIR_BASIC_HEIGHT ) cs = ChangeSet() cs.append(reservoir_pattern) set_pattern(name_c, cs) if project_scada.tanks: # 修改tank初始液位 tank_SCADA_data_dict = TimescaleInternalQueries.query_scada_by_ids_time( device_ids=list(project_scada.tanks.values()), query_time=modify_pattern_start_time, db_name=name, ) tank_dict = { key: tank_SCADA_data_dict[value] for key, value in project_scada.tanks.items() } for tank_name, value in tank_dict.items(): if value and float(value) != 0: tank = get_tank(name_c, tank_name) tank["init_level"] = float(value) cs = ChangeSet() cs.append(tank) set_tank(name_c, cs) if project_scada.fixed_pumps: # 修改工频泵的pattern fixed_pump_SCADA_data_dict = TimescaleInternalQueries.query_scada_by_ids_time( device_ids=list(project_scada.fixed_pumps.values()), query_time=modify_pattern_start_time, db_name=name, ) # print(fixed_pump_SCADA_data_dict) fixed_pump_dict = { key: fixed_pump_SCADA_data_dict[value] for key, value in project_scada.fixed_pumps.items() } # print(fixed_pump_dict) for fixed_pump_name, value in fixed_pump_dict.items(): if value: pump_pattern = get_pattern( name_c, get_pump(name_c, fixed_pump_name)["pattern"] ) # print(pump_pattern) pump_pattern["factors"][modify_index] = float(value) # print(pump_pattern['factors'][modify_index]) cs = ChangeSet() cs.append(pump_pattern) set_pattern(name_c, cs) if project_scada.variable_pumps: # 修改变频泵的pattern variable_pump_SCADA_data_dict = ( TimescaleInternalQueries.query_scada_by_ids_time( device_ids=list(project_scada.variable_pumps.values()), query_time=modify_pattern_start_time, db_name=name, ) ) variable_pump_dict = { key: variable_pump_SCADA_data_dict[value] for key, value in project_scada.variable_pumps.items() } for variable_pump_name, value in variable_pump_dict.items(): if value: pump_pattern = get_pattern( name_c, get_pump(name_c, variable_pump_name)["pattern"] ) pump_pattern["factors"][modify_index] = float(value) / 50 cs = ChangeSet() cs.append(pump_pattern) set_pattern(name_c, cs) if project_scada.demand: # 基于实时数据，修改大用户节点的pattern demand_SCADA_data_dict = TimescaleInternalQueries.query_scada_by_ids_time( device_ids=list(project_scada.demand.values()), query_time=modify_pattern_start_time, db_name=name, ) demand_dict = { key: demand_SCADA_data_dict[value] for key, value in project_scada.demand.items() } for demand_name, value in demand_dict.items(): if value is not None and not np.isnan(float(value)): demand_set = get_demand(name_c, demand_name) demand_pattern = get_pattern( name_c, _primary_demand_pattern(demand_set) ) if get_option(name_c)["UNITS"] == "LPS": demand_pattern["factors"][modify_index] = ( float(value) / 3.6 ) # 默认SCADA数据获取的是流量单位是m3/h， 转换为 L/s elif get_option(name_c)["UNITS"] == "CMH": demand_pattern["factors"][modify_index] = float(value) cs = ChangeSet() cs.append(demand_pattern) set_pattern(name_c, cs) # 水质、压力实时数据使用方法待补充 ############################# # 显式请求参数覆盖实时设备数据，用于扩展模拟。 # 修改清水池(reservoir)液位的pattern if modify_reservoir_head_pattern: for reservoir_name in modify_reservoir_head_pattern.keys(): # 这句代码的作用是判断modify_reservoir_head_pattern[reservoir_name][0]是否不是NaN。 # 如果modify_reservoir_head_pattern[reservoir_name][0]不是NaN，则条件成立，代码块会执行 if not np.isnan(modify_reservoir_head_pattern[reservoir_name][0]): # 给 list 中的所有元素加上 RESERVOIR_BASIC_HEIGHT modified_values = [ value + RESERVOIR_BASIC_HEIGHT for value in modify_reservoir_head_pattern[reservoir_name] ] reservoir_pattern = get_pattern( name_c, get_reservoir(name_c, reservoir_name)["pattern"] ) reservoir_pattern["factors"][ modify_index : modify_index + len(modified_values) ] = modified_values cs = ChangeSet() cs.append(reservoir_pattern) set_pattern(name_c, cs) # 修改调节池(tank)初始液位 if modify_tank_initial_level: for tank_name in modify_tank_initial_level.keys(): if (not np.isnan(modify_tank_initial_level[tank_name])) and ( modify_tank_initial_level[tank_name] != 0 ): tank = get_tank(name_c, tank_name) tank["init_level"] = modify_tank_initial_level[tank_name] cs = ChangeSet() cs.append(tank) set_tank(name_c, cs) # 修改节点（junction）基础水量（demand） if modify_junction_base_demand: for junction_name in modify_junction_base_demand.keys(): if not np.isnan(modify_junction_base_demand[junction_name]): junction = get_demand(name_c, junction_name) _primary_demand(junction)["demand"] = ( modify_junction_base_demand[junction_name] ) cs = ChangeSet() cs.append(junction) set_demand(name_c, cs) # 修改节点（junction）的水量模式（pattern） if modify_junction_damand_pattern: for pattern_name in modify_junction_damand_pattern.keys(): if not np.isnan(modify_junction_damand_pattern[pattern_name][0]): junction_pattern = get_pattern(name_c, pattern_name) junction_pattern["factors"][ modify_index : modify_index + len(modify_junction_damand_pattern[pattern_name]) ] = modify_junction_damand_pattern[pattern_name] cs = ChangeSet() cs.append(junction_pattern) set_pattern(name_c, cs) # 修改工频水泵（fixed_pump）的pattern if modify_fixed_pump_pattern: for pump_name in modify_fixed_pump_pattern.keys(): if not np.isnan(modify_fixed_pump_pattern[pump_name][0]): pump_pattern = get_pattern( name_c, get_pump(name_c, pump_name)["pattern"] ) pump_pattern["factors"][ modify_index : modify_index + len(modify_fixed_pump_pattern[pump_name]) ] = modify_fixed_pump_pattern[pump_name] cs = ChangeSet() cs.append(pump_pattern) set_pattern(name_c, cs) # 修改变频水泵（variable_pump）的pattern if modify_variable_pump_pattern: for pump_name in modify_variable_pump_pattern.keys(): if not np.isnan(modify_variable_pump_pattern[pump_name][0]): # 给 list 中的所有元素除以 50Hz modified_values = [ value / 50 for value in modify_variable_pump_pattern[pump_name] ] pump_pattern = get_pattern( name_c, get_pump(name_c, pump_name)["pattern"] ) pump_pattern["factors"][ modify_index : modify_index + len(modified_values) ] = modified_values cs = ChangeSet() cs.append(pump_pattern) set_pattern(name_c, cs) # 显式阀门控制优先于旧的开度参数。 if valve_control is not None: _apply_valve_control(name_c, valve_control) # 保留原开度参数逻辑，兼容现有方案调用。 elif modify_valve_opening: for valve_name in modify_valve_opening.keys(): if not np.isnan(modify_valve_opening[valve_name]): valve_status = get_status(name_c, valve_name) if modify_valve_opening[valve_name] == 0: valve_status["status"] = "CLOSED" valve_status["setting"] = 0 elif modify_valve_opening[valve_name] < 1: valve_status["status"] = "OPEN" valve_status["setting"] = 0.1036 * pow( modify_valve_opening[valve_name], -3.105 ) elif modify_valve_opening[valve_name] == 1: valve_status["status"] = "OPEN" valve_status["setting"] = 0 cs = ChangeSet() cs.append(valve_status) set_status(name_c, cs) # 运行并返回结果 result_data = json.loads(run_project(name_c)) refresh_materialized_views_after_commit(name_c) time_cost_end = time.perf_counter() print( "{} -- Hydraulic simulation finished, cost time: {:.2f} s.".format( datetime.now(pytz.timezone("Asia/Shanghai")).strftime("%Y-%m-%d %H:%M:%S"), time_cost_end - time_cost_start, ) ) output_data = result_data.get("output") if not isinstance(output_data, dict): raise RuntimeError("run_project did not return JSON output content") node_result = output_data.get("node_results") link_result = output_data.get("link_results") if node_result is None or link_result is None: raise RuntimeError("run_project output missing node_results or link_results") # link_flow = [] # for link in link_result: # link_flow.append(link['result'][-1]['flow']) # print(link_flow) times_info = output_data.get("times") or {} num_periods_result = times_info.get("num_periods") if num_periods_result is None: raise RuntimeError("run_project output missing times.num_periods") print("simulation_type", simulation_type) print("before store result") # print(num_periods_result) # print(node_result) # 存储 starttime = time.time() # 临时处理输入的name问题，后续需要优化，需要传递/获取iot数据库的名字 db_name = result_db_name or name if simulation_type.upper() == "REALTIME": TimescaleInternalStorage.store_realtime_simulation( node_result, link_result, modify_pattern_start_time, db_name=db_name ) elif simulation_type.upper() == "EXTENDED": result_timestep_seconds = times_info.get("report_step") if result_timestep_seconds is None: raise RuntimeError("run_project output missing times.report_step") if not scheme_type or not scheme_name: raise ValueError("extended simulation requires analysis run type and name") detail = dict(scheme_detail or {}) result_units = output_data.get("units") if isinstance(result_units, dict): detail["result_units"] = { key: str(result_units[key]).strip() for key in ("flow", "pressure") if result_units.get(key) } run_id = create_analysis_run( name=db_name, scheme_name=scheme_name, scheme_type=scheme_type, username=scheme_username, scheme_start_time=modify_pattern_start_time, scheme_detail=detail, ) try: TimescaleInternalStorage.store_analysis_simulation( run_id, node_result, link_result, modify_pattern_start_time, num_periods_result, result_timestep_seconds, db_name=db_name, ) except Exception: try: update_analysis_run( db_name, run_id, status="failed", username=scheme_username, scheme_detail=detail, ) except Exception: logging.exception("failed to mark analysis run %s as failed", run_id) raise update_analysis_run( db_name, run_id, status="completed", username=scheme_username, scheme_detail=detail, ) endtime = time.time() logging.info("store time: %f", endtime - starttime) # 暂不需要再次存储 SCADA 模拟信息 # TimescaleInternalQueries.fill_scheme_simulation_result_to_SCADA(scheme_type=scheme_type, scheme_name=scheme_name) print("after store result") return run_id if simulation_type.upper() == "EXTENDED" else None