数学建模实战-城市供水管网片区用水异常模式识别(python实现)

文章目录问题解法1数据分析1.1数据可视化1.1.1读取所有数据并且画图1.1.2读取某一天数据并画图1.2用水流量2典型用水模式2.1独热编码2.1.1映射日期2.1.2编码向量2.2聚类评估2.3 K均值分类问题解法作为小白的我，自己实在没找出什么合适的办法，便只能查阅大佬们的论文，发现了一篇非常棒的，就复盘了一下。原文链接如下http://www.yndxxb.ynu.edu.cn...

喜欢coding的谢同学

7558人浏览 · 2019-04-13 14:43:20

喜欢coding的谢同学 · 2019-04-13 14:43:20 发布

问题

在这里插入图片描述

解法

作为小白的我，自己实在没找出什么合适的办法，便只能查阅大佬们的论文，发现了一篇非常棒的，就复盘了一下。
原文链接如下
http://www.yndxxb.ynu.edu.cn/yndxxb/article/2018/0258-7971-40-5-879.html#outline_anchor_7

1数据分析

1.1数据可视化

这一部分难度不大，就是读数据画图

1.1.1读取所有数据并且画图

读所有数据只需打开excel，取出对应字段的内容就行了，需要注意时间类型的处理

data = xlrd.open_workbook("2018年东北三省数学建模联赛F题附件.xls")  # 打开excel
    table = data.sheet_by_name("朗贝20140415-0612")  # 读sheet
    times = []
    flow1 = []
    flow2 = []
    for i in range(1, table.nrows):
        row = table.row_values(i)  # 行的数据放在数组里
        date = xldate_as_tuple(row[0], 0)		# 日期元组
        times.append(datetime.datetime(*date))
        flow1.append(row[1])
        flow2.append(row[2])
    draw_all(times, flow1, flow2)

根据读取的数据直接绘图，依旧需要注意时间类型的处理

def draw_all(times, flow1, flow2):
    # 配置横坐标
    length = len(times)
    coordinate = [times[0], times[int(length / 4)], times[int(2 * length / 4)], times[int(3 * length / 4)],
                  times[length - 1]]
    plt.xticks(coordinate)
    plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))
    # 配置标题、轴名称
    plt.title("ALL DAY")
    plt.xlabel("time")
    plt.ylabel("flow (L/s)")
    # Plot
    plt.plot(times, flow1, label='DMA1')
    plt.plot(times, flow2, label='DMA2')
    plt.legend(loc=4)  # 指定legend的位置,读者可以自己help它的用法
    plt.gcf().autofmt_xdate()  # 自动旋转日期标记
    plt.show()

结果如下
在这里插入图片描述

1.1.2读取某一天数据并画图

这里需要注意的是如何判断读取的数据为指定日期

def read_oneday(month, day):
    '''
    :param month: 指定月份
    :param day: 指定天数
    :return:    无返回值，而是接着去调用绘图方法。也可将times、flow返回去主函数调用
    '''
    data = xlrd.open_workbook("2018年东北三省数学建模联赛F题附件.xls")  # 打开excel
    table = data.sheet_by_name("朗贝20140415-0612")  # 读sheet
    times = []
    flow1 = []
    flow2 = []
    for i in range(1, table.nrows):
        row = table.row_values(i)  # 行的数据放在数组里
        date = xldate_as_tuple(row[0], 0)
        if date[1] == month and date[2] == day:
            times.append(datetime.datetime(*date))
            flow1.append(row[1])
            flow2.append(row[2])
    draw_oneday(times, flow1, flow2)

与之前的画法差别不大，坐标轴、标题需要修改，依旧需要注意日期的处理

def draw_oneday(times, flow1, flow2):
    # 配置横坐标
    length = len(times)
    coordinate = [times[0], times[int(length / 4)], times[int(2 * length / 4)], times[int(3 * length / 4)],
                  times[length - 1]]
    plt.xticks(coordinate)
    plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))
    # 配置标题、轴名称
    plt.title(times[0].strftime('%Y/%m/%d'))
    plt.xlabel("time")
    plt.ylabel("flow (L/s)")
    # Plot
    plt.plot(times, flow1, label='DMA1')
    plt.plot(times, flow2,  label='DMA2')
    plt.legend(loc=4)  # 指定legend的位置,读者可以自己help它的用法
    plt.gcf().autofmt_xdate()  # 自动旋转日期标记
    plt.show()

