import os import numpy as np from pandas import DataFrame from utils.draw.draw_file import scatter from utils.file.trans_methods import read_file_to_df from utils.log.trans_log import trans_print class ClassIdentifier(object): """ 分类标识 -1:停机 0:好点 1:欠发功率点；2:超发功率点；3:额定风速以上的超发功率点 4: 限电 """ def __init__(self, wind_turbine_number=None, origin_df: DataFrame = None, wind_velocity='wind_velocity', active_power='active_power', pitch_angle_blade='pitch_angle_blade_1', rated_power=1500, file_path: str = None): """ :param file_path: The file path of the input data. :param origin_df: The pandas DataFrame containing the input data. :param wind_velocity: 风速字段 :param active_power: 有功功率字段 :param pitch_angle_blade: 桨距角 :param rated_power: 额定功率 """ self.wind_velocity = wind_velocity self.active_power = active_power self.pitch_angle_blade = pitch_angle_blade self.rated_power = rated_power # 额定功率1500kw,可改为2000kw if self.rated_power is None: trans_print(wind_turbine_number, "WARNING:rated_power配置为空的") self.rated_power = 1500 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 def identifier(self): # 风速 和 有功功率 df # wind_and_power_df = self.df[[self.wind_velocity, self.active_power, "pitch_angle_blade_1"]] wind_and_power_df = self.df wind_and_power_df.reset_index(inplace=True) wind_and_power_df_count = wind_and_power_df.shape[0] power_max = wind_and_power_df[self.active_power].max() power_rated = np.ceil(power_max / 100) * 100 v_cut_out = 25 # 网格法确定风速风向分区数量，功率方向分区数量， p_num = int(np.ceil(power_rated / 25)) # 功率分区间隔25kW v_num = int(np.ceil(v_cut_out / 0.25)) # 风速分区间隔0.25m/s # 存储功率大于零的运行数据 dz_march = np.zeros([wind_and_power_df_count, 2], dtype=float) n_counter1 = 0 for i in range(wind_and_power_df_count): if wind_and_power_df.loc[i, self.active_power] > 0: dz_march[n_counter1, 0] = wind_and_power_df.loc[i, self.wind_velocity] dz_march[n_counter1, 1] = wind_and_power_df.loc[i, self.active_power] n_counter1 = n_counter1 + 1 # 统计各网格落入的散点个数 if v_num == 1: x_box_number = np.ones([p_num], dtype=int) else: x_box_number = np.ones([p_num, v_num], dtype=int) n_which_p = -1 n_which_v = -1 for i in range(n_counter1): for m in range(p_num): if m * 25 < dz_march[i, 1] <= (m + 1) * 25: n_which_p = m break for n in range(v_num): if ((n + 1) * 0.25 - 0.125) < dz_march[i, 0] <= ((n + 1) * 0.25 + 0.125): n_which_v = n break if n_which_p > -1 and n_which_v > -1: x_box_number[n_which_p, n_which_v] = x_box_number[n_which_p, n_which_v] + 1 for m in range(p_num): for n in range(v_num): x_box_number[m, n] = x_box_number[m, n] - 1 # 在功率方向将网格内散点绝对个数转换为相对百分比，备用 p_box_percent = np.zeros([p_num, v_num], dtype=float) p_bin_sum = np.zeros(p_num, dtype=int) for i in range(p_num): for m in range(v_num): p_bin_sum[i] = p_bin_sum[i] + x_box_number[i, m] for m in range(v_num): if p_bin_sum[i] > 0: p_box_percent[i, m] = x_box_number[i, m] / p_bin_sum[i] * 100 # 在风速方向将网格内散点绝对个数转换为相对百分比，备用 v_box_percent = np.zeros([p_num, v_num], dtype=float) v_bin_sum = np.zeros(v_num, dtype=int) for i in range(v_num): for m in range(p_num): v_bin_sum[i] = v_bin_sum[i] + x_box_number[m, i] for m in range(p_num): if v_bin_sum[i] > 0: v_box_percent[m, i] = x_box_number[m, i] / v_bin_sum[i] * 100 # 以水平功率带方向为准，分析每个水平功率带中，功率主带中心，即找百分比最大的网格位置。 p_box_max_index = np.zeros(p_num, dtype=int) # 水平功率带最大网格位置索引 p_box_max_p = np.zeros(p_num, dtype=int) # 水平功率带最大网格百分比 for m in range(p_num): # 确定每一水平功率带的最大网格位置索引即百分比值 p_box_max_p[m], p_box_max_index[m] = p_box_percent[m, :].max(), p_box_percent[m, :].argmax() # 以垂直风速方向为准，分析每个垂直风速带中，功率主带中心，即找百分比最大的网格位置。 v_box_max_index = np.zeros(v_num, dtype=int) v_box_max_v = np.zeros(v_num, dtype=int) for m in range(v_num): [v_box_max_v[m], v_box_max_index[m]] = v_box_percent[:, m].max(), v_box_percent[:, m].argmax() # 切入风速特殊处理，如果切入风速过于偏右，向左拉回 if p_box_max_index[0] > 14: p_box_max_index[0] = 9 # 以水平功率带方向为基准，进行分析 dot_dense = np.zeros(p_num, dtype=int) # 每一水平功率带的功率主带包含的网格数 dot_dense_left_right = np.zeros([p_num, 2], dtype=int) # 存储每一水平功率带的功率主带以最大网格为中心，向向左，向右扩展的网格数 dot_valve = 90 # 从中心向左右对称扩展网格的散点百分比和的阈值。 for i in range(p_num - 6): # 从最下层水平功率带1开始，向上到第PNum-6个水平功率带（额定功率一下水平功率带），逐一分析 p_dot_dense_sum = p_box_max_p[i] # 以中心最大水平功率带为基准，向左向右对称扩展网格，累加各网格散点百分比 i_spread_right = 1 i_spread_left = 1 while p_dot_dense_sum < dot_valve: if (p_box_max_index[i] + i_spread_right) < v_num - 1: p_dot_dense_sum = p_dot_dense_sum + p_box_percent[i, p_box_max_index[i] + i_spread_right] # 向右侧扩展 i_spread_right = i_spread_right + 1 if (p_box_max_index[i] + i_spread_right) > v_num - 1: break if (p_box_max_index[i] - i_spread_left) > 0: p_dot_dense_sum = p_dot_dense_sum + p_box_percent[i, p_box_max_index[i] - i_spread_left] # 向左侧扩展 i_spread_left = i_spread_left + 1 if (p_box_max_index[i] - i_spread_left) <= 0: break i_spread_right = i_spread_right - 1 i_spread_left = i_spread_left - 1 # 向左右对称扩展完毕 dot_dense_left_right[i, 0] = i_spread_left dot_dense_left_right[i, 1] = i_spread_right dot_dense[i] = i_spread_left + i_spread_right + 1 # 各行功率主带右侧宽度的中位数最具有代表性 dot_dense_width_left = np.zeros([p_num - 6, 1], dtype=int) for i in range(p_num - 6): dot_dense_width_left[i] = dot_dense_left_right[i, 1] main_band_right = np.median(dot_dense_width_left) # 散点向右显著延展分布的水平功率带为限功率水平带 power_limit = np.zeros([p_num, 1], dtype=int) # 各水平功率带是否为限功率标识，==1：是；==0：不是 width_average = 0 # 功率主带平均宽度 width_var = 0 # 功率主带方差 # power_limit_valve = 6 #限功率主带判别阈值 power_limit_valve = np.ceil(main_band_right) + 3 # 限功率主带判别阈值 n_counter_limit = 0 n_counter = 0 for i in range(p_num - 6): if dot_dense_left_right[i, 1] > power_limit_valve and p_bin_sum[i] > 20: # 如果向右扩展网格数大于阈值，且该水平功率带点总数>20，是 power_limit[i] = 1 n_counter_limit = n_counter_limit + 1 if dot_dense_left_right[i, 1] <= power_limit_valve: width_average = width_average + dot_dense_left_right[i, 1] # 统计正常水平功率带右侧宽度 n_counter = n_counter + 1 width_average = width_average / n_counter # 功率主带平均宽度 # 各水平功率带的功率主带宽度的方差，反映从下到上宽度是否一致，或是否下宽上窄等异常情况 for i in range(p_num - 6): if dot_dense_left_right[i, 1] <= power_limit_valve: width_var = width_var + (dot_dense_left_right[i, 1] - width_average) * ( dot_dense_left_right[i, 1] - width_average) # 对限负荷水平功率带的最大网格较下面相邻层显著偏右，拉回 for i in range(1, p_num - 6): if power_limit[i] == 1 and abs(p_box_max_index[i] - p_box_max_index[i - 1]) > 5: p_box_max_index[i] = p_box_max_index[i - 1] + 1 # 输出各层功率主带的左右边界网格索引 dot_dense_inverse = np.zeros([p_num, 2], dtype=int) for i in range(p_num): dot_dense_inverse[i, :] = dot_dense_left_right[p_num - i - 1, :] # 功率主带的右边界 curve_width_r = int(np.ceil(width_average) + 2) # curve_width_l = 6 #功率主带的左边界 curve_width_l = curve_width_r b_box_limit = np.zeros([p_num, v_num], dtype=int) # 网格是否为限功率网格的标识，如果为限功率水平功率带，从功率主带右侧边缘向右的网格为限功率网格 for i in range(2, p_num - 6): if power_limit[i] == 1: for j in range(p_box_max_index[i] + curve_width_r, v_num): b_box_limit[i, j] = 1 b_box_remove = np.zeros([p_num, v_num], dtype=int) # 数据异常需要剔除的网格标识，标识==1：功率主带右侧的欠发网格；==2：功率主带左侧的超发网格 for m in range(p_num - 6): for n in range(p_box_max_index[m] + curve_width_r, v_num): b_box_remove[m, n] = 1 for n in range(p_box_max_index[m] - curve_width_l, -1, -1): b_box_remove[m, n] = 2 # 确定功率主带的左上拐点，即额定风速位置的网格索引 curve_top = np.zeros(2, dtype=int) curve_top_valve = 3 # 网格的百分比阈值 b_top_find = 0 for m in range(p_num - 4 - 1, -1, -1): for n in range(v_num): if v_box_percent[m, n] > curve_top_valve and x_box_number[m, n] >= 10: # 如左上角网格的百分比和散点个数大于阈值。 curve_top[0] = m curve_top[1] = n b_top_find = 1 break if b_top_find == 1: break isolate_valve = 3 for m in range(p_num - 6): for n in range(p_box_max_index[m] + curve_width_r, v_num): if p_box_percent[m, n] < isolate_valve: b_box_remove[m, n] = 1 # 功率主带顶部宽度 curve_width_t = 2 for m in range(p_num - curve_width_t - 1, p_num): for n in range(v_num): b_box_remove[m, n] = 3 # 网格为额定功率以上的超发点 # 功率主带拐点左侧的欠发网格标识 for m in range(p_num - 5 - 1, p_num): for n in range(curve_top[1] - 1): b_box_remove[m, n] = 2 # 以网格的标识，决定该网格内数据的标识。dzwind_and_power_sel。散点在哪个网格，此网格的标识即为该点的标识 dzwind_and_power_sel = np.zeros(n_counter1, dtype=int) # -1:停机 0:好点 1:欠发功率点；2:超发功率点；3:额定风速以上的超发功率点 4: 限电 n_which_p = -1 n_which_v = -1 n_bad_a = 0 for i in range(n_counter1): for m in range(p_num): if m * 25 < dz_march[i, 1] <= (m + 1) * 25: n_which_p = m break for n in range(v_num): if ((n + 1) * 0.25 - 0.125) < dz_march[i, 0] <= ((n + 1) * 0.25 + 0.125): n_which_v = n break if n_which_p > -1 and n_which_v > -1: if b_box_remove[n_which_p, n_which_v] == 1: dzwind_and_power_sel[i] = 1 n_bad_a = n_bad_a + 1 if b_box_remove[n_which_p, n_which_v] == 2: dzwind_and_power_sel[i] = 2 if b_box_remove[n_which_p, n_which_v] == 3: dzwind_and_power_sel[i] = 0 # 3 # 额定风速以上的超发功率点认为是正常点，不再标识。 # 限负荷数据标识方法2：把数据切割为若干个窗口。对每一窗口，以第一个点为基准，连续nWindowLength个数据的功率在方差范围内，呈现显著水平分布的点 n_window_length = 3 limit_window = np.zeros(n_window_length, dtype=float) power_std = 15 # 功率波动方差 n_window_num = int(np.floor(n_counter1 / n_window_length)) power_limit_up = self.rated_power - 300 power_limit_low = 200 for i in range(n_window_num): for j in range(n_window_length): limit_window[j] = dz_march[i * n_window_length + j, 1] b_all_in_areas = 1 for j in range(n_window_length): if limit_window[j] < power_limit_low or limit_window[j] > power_limit_up: b_all_in_areas = 0 if b_all_in_areas == 0: continue up_limit = limit_window[0] + power_std low_limit = limit_window[0] - power_std b_all_in_up_low = 1 for j in range(1, n_window_length): if limit_window[j] < low_limit or limit_window[j] > up_limit: b_all_in_up_low = 0 if b_all_in_up_low == 1: for j in range(n_window_length): dzwind_and_power_sel[i * n_window_length + j] = 4 # 标识窗口内的数据为限负荷数据 for i in range(p_num - 6): pv_left_down = np.zeros(2, dtype=float) pv_right_up = np.zeros(2, dtype=float) if (p_box_max_index[i + 1] - p_box_max_index[i]) >= 1: pv_left_down[0] = (p_box_max_index[i] + 1 + curve_width_r) * 0.25 - 0.125 pv_left_down[1] = i * 25 pv_right_up[0] = (p_box_max_index[i + 1] + 1 + curve_width_r) * 0.25 - 0.125 pv_right_up[1] = (i + 1) * 25 for m in range(n_counter1): if pv_left_down[0] < dz_march[m, 0] < pv_right_up[0] and pv_left_down[1] < \ dz_march[m, 1] < pv_right_up[1]: # 在该锯齿中 if (dz_march[m, 1] - pv_left_down[1]) / (dz_march[m, 0] - pv_left_down[0]) > ( pv_right_up[1] - pv_left_down[1]) / ( pv_right_up[0] - pv_left_down[0]): # 斜率大于对角连线，则在锯齿左上三角形中，选中 dzwind_and_power_sel[m] = 0 wind_and_power_df.loc[:, 'lab'] = -1 wind_and_power_df.loc[ wind_and_power_df[wind_and_power_df[self.active_power] > 0].index, 'lab'] = dzwind_and_power_sel # 把部分欠发的优化为限电 # 构建条件表达式 cond1 = (wind_and_power_df['lab'] == 1) & ( (wind_and_power_df[self.active_power] < self.rated_power * 0.75) & (wind_and_power_df[self.pitch_angle_blade] > 0.5) ) cond2 = (wind_and_power_df['lab'] == 1) & ( (wind_and_power_df[self.active_power] < self.rated_power * 0.85) & (wind_and_power_df[self.pitch_angle_blade] > 1.5) ) cond3 = (wind_and_power_df['lab'] == 1) & ( (wind_and_power_df[self.active_power] < self.rated_power * 0.9) & (wind_and_power_df[self.pitch_angle_blade] > 2.5) ) # 使用逻辑或操作符|合并条件 combined_condition = cond1 | cond2 | cond3 wind_and_power_df.loc[combined_condition, 'lab'] = 4 wind_and_power_df.reset_index(drop=True, inplace=True) if 'index' in wind_and_power_df.columns: del wind_and_power_df['index'] return wind_and_power_df def run(self): # Implement your class identification logic here return self.identifier() if __name__ == '__main__': read_dir = r"D:\data\清理数据\和风元宝山\WOF035200003-WOB000005111_MM14号机组0719\minute" files = [read_dir + os.sep + i for i in os.listdir(read_dir)] for file in files: # test = ClassIdentifier(file_path=file, # wind_velocity='wind_velocity', # active_power='active_power', # pitch_angle_blade='pitch_angle_blade_1', # rated_power=1500 # ) # # df = test.run() name = os.path.basename(file).split('.')[0] df = read_file_to_df(file) color_map = {-1: 'red', 0: 'green', 1: 'blue', 2: 'black', 3: 'orange', 4: 'magenta'} c = df['lab'].map(color_map) # -1:停机 0:好点 1:欠发功率点；2:超发功率点；3:额定风速以上的超发功率点 4: 限电 legend_map = {"停机": 'red', "好点": 'green', "欠发": 'blue', "超发": 'black', "额定风速以上的超发": 'orange', "限电": 'magenta'} scatter(name, x_label='风速', y_label='有功功率', x_values=df['wind_velocity'].values, y_values=df['active_power'].values, color=c, col_map=legend_map, save_file_path=os.path.dirname( os.path.dirname( os.path.dirname(__file__))) + os.sep + "tmp_file" + os.sep + "和风元宝山" + os.sep + name + '结果.png')