import pandas as pd import os import numpy as np np.random.seed(42) import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt 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 def segment_statis(self): """ 对总辐射-实际功率进行分段处理，获取分度的中位点，四分位间距和斜率 :return: glob_rp 总辐射分段 """ segs = [x for x in range(50, 2000, 100)] # 对辐照度以100为间隔进行分段 xs = [segs[i-1]+x if i>0 else 25 for i, x in enumerate([50 for _ in segs])] # 分段的中间点 glob_rp = {} # dict: key 辐照度分段中间点 value 分段内的实际功率 for index, row in self.weather_power.iterrows(): glob_ = row[self.opt.usable_power["env"]] rp = row['C_REAL_VALUE'] for i, seg in enumerate(segs): if glob_ <= seg and not (i > 0 and rp < 1): glob_rp.setdefault(xs[i], []).append(rp) break for i, x in enumerate(xs): rps = glob_rp.get(x) if rps is None: glob_rp = {k: v for k, v in glob_rp.items() if k not in xs[xs.index(x):]} break x_l = xs[i-1] if i > 0 else 0 q2_l = glob_rp[xs[i-1]][0] if i > 0 else 0 q1 = np.percentile(rps, self.opt.usable_power['down_fractile']) # 实际功率下四分位点 q2 = np.percentile(rps, 50) # 实际功率中位点 q3 = np.percentile(rps, self.opt.usable_power['up_fractile']) # 实际功率上四分位点 iqr = q3 -q1 # 四分位间距 k1 = round(q2/x, 5) # 整体斜率 k2 = round((q2-q2_l)/(x-x_l), 5) # 趋势斜率，相对上一个中位点 glob_rp[x] = [q2, iqr, k1, k2] # 更新dict return glob_rp def mapping_relation(self, glob_rp): """ 拟合分段处理后的斜率和偏移量 :param glob_rp: 总辐射分段 :return: k_final 斜率 bias 实际功率的分布宽度， glob_rp 总辐射分段 """ ks, iqrs, delete_x, tag_x = [], [], [], [] # ks所有分段斜率集合，iqrs所有分段间距集合,delete_x删除的x坐标集合 for x, values in glob_rp.items(): k1 = values[-2] k2 = values[-1] iqrs.append(values[-3]) if k1 > 0 and k2 > 0: # 清除趋势小于等于0的斜率 ks.append(k1) tag_x.append(x) else: delete_x.append(x) # print("删除的斜率：", k1, k2) bias = round(np.median(iqrs), 3) # 中位点 # print("++++1", ks) mean = np.mean(ks) # 均值 std = np.std(ks) # 标准差 ks = np.array(ks) z_score = (ks-mean)/std # z均值 # print("----", z_score) outliers = np.abs(z_score) > self.opt.usable_power['outliers_threshold'] # 超过阈值为离群点 ks = ks[~outliers] # 消除离群点 delete_x1 = list(np.array(tag_x)[outliers]) # 清除大于阈值的离群点 k_final = round(np.mean(ks), 5) # 对清洗后的斜率做平均 # print("++++2：", ks) delete_x.extend(delete_x1) self.logger.info("拟合可用功率，删除的斜率：" + ' '.join([str(x) for x in delete_x])) glob_rp = {k: v for k, v in glob_rp.items() if k not in delete_x} # 清洗后剩下的分段点位 return k_final, bias, glob_rp def filter_unlimited_power(self, zfs, real_power, k, b): """ 预测可用功主方法 :param zfs: 要预测可用功率的总辐射 :param k: 斜率 :param b: 偏移量 :return: 预测的可用功率 """ high = k*zfs+b/2 if k*zfs+b/2 < self.opt.cap else self.opt.cap low = k*zfs-b/2 if k*zfs-b/2 > 0 else 0 if low <= real_power <= high: return True else: return False def clean_limited_power(self, name, is_repair=False): if is_repair is True: glob_rp = self.segment_statis() k_final, bias, glob_rp = self.mapping_relation(glob_rp) self.opt.usable_power['k'] = float(k_final) self.opt.usable_power['bias'] = float(bias) new_weather_power = [] for index, row in self.weather_power.iterrows(): zfs = row[self.opt.usable_power["env"]] rp = row['C_REAL_VALUE'] if self.filter_unlimited_power(zfs, rp, self.opt.usable_power['k'], self.opt.usable_power['bias'] * self.opt.usable_power['coe']): 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']] 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']] 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__': power = pd.read_csv('2023-12-01至2023-12-23实际功率导出文件.csv', date_parser=['时间']) weather = pd.read_csv('2023-12-01至2023-12-23气象站数据导出文件.csv', date_parser=['时间']) weather_power = pd.merge(weather, power, on='时间') # 联立数据 # glob_rp = segment_statis(weather_power) # k_final, bias, glob_rp = mapping_relation(glob_rp)