Sfoglia il codice sorgente

feat:1. 更新了带异常接口的反渗透脚本,当出现异常时,脚本返回调用脚本的时间及异常变量名 2. 更新了反渗透脚本的相关说明

junc_WHU 21 ore fa
parent
commit
f0c5ad1251

+ 139 - 0
models/ro_mechanism_predict/cip_predict/README.md

@@ -0,0 +1,139 @@
+下面是**精简版 README(交接用)**,重点只保留“项目是什么 + 怎么用”。
+
+---
+
+# CIP Prediction Model(超滤CIP预测模型)
+
+## 1. 项目简介
+
+本项目用于超滤系统运行过程中的膜污染建模与CIP(化学清洗)时间预测。
+
+核心功能包括:
+
+* 基于运行数据计算膜阻力增长(R)
+* 建立污染增长模型(R = a · V^b)
+* 预测剩余运行时间(remaining_days)
+* 支持不同压差通道(DPT_1 / DPT_2)
+* 引入历史周期数据进行参数修正
+* 输出异常状态(ERROR)用于运行监控
+
+---
+
+## 2. 输入数据
+
+---
+
+### 2.1 实时运行数据
+
+文件:
+
+```
+use_data/test_cycle_sample.csv
+```
+
+包含:
+
+* time(时间,5分钟采样)
+* RO1_FluxF ~ RO4_FluxF(各RO段通量)
+* C.M.RO1_FT_JS@out ~ C.M.RO4_FT_JS@out(进水流量)
+* C.M.RO1_FT_NS@out ~ C.M.RO4_FT_NS@out(浓水流量)
+* C.M.RO_TT_ZJS@out(温度)
+* C.M.RO_Cond_ZJS@out(电导率)
+* C.M.RO_PH_ZJS@out(pH)
+* C.M.RO_ORP_ZJS@out(ORP)
+* C.M.RO_PT_ZCS@out(进水压力)
+* C.M.RO1_DB@DPT_1 ~ C.M.RO4_DB@DPT_1(一段压差)
+* C.M.RO1_DB@DPT_2 ~ C.M.RO4_DB@DPT_2(二段压差)
+
+说明:
+
+* 数据为当前CIP周期的连续运行记录
+* 必须为进水运行阶段数据(控制字=24.0)
+* 不允许混入清洗、停机或切换状态数据
+* 可输入全机组数据,模型内部会自动筛选目标机组数据
+
+---
+
+### 2.2 历史周期数据
+
+```
+use_data/cycle_rv_seg1_results.csv
+use_data/cycle_rv_seg2_results.csv
+```
+
+用于提供历史污染参数(a, b)进行修正。
+
+说明:
+
+* seg1 对应 DPT_1 压差模型
+* seg2 对应 DPT_2 压差模型
+* 每条记录代表一个完整历史CIP周期
+* 包含污染模型参数(a, b)及对应运行工况统计特征(flux、cond、pH、ORP、temp)
+* 数据已完成预处理,可直接用于模型匹配,无需额外清洗或转换
+
+---
+
+## 3. 使用方法
+
+### 3.1 基本调用
+
+```python
+out = main(["RO4"], TMP_limit=0.19, dpt_type="DPT_1")
+```
+
+---
+
+### 3.2 参数说明
+
+* `units`:机组列表(如 ["RO4"])
+* `TMP_limit`:最大允许压差(用于CIP判断)
+* `dpt_type`:
+
+  * "DPT_1"(默认)
+  * "DPT_2"
+
+---
+
+## 4. 输出结果
+
+### 正常情况
+
+```text
+remaining_days:预计距离CIP剩余天数
+R_now:当前膜阻力
+V_now:当前产水量
+a_final / b_final:模型参数
+status = OK
+```
+
+---
+
+### 异常情况
+
+当系统判断预测失效时返回:
+
+```text
+status = ERROR
+error_time:异常发生时间
+error_var:使用的压差变量
+```
+
+---
+
+## 5. 运行示例
+
+```python
+if __name__ == "__main__":
+
+    out = main(["RO4"], TMP_limit=0.19, dpt_type="DPT_2")
+
+    row = out.iloc[0]
+
+    if row["status"] == "ERROR":
+        print("系统异常")
+        print(row["error_time"])
+        print(row["error_var"])
+    else:
+        print(f"预计CIP时间:{row['remaining_days']:.2f} 天")
+```
+

