import copy import os import pandas as pd import numpy as np import matplotlib.pyplot as plt import pickle current_path = os.path.dirname(__file__) class LimitPower(object): def __init__(self, logger, args, weather_power): self.logger = logger self.args = args self.opt = self.args.parse_args_and_yaml() self.weather_power = weather_power self.step = self.opt.usable_power['step'] self.segs = np.array([x * self.step for x in range(1, 50)]) # 对风速以50为间隔进行分段 self.xs = np.array([self.segs[i - 1] + x if i > 0 else self.step / 2 for i, x in enumerate([self.step / 2 for _ in self.segs])]) # 分段的中间点 self.polynomial = None self.width = 0 self.max_ws = 50 def segment_statis(self): """ 对机头风速-实际功率进行分段处理，获取分度的中位点，四分位间距和斜率 :return: glob_rp 总辐射分段 """ glob_rp = {} # dict: key 辐照度分段中间点 value 分段内的实际功率 for index, row in self.weather_power.iterrows(): ws_ = row[self.opt.usable_power["env"]] rp = row['C_REAL_VALUE'] for i, seg in enumerate(self.segs): if ws_ <= seg: glob_rp.setdefault(self.xs[i], []).append(rp) break rms = [] for i, x in enumerate(self.xs): rps = glob_rp.get(x) if rps is None: continue rps = np.array(rps) up = self.opt.usable_power['up_fractile'] down = self.opt.usable_power['down_fractile'] offset = 0 while True: index = i-1 if i > 0 else 0 while (self.xs[index] in rms or self.xs[index] not in glob_rp) and index >0: index -= 1 x_l = self.xs[index] if index > 0 else 0 q2_l = glob_rp[self.xs[index]][0] if index > 0 else 0 down = down + offset if down > up: rms.append(x) self.logger.info("删除的坐标点为：{}".format(x)) break q1 = np.percentile(rps, down) # 下四分位点 q2 = round(np.percentile(rps, down+(up-down)/2), 3) # 中位点 q3 = np.percentile(rps, up) # 上四分为点 # q2 = np.around(np.mean(rps[(rps >= q1) & (rps <= q3)]), 3) iqr = q3 - q1 # 四分位间距 k2 = round((q2-q2_l)/(x-x_l), 3) # 趋势斜率 # std = np.around(np.std(rps[(rps >= q1) & (rps <= q3)]), 3) # mean = np.around(np.mean(rps), 3) # 实际功率均值 # std = np.around(np.std(rps), 3) # 实际功率标准差 # print("看看q2={}，mean={}".format(q2, mean)) if k2 >= 0: glob_rp[x] = [q2, iqr] # 更新dict break else: offset += 1 glob_rp = {k: glob_rp[k] for k in glob_rp.keys() if k not in rms} glob_rp = {k: glob_rp[k] for k in sorted(glob_rp.keys())} return glob_rp def mapping_relation(self, glob_rp): degree = self.opt.usable_power['degree'] xs = list(glob_rp.keys()) ys = [y[0] for y in glob_rp.values()] self.width = np.median(np.array([y[1] for y in glob_rp.values()])) coefficients = np.polyfit(xs, ys, degree) self.polynomial = np.poly1d(coefficients) self.max_ws = max(xs) # y_fit = self.polynomial(xs) # plt.scatter(xs, ys, label='Data', color='red') # plt.plot(xs, y_fit, label='Fitted polynomial', color='blue') # plt.plot(xs, y_fit+self.width/2, label='up polynomial', color='purple') # plt.plot(xs, y_fit-self.width/2, label='down polynomial', color='green') # plt.legend() # plt.xlabel('x') # plt.ylabel('y') # plt.title(f'Polynomial Fit (degree {degree})') # plt.show() def saveVar(self, path, data): os.makedirs(os.path.dirname(path), exist_ok=True) with open(path, 'wb') as file: pickle.dump(data, file) def filter_unlimited_power(self, ws, real_power): """ 预测可用功主方法 :param zfs: 要预测可用功率的总辐射 :param k: 斜率 :param b: 偏移量 :return: 预测的可用功率 """ # coe = self.opt.usable_power['outliers_threshold'] # seg = self.xs[np.argmax(self.segs >= ws)] up_offset = self.opt.usable_power['up_offset'] down_offset = self.opt.usable_power['down_offset'] high = self.polynomial(ws) + self.width/up_offset if self.polynomial(ws) + self.width/up_offset < self.opt.cap else self.opt.cap low = self.polynomial(ws) - self.width/down_offset if self.polynomial(ws) - self.width/down_offset > 0 else 0 if low <= real_power <= high: return True else: return False def clean_limited_power(self, name, cluster=False): if cluster is True: glob_rp = self.segment_statis() self.saveVar(os.path.dirname(current_path) + '/var/glob_rp.pickle', glob_rp) self.mapping_relation(glob_rp) new_weather_power, number = [], 0 # fig, ax = plt.subplots() for index, row in self.weather_power.iterrows(): zfs = row[self.opt.usable_power["env"]] rp = row['C_REAL_VALUE'] if zfs < 0 or rp < 0: continue if self.filter_unlimited_power(zfs, rp) and zfs <= self.max_ws: row['c'] = 'red' new_weather_power.append(row) else: row['c'] = 'blue' new_weather_power.append(row) new_weather_power = pd.concat(new_weather_power, axis=1).T new_weather_power.plot.scatter(x=self.opt.usable_power["env"], y='C_REAL_VALUE', c='c') plt.savefig(current_path + '/figs/测风法{}.png'.format(name)) new_weather_power = new_weather_power[new_weather_power['c'] == 'red'] number = len(new_weather_power) self.logger.info("未清洗限电前，总共有：{}条数据".format(len(self.weather_power))) self.logger.info("清除限电后保留的点有：" + str(number) + " 占比：" + str(round(number / len(self.weather_power), 2))) return new_weather_power.loc[:, ['C_TIME', 'C_REAL_VALUE', 'C_ABLE_VALUE', self.opt.usable_power['env']]] def clean_limited_power_by_signal(self, name): weather_power1 = self.weather_power.copy() weather_power1["signal"] = weather_power1.apply( lambda x: self.signal_result(x["C_IS_RATIONING_BY_MANUAL_CONTROL"], x["C_IS_RATIONING_BY_AUTO_CONTROL"]), axis=1) weather_power1['c'] = weather_power1.apply(lambda x: 'cornflowerblue' if bool(x["signal"]) is True else 'pink', axis=1) weather_power1.plot.scatter(x=self.opt.usable_power["env"], y='C_REAL_VALUE', c='c') plt.savefig(current_path + '/figs/信号法{}.png'.format(name)) weather_power1 = weather_power1[weather_power1['signal'] == False] self.logger.info("信号法-未清洗限电前，总共有：{}条数据".format(len(self.weather_power))) self.logger.info("信号法-清除限电后保留的点有：" + str(len(weather_power1)) + " 占比：" + str( round(len(weather_power1) / len(self.weather_power), 2))) return weather_power1.loc[:, ['C_TIME', 'C_REAL_VALUE', 'C_ABLE_VALUE', self.opt.usable_power['env']]] def signal_result(self, manual, auto): if int(manual) == 0: if int(auto) == 0: return False else: return True else: if int(auto) == 1: return True else: return False if __name__ == '__main__': from logs import Log from config import myargparse log = Log().logger # 实例化配置类 args = myargparse(discription="场站端配置", add_help=False) power = pd.read_csv('./data/power.csv') weather = pd.read_csv('./data/tower-1-process.csv') weather_power = pd.merge(weather, power, on='C_TIME') # 联立数据 lp = LimitPower(log, args, weather_power) # glob_rp = lp.segment_statis() # lp.mapping_relation(glob_rp) lp.clean_limited_power('测试1') # glob_rp = {k: glob_rp[k] for k in sorted(glob_rp.keys())} # keys = list(glob_rp.keys()) # values = [v[0] for v in glob_rp.values()] # import matplotlib.pyplot as plt # fig, ax = plt.subplots() # ax.plot(keys, values) # plt.show()