energy-data-trans/tmp_file/ClassIdentifier_2.py

import os

import numpy as np
from pandas import DataFrame

from service.plt_service import get_base_wind_and_power
from utils.draw.draw_file import scatter
from utils.file.trans_methods import read_file_to_df


class ClassIdentifier(object):

    def __init__(self, wind_turbine_number, file_path: str = None, origin_df: DataFrame = None, index='time_stamp',
                 wind_velocity='wind_velocity',
                 active_power='active_power'):
        """
        :param wind_turbine_number: The wind turbine number.
        :param file_path: The file path of the input data.
        :param origin_df: The pandas DataFrame containing the input data.
        :param index: 索引字段
        :param wind_velocity: 风速字段
        :param active_power: 有功功率字段
        """
        self.wind_turbine_number = wind_turbine_number
        self.index = index
        self.wind_velocity = wind_velocity
        self.active_power = active_power

        self.rated_wind_speed = 'rated_wind_speed'
        self.rated_capacity = 'rated_capacity'

        if file_path is None and origin_df is None:
            raise ValueError("Either file_path or origin_df should be provided.")

        if file_path:
            self.df = read_file_to_df(file_path)
        else:
            self.df = origin_df

        self.df = self.df.set_index(keys=self.index)

    def identifier(self):
        # 风速 和 有功功率 df
        wind_and_power_df = self.df[[self.wind_velocity, self.active_power]]
        wind_and_power_df.reset_index(inplace=True)
        wind_and_power_df_count = wind_and_power_df.shape[0]
        PowerMax = wind_and_power_df[self.active_power].max()
        PowerRated = np.ceil(PowerMax / 100) * 100
        PRated = 1500  # 额定功率1500kw,可改为2000kw
        VCutOut = 25
        # 网格法确定风速风向分区数量，功率方向分区数量，
        # PNum = (PRated+100)/25  #功率分区间隔25kW
        PNum = int(np.ceil(PowerRated / 25))  # 功率分区间隔25kW
        VNum = int(np.ceil(VCutOut / 0.25))  # 风速分区间隔0.25m/s

        # 存储功率大于零的运行数据
        DzMarch809 = np.zeros([wind_and_power_df_count, 2], dtype=float)
        nCounter1 = 0
        for i in range(wind_and_power_df_count):
            if wind_and_power_df.loc[i, self.active_power] > 0:
                DzMarch809[nCounter1, 0] = wind_and_power_df.loc[i, self.wind_velocity]
                DzMarch809[nCounter1, 1] = wind_and_power_df.loc[i, self.active_power]

                nCounter1 = nCounter1 + 1

        # 统计各网格落入的散点个数
        if VNum == 1:
            XBoxNumber = np.ones([PNum], dtype=int)
        else:
            XBoxNumber = np.ones([PNum, VNum], dtype=int)
        nWhichP = -1
        nWhichV = -1
        for i in range(nCounter1):
            for m in range(PNum):
                if m * 25 < DzMarch809[i, 1] <= (m + 1) * 25:
                    nWhichP = m
                    break
            for n in range(VNum):
                if ((n + 1) * 0.25 - 0.125) < DzMarch809[i, 0] <= ((n + 1) * 0.25 + 0.125):
                    nWhichV = n
                    break

            if nWhichP > -1 and nWhichV > -1:
                XBoxNumber[nWhichP, nWhichV] = XBoxNumber[nWhichP, nWhichV] + 1

        for m in range(PNum):
            for n in range(VNum):
                XBoxNumber[m, n] = XBoxNumber[m, n] - 1

        # 在功率方向将网格内散点绝对个数转换为相对百分比，备用
        PBoxPercent = np.zeros([PNum, VNum], dtype=float)
        PBinSum = np.zeros(PNum, dtype=int)

        for i in range(PNum):
            for m in range(VNum):
                PBinSum[i] = PBinSum[i] + XBoxNumber[i, m]

            for m in range(VNum):
                if PBinSum[i] > 0:
                    PBoxPercent[i, m] = XBoxNumber[i, m] / PBinSum[i] * 100

        # 在风速方向将网格内散点绝对个数转换为相对百分比，备用
        VBoxPercent = np.zeros([PNum, VNum], dtype=float)
        VBinSum = np.zeros(VNum, dtype=int)

        for i in range(VNum):
            for m in range(PNum):
                VBinSum[i] = VBinSum[i] + XBoxNumber[m, i]

            for m in range(PNum):
                if VBinSum[i] > 0:
                    VBoxPercent[m, i] = XBoxNumber[m, i] / VBinSum[i] * 100

        # 以水平功率带方向为准，分析每个水平功率带中，功率主带中心，即找百分比最大的网格位置。
        PBoxMaxIndex = np.zeros(PNum, dtype=int)  # 水平功率带最大网格位置索引
        PBoxMaxP = np.zeros(PNum, dtype=int)  # 水平功率带最大网格百分比

        for m in range(PNum):
            # 确定每一水平功率带的最大网格位置索引即百分比值
            PBoxMaxP[m], PBoxMaxIndex[m] = PBoxPercent[m, :].max(), PBoxPercent[m, :].argmax()

        # 以垂直风速方向为准，分析每个垂直风速带中，功率主带中心，即找百分比最大的网格位置。
        VBoxMaxIndex = np.zeros(VNum, dtype=int)
        VBoxMaxV = np.zeros(VNum, dtype=int)

        for m in range(VNum):
            [VBoxMaxV[m], VBoxMaxIndex[m]] = VBoxPercent[:, m].max(), VBoxPercent[:, m].argmax()

        # 切入风速特殊处理，如果切入风速过于偏右，向左拉回
        if PBoxMaxIndex[0] > 14:
            PBoxMaxIndex[0] = 9

        # 以水平功率带方向为基准，进行分析
        DotDense = np.zeros(PNum, dtype=int)  # 每一水平功率带的功率主带包含的网格数
        DotDenseLeftRight = np.zeros([PNum, 2], dtype=int)  # 存储每一水平功率带的功率主带以最大网格为中心，向向左，向右扩展的网格数
        DotValve = 90  # 从中心向左右对称扩展网格的散点百分比和的阈值。

        for i in range(PNum - 6):  # 从最下层水平功率带1开始，向上到第PNum-6个水平功率带（额定功率一下水平功率带），逐一分析
            PDotDenseSum = PBoxMaxP[i]  # 以中心最大水平功率带为基准，向左向右对称扩展网格，累加各网格散点百分比
            iSpreadRight = 1
            iSpreadLeft = 1
            while PDotDenseSum < DotValve:

                if (PBoxMaxIndex[i] + iSpreadRight) < VNum - 1:
                    PDotDenseSum = PDotDenseSum + PBoxPercent[i, PBoxMaxIndex[i] + iSpreadRight]  # 向右侧扩展
                    iSpreadRight = iSpreadRight + 1

                if (PBoxMaxIndex[i] + iSpreadRight) > VNum - 1:
                    break

                if (PBoxMaxIndex[i] - iSpreadLeft) > 0:
                    PDotDenseSum = PDotDenseSum + PBoxPercent[i, PBoxMaxIndex[i] - iSpreadLeft]  # 向左侧扩展
                    iSpreadLeft = iSpreadLeft + 1

                if (PBoxMaxIndex[i] - iSpreadLeft) <= 0:
                    break

            iSpreadRight = iSpreadRight - 1
            iSpreadLeft = iSpreadLeft - 1
            # 向左右对称扩展完毕

            DotDenseLeftRight[i, 0] = iSpreadLeft
            DotDenseLeftRight[i, 1] = iSpreadRight
            DotDense[i] = iSpreadLeft + iSpreadRight + 1

        # 各行功率主带右侧宽度的中位数最具有代表性
        DotDenseWidthLeft = np.zeros([PNum - 6, 1], dtype=int)
        for i in range(PNum - 6):
            DotDenseWidthLeft[i] = DotDenseLeftRight[i, 1]

        MainBandRight = np.median(DotDenseWidthLeft)

        # 散点向右显著延展分布的水平功率带为限功率水平带
        PowerLimit = np.zeros([PNum, 1], dtype=int)  # 各水平功率带是否为限功率标识，==1：是；==0：不是
        WidthAverage = 0  # 功率主带平均宽度
        WidthVar = 0  # 功率主带方差
        # PowerLimitValve = 6    #限功率主带判别阈值
        PowerLimitValve = np.ceil(MainBandRight) + 3  # 限功率主带判别阈值

        nCounterLimit = 0
        nCounter = 0

        for i in range(PNum - 6):
            if DotDenseLeftRight[i, 1] > PowerLimitValve and PBinSum[i] > 20:  # 如果向右扩展网格数大于阈值，且该水平功率带点总数>20，是
                PowerLimit[i] = 1
                nCounterLimit = nCounterLimit + 1

            if DotDenseLeftRight[i, 1] <= PowerLimitValve:
                WidthAverage = WidthAverage + DotDenseLeftRight[i, 1]  # 统计正常水平功率带右侧宽度
                nCounter = nCounter + 1

        WidthAverage = WidthAverage / nCounter  # 功率主带平均宽度

        # 各水平功率带的功率主带宽度的方差，反映从下到上宽度是否一致，或是否下宽上窄等异常情况
        for i in range(PNum - 6):
            if DotDenseLeftRight[i, 1] <= PowerLimitValve:
                WidthVar = WidthVar + (DotDenseLeftRight[i, 1] - WidthAverage) * (
                        DotDenseLeftRight[i, 1] - WidthAverage)

        # 对限负荷水平功率带的最大网格较下面相邻层显著偏右，拉回
        for i in range(1, PNum - 6):
            if PowerLimit[i] == 1 and abs(PBoxMaxIndex[i] - PBoxMaxIndex[i - 1]) > 5:
                PBoxMaxIndex[i] = PBoxMaxIndex[i - 1] + 1

        # 输出各层功率主带的左右边界网格索引
        DotDenseInverse = np.zeros([PNum, 2], dtype=int)

        for i in range(PNum):
            DotDenseInverse[i, :] = DotDenseLeftRight[PNum - i - 1, :]

        # 功率主带的右边界
        CurveWidthR = int(np.ceil(WidthAverage) + 2)

        # CurveWidthL = 6    #功率主带的左边界
        CurveWidthL = CurveWidthR

        BBoxLimit = np.zeros([PNum, VNum], dtype=int)  # 网格是否为限功率网格的标识，如果为限功率水平功率带，从功率主带右侧边缘向右的网格为限功率网格
        for i in range(2, PNum - 6):
            if PowerLimit[i] == 1:
                for j in range(PBoxMaxIndex[i] + CurveWidthR, VNum):
                    BBoxLimit[i, j] = 1

        BBoxRemove = np.zeros([PNum, VNum], dtype=int)  # 数据异常需要剔除的网格标识，标识==1：功率主带右侧的欠发网格；==2：功率主带左侧的超发网格
        for m in range(PNum - 6):
            for n in range(PBoxMaxIndex[m] + CurveWidthR, VNum):
                BBoxRemove[m, n] = 1

            for n in range(PBoxMaxIndex[m] - CurveWidthL, -1, -1):
                BBoxRemove[m, n] = 2

        # 确定功率主带的左上拐点，即额定风速位置的网格索引
        CurveTop = np.zeros(2, dtype=int)
        CurveTopValve = 3  # 网格的百分比阈值
        BTopFind = 0
        for m in range(PNum - 4 - 1, -1, -1):
            for n in range(VNum):
                if VBoxPercent[m, n] > CurveTopValve and XBoxNumber[m, n] >= 10:  # 如左上角网格的百分比和散点个数大于阈值。
                    CurveTop[0] = m
                    CurveTop[1] = n
                    BTopFind = 1
                    break

            if BTopFind == 1:
                break

        IsolateValve = 3
        for m in range(PNum - 6):
            for n in range(PBoxMaxIndex[m] + CurveWidthR, VNum):
                if PBoxPercent[m, n] < IsolateValve:
                    BBoxRemove[m, n] = 1

        # 功率主带顶部宽度
        CurveWidthT = 2
        for m in range(PNum - CurveWidthT - 1, PNum):
            for n in range(VNum):
                BBoxRemove[m, n] = 3  # 网格为额定功率以上的超发点

        # 功率主带拐点左侧的欠发网格标识
        for m in range(PNum - 5 - 1, PNum):
            for n in range(CurveTop[1] - 1):
                BBoxRemove[m, n] = 2

        # 以网格的标识，决定该网格内数据的标识。Dzwind_and_power_dfSel功率非零数据的标识位。散点在哪个网格，此网格的标识即为该点的标识
        Dzwind_and_power_dfSel = np.zeros(nCounter1, dtype=int)  # -1:停机 0:好点  1:欠发功率点；2:超发功率点；3:额定风速以上的超发功率点 4: 限电
        nWhichP = -1
        nWhichV = -1
        nBadA = 0

        for i in range(nCounter1):
            for m in range(PNum):
                if m * 25 < DzMarch809[i, 1] <= (m + 1) * 25:
                    nWhichP = m
                    break

            for n in range(VNum):
                if ((n + 1) * 0.25 - 0.125) < DzMarch809[i, 0] <= ((n + 1) * 0.25 + 0.125):
                    nWhichV = n
                    break

            if nWhichP > -1 and nWhichV > -1:

                if BBoxRemove[nWhichP, nWhichV] == 1:
                    Dzwind_and_power_dfSel[i] = 1
                    nBadA = nBadA + 1

                if BBoxRemove[nWhichP, nWhichV] == 2:
                    Dzwind_and_power_dfSel[i] = 2

                if BBoxRemove[nWhichP, nWhichV] == 3:
                    Dzwind_and_power_dfSel[i] = 0  # 3  # 额定风速以上的超发功率点认为是正常点，不再标识。

        # 限负荷数据标识方法2：把数据切割为若干个窗口。对每一窗口，以第一个点为基准，连续nWindowLength个数据的功率在方差范围内，呈现显著水平分布的点
        nWindowLength = 3
        LimitWindow = np.zeros(nWindowLength, dtype=float)
        PowerStd = 15  # 功率波动方差
        nWindowNum = int(np.floor(nCounter1 / nWindowLength))
        PowerLimitUp = PRated - 300
        PowerLimitLow = 200
        for i in range(nWindowNum):
            for j in range(nWindowLength):
                LimitWindow[j] = DzMarch809[i * nWindowLength + j, 1]

            bAllInAreas = 1
            for j in range(nWindowLength):
                if LimitWindow[j] < PowerLimitLow or LimitWindow[j] > PowerLimitUp:
                    bAllInAreas = 0

            if bAllInAreas == 0:
                continue

            UpLimit = LimitWindow[0] + PowerStd
            LowLimit = LimitWindow[0] - PowerStd
            bAllInUpLow = 1
            for j in range(1, nWindowLength):
                if LimitWindow[j] < LowLimit or LimitWindow[j] > UpLimit:
                    bAllInUpLow = 0

            if bAllInUpLow == 1:
                for j in range(nWindowLength):
                    Dzwind_and_power_dfSel[i * nWindowLength + j] = 4  # 标识窗口内的数据为限负荷数据

        nSmooth = 0
        for i in range(PNum - 6):
            PVLeftDown = np.zeros(2, dtype=float)
            PVRightUp = np.zeros(2, dtype=float)

            if (PBoxMaxIndex[i + 1] - PBoxMaxIndex[i]) >= 1:
                PVLeftDown[0] = (PBoxMaxIndex[i] + 1 + CurveWidthR) * 0.25 - 0.125
                PVLeftDown[1] = i * 25

                PVRightUp[0] = (PBoxMaxIndex[i + 1] + 1 + CurveWidthR) * 0.25 - 0.125
                PVRightUp[1] = (i + 1) * 25

                for m in range(nCounter1):
                    if DzMarch809[m, 0] > PVLeftDown[0] and DzMarch809[m, 0] < PVRightUp[0] and PVLeftDown[1] < \
                            DzMarch809[m, 1] < PVRightUp[1]:  # 在该锯齿中
                        if (DzMarch809[m, 1] - PVLeftDown[1]) / (DzMarch809[m, 0] - PVLeftDown[0]) > (
                                PVRightUp[1] - PVLeftDown[1]) / (
                                PVRightUp[0] - PVLeftDown[0]):  # 斜率大于对角连线，则在锯齿左上三角形中，选中
                            Dzwind_and_power_dfSel[m] = 0
                            nSmooth = nSmooth + 1

        print("nSmooth", nSmooth)

        wind_and_power_df.loc[:, 'marker'] = -1
        wind_and_power_df.loc[
            wind_and_power_df[wind_and_power_df[self.active_power] > 0].index, 'marker'] = Dzwind_and_power_dfSel
        wind_and_power_df.to_csv("test.csv", index=False, encoding='utf-8')

        # wind_and_power_df = wind_and_power_df[wind_and_power_df['marker'] == 0]
        color_map = {-1: 'red', 0: 'green', 1: 'blue', 2: 'black', 3: 'orange', 4: 'magenta'}
        c = wind_and_power_df['marker'].map(color_map)

        # -1:停机 0:好点  1:欠发功率点；2:超发功率点；3:额定风速以上的超发功率点 4: 限电
        legend_map = {"停机": 'red', "好点": 'green', "欠发": 'blue', "超发": 'black', "额定风速以上的超发": 'orange', "限电": 'magenta'}

        scatter("测试matlab结果", x_label='风速', y_label='有功功率', x_values=wind_and_power_df[self.wind_velocity].values,
                y_values=wind_and_power_df[self.active_power].values, color=c, col_map=legend_map,
                save_file_path=os.path.dirname(__file__) + os.sep + '测试matlab结果均值.png')

    def run(self):
        # Implement your class identification logic here
        self.identifier()


if __name__ == '__main__':
    test = ClassIdentifier('test', r"D:\中能智能\matlib计算相关\好点坏点matlib计算\A01.csv", index='时间',
                           wind_velocity='风速',
                           active_power='功率')

    test.run()