ClassIdentifier_2.py 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371
  1. import os
  2. import numpy as np
  3. from pandas import DataFrame
  4. from service.plt_service import get_base_wind_and_power
  5. from utils.draw.draw_file import scatter
  6. from utils.file.trans_methods import read_file_to_df
  7. class ClassIdentifier(object):
  8. def __init__(self, wind_turbine_number, file_path: str = None, origin_df: DataFrame = None, index='time_stamp',
  9. wind_velocity='wind_velocity',
  10. active_power='active_power'):
  11. """
  12. :param wind_turbine_number: The wind turbine number.
  13. :param file_path: The file path of the input data.
  14. :param origin_df: The pandas DataFrame containing the input data.
  15. :param index: 索引字段
  16. :param wind_velocity: 风速字段
  17. :param active_power: 有功功率字段
  18. """
  19. self.wind_turbine_number = wind_turbine_number
  20. self.index = index
  21. self.wind_velocity = wind_velocity
  22. self.active_power = active_power
  23. self.rated_wind_speed = 'rated_wind_speed'
  24. self.rated_capacity = 'rated_capacity'
  25. if file_path is None and origin_df is None:
  26. raise ValueError("Either file_path or origin_df should be provided.")
  27. if file_path:
  28. self.df = read_file_to_df(file_path)
  29. else:
  30. self.df = origin_df
  31. self.df = self.df.set_index(keys=self.index)
  32. def identifier(self):
  33. # 风速 和 有功功率 df
  34. wind_and_power_df = self.df[[self.wind_velocity, self.active_power]]
  35. wind_and_power_df.reset_index(inplace=True)
  36. wind_and_power_df_count = wind_and_power_df.shape[0]
  37. PowerMax = wind_and_power_df[self.active_power].max()
  38. PowerRated = np.ceil(PowerMax / 100) * 100
  39. PRated = 1500 # 额定功率1500kw,可改为2000kw
  40. VCutOut = 25
  41. # 网格法确定风速风向分区数量,功率方向分区数量,
  42. # PNum = (PRated+100)/25 #功率分区间隔25kW
  43. PNum = int(np.ceil(PowerRated / 25)) # 功率分区间隔25kW
  44. VNum = int(np.ceil(VCutOut / 0.25)) # 风速分区间隔0.25m/s
  45. # 存储功率大于零的运行数据
  46. DzMarch809 = np.zeros([wind_and_power_df_count, 2], dtype=float)
  47. nCounter1 = 0
  48. for i in range(wind_and_power_df_count):
  49. if wind_and_power_df.loc[i, self.active_power] > 0:
  50. DzMarch809[nCounter1, 0] = wind_and_power_df.loc[i, self.wind_velocity]
  51. DzMarch809[nCounter1, 1] = wind_and_power_df.loc[i, self.active_power]
  52. nCounter1 = nCounter1 + 1
  53. # 统计各网格落入的散点个数
  54. if VNum == 1:
  55. XBoxNumber = np.ones([PNum], dtype=int)
  56. else:
  57. XBoxNumber = np.ones([PNum, VNum], dtype=int)
  58. nWhichP = -1
  59. nWhichV = -1
  60. for i in range(nCounter1):
  61. for m in range(PNum):
  62. if m * 25 < DzMarch809[i, 1] <= (m + 1) * 25:
  63. nWhichP = m
  64. break
  65. for n in range(VNum):
  66. if ((n + 1) * 0.25 - 0.125) < DzMarch809[i, 0] <= ((n + 1) * 0.25 + 0.125):
  67. nWhichV = n
  68. break
  69. if nWhichP > -1 and nWhichV > -1:
  70. XBoxNumber[nWhichP, nWhichV] = XBoxNumber[nWhichP, nWhichV] + 1
  71. for m in range(PNum):
  72. for n in range(VNum):
  73. XBoxNumber[m, n] = XBoxNumber[m, n] - 1
  74. # 在功率方向将网格内散点绝对个数转换为相对百分比,备用
  75. PBoxPercent = np.zeros([PNum, VNum], dtype=float)
  76. PBinSum = np.zeros(PNum, dtype=int)
  77. for i in range(PNum):
  78. for m in range(VNum):
  79. PBinSum[i] = PBinSum[i] + XBoxNumber[i, m]
  80. for m in range(VNum):
  81. if PBinSum[i] > 0:
  82. PBoxPercent[i, m] = XBoxNumber[i, m] / PBinSum[i] * 100
  83. # 在风速方向将网格内散点绝对个数转换为相对百分比,备用
  84. VBoxPercent = np.zeros([PNum, VNum], dtype=float)
  85. VBinSum = np.zeros(VNum, dtype=int)
  86. for i in range(VNum):
  87. for m in range(PNum):
  88. VBinSum[i] = VBinSum[i] + XBoxNumber[m, i]
  89. for m in range(PNum):
  90. if VBinSum[i] > 0:
  91. VBoxPercent[m, i] = XBoxNumber[m, i] / VBinSum[i] * 100
  92. # 以水平功率带方向为准,分析每个水平功率带中,功率主带中心,即找百分比最大的网格位置。
  93. PBoxMaxIndex = np.zeros(PNum, dtype=int) # 水平功率带最大网格位置索引
  94. PBoxMaxP = np.zeros(PNum, dtype=int) # 水平功率带最大网格百分比
  95. for m in range(PNum):
  96. # 确定每一水平功率带的最大网格位置索引即百分比值
  97. PBoxMaxP[m], PBoxMaxIndex[m] = PBoxPercent[m, :].max(), PBoxPercent[m, :].argmax()
  98. # 以垂直风速方向为准,分析每个垂直风速带中,功率主带中心,即找百分比最大的网格位置。
  99. VBoxMaxIndex = np.zeros(VNum, dtype=int)
  100. VBoxMaxV = np.zeros(VNum, dtype=int)
  101. for m in range(VNum):
  102. [VBoxMaxV[m], VBoxMaxIndex[m]] = VBoxPercent[:, m].max(), VBoxPercent[:, m].argmax()
  103. # 切入风速特殊处理,如果切入风速过于偏右,向左拉回
  104. if PBoxMaxIndex[0] > 14:
  105. PBoxMaxIndex[0] = 9
  106. # 以水平功率带方向为基准,进行分析
  107. DotDense = np.zeros(PNum, dtype=int) # 每一水平功率带的功率主带包含的网格数
  108. DotDenseLeftRight = np.zeros([PNum, 2], dtype=int) # 存储每一水平功率带的功率主带以最大网格为中心,向向左,向右扩展的网格数
  109. DotValve = 90 # 从中心向左右对称扩展网格的散点百分比和的阈值。
  110. for i in range(PNum - 6): # 从最下层水平功率带1开始,向上到第PNum-6个水平功率带(额定功率一下水平功率带),逐一分析
  111. PDotDenseSum = PBoxMaxP[i] # 以中心最大水平功率带为基准,向左向右对称扩展网格,累加各网格散点百分比
  112. iSpreadRight = 1
  113. iSpreadLeft = 1
  114. while PDotDenseSum < DotValve:
  115. if (PBoxMaxIndex[i] + iSpreadRight) < VNum - 1:
  116. PDotDenseSum = PDotDenseSum + PBoxPercent[i, PBoxMaxIndex[i] + iSpreadRight] # 向右侧扩展
  117. iSpreadRight = iSpreadRight + 1
  118. if (PBoxMaxIndex[i] + iSpreadRight) > VNum - 1:
  119. break
  120. if (PBoxMaxIndex[i] - iSpreadLeft) > 0:
  121. PDotDenseSum = PDotDenseSum + PBoxPercent[i, PBoxMaxIndex[i] - iSpreadLeft] # 向左侧扩展
  122. iSpreadLeft = iSpreadLeft + 1
  123. if (PBoxMaxIndex[i] - iSpreadLeft) <= 0:
  124. break
  125. iSpreadRight = iSpreadRight - 1
  126. iSpreadLeft = iSpreadLeft - 1
  127. # 向左右对称扩展完毕
  128. DotDenseLeftRight[i, 0] = iSpreadLeft
  129. DotDenseLeftRight[i, 1] = iSpreadRight
  130. DotDense[i] = iSpreadLeft + iSpreadRight + 1
  131. # 各行功率主带右侧宽度的中位数最具有代表性
  132. DotDenseWidthLeft = np.zeros([PNum - 6, 1], dtype=int)
  133. for i in range(PNum - 6):
  134. DotDenseWidthLeft[i] = DotDenseLeftRight[i, 1]
  135. MainBandRight = np.median(DotDenseWidthLeft)
  136. # 散点向右显著延展分布的水平功率带为限功率水平带
  137. PowerLimit = np.zeros([PNum, 1], dtype=int) # 各水平功率带是否为限功率标识,==1:是;==0:不是
  138. WidthAverage = 0 # 功率主带平均宽度
  139. WidthVar = 0 # 功率主带方差
  140. # PowerLimitValve = 6 #限功率主带判别阈值
  141. PowerLimitValve = np.ceil(MainBandRight) + 3 # 限功率主带判别阈值
  142. nCounterLimit = 0
  143. nCounter = 0
  144. for i in range(PNum - 6):
  145. if DotDenseLeftRight[i, 1] > PowerLimitValve and PBinSum[i] > 20: # 如果向右扩展网格数大于阈值,且该水平功率带点总数>20,是
  146. PowerLimit[i] = 1
  147. nCounterLimit = nCounterLimit + 1
  148. if DotDenseLeftRight[i, 1] <= PowerLimitValve:
  149. WidthAverage = WidthAverage + DotDenseLeftRight[i, 1] # 统计正常水平功率带右侧宽度
  150. nCounter = nCounter + 1
  151. WidthAverage = WidthAverage / nCounter # 功率主带平均宽度
  152. # 各水平功率带的功率主带宽度的方差,反映从下到上宽度是否一致,或是否下宽上窄等异常情况
  153. for i in range(PNum - 6):
  154. if DotDenseLeftRight[i, 1] <= PowerLimitValve:
  155. WidthVar = WidthVar + (DotDenseLeftRight[i, 1] - WidthAverage) * (
  156. DotDenseLeftRight[i, 1] - WidthAverage)
  157. # 对限负荷水平功率带的最大网格较下面相邻层显著偏右,拉回
  158. for i in range(1, PNum - 6):
  159. if PowerLimit[i] == 1 and abs(PBoxMaxIndex[i] - PBoxMaxIndex[i - 1]) > 5:
  160. PBoxMaxIndex[i] = PBoxMaxIndex[i - 1] + 1
  161. # 输出各层功率主带的左右边界网格索引
  162. DotDenseInverse = np.zeros([PNum, 2], dtype=int)
  163. for i in range(PNum):
  164. DotDenseInverse[i, :] = DotDenseLeftRight[PNum - i - 1, :]
  165. # 功率主带的右边界
  166. CurveWidthR = int(np.ceil(WidthAverage) + 2)
  167. # CurveWidthL = 6 #功率主带的左边界
  168. CurveWidthL = CurveWidthR
  169. BBoxLimit = np.zeros([PNum, VNum], dtype=int) # 网格是否为限功率网格的标识,如果为限功率水平功率带,从功率主带右侧边缘向右的网格为限功率网格
  170. for i in range(2, PNum - 6):
  171. if PowerLimit[i] == 1:
  172. for j in range(PBoxMaxIndex[i] + CurveWidthR, VNum):
  173. BBoxLimit[i, j] = 1
  174. BBoxRemove = np.zeros([PNum, VNum], dtype=int) # 数据异常需要剔除的网格标识,标识==1:功率主带右侧的欠发网格;==2:功率主带左侧的超发网格
  175. for m in range(PNum - 6):
  176. for n in range(PBoxMaxIndex[m] + CurveWidthR, VNum):
  177. BBoxRemove[m, n] = 1
  178. for n in range(PBoxMaxIndex[m] - CurveWidthL, -1, -1):
  179. BBoxRemove[m, n] = 2
  180. # 确定功率主带的左上拐点,即额定风速位置的网格索引
  181. CurveTop = np.zeros(2, dtype=int)
  182. CurveTopValve = 3 # 网格的百分比阈值
  183. BTopFind = 0
  184. for m in range(PNum - 4 - 1, -1, -1):
  185. for n in range(VNum):
  186. if VBoxPercent[m, n] > CurveTopValve and XBoxNumber[m, n] >= 10: # 如左上角网格的百分比和散点个数大于阈值。
  187. CurveTop[0] = m
  188. CurveTop[1] = n
  189. BTopFind = 1
  190. break
  191. if BTopFind == 1:
  192. break
  193. IsolateValve = 3
  194. for m in range(PNum - 6):
  195. for n in range(PBoxMaxIndex[m] + CurveWidthR, VNum):
  196. if PBoxPercent[m, n] < IsolateValve:
  197. BBoxRemove[m, n] = 1
  198. # 功率主带顶部宽度
  199. CurveWidthT = 2
  200. for m in range(PNum - CurveWidthT - 1, PNum):
  201. for n in range(VNum):
  202. BBoxRemove[m, n] = 3 # 网格为额定功率以上的超发点
  203. # 功率主带拐点左侧的欠发网格标识
  204. for m in range(PNum - 5 - 1, PNum):
  205. for n in range(CurveTop[1] - 1):
  206. BBoxRemove[m, n] = 2
  207. # 以网格的标识,决定该网格内数据的标识。Dzwind_and_power_dfSel功率非零数据的标识位。散点在哪个网格,此网格的标识即为该点的标识
  208. Dzwind_and_power_dfSel = np.zeros(nCounter1, dtype=int) # -1:停机 0:好点 1:欠发功率点;2:超发功率点;3:额定风速以上的超发功率点 4: 限电
  209. nWhichP = -1
  210. nWhichV = -1
  211. nBadA = 0
  212. for i in range(nCounter1):
  213. for m in range(PNum):
  214. if m * 25 < DzMarch809[i, 1] <= (m + 1) * 25:
  215. nWhichP = m
  216. break
  217. for n in range(VNum):
  218. if ((n + 1) * 0.25 - 0.125) < DzMarch809[i, 0] <= ((n + 1) * 0.25 + 0.125):
  219. nWhichV = n
  220. break
  221. if nWhichP > -1 and nWhichV > -1:
  222. if BBoxRemove[nWhichP, nWhichV] == 1:
  223. Dzwind_and_power_dfSel[i] = 1
  224. nBadA = nBadA + 1
  225. if BBoxRemove[nWhichP, nWhichV] == 2:
  226. Dzwind_and_power_dfSel[i] = 2
  227. if BBoxRemove[nWhichP, nWhichV] == 3:
  228. Dzwind_and_power_dfSel[i] = 0 # 3 # 额定风速以上的超发功率点认为是正常点,不再标识。
  229. # 限负荷数据标识方法2:把数据切割为若干个窗口。对每一窗口,以第一个点为基准,连续nWindowLength个数据的功率在方差范围内,呈现显著水平分布的点
  230. nWindowLength = 3
  231. LimitWindow = np.zeros(nWindowLength, dtype=float)
  232. PowerStd = 15 # 功率波动方差
  233. nWindowNum = int(np.floor(nCounter1 / nWindowLength))
  234. PowerLimitUp = PRated - 300
  235. PowerLimitLow = 200
  236. for i in range(nWindowNum):
  237. for j in range(nWindowLength):
  238. LimitWindow[j] = DzMarch809[i * nWindowLength + j, 1]
  239. bAllInAreas = 1
  240. for j in range(nWindowLength):
  241. if LimitWindow[j] < PowerLimitLow or LimitWindow[j] > PowerLimitUp:
  242. bAllInAreas = 0
  243. if bAllInAreas == 0:
  244. continue
  245. UpLimit = LimitWindow[0] + PowerStd
  246. LowLimit = LimitWindow[0] - PowerStd
  247. bAllInUpLow = 1
  248. for j in range(1, nWindowLength):
  249. if LimitWindow[j] < LowLimit or LimitWindow[j] > UpLimit:
  250. bAllInUpLow = 0
  251. if bAllInUpLow == 1:
  252. for j in range(nWindowLength):
  253. Dzwind_and_power_dfSel[i * nWindowLength + j] = 4 # 标识窗口内的数据为限负荷数据
  254. nSmooth = 0
  255. for i in range(PNum - 6):
  256. PVLeftDown = np.zeros(2, dtype=float)
  257. PVRightUp = np.zeros(2, dtype=float)
  258. if (PBoxMaxIndex[i + 1] - PBoxMaxIndex[i]) >= 1:
  259. PVLeftDown[0] = (PBoxMaxIndex[i] + 1 + CurveWidthR) * 0.25 - 0.125
  260. PVLeftDown[1] = i * 25
  261. PVRightUp[0] = (PBoxMaxIndex[i + 1] + 1 + CurveWidthR) * 0.25 - 0.125
  262. PVRightUp[1] = (i + 1) * 25
  263. for m in range(nCounter1):
  264. if DzMarch809[m, 0] > PVLeftDown[0] and DzMarch809[m, 0] < PVRightUp[0] and PVLeftDown[1] < \
  265. DzMarch809[m, 1] < PVRightUp[1]: # 在该锯齿中
  266. if (DzMarch809[m, 1] - PVLeftDown[1]) / (DzMarch809[m, 0] - PVLeftDown[0]) > (
  267. PVRightUp[1] - PVLeftDown[1]) / (
  268. PVRightUp[0] - PVLeftDown[0]): # 斜率大于对角连线,则在锯齿左上三角形中,选中
  269. Dzwind_and_power_dfSel[m] = 0
  270. nSmooth = nSmooth + 1
  271. print("nSmooth", nSmooth)
  272. wind_and_power_df.loc[:, 'marker'] = -1
  273. wind_and_power_df.loc[
  274. wind_and_power_df[wind_and_power_df[self.active_power] > 0].index, 'marker'] = Dzwind_and_power_dfSel
  275. wind_and_power_df.to_csv("test.csv", index=False, encoding='utf-8')
  276. # wind_and_power_df = wind_and_power_df[wind_and_power_df['marker'] == 0]
  277. color_map = {-1: 'red', 0: 'green', 1: 'blue', 2: 'black', 3: 'orange', 4: 'magenta'}
  278. c = wind_and_power_df['marker'].map(color_map)
  279. # -1:停机 0:好点 1:欠发功率点;2:超发功率点;3:额定风速以上的超发功率点 4: 限电
  280. legend_map = {"停机": 'red', "好点": 'green', "欠发": 'blue', "超发": 'black', "额定风速以上的超发": 'orange', "限电": 'magenta'}
  281. scatter("测试matlab结果", x_label='风速', y_label='有功功率', x_values=wind_and_power_df[self.wind_velocity].values,
  282. y_values=wind_and_power_df[self.active_power].values, color=c, col_map=legend_map,
  283. save_file_path=os.path.dirname(__file__) + os.sep + '测试matlab结果均值.png')
  284. def run(self):
  285. # Implement your class identification logic here
  286. self.identifier()
  287. if __name__ == '__main__':
  288. test = ClassIdentifier('test', r"D:\中能智能\matlib计算相关\好点坏点matlib计算\A01.csv", index='时间',
  289. wind_velocity='风速',
  290. active_power='功率')
  291. test.run()