WTOAAM/dataAnalysisBehavior/behavior/dataMarkerOfScada3.py

import numpy as np
import pandas as pd
import os
import plotly.graph_objects as go

# 声明全局变量
fieldActivePower = "功率"
fieldWindSpeed = "风速"
fieldPitchAngle = "叶片角度"
PRated = 1500  # 风机额定功率
VCutOut = 25  # 切出风速
VCutIn = 3  # 切入风速
VRated = 10  # 额定风速
VNum = int(VCutOut / 0.25)  # 风速分区数量

# 读取数据函数


def read_data():
    # 读取 CSV 文件
    March809 = pd.read_csv(
        "./data/DataClassificationIdentification/A01-G.csv", encoding="utf-8")
    IdealCurve = pd.read_csv(
        "./data/DataClassificationIdentification/A型风机设计功率曲线.csv", encoding="utf-8")
    return March809, IdealCurve

# 计算统计信息函数


def calculate_statistics(March809, IdealCurve):
    AA = len(March809)
    BB = len(IdealCurve)
    PowerMax = March809[fieldActivePower].max()
    PowerRated = int(np.ceil(PowerMax / 100) * 100)
    PNum = PowerRated // 25  # 功率分区数量

    # 计算实际发电量
    EPActualTotal = March809[March809[fieldActivePower]
                             >= 0][fieldActivePower].sum() / 6
    WindSpeedAvr = March809[fieldWindSpeed].mean()

    # 计算风机可利用率
    nShouldGP = np.sum(March809[fieldWindSpeed] >= VCutIn)
    nRealGP = np.sum((March809[fieldWindSpeed] >= VCutIn)
                     & (March809[fieldActivePower] > 0))
    TurbineRunRate = (nRealGP / nShouldGP * 100) if nShouldGP > 0 else 0

    # 计算理论发电量
    EPIdealTotalAAA = 0
    for i in range(AA):
        nWhichBin = 0
        for m in range(BB - 1):
            if IdealCurve.iloc[m][fieldWindSpeed] < March809.iloc[i][fieldWindSpeed] <= IdealCurve.iloc[m + 1][fieldWindSpeed]:
                nWhichBin = m
                break
        if nWhichBin > 0:
            IdealPower = (March809.iloc[i][fieldWindSpeed] - IdealCurve.iloc[nWhichBin][fieldWindSpeed]) / (IdealCurve.iloc[nWhichBin + 1][fieldWindSpeed] - IdealCurve.iloc[nWhichBin]
                                                                                                            [fieldWindSpeed]) * (IdealCurve.iloc[nWhichBin + 1][fieldActivePower] - IdealCurve.iloc[nWhichBin][fieldActivePower]) + IdealCurve.iloc[nWhichBin][fieldActivePower]
            EPIdealTotalAAA += IdealPower / 6

    return AA, BB, PNum, EPActualTotal, WindSpeedAvr, TurbineRunRate, EPIdealTotalAAA

# 分类数据函数


def classify_data(March809, PNum, VNum):
    DzMarch809 = March809[March809[fieldActivePower] > 0]
    nCounter1 = len(DzMarch809)

    XBoxNumber = np.ones((PNum, VNum), dtype=int)
    for i in range(nCounter1):
        nWhichP = np.digitize(
            DzMarch809.iloc[i][fieldActivePower], np.arange(0, PNum * 25, 25)) - 1
        nWhichV = np.digitize(
            DzMarch809.iloc[i][fieldWindSpeed], np.arange(0.125, VNum * 0.25, 0.25))
        if nWhichP < PNum and nWhichV < VNum:
            XBoxNumber[nWhichP, nWhichV] += 1
    XBoxNumber -= 1

    return DzMarch809, XBoxNumber

# 计算百分比函数


