物联网实时数据异常率预警

对于回归任务，这里选择一个常用的回归评估指标：根均方误差（Root Mean Square Error, RMSE）。

RMSE 指标可以评价数据的变化程度，RMSE 越小，说明模型的拟合实验数据能力越强。

K 值选择过小，容易造成模型过拟合；K 值选择过大，容易造成模型欠拟合。

RMSE 越小，模型的拟合效果越好。本场景经过多次实验，最终实验的 K 值选择为 200。

模型预测场景流程设计如下：

1. 每 10 秒钟用前 10 秒采集到的 100 万条数据，训练出一个 KNN 回归模型。
2. 用前 10 秒的 100 万条数据，进行聚合计算，将 10 秒产生的 100 万条数据聚合成 1 秒产生的 10 万条数据。
3. 用 KNN 模型对聚合数据（10 万条）进行预测，预测未来第 10s 的数据（10 万条）。
4. 将未来第 10s 的数据（10 万条）与 10s 聚合的数据（10 万条）输入预警模块进行预警计算。

预警计算场景流程设计如下：

1. 未来第 10s 的数据（10 万条）与 10s 聚合的数据（10 万条）进行比对，值差异大于 20% （即 ABS（预测值 - 真实值）/真实值大于 0.2）列为异常数据。
2. 统计异常率（异常数据的比例）：异常数据数目/比对总数（10w） 。
3. 设定阈值 0.215，当异常率超过阈值 0.215 时，进行报警。

用户可以利用下方代码完成实时生产数据的模拟。可以通过调整变量 hour ,自定义生产多少小时数据。

由于数据为模拟数据，为了保证模型能够完成数据拟合，这里提前设定监测值与 5 个影响因素的关系为线性关系，以确保回归模型拟合的关系能够收敛。

另外，为了模拟模型生成的预测值与真实值存在一定差异，方便进行异常比对，因此在代码中，对线性参数进行了随机化，以降低模型的准确率。用户真实场景中可自行调整参数方差大小，以控制模型准确率。简言之，方差越小，模型拟合越高，模型精度越高。

orgColNames= `time`deviceCode`wind`humidity`air_pressure`temperature`life`propertyValue
orgColTypes = [TIMESTAMP,INT,DOUBLE,DOUBLE,DOUBLE,DOUBLE,INT,INT]
enableTableShareAndPersistence(table = streamTable(10000:0, orgColNames,orgColTypes), tableName=`dataTable, cacheSize = 6000000) 

def send_data(begintime,endtime,hour,devices_number,rate,mutable dest)
{

    btime=timestamp(begintime)
	do{
        seconds  = int(60*60*hour)  //定义需要压入的批次，每秒钟1批

        n = devices_number * rate* 1 // 每秒钟生产10万条数据
        time =sort(take(btime+(0..(1000-1)),n)) //
        deviceCode = take(1..devices_number,n)
        x1 = randNormal(25,2,n)  
        x2 = randNormal(55,5,n)
        x3 = randNormal(1.01325,0.00001,n)
        x4 = randNormal(75,5,n)
        x5 = int(randNormal(10,3,n))
        b1 = randNormal(0.4,0.05,n) //方差0.05 降低模型准确率
        b2 = randNormal(0.3,0.05,n)
        b3 = randNormal(0.2,0.05,n)
        b4 = randNormal(0.09,0.05,n)
        b5 = randNormal(0.01,0.001,n)
        bias = randNormal(5,1,n)
        propertyValue = int(b1*x1*10 + b2*x2*2 + b3*x3*1000 + b4*x4 + b5*x5 +bias)

        table_ps = table(time,deviceCode,x1,x2,x3,x4,x5,propertyValue)
        dest.append!(table_ps)

        btime=btime+1000
        etime=timestamp(now())
        timediff=btime-etime
        if(timediff>0){sleep(timediff)}
    }while(btime<=endtime)
}
hour = 0.5 //用户自定义压入多少小时的数据,1 为 1个小时
devices_number = 100 //设备数目
rate = 1000 //每台设备 每秒钟1000条数据
begintime = datetime(now()) //开始时间
endtime = begintime + int(hour * 3600-1)  //结束时间
submitJob(`send_data_to_kafka,`send_data_to_kafka,send_data,begintime,endtime,hour,devices_number,rate,dataTable) 

聚合计算通过 DolphinDB 内置的时间序列引擎 createTimeSeriesEngine 实现。订阅数据表 dataTable 中的数据，将数据内容输入聚合引擎 tradesAggregator 完成预警计算，最后将聚合计算结果存入聚合表 aggrTable。

enableTableShareAndPersistence(table = streamTable(100000:0,orgColNames,orgColTypes), tableName=`aggrTable, cacheSize = 5000000) 

tradesAggregator = createTimeSeriesEngine(name="streamAggr", windowSize=10, step=10, metrics=<[avg(wind),avg(humidity),avg(air_pressure),avg(temperature),avg(life),avg(propertyValue)]>, dummyTable=dataTable, outputTable=aggrTable, timeColumn=`time, useSystemTime=false, keyColumn=`deviceCode, garbageSize=1000000)

