# -*- coding: utf-8 -*- """ Created on Mon Apr 8 15:01:43 2024 @author: LDDN """ import math import pandas as pd import numpy as np import matplotlib.pyplot as plt from matplotlib.pyplot import MultipleLocator#设定固定刻度 scada_10min = pd.read_csv(r'E:\BaiduNetdiskDownload\test\min_scada_LuoTuoGou\72\82.csv',encoding="utf-8") #.value是将单元格 turbine_info = pd.read_csv(r'E:\BaiduNetdiskDownload\test\min_scada_LuoTuoGou\72\info.csv') #.value是将单元格 PRated = turbine_info.loc[:,["额定功率"]] #2000 PRated = PRated.values VCutOut = turbine_info.loc[:,["切出风速"]] #25 VCutOut = VCutOut.values VCutIn = turbine_info.loc[:,["切入风速"]] #3 VCutIn = VCutIn.values VRated = turbine_info.loc[:,["额定风速"]] #10 VRated = VRated.values time_stamp = scada_10min.loc[:,['时间']] #dataframe active_power = scada_10min.loc[:,['变频器电网侧有功功率']] wind_speed = scada_10min.loc[:,['风速']] LM = pd.concat([time_stamp,active_power,wind_speed],axis=1) #dataframe Labeled_March809 = LM APower = Labeled_March809["变频器电网侧有功功率"] #series读入有功功率 WSpeed = Labeled_March809["风速"] #读入风速 maxP=np.max(APower) intervalP=25 #ceil(PRated*0.01)#功率分区间隔为额定功率的1% intervalwindspeed=0.25 #风速分区间隔0.25m/s #初始化 PNum = 0 TopP = 0 # 根据条件计算PNum和TopP if maxP >= PRated: PNum = math.floor(maxP / intervalP) + 1 TopP = math.floor((maxP - PRated) / intervalP) + 1 else: PNum = math.floor(PRated / intervalP) TopP = 0 VNum = math.ceil(VCutOut / intervalwindspeed) SM1 = Labeled_March809.shape AA1 = SM1[0] lab = [[0] for _ in range(AA1)] lab = pd.DataFrame(lab,columns=['lab']) Labeled_March809 = pd.concat([Labeled_March809,lab],axis=1) #在tpv后加一列标签列 Labeled_March809 = Labeled_March809.values SM = Labeled_March809.shape #(52561,4) AA = SM[0] #存储功率大于0的运行数据 #标识功率为0的点，标识-1 DzMarch809_0 = np.zeros((AA, 3)) # 初始化数组来存储功率大于零的运行数据 nCounter1 = 0 Point_line = np.zeros(AA, dtype=int) #考虑到很多功率小于10的数据存在，将<10的功率视为0 for i in range(AA): if (APower[i] > 10) & (WSpeed[i] > 0): nCounter1 += 1 #共有nCounter1个功率大于0的正常数据 DzMarch809_0[nCounter1-1, 0] = WSpeed[i] DzMarch809_0[nCounter1-1, 1] = APower[i] Point_line[nCounter1-1] = i+1 # 记录nCounter1记下的数据在原始数据中的位置 if APower[i] <= 10: Labeled_March809[i,SM[1]-1] = -1 # 功率为0标识为-1 array类型 # 截取DzMarch809_0中实际存储的数据 其他全为0 DzMarch809 = DzMarch809_0[:nCounter1, :] #统计各网格落入的散点个数 XBoxNumber = np.ones((PNum, VNum),dtype=int) #(86 100) nWhichP = 0 nWhichV = 0 # 循环遍历DzMarch809中的有效数据 for i in range(nCounter1): # 查找功率所在的区间 for m in range(1, PNum + 1): # 注意Python的range是左闭右开的，所以需要+1 if (DzMarch809[i,1] > (m - 1) * intervalP) and (DzMarch809[i,1] <= m * intervalP): nWhichP = m break # 查找风速所在的区间 for n in range(1, VNum + 1): # 同样需要+1 if (DzMarch809[i, 0] > (n - 1)*intervalwindspeed) and (DzMarch809[i, 0] <= n*intervalwindspeed): nWhichV = n break # 如果功率和风速都在有效区间内，增加对应网格的计数 if (nWhichP > 0) and (nWhichV > 0): XBoxNumber[nWhichP - 1, nWhichV - 1] += 1 # 注意Python的索引是从0开始的，所以需要减1 # XBoxNumber现在包含了每个网格的计数[PNum行, VNum列] for m in range(1,PNum+1): for n in range(1,VNum+1): XBoxNumber[m-1,n-1] = XBoxNumber[m-1,n-1] - 1 #在功率方向将网格内散点绝对个数转换为相对百分比，备用 PBoxPercent = np.zeros((PNum, VNum),dtype = float) #(86 100) #计算后会出现浮点型，所以不能定义int类型 PBinSum = np.zeros((PNum,1),dtype=int) for i in range(1,PNum+1): for m in range(1,VNum+1): PBinSum[i-1] = PBinSum[i-1] + XBoxNumber[i-1,m-1] for m in range(1,VNum+1): if PBinSum[i-1]>0: PBoxPercent[i-1,m-1] = (XBoxNumber[i-1,m-1] / PBinSum[i-1])*100 #在风速方向将网格内散点绝对个数转换为相对百分比，备用 VBoxPercent = np.zeros((PNum, VNum)) #(86 100) #计算后会出现浮点型，所以不能定义int类型 VBinSum = np.zeros((VNum,1),dtype=int) for i in range(1,VNum+1): for m in range(1,PNum+1): VBinSum[i-1] = VBinSum[i-1] + XBoxNumber[m-1,i-1] for m in range(1,PNum+1): if VBinSum[i-1]>0: VBoxPercent[m-1,i-1] = (XBoxNumber[m-1,i-1] / VBinSum[i-1])*100 # VBoxPercent PBoxPercent 左上-右下 # 将数据颠倒一下 左下-右上 第一行换为倒数第一行 方便可视化 InvXBoxNumber = np.zeros((PNum,VNum),dtype = int) InvPBoxPercent = np.zeros((PNum,VNum),dtype = float) InvVBoxPercent = np.zeros((PNum,VNum),dtype = float) for m in range(1,PNum+1): for n in range(1,VNum+1): InvXBoxNumber[m-1,n-1] = XBoxNumber[PNum-(m-1)-1,n-1] InvPBoxPercent[m-1,n-1] = PBoxPercent[PNum-(m-1)-1,n-1] InvVBoxPercent[m-1,n-1] = VBoxPercent[PNum-(m-1)-1,n-1] #以水平功率带方向为准，分析每个水平功率带中，功率主带中心，即找百分比最大的网格位置。 PBoxMaxIndex = np.zeros((PNum,1),dtype = int) #水平功率带最大网格位置索引 PBoxMaxP = np.zeros((PNum,1),dtype = float) #水平功率带最大网格百分比 for m in range(1,PNum+1): PBoxMaxIndex[m-1] = np.argmax(PBoxPercent[m-1, :]) #argmax返回最大值的索引 PBoxMaxP[m-1] = np.max(PBoxPercent[m-1, :]) #以垂直风速方向为准，分析每个垂直风速带中，功率主带中心，即找百分比最大的网格位置。 VBoxMaxIndex = np.zeros((VNum,1),dtype = int) VBoxMaxV = np.zeros((VNum,1),dtype = float) for m in range(1,VNum+1): VBoxMaxIndex[m-1] = np.argmax(VBoxPercent[:, m-1]) VBoxMaxV[m-1] = np.max(VBoxPercent[:, m-1]) #切入风速特殊处理，如果切入风速过于偏右，向左拉回 if PBoxMaxIndex[0]>14: #第一个值对应的是风速最小处 即切入风速 PBoxMaxIndex[0] = 9 #以水平功率带方向为基准，进行分析 DotDense = np.zeros(PNum) #每一水平功率带的功率主带包含的网格数 DotDenseLeftRight = np.zeros((PNum,2)) #存储每一水平功率带的功率主带以最大网格为中心，向左，向右扩展的网格数 DotValve = 90 #从中心向左右对称扩展网格的散点百分比和的阈值。 PDotDenseSum = 0 for i in range(PNum - TopP): # 从最下层水平功率带开始，向上分析到特定的功率带 PDotDenseSum = PBoxMaxP[i] # 以中心最大水平功率带为基准，向左向右对称扩展网格，累加各网格散点百分比 iSpreadRight = 1 iSpreadLeft = 1 while PDotDenseSum < DotValve: if (PBoxMaxIndex[i] + iSpreadRight) < VNum-1-1: PDotDenseSum += PBoxPercent[i, PBoxMaxIndex[i] + iSpreadRight] # 向右侧扩展 iSpreadRight += 1 else: break if (PBoxMaxIndex[i] - iSpreadLeft) > 0: PDotDenseSum += PBoxPercent[i, PBoxMaxIndex[i] - iSpreadLeft] # 向左侧扩展 iSpreadLeft += 1 else: break iSpreadRight = iSpreadRight-1 iSpreadLeft = iSpreadLeft-1 #向左右扩展完毕 DotDenseLeftRight[i, 0] = iSpreadLeft # 左 DotDenseLeftRight[i, 1] = iSpreadRight # 右 DotDense[i] = iSpreadLeft + iSpreadRight + 1 # 记录向左向右扩展的个数及每个功率仓内网格的个数 # 此时DotDense和DotDenseLeftRight数组已经包含了所需信息 #各行功率主带右侧宽度的中位数最具有代表性（因为先右后左） DotDenseWidthLeft = np.zeros((PNum-TopP)) for i in range(PNum-TopP): DotDenseWidthLeft[i] = DotDenseLeftRight[i,1] #DotDenseLeftRight[i,1]:向右延伸个数 MainBandRight = np.median(DotDenseWidthLeft) #计算中位数 # 初始化变量 PowerLimit = np.zeros(PNum, dtype=int) # 各水平功率带是否为限功率标识，1：是；0：不是 WidthAverage = 0 # 功率主带右侧平均宽度 WidthAverage_L = 0 # 功率主带左侧平均宽度 WidthVar = 0 # 功率主带方差（此变量在提供的代码中并未使用） PowerLimitValve = 6 # 限功率主带判别阈值 N_Pcount = 20 # 阈值 nCounterLimit = 0 # 限功率的个数 nCounter = 0 # 正常水平功率带的个数 # 循环遍历水平功率带，从第1个到第PNum-TopP个 for i in range(PNum - TopP): # 如果向右扩展网格数大于阈值，且该水平功率带点总数大于20，则标记为限功率带 if (DotDenseLeftRight[i, 1] > PowerLimitValve) and (PBinSum[i] > N_Pcount): PowerLimit[i] = 1 nCounterLimit += 1 #限功率的个数 # 如果向右扩展网格数小于等于阈值，则累加右侧宽度（左侧宽度在代码中似乎有误） if DotDenseLeftRight[i, 1] <= PowerLimitValve: WidthAverage += DotDenseLeftRight[i, 1] # 统计正常水平功率带右侧宽度 WidthAverage_L += DotDenseLeftRight[i,1] #统计正常水平功率带左侧宽度 nCounter += 1 # 计算平均宽度 WidthAverage /= nCounter if nCounter > 0 else 1 # 避免除以0的情况 WidthAverage_L /= nCounter if nCounter > 0 else 1 #计算正常（即非限功率）水平功率带的功率主带宽度的方差，以此来反映从下到上宽度是否一致 WidthVar = 0 # 功率主带宽度的方差 for i in range(PNum - TopP): # 如果向右扩展网格数小于等于阈值，则计算当前宽度与平均宽度的差值平方 if DotDenseLeftRight[i, 1] <= PowerLimitValve: WidthVar += (DotDenseLeftRight[i, 1] - WidthAverage) ** 2 # 计算方差（注意：除以nCounter-1是为了得到样本方差） WidthVar = np.sqrt(WidthVar / (nCounter - 1) if nCounter > 1 else 0) # 避免除以0的情况 #各水平功率带，功率主带的风速范围，右侧扩展网格数*2*0.25 PowerBandWidth = WidthAverage*intervalwindspeed+WidthAverage_L*intervalwindspeed # 对限负荷水平功率带的最大网格进行修正 for i in range(1, PNum - TopP+1): if (PowerLimit[i] == 1) and (abs(PBoxMaxIndex[i] - PBoxMaxIndex[i - 1]) > 5): PBoxMaxIndex[i] = PBoxMaxIndex[i - 1] + 1 # 输出各层功率主带的左右边界网格索引 DotDenseInverse = np.flipud(DotDenseLeftRight) # 上下翻转数组以得到反向顺序 # 计算功率主带的左右边界 CurveWidthR = np.ceil(WidthAverage) + 2 # 功率主带的右边界 + 2 CurveWidthL = np.ceil(WidthAverage_L) + 2 # 功率主带的左边界 + 2 # 网格是否为限功率网格的标识数组 BBoxLimit = np.zeros((PNum, VNum), dtype=int) # 标记限功率网格 for i in range(2, PNum - TopP): if PowerLimit[i] == 1: BBoxLimit[i, int(PBoxMaxIndex[i] + CurveWidthR + 1):VNum] = 1 # 初始化数据异常需要剔除的网格标识数组 BBoxRemove = np.zeros((PNum, VNum), dtype=int) # 标记需要剔除的网格 for m in range(PNum - TopP): for n in range(int(PBoxMaxIndex[m]) + int(CurveWidthR), VNum): # 注意Python中的索引从0开始，因此需要减去1 BBoxRemove[m, n] = 1 # 功率主带左侧的超发网格，从最大索引向左直到第一个网格 for n in range(int(PBoxMaxIndex[m]) - int(CurveWidthL)+1, 0, -1): # 使用range的步长参数来实现从右向左的迭代 BBoxRemove[m, n-1] = 2 # 注意Python中的索引从0开始，因此需要减去1 # 初始化变量 CurveTop = np.zeros((2, 1), dtype=int) CurveTopValve = 1 # 网格的百分比阈值 BTopFind = 0 mm = 0 #确定功率主带的左上拐点，即额定风速位置的网格索引 CurveTop = np.zeros((2, 1), dtype=int) CurveTopValve = 1 # 网格的百分比阈值 BTopFind = 0 mm = 0 for m in range(PNum - TopP, 0, -1): # 注意Python的range是左闭右开区间，所以这里从PNum-TopP开始到1（不包括0） for n in range(int(np.floor(int(VCutIn) / intervalwindspeed)), VNum - 1): # 使用floor函数来向下取整 if (VBoxPercent[m, n - 1] < VBoxPercent[m, n]) and (VBoxPercent[m, n] <= VBoxPercent[m, n + 1]) and (XBoxNumber[m, n] >= 3): CurveTop[0] = m CurveTop[1] = n #[第80个,第40个] BTopFind = 1 mm = m # mm是拐点所在功率仓，对应其index break # 找到后退出内层循环 if BTopFind == 1: break # 找到后退出外层循环 IsolateValve = 3 #功率主带右侧孤立点占比功率仓阈值 3% # 遍历功率仓和网格 for m in range(PNum - TopP): for n in range(int(PBoxMaxIndex[m]) + int(CurveWidthR), VNum): # 检查PBoxPercent是否小于阈值，如果是，则标记BBoxRemove为1 if PBoxPercent[m, n] < IsolateValve: BBoxRemove[m, n] = 1 #功率主带顶部宽度 CurveWidthT = np.floor((maxP - PRated) / intervalP) + 1 # 标记额定功率以上的超发点(PNum-PTop之间) for m in range(PNum - TopP, PNum): for n in range(VNum): BBoxRemove[m, n] = 3 # 标记功率主带拐点左侧的欠发网格 for m in range(mm-1, PNum - TopP): for n in range(int(CurveTop[1]) - 2): BBoxRemove[m, n] = 2 # BBoxRemove数组现在包含了根据条件标记的超发点和欠发网格的信息 #以网格的标识，决定该网格内数据的标识。 # DzMarch809Sel数组现在包含了每个数据点的标识 DzMarch809Sel = np.zeros(nCounter1, dtype=int) # 初始化标识数组 nWhichP = 0 nWhichV = 0 nBadA = 0 for i in range(nCounter1): for m in range( PNum ): if (DzMarch809[i, 1] > m * intervalP) and (DzMarch809[i, 1] <= (m+1) * intervalP): nWhichP = m #m记录的是index break for n in range( VNum ): # 注意Python的range是左闭右开区间，所以这里到VNum+1 if DzMarch809[i, 0] > ((n+1) * intervalwindspeed - intervalwindspeed/2) and DzMarch809[i, 0] <= ((n+1) * intervalwindspeed + intervalwindspeed / 2): nWhichV = n #index break if nWhichP >= 0 and nWhichV >= 0: if BBoxRemove[nWhichP, nWhichV] == 1: DzMarch809Sel[i] = 1 nBadA += 1 elif BBoxRemove[nWhichP, nWhichV] == 2: DzMarch809Sel[i] = 2 elif BBoxRemove[nWhichP , nWhichV] == 3: DzMarch809Sel[i] = 0 # 额定风速以上的超发功率点认为是正常点，不再标识 # DzMarch809Sel数组现在包含了每个数据点的标识 ##############################滑动窗口方法 # 存储限负荷数据 PVLimit = np.zeros((nCounter1, 3)) #存储限负荷数据 %第3列用于存储限电的点所在的行数 nLimitTotal = 0 nWindowLength = 6 #滑动窗口长度设置为6 LimitWindow = np.zeros(nWindowLength) #滑动窗口空列表 UpLimit = 0 #上限 LowLimit = 0 #下限 PowerStd = 30 # 功率波动方差 nWindowNum = np.floor(nCounter1/nWindowLength) #6587 PowerLimitUp = PRated - 100 PowerLimitLow = 100 # 循环遍历每个窗口 for i in range(int(nWindowNum)): start_idx = i * nWindowLength end_idx = start_idx + nWindowLength LimitWindow = DzMarch809[start_idx:end_idx, 1] # 提取当前窗口的数据 # 检查窗口内所有数据是否在功率范围内 bAllInAreas = np.all(LimitWindow >= PowerLimitLow) and np.all(LimitWindow <= PowerLimitUp) if not bAllInAreas: continue # 计算方差上下限 UpLimit = LimitWindow[0] + PowerStd LowLimit = LimitWindow[0] - PowerStd # 检查窗口内数据是否在方差范围内 bAllInUpLow = np.all(LimitWindow >= LowLimit) and np.all(LimitWindow <= UpLimit) if bAllInUpLow: # 标识窗口内的数据为限负荷数据 DzMarch809Sel[start_idx:end_idx] = 4 # 存储限负荷数据 for j in range(nWindowLength): PVLimit[nLimitTotal, :2] = DzMarch809[start_idx + j, :2] PVLimit[nLimitTotal, 2] = Point_line[start_idx + j] # 对数据进行标识 nLimitTotal += 1 # PVLimit现在包含了限负荷数据，nLimitTotal是限负荷数据的总数 #将功率滑动窗口主带平滑化 # 初始化锯齿平滑的计数器 nSmooth = 0 # 遍历除了最后 TopP+1 个元素之外的所有 PBoxMaxIndex 元素 for i in range(PNum - TopP - 1): PVLeftDown = np.zeros(2) PVRightUp = np.zeros(2) # 检查当前与下一个 PBoxMaxIndex 之间的距离是否大于等于1 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 # 遍历 DzMarch809 数组 for m in range(nCounter1): # 检查当前点是否在锯齿区域内 if (DzMarch809[m, 0] > PVLeftDown[0]) and (DzMarch809[m, 0] < PVRightUp[0]) and (DzMarch809[m, 1] > PVLeftDown[1]) and (DzMarch809[m, 1] < PVRightUp[1]): # 检查斜率是否大于对角连线 if (DzMarch809[m, 1] - PVLeftDown[1]) / (DzMarch809[m, 0] - PVLeftDown[0]) > (PVRightUp[1] - PVLeftDown[1]) / (PVRightUp[0] - PVLeftDown[0]): # 如果在锯齿左上三角形中，则选中并增加锯齿平滑计数器 DzMarch809Sel[m] = 0 nSmooth += 1 # DzMarch809Sel 数组现在＋包含了锯齿平滑的选择结果，nSmooth 是选中的点数 ###################################存储数据 # 存储好点 nCounterPV = 0 # 初始化计数器 PVDot = np.zeros((nCounter1, 3)) # 初始化存储好点的数组 nCounter1是p>0的数 for i in range(nCounter1): if DzMarch809Sel[i] == 0: nCounterPV += 1 PVDot[nCounterPV-1, :2] = DzMarch809[i, :2] PVDot[nCounterPV-1, 2] = Point_line[i] # 好点 Point_line记录nCounter1在原始数据中的位置 nCounterVP = nCounterPV # 对所有数据中的好点进行标注 for i in range(nCounterVP): Labeled_March809[int(PVDot[i, 2] - 1), (SM[1]-1)] = 1 # 注意Python的索引从0开始，并且需要转换为整数索引 # 存储坏点 nCounterBad = 0 # 初始化计数器 PVBad = np.zeros((nCounter1, 3)) # 初始化存储坏点的数组 for i in range(nCounter1): if DzMarch809Sel[i] in [1, 2, 3]: nCounterBad += 1 PVBad[nCounterBad-1, :2] = DzMarch809[i, :2] PVBad[nCounterBad-1, 2] = Point_line[i] # 对所有数据中的坏点进行标注 for i in range(nCounterBad): Labeled_March809[int(PVBad[i, 2] - 1),(SM[1]-1)] = 5 # 坏点标识 # 对所有数据中的限电点进行标注 for i in range(nLimitTotal): Labeled_March809[int(PVLimit[i, 2] - 1),(SM[1]-1)] = 4 # 限电点标识 # 对所有的数据点进行标注 # Labeled_March809是array，提取所第四列的值保存为dataframe A = Labeled_March809[:,3] A=pd.DataFrame(A,columns=['lab']) mergedTable = pd.concat([scada_10min,A],axis=1)#合并dataframe D = mergedTable[mergedTable['lab'] == 1]#选择为1的行 ws = D["风速"].values #array ap = D["变频器电网侧有功功率"] # fig=plt.figure(figsize=(10,6),dpi=500) #figsize是图形大小，dpi像素 fig=plt.figure() #figsize是图形大小，dpi像素 plt.scatter(ws,ap,s=1,c='black',marker='.') #'.'比'o'要更小 # plt.scatter(x2,y2,s=10,c='b',marker='.',label='5.8-6.5建模噪声点') x_major_locator=MultipleLocator(5) y_major_locator=MultipleLocator(500) ax=plt.gca() ax.xaxis.set_major_locator(x_major_locator) ax.yaxis.set_major_locator(y_major_locator) plt.xlim((0,30)) plt.ylim((0,2200)) plt.tick_params(labelsize=8) # plt.grid(c='dimgray',alpha=0.2) plt.xlabel("V/(m$·$s$^{-1}$)",fontsize=8) plt.ylabel("P/kW",fontsize=8) # plt.savefig(r'D:\赵雅丽\研究生学习资料\学习资料\劣化度健康度\spyder\大论文\图\风速-功率.jpg',bbox_inches='tight') plt.show()