def compute_percentages(XBoxNumber, PNum, VNum):
    PBoxPercent = np.zeros((PNum, VNum))
    PBinSum = XBoxNumber.sum(axis=1)
    for i in range(PNum):
        if PBinSum[i] > 0:
            PBoxPercent[i, :] = XBoxNumber[i, :] / PBinSum[i] * 100

    VBoxPercent = np.zeros((PNum, VNum))
    VBinSum = XBoxNumber.sum(axis=0)
    for i in range(VNum):
        if VBinSum[i] > 0:
            VBoxPercent[:, i] = XBoxNumber[:, i] / VBinSum[i] * 100

    return PBoxPercent, VBoxPercent

# 查找主带函数


def find_main_band(PBoxPercent, PNum, VNum, XBoxNumber):
    PBoxMaxIndex = np.argmax(PBoxPercent, axis=1)
    PBoxMaxP = np.max(PBoxPercent, axis=1)

    DotDense = np.zeros(PNum)
    DotDenseLeftRight = np.zeros((PNum, 2), dtype=int)
    DotValve = 90
    for i in range(PNum - 6):
        PDotDenseSum = PBoxMaxP[i]
        iSpreadRight = iSpreadLeft = 1
        while PDotDenseSum < DotValve:
            if (PBoxMaxIndex[i] + iSpreadRight) < VNum - 1:
                PDotDenseSum += PBoxPercent[i, PBoxMaxIndex[i] + iSpreadRight]
                iSpreadRight += 1
            if (PBoxMaxIndex[i] + iSpreadRight) > VNum - 1:
                break
            if (PBoxMaxIndex[i] - iSpreadLeft) > 0:
                PDotDenseSum += PBoxPercent[i, PBoxMaxIndex[i] - iSpreadLeft]
                iSpreadLeft += 1
            if (PBoxMaxIndex[i] - iSpreadLeft) <= 0:
                break
        iSpreadRight -= 1
        iSpreadLeft -= 1
        DotDenseLeftRight[i, :] = [iSpreadLeft, iSpreadRight]
        DotDense[i] = iSpreadLeft + iSpreadRight + 1

    DotDenseWidthLeft = DotDenseLeftRight[:, 1]
    MainBandRight = np.median(DotDenseWidthLeft)
    PowerLimitValve = int(np.ceil(MainBandRight)) + 3

    PowerLimit = np.zeros(PNum)
    nCounterLimit = nCounter = 0
    WidthAverage = 0
    for i in range(PNum - 6):
        if DotDenseLeftRight[i, 1] > PowerLimitValve and XBoxNumber[i, :].sum() > 20:
            PowerLimit[i] = 1
            nCounterLimit += 1
        else:
            WidthAverage += DotDenseLeftRight[i, 1]
            nCounter += 1
    WidthAverage /= nCounter

    WidthVar = np.sqrt(np.mean((DotDenseLeftRight[:, 1] - WidthAverage) ** 2))
    PowerBandWidth = WidthAverage * 2 * 0.25

    return DotDense, DotDenseLeftRight, PowerLimit, WidthAverage, WidthVar, PowerBandWidth, PBoxMaxIndex

# 标记坏点函数