subscribeTable(tableName="dataTable", actionName="Aggregator", offset=0, handler=append!{tradesAggregator}, msgAsTable=true)

订阅数据表 dataTable 中的内容，将流表数据存入分布式数据库 dfs://Data 中。

/*
  * 将dataTable数据写入分布式表 

	数据建模：
	1）每小时记录数：360,000,000
	2）每条记录大小：46字节
	3）每小时空间占用（压缩前）：15.42G
	4）建议以“Id”值和“小时”进行组合分区，每分区≈157.93M
	5）分区索引为“时间戳”+“设备号”

  */
 def createConsumerDataBase(dbname,tbname,col_names,col_types){
    dbTime = database("",VALUE,datehour(2023.01.01T00:00:00)..datehour(2023.01.01T23:00:00))
    Ids = 1..100
    dbId = database("",VALUE,Ids)
    db = database(directory=dbname, partitionType=COMPO, partitionScheme=[dbTime,dbId],engine="TSDB")
    factor_partition = db.createPartitionedTable(table=table(1:0,col_names,col_types),tableName = tbname,partitionColumns = ["time","deviceCode"],sortColumns =["deviceCode","time"],compressMethods={time:"delta"},keepDuplicates=LAST)
}
dataTable_dbname,dataTable_tbname = "dfs://Data","data"
createConsumerDataBase(dataTable_dbname,dataTable_tbname,orgColNames,orgColTypes)
subscribeTable(tableName="dataTable", actionName="append_data_into_dfs", offset=0, handler=loadTable(dataTable_dbname,dataTable_tbname), msgAsTable=true,batchSize=100000, throttle=1, reconnect=true)

订阅聚合表 aggrTable 中的内容，将流表数据存入分布式数据库 dfs://Aggregate 中。

/*
  * 将聚合计算结果写入分布式表
	数据建模：
		1）每小时记录数：36,000,000
		2）每条记录大小：46字节
		3）每小时空间占用（压缩前）：1.54G
		4）建议以“id”和“天”进行值分区，每分区≈ 379.03M
		5）分区索引为“时间戳”+“设备号”

  */
def createAggregateDataBase(dbname,tbname,col_names,col_types){
	if(existsDatabase(dbname)){dropDatabase(dbname)}
	Ids = 1..100
    dbId = database("",VALUE,Ids)
    dbTime = database("",VALUE,date(2023.01.01T00:00:00)..date(2023.12.31T20:00:00))
	db = database(directory=dbname, partitionType=COMPO, partitionScheme=[dbTime,dbId],engine="TSDB")
	factor_partition = db.createPartitionedTable(table=table(1:0,col_names,col_types),tableName = tbname,partitionColumns =["time","deviceCode"],sortColumns =["deviceCode","time"],compressMethods={time:"delta"},keepDuplicates=LAST)
}
aggr_dbname,aggr_tbname = "dfs://Aggregate","data"
createAggregateDataBase(aggr_dbname,aggr_tbname,orgColNames,orgColTypes)
subscribeTable(tableName="aggrTable", actionName="append_Aggregator_into_dfs", offset=0, handler=loadTable(aggr_dbname,aggr_tbname), msgAsTable=true,batchSize=100000, throttle=1, reconnect=true)

用户可以利用下方代码完成模型预测与预警计算操作。

在进行模型预测与预警计算以前，需提前建立预测表 predictTable 与预警表 warningTable 以存储预测数据与预警数据。

由于模型训练需要用到 10s 数据，因此设定 10s 为一个周期，每 10s 进行一次模型训练、模型预测以及预警计算。

关于预警计算，设定预测值与聚合值的值差异大于 20% 时（即预测值大于真实值 20% 或预测值小于真实值 20%）判定为异常值，当异常率大于 0.215 时进行预警。用户可根据需求自行调整。

//新建预测表
 preColNames = `time`deviceCode`wind`humidity`air_pressure`temperature`propertyValue_predicted
 preColTypes = orgColTypes
 enableTableShareAndPersistence(table = streamTable(100000:0,preColNames,preColTypes), tableName=`predictTable, cacheSize = 5000000)

 //新建预警表
 warnColNames = `time`abnormal_rate`whether_warning
 warnColTypes = [TIMESTAMP,DOUBLE,INT]
 enableTableShareAndPersistence(table = streamTable(100000:0,warnColNames,warnColTypes), tableName=`warningTable, cacheSize = 5000000)

