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()