def mark_bad_points(DzMarch809, DotDenseLeftRight, PBoxMaxIndex, PowerLimit, PNum, VNum, XBoxNumber, PBoxPercent):
    CurveWidthR = int(np.ceil(DotDenseLeftRight[:, 1].mean())) + 2
    CurveWidthL = CurveWidthR
    BBoxLimit = np.zeros((PNum, VNum))
    for i in range(3, PNum - 6):
        if PowerLimit[i] == 1:
            BBoxLimit[i, PBoxMaxIndex[i] + CurveWidthR + 1: VNum] = 1

    BBoxRemove = np.zeros((PNum, VNum))
    for m in range(PNum - 6):
        BBoxRemove[m, PBoxMaxIndex[m] + CurveWidthR: VNum] = 1
        BBoxRemove[m, :PBoxMaxIndex[m] - CurveWidthL + 1] = 2

    CurveTop = [0, 0]
    CurveTopValve = 3
    BTopFind = False
    for m in range(PNum - 4, 0, -1):
        for n in range(VNum):
            if PBoxPercent[m, n] > CurveTopValve and XBoxNumber[m, n] >= 10:
                CurveTop = [m, n]
                BTopFind = True
                break
        if BTopFind:
            break

    IsolateValve = 3
    for m in range(PNum - 6):
        for n in range(PBoxMaxIndex[m] + CurveWidthR, VNum):
            if PBoxPercent[m, n] < IsolateValve:
                BBoxRemove[m, n] = 1

    for m in range(PNum - 6, PNum):
        BBoxRemove[m, :] = 3
    for m in range(PNum - 5, PNum):
        BBoxRemove[m, :CurveTop[1] - 2] = 2

    DzMarch809Sel = np.zeros(len(DzMarch809), dtype=int)
    for i in range(len(DzMarch809)):
        nWhichP = np.digitize(
            DzMarch809.iloc[i][fieldActivePower], np.arange(0, PNum * 25, 25)) - 1
        nWhichV = np.digitize(
            DzMarch809.iloc[i][fieldWindSpeed], np.arange(0.125, VNum * 0.25, 0.25))
        
        # if (
        #     (DzMarch809.iloc[i][fieldActivePower] < PRated * 0.75) and
        #     (DzMarch809.iloc[i][fieldPitchAngle] > 0.5) or
        #     (DzMarch809.iloc[i][fieldActivePower] < PRated * 0.85) and
        #     (DzMarch809.iloc[i][fieldPitchAngle] > 1.5) or
        #     (DzMarch809.iloc[i][fieldActivePower] < PRated * 0.9) and
        #     (DzMarch809.iloc[i][fieldPitchAngle] > 2.5)
        # ):
        #     continue

        if nWhichP < PNum and nWhichV < VNum:
            if BBoxRemove[nWhichP, nWhichV] == 1:
                DzMarch809Sel[i] = 1
            elif BBoxRemove[nWhichP, nWhichV] == 2:
                DzMarch809Sel[i] = 2
            elif BBoxRemove[nWhichP, nWhichV] == 3:
                DzMarch809Sel[i] = 0

    return DzMarch809Sel


def identify_limit_load_data(DzMarch809: pd.DataFrame, DzMarch809Sel, PRated):
    nCounter1 = len(DzMarch809)
    PVLimit = np.zeros((nCounter1, 2))
    nLimitTotal = 0
    nWindowLength = 3
    PowerStd = 15  # 功率波动方差
    PowerLimitUp =  PRated - 300
    PowerLimitLow = 5  # 200kW 

    nWindowNum = nCounter1 // nWindowLength

    for i in range(nWindowNum):
        LimitWindow = DzMarch809.iloc[i * nWindowLength:(i + 1) * nWindowLength][fieldActivePower].values
               
        # 检查所有数据是否在 PowerLimitLow值~PowerLimitUp值范围内
        if not ((LimitWindow >= PowerLimitLow) & (LimitWindow <= PowerLimitUp)).all():
            continue
        
        """
        限功率识别策略（参考）：
        1. 功率<额定功率(PRated)*0.75 and  叶片角度>0.5 ; 2. 功率<额定功率*0.85 and 叶片角度>1.5 ; 3. 功率<额定功率*0.9 and 叶片角度>2.5；
        示例：
        1. 功率<1100，and  叶片角度>0.5 ; 2. 功率<1250 and 叶片角度>1.5 ; 3. 功率<1400 and 叶片角度>2.5 ;
        """        
        # 额外的限功率识别策略
        pitch_angle_window = DzMarch809.iloc[i * nWindowLength:(i + 1) * nWindowLength][fieldPitchAngle].values

        if (
            (DzMarch809.iloc[i * nWindowLength:(i + 1) * nWindowLength][fieldActivePower] < PRated * 0.75).any() and
            (pitch_angle_window > 0.5).any() or
            (DzMarch809.iloc[i * nWindowLength:(i + 1) * nWindowLength][fieldActivePower] < PRated * 0.85).any() and
            (pitch_angle_window > 1.5).any() or
            (DzMarch809.iloc[i * nWindowLength:(i + 1) * nWindowLength][fieldActivePower] < PRated * 0.9).any() and
            (pitch_angle_window > 2.5).any()
        ):
            # 标识限负荷数据
            DzMarch809Sel[i * nWindowLength:(i + 1) * nWindowLength] = 4

        UpLimit = LimitWindow[0] + PowerStd
        LowLimit = LimitWindow[0] - PowerStd

        # 检查所有数据是否在方差范围内
        if not ((LimitWindow[1:] >= LowLimit) & (LimitWindow[1:] <= UpLimit)).all():
            continue

        # 标识限负荷数据
        DzMarch809Sel[i * nWindowLength:(i + 1) * nWindowLength] = 4

        for j in range(nWindowLength):
            # 只提取功率和风速数据
            PVLimit[nLimitTotal] = DzMarch809.iloc[i * nWindowLength + j][[fieldWindSpeed, fieldActivePower]]
            nLimitTotal += 1

    PVLimit = PVLimit[:nLimitTotal]  # 截取实际数据部分

    return PVLimit, DzMarch809Sel