结果如下
在这里插入图片描述

1.2用水流量

先求出每天的漏水量

def find_leak():
    data = xlrd.open_workbook("2018年东北三省数学建模联赛F题附件.xls")  # 打开excel
    table = data.sheet_by_name("朗贝20140415-0612")  # 读sheet
    # 漏水量，字典类型，键为日期，值为当天的漏水量
    leak1 = {}
    leak2 = {}
    # 当天2-5点的最小流量作为漏水量
    minFlow1 = 999
    minFlow2 = 999
    for i in range(1, table.nrows):
        row = table.row_values(i)  # 行的数据放在数组里
        date = xldate_as_tuple(row[0], 0)
        if i != 1 and date[2] != oldDate[2]:  # 比较日期变化
            key = str(oldDate[1])+"-"+str(oldDate[2])
            leak1[key] = minFlow1
            leak2[key] = minFlow2
            minFlow1 = 999
            minFlow2 = 999
        if date[3] >= 2 and date[3] <= 5:
            minFlow1 = min(minFlow1, row[1])
            minFlow2 = min(minFlow2, row[2])
        oldDate = date  # 记录上一个日期
        if i == table.nrows - 1:                        ##单独处理最后一天
            key = str(date[1]) + "-" + str(date[2])
            leak1[key] = minFlow1
            leak2[key] = minFlow2
    return leak1, leak2

根据漏水量计算用水量，注意将负数替换为0

def cal_usewater():
    data = xlrd.open_workbook("result.xls")  # 打开excel
    table = data.sheet_by_name("朗贝20140415-0612")  # 读sheet
    from xlutils.copy import copy
    wb = xlutils.copy.copy(data)
    sheet = wb.get_sheet(0)
    sheet.write(0, 3, 'use1')
    sheet.write(0, 4, 'use2')
    leak1, leak2 = find_leak()
    for i in range(1, table.nrows):
        row = table.row_values(i)  # 行的数据放在数组里
        date = xldate_as_tuple(row[0], 0)
        key = str(date[1]) + "-" + str(date[2])
        sheet.write(i, 3, str(max(row[1] - leak1[key], 0)))
        sheet.write(i, 4, str(max(row[2] - leak2[key], 0)))
    wb.save("result.xls")

由于结果太长，只展示部分。记得复制一下原来的文档，直接将结果追加到文件后，便于再次使用
在这里插入图片描述

2典型用水模式

2.1独热编码

这一步是将样本集映射到欧式空间，便于算距离，其中时间特征采用独热编码，用水流量不变。

2.1.1映射日期

def get_timeMap():
    data = xlrd.open_workbook("2018年东北三省数学建模联赛F题附件.xls")  # 打开excel
    table = data.sheet_by_name("朗贝20140415-0612")  # 读sheet
    timeMap = {}
    for i in range(1, table.nrows):
        row = table.row_values(i)  # 行的数据放在数组里
        date = xldate_as_tuple(row[0], 0)
        if date[1] == 4 and date[2] == 15:
            key = str(date[3]) + ":" + str(date[4]) + ":" +str(date[5])
            timeMap[key] = i
    enc = preprocessing.OneHotEncoder()
    studyData = []
    for item in timeMap.values():
        studyData.append([item])
    enc.fit(studyData)  # fit来学习编码
    endcodeData = enc.transform(studyData).toarray().tolist()   # 进行编码
    count = 0
    for key in timeMap.keys():
        timeMap[key] = endcodeData[count]
        count = count + 1
    return timeMap

2.1.2编码向量

将映射的日期与用水量组合构成样本向量

