test/processing_limit_power/processing_limit_power_by_wind.py

#!/usr/bin/env python
# -*- coding:utf-8 -*-
# @FileName  :processing_limit_power_by_wind.py
# @Time      :2024/12/3 14:32
# @Author    :David
# @Company: shenyang JY

import os, time
import pandas as pd
import numpy as np
from pymongo import MongoClient
from flask import request, app
from logs import Log
import matplotlib.pyplot as plt
import pickle, traceback
current_path = os.path.dirname(__file__)
parent_path = os.path.dirname(current_path)

def get_data_from_mongo(args):
    mongodb_connection,mongodb_database,mongodb_read_table = "mongodb://root:sdhjfREWFWEF23e@192.168.1.43:30000/",args['mongodb_database'],args['mongodb_read_table']
    client = MongoClient(mongodb_connection)
    # 选择数据库（如果数据库不存在，MongoDB 会自动创建）
    db = client[mongodb_database]
    collection = db[mongodb_read_table]  # 集合名称
    data_from_db = collection.find()  # 这会返回一个游标（cursor）
    # 将游标转换为列表，并创建 pandas DataFrame
    df = pd.DataFrame(list(data_from_db))
    client.close()
    return df

def insert_data_into_mongo(res_df,args):
    mongodb_connection,mongodb_database,mongodb_write_table = "mongodb://root:sdhjfREWFWEF23e@192.168.1.43:30000/",args['mongodb_database'],args['mongodb_write_table']
    client = MongoClient(mongodb_connection)
    db = client[mongodb_database]
    if mongodb_write_table in db.list_collection_names():
        db[mongodb_write_table].drop()
        print(f"Collection '{mongodb_write_table} already exist, deleted successfully!")
    collection = db[mongodb_write_table]  # 集合名称
    # 将 DataFrame 转为字典格式
    data_dict = res_df.to_dict("records")  # 每一行作为一个字典
    # 插入到 MongoDB
    collection.insert_many(data_dict)
    print("data inserted successfully!")

@app.route('/processing_limit_power_by_wind', methods=['POST', 'GET'])
def processing_limit_power_by_agcavc():
    # 获取程序开始时间
    start_time = time.time()
    result = {}
    success = 0
    print("Program starts execution!")
    try:
        logger = Log().logger
        args = request.values.to_dict()
        weather_power = get_data_from_mongo(args)
        lp = LimitPower(logger, args, weather_power)
        weather_power = lp.clean_limited_power('')
        print('args', args)
        logger.info(args)
        insert_data_into_mongo(weather_power, args)
        success = 1
    except Exception as e:
        my_exception = traceback.format_exc()
        my_exception.replace("\n", "\t")
        result['msg'] = my_exception
    end_time = time.time()

    result['success'] = success
    result['args'] = args
    result['start_time'] = time.strftime('%Y-%m-%d %H:%M:%S', time.localtime(start_time))
    result['end_time'] = time.strftime('%Y-%m-%d %H:%M:%S', time.localtime(end_time))
    print("Program execution ends!")
    return result

class LimitPower(object):
    def __init__(self, logger, args, weather_power):
        self.logger = logger
        self.args = args
        self.weather_power = weather_power
        self.step = self.args.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.args.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.args.usable_power['up_offset']
        down_offset = self.args.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):
        glob_rp = self.segment_statis()
        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.args.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.args.usable_power["env"], y='C_REAL_VALUE', c='c')
        plt.savefig(parent_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.args.usable_power['env']]]


if __name__ == '__main__':
    from logs import Log
    log = Log().logger
    args = {}
    # 实例化配置类
    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()