# 计算能量损失函数
def calculate_energy_loss(DzMarch809Sel, DzMarch809, IdealCurve, PNum, BB, PRated, EPIdealTotalAAA, EPActualTotal):
    EPLostStopTotal = 0
    nStopTotal = 0
    for i in range(len(DzMarch809)):
        if DzMarch809.iloc[i][fieldActivePower] <= 0:
            nWhichBin = 0
            for m in range(BB - 1):
                if IdealCurve.iloc[m][fieldWindSpeed] < DzMarch809.iloc[i][fieldWindSpeed] <= IdealCurve.iloc[m + 1][fieldWindSpeed]:
                    nWhichBin = m
                    break
            if nWhichBin > 0:
                IdealPower = (DzMarch809.iloc[i][fieldWindSpeed] - IdealCurve.iloc[nWhichBin][fieldWindSpeed]) / (IdealCurve.iloc[nWhichBin + 1][fieldWindSpeed] - IdealCurve.iloc[nWhichBin]
                                                                                                                  [fieldWindSpeed]) * (IdealCurve.iloc[nWhichBin + 1][fieldActivePower] - IdealCurve.iloc[nWhichBin][fieldActivePower]) + IdealCurve.iloc[nWhichBin][fieldActivePower]
                EPLostStopTotal += IdealPower / 6
                nStopTotal += 1

    EPLostBadTotal = 0
    EPLost = 0
    nBadTotal = 0
    for i in range(len(DzMarch809)):
        if DzMarch809Sel[i] == 1:
            nWhichBin = 0
            for m in range(BB - 1):
                if IdealCurve.iloc[m][fieldWindSpeed] < DzMarch809.iloc[i][fieldWindSpeed] <= IdealCurve.iloc[m + 1][fieldWindSpeed]:
                    nWhichBin = m
                    break
            if nWhichBin > 0:
                IdealPower = (DzMarch809.iloc[i][fieldWindSpeed] - IdealCurve.iloc[nWhichBin][fieldWindSpeed]) / (IdealCurve.iloc[nWhichBin + 1][fieldWindSpeed] - IdealCurve.iloc[nWhichBin]
                                                                                                                  [fieldWindSpeed]) * (IdealCurve.iloc[nWhichBin + 1][fieldActivePower] - IdealCurve.iloc[nWhichBin][fieldActivePower]) + IdealCurve.iloc[nWhichBin][fieldActivePower]
                EPLost += abs(IdealPower -
                              DzMarch809.iloc[i][fieldActivePower]) / 6
                EPLostBadTotal += EPLost
                nBadTotal += 1

    EPOverTotal = 0
    nOverTotal = 0
    for i in range(len(DzMarch809)):
        if DzMarch809Sel[i] == 3:
            EPOver = (DzMarch809.iloc[i][fieldActivePower] - PRated) / 6
            EPOverTotal += EPOver
            nOverTotal += 1

    EPLostPerformTotal = 0
    for i in range(len(DzMarch809)):
        nWhichBinI = 0
        for m in range(BB - 1):
            if IdealCurve.iloc[m][fieldWindSpeed] < DzMarch809.iloc[i][fieldWindSpeed] <= IdealCurve.iloc[m + 1][fieldWindSpeed]:
                nWhichBinI = m
                break
        if nWhichBinI > 0:
            IdealPower = (DzMarch809.iloc[i][fieldWindSpeed] - IdealCurve.iloc[nWhichBinI][fieldWindSpeed]) / (IdealCurve.iloc[nWhichBinI + 1][fieldWindSpeed] - IdealCurve.iloc[nWhichBinI]
                                                                                                               [fieldWindSpeed]) * (IdealCurve.iloc[nWhichBinI + 1][fieldActivePower] - IdealCurve.iloc[nWhichBinI][fieldActivePower]) + IdealCurve.iloc[nWhichBinI][fieldActivePower]
            EPLostPerformTotal += (IdealPower -
                                   DzMarch809.iloc[i][fieldActivePower]) / 6

    EPIdealTotal = EPActualTotal + EPLostStopTotal + EPLostBadTotal + \
        EPLostPerformTotal if EPLostPerformTotal >= 0 else EPActualTotal + \
        EPLostStopTotal + EPLostBadTotal

    RemoveOverEP = 0
    for i in range(len(DzMarch809)):
        if DzMarch809Sel[i] == 2:
            nWhichBin = 0
            for m in range(BB - 1):
                if IdealCurve.iloc[m][fieldWindSpeed] < DzMarch809.iloc[i][fieldWindSpeed] <= IdealCurve.iloc[m + 1][fieldWindSpeed]:
                    nWhichBin = m
                    break
            if nWhichBin > 0:
                IdealPower = (DzMarch809.iloc[i][fieldWindSpeed] - IdealCurve.iloc[nWhichBin][fieldWindSpeed]) / (IdealCurve.iloc[nWhichBin + 1][fieldWindSpeed] - IdealCurve.iloc[nWhichBin]
                                                                                                                  [fieldWindSpeed]) * (IdealCurve.iloc[nWhichBin + 1][fieldActivePower] - IdealCurve.iloc[nWhichBin][fieldActivePower]) + IdealCurve.iloc[nWhichBin][fieldActivePower]
                RemoveOverEP += (DzMarch809.iloc[i]
                                 [fieldActivePower] - IdealPower) / 6

    for i in range(len(DzMarch809)):
        if DzMarch809.iloc[i][fieldActivePower] > PRated:
            RemoveOverEP += (DzMarch809.iloc[i][fieldActivePower] - PRated) / 6

    return EPLostStopTotal, EPLostBadTotal, EPOverTotal, EPLostPerformTotal, EPIdealTotal - RemoveOverEP