timeMap = get_timeMap()
    data = xlrd.open_workbook("2018年东北三省数学建模联赛F题附件.xls")  # 打开excel
    table = data.sheet_by_name("朗贝20140415-0612")  # 读sheet
    vectors1 = []
    vectors2 = []
    for i in range(1, table.nrows):
        row = table.row_values(i)  # 行的数据放在数组里
        date = xldate_as_tuple(row[0], 0)
        key = str(date[3]) + ":" + str(date[4]) + ":" + str(date[5])
        item1 = copy.deepcopy(timeMap[key])         # 注意要用深层复制
        item1.append(use1[i-1])
        vectors1.append(item1)
        item2 = copy.deepcopy(timeMap[key])
        item2.append(use2[i - 1])
        vectors2.append(item2)
    return np.array(vectors1), np.array(vectors2)

结果如下
在这里插入图片描述

2.2聚类评估

采取轮廓系数

def choose_by_silhouette(X):
    tests = [2, 3, 4, 5]
    sc_scores = []
    for k in tests:
        kmeans_model = KMeans(n_clusters=k).fit(X)
        sc_score = metrics.silhouette_score(X, kmeans_model.labels_, metric='euclidean')
        sc_scores.append(sc_score)
return sc_scores

k取2、3、4、5的结果如下
在这里插入图片描述

2.3 K均值分类

在分类之前需要再获取时间列表，用之前read_all()的方法即可，只是需要将times中的日期变为同一天

 for count in range(0, len(times)):
         times[count] = times[count].replace(2014, 4, 15)

传入之前获得的参数即可分类

def cluster_by_K(x1, x2, X, k, name):
    '''
    :param x1:   时间列表
    :param x2:   用水量列表
    :param X:    样本向量构成的矩阵
    :param k:    k值
    :param name: 图片名称
    :return:
    '''
    data = xlrd.open_workbook("result.xls")  # 打开excel
    from xlutils.copy import copy
    wb = xlutils.copy.copy(data)
    sheet = wb.get_sheet(0)
    if name == 'DMA1':
        cloumn = 3
    else:
        cloumn = 4
    sheet.write(0, cloumn, name)
    plt.figure(figsize=(8, 6))
    colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'b']
    markers = ['o', 's', 'D', 'v', '^', 'p', '*', '+']
    kmeans_model = KMeans(n_clusters=k).fit(X)
    for i, l in enumerate(kmeans_model.labels_):
        plt.subplot(k, 1, 1 + l)
        sheet.write(i+1, cloumn, str(l))
        print(i)
        plt.plot(x1[i], x2[i], color=colors[l], marker=markers[l])
    wb.save("result.xls")
    for t in range(1, k+1):
        plt.subplot(k, 1, t)
        plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%H:%M:%S'))
        plt.gcf().autofmt_xdate()  # 自动旋转日期标记
        plt.title('第'+str(t)+'种', fontproperties=font)
    plt.savefig(name+'.png')
    return kmeans_model.labels_

DMA1结果如下
在这里插入图片描述
DMA2结果如下

同时将结果追加到文档后面便于再次使用

3异常模式识别

3.1提取特征

3.1.1计算变化量

def cal_change():
    data = xlrd.open_workbook("result.xls")  # 打开excel
    table = data.sheet_by_name("朗贝20140415-0612")  # 读sheet
    from xlutils.copy import copy
    wb = xlutils.copy.copy(data)
    sheet = wb.get_sheet(0)
    sheet.write(0, 7, 'change1')
    sheet.write(0, 8, 'change2')
    oldRow = table.row_values(1)
    for i in range(2, table.nrows):
        row = table.row_values(i)  # 行的数据放在数组里
        sheet.write(i, 7, abs(float(row[3]) - float(oldRow[3])))
        sheet.write(i, 8, abs(float(row[4]) - float(oldRow[4])))
        oldRow = row
    wb.save("result.xls")

将计算出来的变化量追加到文档后面
在这里插入图片描述

3.1.2构造矩阵

对不同的模式构造不同的矩阵用于之后的聚类

