一、【时序预测】之水质净化厂工艺控制-曝气量预测
竞赛地址： www.datafountain.cn/competition…
1.大赛背景
为推进智慧水务建设，激发数字化创新能力，助力创新应用挖掘与落地，加快水务行业现代化进程，深圳市环境水务集团有限公司发起主办首届“深水云脑杯”智慧水务数据创新大赛。深水渠成，群智创新！本届大赛将以数字化创新模式为抓手，把握发展脉搏、汇聚产学研力量，成为汇聚大数据、人工智能等数据智能相关信息技术在水务领域应用创新的擂台，揭榜挂帅，推动赛事成果转化，解决社会和行业聚焦的难题，促进水务行业实现智慧升级及高质量发展。
2. 赛题背景
当今社会，水资源匮乏的问题越来越突出，污水处理是解决水资源匮乏的有效手段之一。膜生物反应器（MBR）工艺作为近年来的一种新型污水工艺，较传统的活性污泥法来说，具有占地面积小，产水水质高、剩余污泥少、自控程度高等优势，在用地资源日益紧张的今天，MBR工艺在全国各地的污水处理厂均得到了一定的应用。但同时，由于其基础造价较高、膜污染及能耗较高等问题，其进一步应用也得到了了一定的限制。以采用厌氧-缺氧-好氧生物脱氮除磷（A-A-O）耦合膜生物反应器工艺（A2O-MBR）的某污水厂为例，通过对各工艺段能耗进行对比分析，得出其生化段曝气能耗占全厂能耗的比例高达49%。从能量转换的角度来看，其实质是以能耗换取水质。
本次赛题主要是通过对采用A2O-MBR工艺的某污水厂运行过程中的历史数据，利用大数据建模，形成可供推理的智能曝气数理模型，通过算法迭代计算出最优曝气量。结合污水处理厂工艺流程，通过数学建模，建立污水厂精准曝气机理模型，实现生化处理系统运行效果的优化控制，在保障污水厂出水水质满足行标的前提下，采用智能化、自动化手段降低能耗，有效解决实际问题，助力市政污水处理行业低碳发展。
3.赛题任务
水质净化厂运营过程中，曝气量需根据进出水水质等参数实时进行调节，以保障出水水质达标，而在实际生产过程中，由于影响因素多，目前尚不能对曝气量进行精确控制，希望通过机器学习模型的应用预测曝气量，以指导实际生产。
4.数据简介
数据来源于使用A2O-MBR工艺的污水厂，生化池好氧段曝气生产环境。在生化池的好氧工艺段，我们通过控制曝气量，来改变水中的溶解氧（DO），将污水中的氨氮（NH4）转化成硝氮（NO3），同时将水中的有机物（COD）降解。
5.数据说明
- 1.部分特征中的0值并不代表实际值。
- 2.特征在邻近时序内数值未发生改变，通常是因为数据采集频率原因，并非生产环境中未发生显著改变。
- 3.北生化池和南生化池在生产过程中不会互相影响。
| time | Label1 | Label2 |
| 2022/7/18 2:40 | 0 | 0 |
| 2022/7/18 2:42 | 0 | 0 |
| 2022/7/18 2:44 | 0 | 0 |
| 2022/7/18 2:46 | 0 | 0 |
| 2022/7/18 2:48 | 0 | 0 |
| 2022/7/18 2:50 | 0 | 0 |
| … | 0 | 0 |
二、特征提取
1.lightgbm升级
lightgbm aistudio默认为3.1.1，许多方法不支持，建议升级最新版3.3.2.
!pip list | grep lightgbm
lightgbm 3.1.1
!pip install -U -q lightgbm !pip list | grep lightgbm
lightgbm 3.3.2
2.导入常用库
2.1 gc库简介
gc模块即Python中垃圾回收模块，它提供可选的垃圾回收器的接口。同时提供对回收器找到但是无法释放的不可达对象的访问。由于 Python 使用了带有引用计数的回收器，如果你确定你的程序不会产生循环引用，你可以关闭回收器。可以通过调用 gc.disable() 关闭自动垃圾回收。
- enable() --启用自动垃圾回收。
- disable() --禁用自动垃圾回收。
- isenabled() --如果启用了自动收集，则返回true。
- collect() --立即执行完全收集。
- get_count() --返回当前集合计数。
- get_stats() --返回包含每代统计信息的词典列表。
- set_debug() --设置调试标志。
- get_debug() --获取调试标志。
- set_threshold() --设置收集阈值。
- get_threshold() --返回集合阈值的当前值。
- get_objects() --返回收集器跟踪的所有对象的列表。
- is_tracked() --如果跟踪给定对象，则返回true。
- is_finalized() --如果给定对象已定稿，则返回true。
- get_referrers() --返回引用对象的对象列表。
- get_referents() --返回对象引用的对象列表。
- freeze() --冻结所有跟踪对象，并在将来的收集中忽略它们。
- unfreeze() --解冻永久生成中的所有对象。
- get_freeze_count() --返回永久生成中的对象数。
最常用的方法：gc.collect() --立即执行完全收集，释放出不使用的资源，归还内存。可以通过参数generation，单独对0,1,2代进行回收释放。
例如：编写无限循环C程序：
#include<stdio.h> int main(){ while(1){ printf("okok"); } return 0; }
编译备用
cc test.c -o test
!cc test.c -o test
import subprocess, psutil, gc mem1 = psutil.virtual_memory() print(f"某程序前内存已使用：{mem1.used}") print(f"某程序前内存剩余：{mem1.free}") print(f"某程序前内存百分比：{mem1.percent}") app1 = subprocess.Popen('test') app2 = subprocess.Popen(r'test') app3 = subprocess.Popen(r'test') mem2 = psutil.virtual_memory() print(f"某程序后内存已使用：{mem2.used}") print(f"某程序后内存剩余：{mem2.free}") print(f"某程序后内存百分比：{mem2.percent}") app1.kill() app2.kill() app3.kill() gc.collect() mem3 = psutil.virtual_memory() print(f"GC回收后内存已使用：{mem3.used}") print(f"GC回收后内存剩余：{mem3.free}") print(f"GC回收后内存百分比：{mem3.percent}")
某程序前内存已使用：6770204672 某程序前内存剩余：68427513856 某程序前内存百分比：6.5 某程序后内存已使用：6771159040 某程序后内存剩余：68426559488 某程序后内存百分比：6.5 GC回收后内存已使用：6769913856 GC回收后内存剩余：68427763712 GC回收后内存百分比：6.5
需要注意的是：执行收集本身也需要一点的内存代价，所以可能存在收集完成后内存反而增加的情况。
2.2 tqdm介绍
2.2.1 不带参数
# 不带参数 from tqdm import tqdm import time for i in tqdm(range(50)): time.sleep(0.1)
100%|██████████| 50/50 [00:05<00:00, 9.92it/s]
2.2.2 带参数
# 带参数 from tqdm import tqdm import time d = {'loss':0.2,'learn':0.8} for i in tqdm(range(50),desc='进行中',ncols=100,postfix=d): #desc设置名称,ncols设置进度条长度.postfix以字典形式传入详细信息 time.sleep(0.1)
进行中: 100%|█████████████████████████████████████| 50/50 [00:05<00:00, 9.89it/s, learn=0.8, loss=0.2]
2.2.3 处理列表
# 用tqdm处理列表中的对象，显示处理进度 from tqdm import tqdm import time bar = tqdm(['p1','p2','p3','p4','p5']) for b in bar: time.sleep(0.5) bar.set_description("处理{0}中".format(b))
处理p5中: 100%|██████████| 5/5 [00:02<00:00, 1.99it/s]
import warnings warnings.simplefilter('ignore') import gc import numpy as np import pandas as pd pd.set_option('display.max_columns', None) from tqdm.auto import tqdm from sklearn.model_selection import KFold from sklearn.metrics import mean_squared_error import lightgbm as lgb
2.3 rolling滑动窗口
为了提升数据的准确性，将某个点的取值扩大到包含这个点的一段区间，用区间来进行判断，这个区间就是窗口。移动窗口就是窗口向一端滑行，默认是从右往左，每次滑行并不是区间整块的滑行，而是一个单位一个单位的滑行。
例如：
import pandas as pd s = [1,2,3,5,6,10,12,14,12,30] pd.Series(s).rolling(window=3).mean()
0 NaN 1 NaN 2 2.000000 3 3.333333 4 4.666667 5 7.000000 6 9.333333 7 12.000000 8 12.666667 9 18.666667 dtype: float64
设置的窗口window=3，也就是3个数取一个均值。index 0,1 为NaN，是因为它们前面都不够3个数，等到index2 的时候，它的值是怎么算的呢，就是（index0+index1+index2 ）/3 index3 的值就是（ index1+index2+index3）/ 3
DataFrame.rolling(window, min_periods=None, center=False, win_type=None, on=None, axis=0, closed=None)
- window： 也可以省略不写。表示时间窗的大小，注意有两种形式（int or offset）。如果使用int，则数值表示计算统计量的观测值的数量即向前几个数据。如果是offset类型，表示时间窗的大小。offset详解
- min_periods：每个窗口最少包含的观测值数量，小于这个值的窗口结果为NA。值可以是int，默认None。offset情况下，默认为1。
- center: 把窗口的标签设置为居中。布尔型，默认False，居右
- win_type: 窗口的类型。截取窗的各种函数。字符串类型，默认为None。各种类型
- on: 可选参数。对于dataframe而言，指定要计算滚动窗口的列。值为列名。
- axis: int、字符串，默认为0，即对列进行计算
- closed：定义区间的开闭，支持int类型的window。对于offset类型默认是左开右闭的即默认为right。可以根据情况指定为left both等。
3.特征提取
train = pd.read_csv('data/data169443/train_dataset.csv')
test = pd.read_csv('data/data169443/evaluation_public.csv')
# 合并数据集 df = pd.concat([train,test])
roll_cols = ['JS_NH3', 'CS_NH3', 'JS_TN', 'CS_TN', 'JS_LL', 'CS_LL', 'MCCS_NH4', 'MCCS_NO3', 'JS_COD', 'CS_COD', 'JS_SW', 'CS_SW', 'B_HYC_NH4', 'B_HYC_XD', 'B_HYC_MLSS', 'B_HYC_JS_DO', 'B_HYC_DO', 'B_CS_MQ_SSLL', 'B_QY_ORP', 'N_HYC_NH4', 'N_HYC_XD', 'N_HYC_MLSS', 'N_HYC_JS_DO', 'N_HYC_DO', 'N_CS_MQ_SSLL', 'N_QY_ORP'] df['time'] = pd.to_datetime(df['time']) for i in range(1,5): df[[ii+f'_roll_{i}_mean_diff' for ii in roll_cols]] = df[roll_cols].rolling(i, min_periods=1).sum().diff() df[[ii+'_roll_8_mean' for ii in roll_cols]] = df[roll_cols].rolling(8, min_periods=1).mean() df[[ii+'_roll_16_mean' for ii in roll_cols]] = df[roll_cols].rolling(16, min_periods=1).mean() df[[ii+'_roll_16_mean_diff' for ii in roll_cols]] = df[[ii+'_roll_16_mean' for ii in roll_cols]].diff() df[[ii+'_roll_8_mean_diff' for ii in roll_cols]] = df[[ii+'_roll_8_mean' for ii in roll_cols]].diff() df[[ii+'_roll_8_std' for ii in roll_cols]] = df[roll_cols].rolling(8, min_periods=1).std()
train = df.iloc[:train.shape[0]] test = df.iloc[train.shape[0]:]
N_col = ['N_HYC_NH4', 'N_HYC_XD', 'N_HYC_MLSS', 'N_HYC_JS_DO', 'N_HYC_DO', 'N_CS_MQ_SSLL', 'N_QY_ORP'] B_col = ['B_HYC_NH4', 'B_HYC_XD', 'B_HYC_MLSS', 'B_HYC_JS_DO', 'B_HYC_DO', 'B_CS_MQ_SSLL', 'B_QY_ORP'] NB_col = ['A_'+ ii[2:] for ii in ['B_HYC_NH4', 'B_HYC_XD', 'B_HYC_MLSS', 'B_HYC_JS_DO', 'B_HYC_DO', 'B_CS_MQ_SSLL', 'B_QY_ORP']] train[NB_col] = train[B_col].values/(train[N_col].values+ 1e-3) test[NB_col] = test[B_col].values/(test[N_col].values+ 1e-3) # NB_col
# 1. 数据说明里表示，北生化池和南生化池在生产过程中不会互相影响, 可以先试下分开两部分 # 2. 只用有 label 的数据 train_B = train[[i for i in train.columns if (i != 'Label2' and not i.startswith('N_'))]].copy() train_N = train[[i for i in train.columns if (i != 'Label1' and not i.startswith('B_'))]].copy() train_B = train_B[train_B['Label1'].notna()].copy().reset_index(drop=True) train_N = train_N[train_N['Label2'].notna()].copy().reset_index(drop=True) test_B = test[[i for i in test.columns if not i.startswith('N_')]].copy() test_N = test[[i for i in test.columns if not i.startswith('B_')]].copy()
# 时间特征 def add_datetime_feats(df): df['time'] = pd.to_datetime(df['time']) df['day'] = df['time'].dt.day df['hour'] = df['time'].dt.hour df['dayofweek'] = df['time'].dt.dayofweek return df train_B = add_datetime_feats(train_B) train_N = add_datetime_feats(train_N) test_B = add_datetime_feats(test_B) test_N = add_datetime_feats(test_N)
# 做点比率数值特征 def add_ratio_feats(df, type_='B'): df['JS_CS_NH3_ratio'] = df['JS_NH3'] / (df['CS_NH3'] + 1e-3) df['JS_CS_TN_ratio'] = df['JS_TN'] / (df['CS_TN'] + 1e-3) df['JS_CS_LL_ratio'] = df['JS_LL'] / (df['CS_LL'] + 1e-3) df['MCCS_NH4_NH3_ratio'] = df['MCCS_NH4'] / (df['CS_NH3'] + 1e-3) df['MCCS_NO3_NH3_ratio'] = df['MCCS_NO3'] / (df['CS_NH3'] + 1e-3) df['JS_CS_COD_ratio'] = df['JS_COD'] / (df['CS_COD'] + 1e-3) df['JS_CS_SW_ratio'] = df['JS_SW'] / (df['CS_SW'] + 1e-3) df['HYC_DO_ratio'] = df[f'{type_}_HYC_JS_DO'] / (df[f'{type_}_HYC_DO'] + 1e-3) df['CS_MQ_LL_ratio'] = df[f'{type_}_CS_MQ_SSLL'] / (df['CS_LL'] + 1e-3) return df train_B = add_ratio_feats(train_B, type_='B') train_N = add_ratio_feats(train_N, type_='N') test_B = add_ratio_feats(test_B, type_='B') test_N = add_ratio_feats(test_N, type_='N')
# target log1p 转换 B_max, B_min = train_B['Label1'].max(), train_B['Label1'].min() N_max, N_min = train_N['Label2'].max(), train_N['Label2'].min() train_B['Label1'] = np.log1p(train_B['Label1']) train_N['Label2'] = np.log1p(train_N['Label2'])
三、模型训练
需安装更新 tscv
----> 1 df_oof_B, pred_B = run_lgb(train_B, test_B, ycol='Label1',n_splits=10) 2 df_oof_N, pred_N = run_lgb(train_N, test_N, ycol='Label2',n_splits=10) /tmp/ipykernel_110/1103088436.py in run_lgb(df_train, df_test, ycol, n_splits, seed) 15 prediction[ycol] = 0 16 df_importance_list = [] ---> 17 from tscv import GapKFold 18 cv = GapKFold(n_splits=n_splits, gap_before=0, gap_after=0) 19 for fold_id, (trn_idx, val_idx) in enumerate(cv.split(df_train[use_feats])): ModuleNotFoundError: No module named 'tscv'
!pip install -U tscv
Looking in indexes: https://pypi.tuna.tsinghua.edu.cn/simple Collecting tscv Downloading https://pypi.tuna.tsinghua.edu.cn/packages/b8/5f/dfdbec6c4441f484e15d1f89ff6f3cbb33009e67a282fa7f6d31d16de13a/tscv-0.1.2-py3-none-any.whl (18 kB) Requirement already satisfied: scikit-learn>=0.22 in /opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages (from tscv) (0.24.2) Requirement already satisfied: numpy>=1.13.3 in /opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages (from tscv) (1.19.5) Requirement already satisfied: scipy>=0.19.1 in /opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages (from scikit-learn>=0.22->tscv) (1.6.3) Requirement already satisfied: threadpoolctl>=2.0.0 in /opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages (from scikit-learn>=0.22->tscv) (2.1.0) Requirement already satisfied: joblib>=0.11 in /opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages (from scikit-learn>=0.22->tscv) (0.14.1) Installing collected packages: tscv Successfully installed tscv-0.1.2 [1m[[0m[34;49mnotice[0m[1;39;49m][0m[39;49m A new release of pip available: [0m[31;49m22.1.2[0m[39;49m -> [0m[32;49m22.2.2[0m [1m[[0m[34;49mnotice[0m[1;39;49m][0m[39;49m To update, run: [0m[32;49mpip install --upgrade pip[0m