+ 229 - 51
models/ro_mechanism_predict/cip_predict/run_this.py

@@ -1,18 +1,28 @@
 import pandas as pd
 import numpy as np
 import os
+from datetime import datetime
 from sklearn.linear_model import LinearRegression
 
+# ============================================================
+# 路径配置
 # ============================================================
 BASE_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
 
 FEATURE_PATH = f"{BASE_DIR}/use_data/test_cycle_sample.csv"
+HIST_PATHS = {
+    "DPT_1": f"{BASE_DIR}/use_data/cycle_rv_seg1_results.csv",
+    "DPT_2": f"{BASE_DIR}/use_data/cycle_rv_seg2_results.csv"
+}
 
 
 # ============================================================
-# viscosity
+# 粘度模型(温度修正)
 # ============================================================
 def xishan_viscosity(T):
+    """
+    根据温度计算水的动力粘度(Pa·s)
+    """
     x = (T + 273.15) / 300
     factor = 890 / (
         280.68 * x ** -1.9 +
@@ -24,59 +34,65 @@ def xishan_viscosity(T):
 
 
 # ============================================================
-# unit mapping
+# 字段映射(适配不同机组)
 # ============================================================
-def get_unit_cols(unit):
+def get_unit_cols(unit, dpt_type="DPT_1"):
+
+    dpt_map = {
+        "DPT_1": "C.M.{}_DB@DPT_1",
+        "DPT_2": "C.M.{}_DB@DPT_2"
+    }
+
     return {
         "flux": f"{unit}_FluxF",
-        "dpt1": f"C.M.{unit}_DB@DPT_1",
-        "dpt2": f"C.M.{unit}_DB@DPT_2"
+        "dpt": dpt_map[dpt_type].format(unit)
     }
 
 
 # ============================================================
-# resistance
+# 膜阻力计算
 # ============================================================
-def compute_resistance(df, unit):
+def compute_resistance(df, unit, dpt_type="DPT_1"):
 
-    cols = get_unit_cols(unit)
+    cols = get_unit_cols(unit, dpt_type)
 
     flux = df[cols["flux"]]
-    dpt1 = df[cols["dpt1"]]
-    dpt2 = df[cols["dpt2"]]
+    dpt = df[cols["dpt"]]
 
     mu = xishan_viscosity(df["C.M.RO_TT_ZJS@out"])
 
-    df["R1"] = dpt1 * 3.6e12 / (mu * flux)
+    df["R1"] = dpt * 3.6e12 / (mu * flux)
+
     return df
 
 
 # ============================================================
-# V accumulation
+# 累积产水量计算
 # ============================================================
-def compute_V(df, unit):
+def compute_V(df, unit, dpt_type="DPT_1"):
 
-    cols = get_unit_cols(unit)
+    cols = get_unit_cols(unit, dpt_type)
     flux = df[cols["flux"]].values
 
     dt = 300
     J = flux * (1e-3 / 3600)
 
     V = np.zeros(len(df))
+
     for i in range(1, len(df)):
-        V[i] = V[i-1] + J[i-1] * dt
+        V[i] = V[i - 1] + J[i - 1] * dt
 
     df["V"] = V
     return df
 
 
 # ============================================================
-# fit model (POWER LAW)
-# R = a * V^b
-# log R = log a + b log V
+# 当前周期拟合(幂律模型)
 # ============================================================
 def fit_model(df, R_col):
-
+    """
+    拟合 R = a * V^b
+    """
     df = df.dropna(subset=[R_col, "V"])
     df = df[(df[R_col] > 0) & (df["V"] > 0)]
 
@@ -96,70 +112,207 @@ def fit_model(df, R_col):
 
 
 # ============================================================
-# solve V limit (analytical)
+# 读取历史周期库
 # ============================================================
-def solve_v_limit(model, R_limit):
+def load_history(dpt_type="DPT_1"):
+    """
+    根据压差类型加载对应历史库
+    """
+    return pd.read_csv(HIST_PATHS[dpt_type])
+
+
+# ============================================================
+# 当前周期特征提取
+# ============================================================
+def extract_features(df, unit, dpt_type="DPT_1"):
+
+    cols = get_unit_cols(unit, dpt_type)
+
+    return {
+        "flux_mean": df[cols["flux"]].mean(),
+        "cond_mean": df["C.M.RO_Cond_ZJS@out"].mean(),
+        "ph_mean": df["C.M.RO_PH_ZJS@out"].mean(),
+        "orp_mean": df["C.M.RO_ORP_ZJS@out"].mean(),
+        "temp_mean": df["C.M.RO_TT_ZJS@out"].mean()
+    }
 
-    if model is None:
-        return None
 
-    b = model.coef_[0]
-    log_a = model.intercept_
+# ============================================================
+# 历史周期匹配(KNN + 加权)
+# ============================================================
+def match_history(hist_df, feat, unit, k=5):
+    """
+    在历史周期中寻找最相似的 k 个周期
+    并加权得到 a_hist, b_hist
+    """
+    df = hist_df[hist_df["unit"] == unit].copy()
+
+    feature_cols = ["flux_mean", "cond_mean", "ph_mean", "orp_mean", "temp_mean"]
+
+    X_hist = df[feature_cols].values
+    x_now = np.array([feat[c] for c in feature_cols]).reshape(1, -1)
+
+    # 标准化
+    mean = X_hist.mean(axis=0)
+    std = X_hist.std(axis=0) + 1e-9
+
+    Xn = (X_hist - mean) / std
+    xn = (x_now - mean) / std
+
+    # 欧氏距离
+    dist = np.linalg.norm(Xn - xn, axis=1)
+
+    df["dist"] = dist
+    df = df.sort_values("dist").head(k)
+
+    # 权重
+    w = 1 / (df["dist"].values + 1e-6)
+    w = w / w.sum()
 
-    target = np.log(R_limit / 1e12)
+    a_hist = np.sum(w * df["a"].values)
+    b_hist = np.sum(w * df["b"].values)
 
+    return a_hist, b_hist
+
+
+# ============================================================
+# 根据 a, b 解析求解 V_limit
+# ============================================================
+def solve_v_limit_ab(a, b, R_limit):
+    """
+    解 R_limit = a * V^b
+    """
     if abs(b) < 1e-12:
         return np.inf
 
-    logV = (target - log_a) / b
-    V_limit = np.exp(logV)
-
+    V_limit = (R_limit / 1e12 / a) ** (1 / b)
     return V_limit
 
 
 # ============================================================
-# main per unit
+# 单机组处理(核心函数)
 # ============================================================
-def process_unit(unit, df_new, TMP_limit=0.23):
+def process_unit(unit, df_new, TMP_limit=0.23, dpt_type="DPT_1"):
 
     df = df_new.copy()
 
-    df = compute_resistance(df, unit)
-    df = compute_V(df, unit)
+    # =========================
+    # 计算 R 和 V
+    # =========================
+    df = compute_resistance(df, unit, dpt_type)
+    df = compute_V(df, unit, dpt_type)
+
+    # =========================
+    # 当前周期时长(天)
+    # =========================
+    days = max(
+        (df["time"].max() - df["time"].min()).total_seconds() / 86400,
+        1e-6
+    )
 
     # =========================
     # 当前状态
     # =========================
     mu_avg = xishan_viscosity(df["C.M.RO_TT_ZJS@out"]).mean()
-
-    cols = get_unit_cols(unit)
+    cols = get_unit_cols(unit, dpt_type)
     flux_avg = df[cols["flux"]].mean()
 
     V_now = df["V"].iloc[-1]
     R_now = df["R1"].iloc[-1]
 
-    days = max((df["time"].max() - df["time"].min()).total_seconds() / 86400, 1e-6)
     V_daily = V_now / days
 
     # =========================
-    # model
+    # 当前周期拟合
     # =========================
     model = fit_model(df, "R1")
-
     if model is None:
         return None
 
+    a_cur = np.exp(model.intercept_)
+    b_cur = model.coef_[0]
+
+    # =========================
+    # 历史参数(仅用于修正)
+    # =========================
+    hist_df = load_history(dpt_type)
+    feat = extract_features(df, unit, dpt_type)
+    a_hist, b_hist= match_history(hist_df, feat, unit)
+
+    # =========================
+    # 时间权重(只依赖周期长度)
+    # 60天后完全关闭历史
+    # =========================
+    w = max(0.0, 1.0 - days / 60.0)
+    w = w ** 2   # 加速衰减(关键)
+
+    # =========================
+    # 数值保护
+    # =========================
+    a_cur_safe = max(a_cur, 1e-12)
+    b_cur_safe = np.sign(b_cur) * max(abs(b_cur), 1e-6)
+
+    # =========================
+    # a 修正(对数空间,更稳定)
+    # 最大允许 ±20%
+    # =========================
+    log_a_cur = np.log(a_cur_safe)
+    log_a_hist = np.log(max(a_hist, 1e-12))
+
+    delta_log_a = w * (log_a_hist - log_a_cur)
+    delta_log_a = np.clip(delta_log_a, -0.2, 0.2)
+
+    a_final = np.exp(log_a_cur + delta_log_a)
+
     # =========================
-    # R limit from TMP
+    # b 修正(比例修正)
+    # 最大允许 ±15%(更严格)
+    # =========================
+    delta_b = w * (b_hist - b_cur) / (abs(b_cur) + 1e-6)
+    delta_b = np.clip(delta_b, -0.15, 0.15)
+
+    b_final = b_cur * (1 + delta_b)
+
+    # =========================
+    # 物理约束(防止异常)
+    # =========================
+    b_final = np.clip(b_final, 0.01, 3.0)
+
+    # =========================
+    # TMP → R_limit
     # =========================
     R_limit = TMP_limit * 3.6e12 / (mu_avg * flux_avg)
 
     # =========================
-    # solve
+    # 求解 V_limit
     # =========================
-    V_limit = solve_v_limit(model, R_limit)
+    V_limit = solve_v_limit_ab(a_final, b_final, R_limit)
+
+    # 数值安全限制
+    if (not np.isfinite(V_limit)) or (V_limit <= 0):
+        V_limit = V_now
 
+    V_limit = max(V_limit, V_now)
+
+    # 剩余时间计算
     remaining_days = (V_limit - V_now) / max(V_daily, 1e-9)
+
+    # 异常检测
+    if remaining_days <= 0 or not np.isfinite(remaining_days):
+        cols = get_unit_cols(unit, dpt_type)
+
+        return {
+            "unit": unit,
+            "status": "ERROR",
+            "error_time": datetime.now().strftime("%Y-%m-%d %H:%M:%S"),
+            "error_var": cols["dpt"],  # 当前使用的压差变量
+            "V_now": V_now,
+            "V_limit": V_limit,
+            "a_final": a_final,
+            "b_final": b_final
+        }
+
+    # 正常情况
     remaining_days = min(remaining_days, 90 - days)
 
     return {
@@ -168,14 +321,15 @@ def process_unit(unit, df_new, TMP_limit=0.23):
         "V_now": V_now,
         "R_limit": R_limit,
         "V_limit": V_limit,
-        "remaining_days": remaining_days
+        "remaining_days": remaining_days,
+        "a_final": a_final,
+        "b_final": b_final,
     }
 
-
 # ============================================================
-# main
+# 主函数
 # ============================================================
-def main(units, TMP_limit=0.23):
+def main(units, TMP_limit=0.21, dpt_type="DPT_1"):
 
     df_new = pd.read_csv(FEATURE_PATH)
     df_new["time"] = pd.to_datetime(df_new["time"])
@@ -183,20 +337,44 @@ def main(units, TMP_limit=0.23):
     results = []
 
     for u in units:
-        res = process_unit(u, df_new, TMP_limit)
-        if res is not None:
+        res = process_unit(u, df_new, TMP_limit, dpt_type)
+
+        if res is None:
+            results.append({
+                "unit": u,
+                "status": "ERROR",
+                "error_time": None,
+                "error_var": None,
+                "error_type": "MODEL_FIT_FAILED"
+            })
+            continue
+
+        if "status" in res and res["status"] == "ERROR":
+            results.append(res)
+        else:
+            res["status"] = "OK"
             results.append(res)
 
     return pd.DataFrame(results)
 
-
 # ============================================================
-# run
+# 运行入口
 # ============================================================
 if __name__ == "__main__":
 
     units_to_run = ["RO4"]
 
-    out = main(units_to_run, TMP_limit=0.20)
+    out = main(["RO4"], TMP_limit=0.19, dpt_type="DPT_2")
+
+    # 防止空结果
+    if out is None or len(out) == 0:
+        print("无有效结果")
+    else:
+        row = out.iloc[0]
 
-    print(f"该机组预计下次CIP为:{out['remaining_days'].iloc[0]:.2f} 天后")
+        if row["status"] == "ERROR":
+            print("系统异常:")
+            print(f"时间: {row['error_time']}")
+            print(f"变量: {row['error_var']}")
+        else:
+            print(f"该机组预计下次CIP为:{row['remaining_days']:.2f} 天后")