def get_change():
    data = xlrd.open_workbook("result.xls")  # 打开excel
    table = data.sheet_by_name("朗贝20140415-0612")  # 读sheet
    DMA10 = []
    DMA11 = []
    DMA20 = []
    DMA21 = []
    DMA22 = []
    oldRow = table.row_values(1)
    for i in range(2, table.nrows):
        row = table.row_values(i)  # 行的数据放在数组里
        item = [float(row[3]), abs(float(row[3]) - float(oldRow[3]))]
        if int(row[5]) == 0:
            DMA10.append(item)
        else:
            DMA11.append(item)
        item = [float(row[4]), abs(float(row[4]) - float(oldRow[4]))]
        if int(row[6]) == 0:
            DMA20.append(item)
        elif int(row[6]) == 1:
            DMA21.append(item)
        else:
            DMA22.append(item)
        oldRow = row
    return np.array(DMA10), np.array(DMA11),np.array(DMA20), np.array(DMA21), np.array(DMA22)

3.2DBSCAN聚类

聚类后根据标签值绘制点的分布图

def find_by_DBSCAN(X, eps, min_samples, name):
    y_pred = DBSCAN(eps=eps, min_samples=min_samples).fit_predict(X)
    plt.scatter(X[:, 0], X[:, 1], c=y_pred)
    # plt.show()
    plt.savefig(name)

记得对不同的类型调整参数，我这里给出eps=10，min_samples=10的DMA1模式0的情况(对应原文中的模式1)
在这里插入图片描述

3.3打印异常点

在聚类之后调用即可

def print_exception(lables, name):
    if name == "DMA10":
        columnIndex = 5
        typeValue = 0
    elif name == "DMA11":
        columnIndex = 5
        typeValue = 1
    elif name == "DMA20":
        columnIndex = 6
        typeValue = 0
    elif name == "DMA21":
        columnIndex = 6
        typeValue = 1
    elif name == "DMA22":
        columnIndex = 6
        typeValue = 2
    print("---"+name+"异常点如下---")
    data = xlrd.open_workbook("result.xls")  # 打开excel
    table = data.sheet_by_name("朗贝20140415-0612")  # 读sheet
    count = 0
    for i in range(2, table.nrows):
        row = table.row_values(i)  # 行的数据放在数组里
        if int(row[columnIndex]) == int(typeValue):
            count = count + 1
            if int(lables[count - 1]) != 0:
                date = xldate_as_tuple(row[0], 0)
                print([date, row[1], row[2]])

可以看到结果与论文中给出的差异不大，需要注意的是我的0与1与论文的模式0与1正好相反
在这里插入图片描述

CSDN学习社区

CSDN联合极客时间，共同打造面向开发者的精品内容学习社区，助力成长！

更多推荐

嵌入式作业（七）：基于Ardunio的STM32串口通信

嵌入式作业（七）0作业要求1Ardunio 完成STM32的串口通信（1）安装Ardunio IDE（2）stm32串口通信2关于 stduino IDE0作业要求安装 Ardunio IDE 和相关软件支持库，在Ardunio 完成STM32板子的串口通信程序：（1）持续向串口输出“Hello world！”；（2）当接收到“stop!”时，停止输出。网上有一个国人版的MCU集成开发平台， st

CSDN学习社区

JDBC详解

JDBC文章目录JDBC什么是JDBC?JDBC驱动程序:Java使用JDBC访问数据库的步骤:设置classpath:Oracle连接字符串的书写格式:简单的例子:常用数据库的驱动程序及JDBC URL:Oracle数据库:SQL Server数据库MySQL数据库Access数据库PreparedStatement接口:JNDI-数据源（Data Source）与连接池（Connection

CSDN学习社区

“模式识别与机器学习”学习笔记no2.再谈感知机

接**上篇：上篇主要进行了PLA，Pocket算法的理论过程分析和在给定数据集上利用pocket算法对数据集进行分类学习，得到错分数量最少的分类面。上篇中pocket算法的过程已经进行了编程和测试，框架已经建立了起来，这一篇主要上篇中没有提到或涉及不深的几个问题。1.数据集的构造。上篇是直接使用了题目给的向量，这次来根据正态分布来产生数据集。np.random.normal函数可以根据均值和方差生