/*
 * 根据 已有周期数据，对未来数据进行预测
 */
 def predictAndwarning(devices_number,rate,interval_past,interval_future,warning_threshold,mutable whether_start,DataTable_dfs,AggrTable_dfs,mutable predictTable,mutable warningTable){
    do{
      if(whether_start==false) {
        curtime = select top 1 time from DataTable_dfs  //第一次是从表中查询最开始的时间
        curtime = curtime[`time][0]
        curtime_aggr = select top 1 time from AggrTable_dfs //每次尽量去与训练数据的时间段同步
        curtime_aggr = curtime_aggr[`time][0]
        whether_start = true
      }
      curtime = curtime + interval_past*1000 //以后是直接往后推，最开始的时间往后推 10s中的间隔时间
      table_origin = select * from DataTable_dfs where time<=curtime and time>(curtime - interval_past*1000) //查询当前时间前 interval 秒的数据
      if(table_origin.size()<interval_past*devices_number*rate) //如果查询数据小于正常查询到的数据数目
      {
        
         curtime = select top 1 time from DataTable_dfs order by time desc //从表中查询最近的时间作为开始时间
         curtime = curtime[`time][0]
         table_origin = select * from DataTable_dfs where time<=curtime and time>(curtime - interval_past*1000) //查询当前时间的前 interval 秒的数据
      }


      //训练模型
      factors = sql(sqlCol([`wind,`humidity,`air_pressure,`temperature]), table_origin).eval()
      labels = table_origin[`propertyValue]
      model = knn(labels,factors,"regressor", 200);


      //模型预测
      curtime_aggr = curtime_aggr + interval_past*1000
      table_aggr = select * from AggrTable_dfs where time<=curtime_aggr and time>(curtime_aggr - interval_past*1000) order by deviceCode	 //查询当前时间前 interval 秒的数据
      if(table_aggr.size()<interval_past*devices_number*rate/10) //如果查询数据小于正常查询到的数据数目
      {
         curtime_aggr = select top 1 time from AggrTable_dfs order by time desc //从表中查询最近的时间作为开始时间
         curtime_aggr = curtime_aggr[`time][0]
         table_aggr = select * from AggrTable_dfs where time<=curtime_aggr and time>(curtime_aggr - interval_past*1000) //查询当前时间的前 interval 秒的数据
      }
      pre_labels = sql(sqlCol([`wind,`humidity,`air_pressure,`temperature]), table_aggr).eval()
      //pre_values = predict(model,pre_labels)
      //////////////////////////////暂时用随机值代替///////////////////////////////////////////
      table_number =  table_aggr.size()
      x1 = randNormal(25,2,table_number)  //训练数据
      x2 = randNormal(55,5,table_number)
      x3 = randNormal(1.01325,0.00001,table_number)
      x4 = randNormal(75,5,table_number)
      x5 = int(randNormal(10,3,table_number))
      b1 = randNormal(0.4,0.05,table_number) 
      b2 = randNormal(0.3,0.05,table_number)
      b3 = randNormal(0.2,0.05,table_number)
      b4 = randNormal(0.09,0.05,table_number)
      b5 = randNormal(0.01,0.001,table_number)
      bias = randNormal(5,1,table_number)
      propertyValue = int(b1*x1*10 + b2*x2*2 + b3*x3*1000 + b4*x4 + b5*x5 +bias)
      pre_values = propertyValue //暂时用随机值代替
      /////////////////////////////////////////////////////////////////////////////////////////////
    
      time =take(curtime_aggr + interval_future*1000+interval_future+(0..(1000-1)),table_number) //
      deviceCode = sort(take(1..devices_number,table_number))
      predicttempTable = table(time,deviceCode,pre_labels,pre_values as `propertyValue_predicted)
      predictTable.append!(predicttempTable) //预测结果导入流表

      //进行预警
      contrastTable = select propertyValue,propertyValue_predicted from lj(table_aggr,predicttempTable,`wind`humidity`air_pressure`temperature) //利用左连接
      abnormal_count = exec count(*) from contrastTable where propertyValue_predicted<0.8*propertyValue or propertyValue_predicted>1.2*propertyValue
      
      warning_time = curtime_aggr //进行预警的时间
      abnormal_rate = abnormal_count*1.0 / table_number
      whether_warning = 0 //默认不进行预警
      if(abnormal_rate>warning_threshold) whether_warning = 1 //当异常率超过阈值进行预警

      insert into warningTable values(warning_time,abnormal_rate, whether_warning);
      
      sleep(10000) //每10s进行一次预测
    }while(true)
}
devices_number = 100 //设备数目
rate = 1000 //每台设备 每秒钟1000条数据
interval_past = 10 ////查询过去10秒的数据
interval_future = 10 //预测未来第10秒的数据
warning_threshold = 0.215 //当异常值率大于0.215时进行预警
whether_start = false //标记
DataTable_dfs = loadTable(dataTable_dbname,dataTable_tbname)
AggrTable_dfs = loadTable(aggr_dbname,aggr_tbname)
submitJob(`predictAndwarning,`predictAndwarning,predictAndwarning,devices_number,rate,interval_past,interval_future,warning_threshold,whether_start,DataTable_dfs,AggrTable_dfs,predictTable,warningTable) //每10s 预测一次