一、【时序预测】之水质净化厂工艺控制-曝气量预测

竞赛地址： https://www.datafountain.cn/competitions/602

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.北生化池和南生化池在生产过程中不会互相影响。

二、特征提取

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 

def run_lgb(df_train, df_test, ycol, n_splits=5, seed=2022):
    use_feats = [col for col in df_test.columns if col not in ['time','Label1','Label2','label']]
    model = lgb.LGBMRegressor(num_leaves=32,objective='mape',
                              max_depth=16,
                              learning_rate=0.1,
                              n_estimators=10000,
                              subsample=0.8,
                              feature_fraction=0.6,
                              reg_alpha=0.5,
                              reg_lambda=0.25,
                              random_state=seed,
                              metric=None)
    oof = []
    prediction = df_test[['time']]
    prediction[ycol] = 0
    df_importance_list = []
    from tscv import GapKFold
    cv = GapKFold(n_splits=n_splits, gap_before=0, gap_after=0)
    for fold_id, (trn_idx, val_idx) in enumerate(cv.split(df_train[use_feats])):
        X_train = df_train.iloc[trn_idx][use_feats]
        Y_train = df_train.iloc[trn_idx][ycol]        
        X_val = df_train.iloc[val_idx][use_feats]
        Y_val = df_train.iloc[val_idx][ycol]
        lgb_model = model.fit(X_train,
                              Y_train,
                              eval_names=['train', 'valid'],
                              eval_set=[(X_train, Y_train), (X_val, Y_val)],
                              verbose=100,
                              eval_metric='rmse',
                              early_stopping_rounds=100)
        pred_val = lgb_model.predict(X_val, num_iteration=lgb_model.best_iteration_)
        df_oof = df_train.iloc[val_idx][['time', ycol]].copy()
        df_oof['pred'] = pred_val
        oof.append(df_oof)
        pred_test = lgb_model.predict(df_test[use_feats], num_iteration=lgb_model.best_iteration_)
        prediction[ycol] += pred_test / n_splits
        df_importance = pd.DataFrame({
            'column': use_feats,
            'importance': lgb_model.feature_importances_,
        })
        df_importance_list.append(df_importance)
        del lgb_model, pred_val, pred_test, X_train, Y_train, X_val, Y_val
        gc.collect()
    df_importance = pd.concat(df_importance_list)
    df_importance = df_importance.groupby(['column'])['importance'].agg(
        'mean').sort_values(ascending=False).reset_index()
    display(df_importance.head(50))
    df_oof = pd.concat(oof).reset_index(drop=True)
    df_oof[ycol] = np.expm1(df_oof[ycol])
    df_oof['pred'] = np.expm1(df_oof['pred'])
    prediction[ycol] = np.expm1(prediction[ycol])
    
    return df_oof, prediction

df_oof_B, pred_B = run_lgb(train_B, test_B, ycol='Label1',n_splits=10)
df_oof_N, pred_N = run_lgb(train_N, test_N, ycol='Label2',n_splits=10)

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

四、预测提交

def calc_score(df1, df2):
    rmse_1 = np.sqrt(mean_squared_error(df1['pred'], (df1['Label1'])))
    rmse_2 = np.sqrt(mean_squared_error(df2['pred'], (df2['Label2'])))
    loss = (rmse_1+rmse_2)/2
    print(rmse_1,rmse_2)
    score = (1 / (1 + loss)) * 1000
    return score

calc_score(df_oof_B, df_oof_N)

3091.5013527627148 2248.255071349608





0.37440868531034793

# 提交
sub = pd.read_csv('data/data169443/sample_submission.csv')
sub['Label1'] = pred_B['Label1'].values
sub['Label2'] = pred_N['Label2'].values
sub

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

sub.to_csv('result.csv', index=False)

下载提交即可。