WTOAAM/wtoaamapi/apps/business/algorithm/yawErrorAnalyst.py

import os
import numpy as np
from plotly.subplots import make_subplots
import plotly.graph_objects as go
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import pandas as pd
from .analyst import Analyst
from .utils.directoryUtil import DirectoryUtil as dir
from confBusiness import ConfBusiness


class YawErrorAnalyst(Analyst):
    """
    风电机组静态偏航误差分析
    """

    def typeAnalyst(self):
        return "yaw_error"

    def turbineAnalysis(self,
                        dataFrame,
                        outputAnalysisDir,
                        outputFilePath,
                        confData: ConfBusiness,
                        turbineName):

        self.yaw_error(dataFrame, outputFilePath,
                       confData.field_turbine_time, confData.field_angle_included, confData.field_power, confData.field_pitch_angle1)

    def yaw_error(self, dataFrame, output_path,  field_time, field_angle_included_wind_dir, field_power_active, pitch_col):
        # Calculate floor values and other transformations
        dataFrame['wind_dir_floor'] = np.floor(dataFrame[field_angle_included_wind_dir]).astype(int)
        # dataFrame['wind_dir_floor'] = dataFrame[field_angle_included_wind_dir].round().astype(int)

        dataFrame['power'] = dataFrame[field_power_active].astype(float)
        dataFrame['time'] = pd.to_datetime(dataFrame[field_time])

        # Calculate aggregated metrics for power
        grouped = dataFrame.groupby('wind_dir_floor').agg({
            'power': ['mean', 'max', 'min', 'median', lambda x: (x > 0).sum(), lambda x: (x > x.median()).sum()]
        }).reset_index()

        # Rename columns for clarity
        grouped.columns = ['wind_dir_floor', 'mean_power', 'max_power',
                           'min_power', 'median_power', 'power_gt_0', 'power_gt_80p']

        # Calculate total sums for conditions
        power_gt_0_sum = grouped['power_gt_0'].sum()
        power_gt_80p_sum = grouped['power_gt_80p'].sum()

        # Calculate ratios
        grouped['ratio_0'] = grouped['power_gt_0'] / power_gt_0_sum
        grouped['ratio_80p'] = grouped['power_gt_80p'] / power_gt_80p_sum

        # Filter out zero ratios and calculate slope
        grouped = grouped[grouped['ratio_0'] > 0]
        grouped['slop'] = grouped['ratio_80p'] / grouped['ratio_0']

        # Sort by wind direction floor
        grouped.sort_values('wind_dir_floor', inplace=True)

        # Write to CSV
        grouped.to_csv(output_path, index=False)

    def turbinesAnalysis(self, dataFrameMerge, outputAnalysisDir, confData: ConfBusiness):
        self.yaw_result(outputAnalysisDir)

    def poly_func(self, x, a, b, c, d, e):
        return a * x**4 + b * x ** 3 + c * x ** 2 + d * x + e

    def yaw_result(self, csvFileDir):
        files = os.listdir(csvFileDir)
        for file in files:
            old_name = os.path.join(csvFileDir, file)
            print(file)
            # if os.path.isdir(old_name):
            #    data_stat(old_name)
            #    continue
            # print(old_name)
            if not file.endswith(".csv"):
                continue
            data_path = old_name
            # limit_path = data_path+".limit"

            df = pd.read_csv(data_path)
            
            # 不具备普适性，视静态偏航误差输出数据结果调整，经2024-3-13验证输出结果不可信
            # df = df[df["wind_dir_floor"] > -12]
            # df = df[df["wind_dir_floor"] < 12]

            df.dropna(inplace=True)
            # drop extreme value, the 3 largest and 3 smallest

            df_ = df[df["ratio_80p"] > 0.000]

            xdata = df_['wind_dir_floor']
            ydata = df_['slop']

            # make ydata smooth
            ydata = ydata.rolling(7).median()[6:]
            xdata = xdata[3:-3]
            
            # Curve fitting
            popt, pcov = curve_fit(self.poly_func, xdata, ydata)
            # popt contains the optimized parameters a, b, and c
            # Generate fitted y-dataFrame using the optimized parameters
            fitted_ydata = self.poly_func(xdata, *popt)

            # get the max value of fitted_ydata and its index
            max_pos = fitted_ydata.idxmax()

            # subplots
            plt.figure(figsize=(10, 6))
            fig = plt.subplot(211)
            plt.title(file)
            plt.scatter(x=df["wind_dir_floor"], y=df["slop"],
                        s=5, label="energy gain")
            plt.plot(xdata, fitted_ydata, 'r-', label="fit")
            plt.scatter(xdata[max_pos], fitted_ydata[max_pos], s=20,
                        label="max pos:"+str(df["wind_dir_floor"][max_pos]))
            plt.legend()
            plt.grid(True)
            fig = plt.subplot(212)
            plt.scatter(x=df["wind_dir_floor"],
                        y=df["ratio_0"], s=5, label="base energy")
            plt.scatter(x=df["wind_dir_floor"],
                        y=df["ratio_80p"], s=5, label="slope energy")
            plt.legend()
            plt.grid(True)
            plt.savefig(data_path+".png")
            plt.close()
            # calc squear error of fitted_ydata and ydata

            print(file, df["wind_dir_floor"][max_pos])

    # def yaw_result(self, csvFileDir):
    #     files = os.listdir(csvFileDir)
    #     for file in files:
    #         if not file.endswith(".csv"):
    #             continue

    #         data_path = os.path.join(csvFileDir, file)
    #         df = pd.read_csv(data_path)
    #         df.dropna(inplace=True)

    #         # Filter and smooth data
    #         df_ = df[df["ratio_80p"] > 0.000]
    #         xdata = df_['wind_dir_floor']
    #         ydata = df_['slop'].rolling(7).median()[6:]
    #         xdata = xdata[3:-3].reset_index(drop=True)
    #         ydata = ydata.reset_index(drop=True)

    #         # Curve fitting
    #         popt, _ = curve_fit(self.poly_func, xdata, ydata)
    #         fitted_ydata = self.poly_func(xdata, *popt)
    #         max_pos = fitted_ydata.argmax()

    #         # Plotting
    #         fig = make_subplots(rows=2, cols=1)

    #         # Scatter plot of energy gain
    #         fig.add_trace(go.Scatter(x=df["wind_dir_floor"], y=df["slop"], mode='markers', name="energy gain", marker=dict(size=5)), row=1, col=1)
    #         # Line plot of fit
    #         fig.add_trace(go.Scatter(x=xdata, y=fitted_ydata, mode='lines', name="fit", line=dict(color='red')), row=1, col=1)
    #         # Marker for max position
    #         fig.add_trace(go.Scatter(x=[xdata[max_pos]], y=[fitted_ydata[max_pos]], mode='markers', name=f"max pos: {xdata[max_pos]}", marker=dict(size=10, color='orange')), row=1, col=1)

    #         # Scatter plot of base and slope energy
    #         fig.add_trace(go.Scatter(x=df["wind_dir_floor"], y=df["ratio_0"], mode='markers', name="base energy", marker=dict(size=5)), row=2, col=1)
    #         fig.add_trace(go.Scatter(x=df["wind_dir_floor"], y=df["ratio_80p"], mode='markers', name="slope energy", marker=dict(size=5)), row=2, col=1)

    #         fig.update_layout(height=600, width=800, title_text=file)
    #         fig.write_image(data_path+".png")

    #         print(file, xdata[max_pos])