# 计算实测功率曲线函数


def calculate_measured_power_curve(PVDot, VRated, PRated):
    XBinNumber = np.ones(50)
    PCurve = np.zeros((50, 2))
    PCurve[:, 0] = np.arange(0.5, 25.5, 0.5)
    XBinSum = np.zeros((50, 2))

    for i in range(len(PVDot)):
        nWhichBin = 0
        for b in range(50):
            if (b * 0.5 - 0.25) < PVDot.iloc[i][fieldWindSpeed] <= (b * 0.5 + 0.25):
                nWhichBin = b
                break

        if nWhichBin > 0:
            XBinSum[nWhichBin, 0] += PVDot.iloc[i][fieldWindSpeed]
            XBinSum[nWhichBin, 1] += PVDot.iloc[i][fieldActivePower]
            XBinNumber[nWhichBin] += 1

    XBinNumber -= 1

    for b in range(50):
        if XBinNumber[b] > 0:
            PCurve[b, 0] = XBinSum[b, 0] / XBinNumber[b]
            PCurve[b, 1] = XBinSum[b, 1] / XBinNumber[b]

    VRatedNum = int(VRated / 0.5)
    for m in range(VRatedNum, 50):
        if PCurve[m, 1] == 0:
            PCurve[m, 1] = PRated

    return PCurve

# 计算标准正则功率曲线函数
def calculate_normalized_power_curve(IdealCurve, VRated, PRated):
    PCurveNorm = np.zeros((50, 2))
    VRatedNum = int(VRated / 0.5)

    # 15m/s以上为额定功率
    high_wind_speeds = np.arange(15, 25.5, 0.5)
    PCurveNorm[VRatedNum:VRatedNum+len(high_wind_speeds), 0] = high_wind_speeds
    PCurveNorm[VRatedNum:VRatedNum+len(high_wind_speeds), 1] = PRated

    # 15m/s以下正则功率曲线
    VSpeed = np.arange(0.5, 15.5, 0.5)
    CurveData = IdealCurve[IdealCurve[fieldWindSpeed]
                           <= 15].to_numpy()  # 提取风速<=15的数据

    for i, v in enumerate(VSpeed):
        if i < len(CurveData) - 1:
            # 插值计算
            x0, y0 = CurveData[i]
            x1, y1 = CurveData[i + 1]
            PCurveNorm[i, 0] = v
            PCurveNorm[i, 1] = y0 + (v - x0) * (y1 - y0) / (x1 - x0)
        else:
            PCurveNorm[i, 0] = v
            PCurveNorm[i, 1] = PRated  # 防止超出范围的情况

    return PCurveNorm

# 保存结果函数


def save_results(PCurve, PCurveNorm, EPKW, EPPer, outputDir='./output/A01'):
    if not os.path.exists(outputDir):
        os.makedirs(outputDir)
    pd.DataFrame(PCurve, columns=[fieldWindSpeed, fieldActivePower]).to_csv(
        os.path.join(outputDir, 'PCurve.csv'), index=False)
    pd.DataFrame(PCurveNorm, columns=[fieldWindSpeed, fieldActivePower]).to_csv(
        os.path.join(outputDir, 'PCurveNorm.csv'), index=False)
    pd.DataFrame([EPKW], columns=['EPIdealTotal', 'EPActualTotal', 'EPLostStopTotal', 'EPLostBadLimitTotal', 'EPLostPerformTotal', 'EPLostBadTotal',
                 'EPLostLimitTotal', 'EPOverTotal', 'WindSpeedAvr', 'TurbineRunRate']).to_csv(os.path.join(outputDir, 'EPKW.csv'), index=False)
    pd.DataFrame([EPPer], columns=['Percent1', 'Percent2', 'Percent3', 'Percent4', 'Percent5', 'Percent6', 'Percent7',
                 'Percent8', 'WindSpeedAvr', 'TurbineRunRate']).to_csv(os.path.join(outputDir, 'EPPer.csv'), index=False)

# 绘制结果函数，使用 Plotly


def plot_results(PVBad, PVLimit, PVDot, PCurve, IdealCurve):
    fig = go.Figure()

    # 添加坏点数据
    if not PVBad.empty:
        fig.add_trace(go.Scatter(x=PVBad[fieldWindSpeed], y=PVBad[fieldActivePower],
                      mode='markers', name='坏点', marker=dict(color='red')))

    # 添加限电点数据
    if not PVLimit.empty:
        fig.add_trace(go.Scatter(x=PVLimit[fieldWindSpeed], y=PVLimit[fieldActivePower],
                      mode='markers', name='限功率', marker=dict(color='blue')))

    # 添加正常点数据
    if not PVDot.empty:
        fig.add_trace(go.Scatter(x=PVDot[fieldWindSpeed], y=PVDot[fieldActivePower],
                      mode='markers', name='好点', marker=dict(color='black')))

    # 添加实测功率曲线
    if PCurve.shape[0] > 0:
        fig.add_trace(go.Scatter(
            x=PCurve[:, 0], y=PCurve[:, 1], mode='lines+markers', name='实测功率曲线', line=dict(color='green')))

    # 添加设计功率曲线
    if IdealCurve.shape[0] > 0:
        fig.add_trace(go.Scatter(x=IdealCurve[:, 0], y=IdealCurve[:, 1],
                      mode='lines+markers', name='合同功率曲线', line=dict(color='yellow')))

    # 更新布局
    fig.update_layout(
        title={'text':'风力机功率散点数据分类标识','x':0.5},
        xaxis_title=fieldWindSpeed,
        yaxis_title=fieldActivePower,
        legend=dict(x=0.01, y=0.99),
        template='plotly_white'
    )

    # fig.show()
    outputDir = './output/A01'
    filePath = os.path.join(outputDir, f'power_scatter.html')
    fig.write_html(filePath)


# 主函数
def main():
    # 读取数据
    March809, IdealCurve = read_data()
    # 计算统计信息
    AA, BB, PNum, EPActualTotal, WindSpeedAvr, TurbineRunRate, EPIdealTotalAAA = calculate_statistics(
        March809, IdealCurve)
    # 分类数据
    DzMarch809, XBoxNumber = classify_data(March809, PNum, VNum)
    # 计算百分比
    PBoxPercent, VBoxPercent = compute_percentages(XBoxNumber, PNum, VNum)
    # 查找主带
    DotDense, DotDenseLeftRight, PowerLimit, WidthAverage, WidthVar, PowerBandWidth, PBoxMaxIndex = find_main_band(
        PBoxPercent, PNum, VNum, XBoxNumber)
    # 标记坏点
    DzMarch809Sel = mark_bad_points(
        DzMarch809, DotDenseLeftRight, PBoxMaxIndex, PowerLimit, PNum, VNum, XBoxNumber, PBoxPercent)
    # 标识限负荷数据点
    PVLimitArray, DzMarch809Sel = identify_limit_load_data(
        DzMarch809, DzMarch809Sel, PRated)

    # 初始化存储好点和坏点的 DataFrame
    PVBad = pd.DataFrame(columns=DzMarch809.columns)
    PVLimit = pd.DataFrame(columns=DzMarch809.columns)
    PVDot = pd.DataFrame(columns=DzMarch809.columns)

    # 存储好点和坏点
    for i in range(len(DzMarch809)):
        if DzMarch809Sel[i] in [1, 2, 3]:
            if not DzMarch809.iloc[[i]].isna().all().all():
                PVBad = pd.concat(
                    [PVBad, DzMarch809.iloc[[i]]], ignore_index=True)
        elif DzMarch809Sel[i] == 4:
            if not DzMarch809.iloc[[i]].isna().all().all():
                PVLimit = pd.concat(
                    [PVLimit, DzMarch809.iloc[[i]]], ignore_index=True)
        else:
            if not DzMarch809.iloc[[i]].isna().all().all():
                PVDot = pd.concat(
                    [PVDot, DzMarch809.iloc[[i]]], ignore_index=True)

    # 计算能量损失
    EPLostStopTotal, EPLostBadTotal, EPOverTotal, EPLostPerformTotal, EPIdealTotal = calculate_energy_loss(
        DzMarch809Sel, DzMarch809, IdealCurve, PNum, BB, PRated, EPIdealTotalAAA, EPActualTotal)

    EPKW = [EPIdealTotal, EPActualTotal, EPLostStopTotal, EPLostBadTotal,
            EPLostPerformTotal, EPLostBadTotal, 0, EPOverTotal, WindSpeedAvr, TurbineRunRate]
    EPPer = [100, EPActualTotal / EPIdealTotal * 100, EPLostStopTotal / EPIdealTotal * 100, EPLostBadTotal / EPIdealTotal * 100,
             EPLostPerformTotal / EPIdealTotal * 100, EPLostBadTotal / EPIdealTotal * 100, 0, EPOverTotal / EPActualTotal * 100, WindSpeedAvr, TurbineRunRate]

    # 计算实测功率曲线
    PCurve = calculate_measured_power_curve(PVDot, VRated, PRated)
    # 计算标准正则功率曲线
    PCurveNorm = calculate_normalized_power_curve(IdealCurve, VRated, PRated)

    # 保存结果
    save_results(PCurve, PCurveNorm, EPKW, EPPer)
    # 绘制结果
    plot_results(PVBad, PVLimit, PVDot, PCurve, IdealCurve.to_numpy())


if __name__ == "__main__":
    main()