WTOAAM/dataAnalysisBusiness/demo/SCADA_10min_category_0.py

# -*- coding: utf-8 -*-
"""
Created on Mon Apr  8 15:01:43 2024

@author: LDDN
"""
import math
import pandas as pd  
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator#设定固定刻度


scada_10min = pd.read_csv(r'E:\BaiduNetdiskDownload\test\min_scada_LuoTuoGou\72\82.csv',encoding="utf-8")  #.value是将单元格
turbine_info = pd.read_csv(r'E:\BaiduNetdiskDownload\test\min_scada_LuoTuoGou\72\info.csv')  #.value是将单元格
PRated = turbine_info.loc[:,["额定功率"]] #2000
PRated = PRated.values
VCutOut = turbine_info.loc[:,["切出风速"]]  #25
VCutOut = VCutOut.values
VCutIn = turbine_info.loc[:,["切入风速"]]  #3
VCutIn = VCutIn.values
VRated = turbine_info.loc[:,["额定风速"]] #10
VRated = VRated.values

time_stamp = scada_10min.loc[:,['时间']] #dataframe
active_power = scada_10min.loc[:,['变频器电网侧有功功率']]
wind_speed = scada_10min.loc[:,['风速']]
LM = pd.concat([time_stamp,active_power,wind_speed],axis=1)  #dataframe


Labeled_March809 = LM
APower = Labeled_March809["变频器电网侧有功功率"]  #series读入有功功率
WSpeed = Labeled_March809["风速"]  #读入风速
maxP=np.max(APower)
intervalP=25  #ceil(PRated*0.01)#功率分区间隔为额定功率的1%
intervalwindspeed=0.25  #风速分区间隔0.25m/s

#初始化
PNum = 0  
TopP = 0   
# 根据条件计算PNum和TopP  
if maxP >= PRated:  
    PNum = math.floor(maxP / intervalP) + 1  
    TopP = math.floor((maxP - PRated) / intervalP) + 1  
else:  
    PNum = math.floor(PRated / intervalP)  
    TopP = 0   
VNum = math.ceil(VCutOut / intervalwindspeed)  
  
SM1 = Labeled_March809.shape
AA1 = SM1[0]  
lab = [[0] for _ in range(AA1)]
lab = pd.DataFrame(lab,columns=['lab'])
Labeled_March809 = pd.concat([Labeled_March809,lab],axis=1)  #在tpv后加一列标签列
Labeled_March809 = Labeled_March809.values
SM = Labeled_March809.shape #(52561,4)
AA = SM[0]  
#存储功率大于0的运行数据
#标识功率为0的点，标识-1
DzMarch809_0 = np.zeros((AA, 3)) # 初始化数组来存储功率大于零的运行数据  
nCounter1 = 0 
Point_line = np.zeros(AA, dtype=int)  
#考虑到很多功率小于10的数据存在，将<10的功率视为0
for i in range(AA):
    if (APower[i] > 10) & (WSpeed[i] > 0):
        nCounter1 += 1   #共有nCounter1个功率大于0的正常数据
        DzMarch809_0[nCounter1-1, 0] = WSpeed[i]  
        DzMarch809_0[nCounter1-1, 1] = APower[i]  
        Point_line[nCounter1-1] = i+1  # 记录nCounter1记下的数据在原始数据中的位置  
    if APower[i] <= 10: 
        Labeled_March809[i,SM[1]-1] = -1  # 功率为0标识为-1  array类型
# 截取DzMarch809_0中实际存储的数据  其他全为0
DzMarch809 = DzMarch809_0[:nCounter1, :]  
#统计各网格落入的散点个数
XBoxNumber = np.ones((PNum, VNum),dtype=int)  #(86 100)
nWhichP = 0
nWhichV = 0

# 循环遍历DzMarch809中的有效数据  
for i in range(nCounter1):  
    
    # 查找功率所在的区间  
    for m in range(1, PNum + 1):  # 注意Python的range是左闭右开的，所以需要+1  
        if (DzMarch809[i,1] > (m - 1) * intervalP) and (DzMarch809[i,1] <= m * intervalP):  
            nWhichP = m  
            break  
      
    # 查找风速所在的区间  
    for n in range(1, VNum + 1):  # 同样需要+1  
        if (DzMarch809[i, 0] > (n - 1)*intervalwindspeed) and (DzMarch809[i, 0] <= n*intervalwindspeed):  
            nWhichV = n  
            break  
      
    # 如果功率和风速都在有效区间内，增加对应网格的计数  
    if (nWhichP > 0) and (nWhichV > 0):  
        XBoxNumber[nWhichP - 1, nWhichV - 1] += 1  # 注意Python的索引是从0开始的，所以需要减1  
# XBoxNumber现在包含了每个网格的计数[PNum行, VNum列]

for m in range(1,PNum+1):
    for n in range(1,VNum+1):
        XBoxNumber[m-1,n-1] = XBoxNumber[m-1,n-1] - 1

#在功率方向将网格内散点绝对个数转换为相对百分比，备用
PBoxPercent = np.zeros((PNum, VNum),dtype = float)  #(86 100) #计算后会出现浮点型，所以不能定义int类型
PBinSum = np.zeros((PNum,1),dtype=int)
for i in range(1,PNum+1):
    for m in range(1,VNum+1):
        PBinSum[i-1] = PBinSum[i-1] + XBoxNumber[i-1,m-1] 
    for m in range(1,VNum+1):
        if PBinSum[i-1]>0:
            PBoxPercent[i-1,m-1] = (XBoxNumber[i-1,m-1] / PBinSum[i-1])*100
#在风速方向将网格内散点绝对个数转换为相对百分比，备用          
VBoxPercent = np.zeros((PNum, VNum))  #(86 100) #计算后会出现浮点型，所以不能定义int类型
VBinSum = np.zeros((VNum,1),dtype=int)
for i in range(1,VNum+1):
    for m in range(1,PNum+1):
        VBinSum[i-1] = VBinSum[i-1] + XBoxNumber[m-1,i-1] 
    for m in range(1,PNum+1):
        if VBinSum[i-1]>0:
            VBoxPercent[m-1,i-1] = (XBoxNumber[m-1,i-1] / VBinSum[i-1])*100
# VBoxPercent PBoxPercent 左上-右下
# 将数据颠倒一下  左下-右上         第一行换为倒数第一行 方便可视化
InvXBoxNumber = np.zeros((PNum,VNum),dtype = int)
InvPBoxPercent = np.zeros((PNum,VNum),dtype = float)
InvVBoxPercent = np.zeros((PNum,VNum),dtype = float)
for m in range(1,PNum+1):
    for n in range(1,VNum+1):
        InvXBoxNumber[m-1,n-1] = XBoxNumber[PNum-(m-1)-1,n-1]
        InvPBoxPercent[m-1,n-1] = PBoxPercent[PNum-(m-1)-1,n-1]
        InvVBoxPercent[m-1,n-1] = VBoxPercent[PNum-(m-1)-1,n-1]

#以水平功率带方向为准，分析每个水平功率带中，功率主带中心，即找百分比最大的网格位置。
PBoxMaxIndex = np.zeros((PNum,1),dtype = int)  #水平功率带最大网格位置索引
PBoxMaxP = np.zeros((PNum,1),dtype = float)       #水平功率带最大网格百分比
for m in range(1,PNum+1):
    PBoxMaxIndex[m-1] = np.argmax(PBoxPercent[m-1, :])   #argmax返回最大值的索引
    PBoxMaxP[m-1] = np.max(PBoxPercent[m-1, :])
#以垂直风速方向为准，分析每个垂直风速带中，功率主带中心，即找百分比最大的网格位置。
VBoxMaxIndex = np.zeros((VNum,1),dtype = int)  
VBoxMaxV = np.zeros((VNum,1),dtype = float)       
for m in range(1,VNum+1):
    VBoxMaxIndex[m-1] = np.argmax(VBoxPercent[:, m-1])   
    VBoxMaxV[m-1] = np.max(VBoxPercent[:, m-1])

#切入风速特殊处理，如果切入风速过于偏右，向左拉回
if PBoxMaxIndex[0]>14:                     #第一个值对应的是风速最小处 即切入风速
    PBoxMaxIndex[0] = 9 
#以水平功率带方向为基准，进行分析
DotDense = np.zeros(PNum)   #每一水平功率带的功率主带包含的网格数
DotDenseLeftRight = np.zeros((PNum,2))  #存储每一水平功率带的功率主带以最大网格为中心，向左，向右扩展的网格数
DotValve = 90  #从中心向左右对称扩展网格的散点百分比和的阈值。
PDotDenseSum = 0
for i in range(PNum - TopP):  # 从最下层水平功率带开始，向上分析到特定的功率带  
    PDotDenseSum = PBoxMaxP[i]  # 以中心最大水平功率带为基准，向左向右对称扩展网格，累加各网格散点百分比  
    iSpreadRight = 1  
    iSpreadLeft = 1  
      
    while PDotDenseSum < DotValve:  
        if (PBoxMaxIndex[i] + iSpreadRight) < VNum-1-1:  
            PDotDenseSum += PBoxPercent[i, PBoxMaxIndex[i] + iSpreadRight]  # 向右侧扩展  
            iSpreadRight += 1  
        else:
            break  
          
        if (PBoxMaxIndex[i] - iSpreadLeft) > 0:  
            PDotDenseSum += PBoxPercent[i, PBoxMaxIndex[i] - iSpreadLeft]  # 向左侧扩展  
            iSpreadLeft += 1  
        else:  
            break  
    iSpreadRight = iSpreadRight-1
    iSpreadLeft = iSpreadLeft-1
    #向左右扩展完毕
    DotDenseLeftRight[i, 0] = iSpreadLeft  # 左  
    DotDenseLeftRight[i, 1] = iSpreadRight  # 右  
    DotDense[i] = iSpreadLeft + iSpreadRight + 1  # 记录向左向右扩展的个数及每个功率仓内网格的个数  
# 此时DotDense和DotDenseLeftRight数组已经包含了所需信息    
#各行功率主带右侧宽度的中位数最具有代表性（因为先右后左）
DotDenseWidthLeft = np.zeros((PNum-TopP))
for i in range(PNum-TopP):
    DotDenseWidthLeft[i] = DotDenseLeftRight[i,1]  #DotDenseLeftRight[i,1]:向右延伸个数
MainBandRight = np.median(DotDenseWidthLeft) #计算中位数

# 初始化变量  
PowerLimit = np.zeros(PNum, dtype=int)  # 各水平功率带是否为限功率标识，1：是；0：不是  
WidthAverage = 0  # 功率主带右侧平均宽度  
WidthAverage_L = 0  # 功率主带左侧平均宽度  
WidthVar = 0  # 功率主带方差（此变量在提供的代码中并未使用）  
PowerLimitValve = 6  # 限功率主带判别阈值  
N_Pcount = 20  # 阈值  
  
nCounterLimit = 0  # 限功率的个数  
nCounter = 0  # 正常水平功率带的个数  
  
# 循环遍历水平功率带，从第1个到第PNum-TopP个  
for i in range(PNum - TopP):  
    # 如果向右扩展网格数大于阈值，且该水平功率带点总数大于20，则标记为限功率带  
    if (DotDenseLeftRight[i, 1] > PowerLimitValve) and (PBinSum[i] > N_Pcount):  
        PowerLimit[i] = 1  
        nCounterLimit += 1  #限功率的个数
      
    # 如果向右扩展网格数小于等于阈值，则累加右侧宽度（左侧宽度在代码中似乎有误）  
    if DotDenseLeftRight[i, 1] <= PowerLimitValve:  
        WidthAverage += DotDenseLeftRight[i, 1]  # 统计正常水平功率带右侧宽度
        WidthAverage_L += DotDenseLeftRight[i,1]   #统计正常水平功率带左侧宽度
        nCounter += 1  
# 计算平均宽度  
WidthAverage /= nCounter if nCounter > 0 else 1  # 避免除以0的情况  
WidthAverage_L /= nCounter if nCounter > 0 else 1   

#计算正常（即非限功率）水平功率带的功率主带宽度的方差，以此来反映从下到上宽度是否一致
WidthVar = 0  # 功率主带宽度的方差   
for i in range(PNum - TopP):  
    # 如果向右扩展网格数小于等于阈值，则计算当前宽度与平均宽度的差值平方  
    if DotDenseLeftRight[i, 1] <= PowerLimitValve:  
        WidthVar += (DotDenseLeftRight[i, 1] - WidthAverage) ** 2  
# 计算方差（注意：除以nCounter-1是为了得到样本方差）  
WidthVar = np.sqrt(WidthVar / (nCounter - 1) if nCounter > 1 else 0)  # 避免除以0的情况

#各水平功率带，功率主带的风速范围，右侧扩展网格数*2*0.25
PowerBandWidth = WidthAverage*intervalwindspeed+WidthAverage_L*intervalwindspeed

# 对限负荷水平功率带的最大网格进行修正  
for i in range(1, PNum - TopP+1):  
    if (PowerLimit[i] == 1) and (abs(PBoxMaxIndex[i] - PBoxMaxIndex[i - 1]) > 5):  
        PBoxMaxIndex[i] = PBoxMaxIndex[i - 1] + 1  
  
# 输出各层功率主带的左右边界网格索引  
DotDenseInverse = np.flipud(DotDenseLeftRight)  # 上下翻转数组以得到反向顺序  
  
# 计算功率主带的左右边界  
CurveWidthR = np.ceil(WidthAverage) + 2  # 功率主带的右边界 + 2  
CurveWidthL = np.ceil(WidthAverage_L) + 2  # 功率主带的左边界 + 2  
  
# 网格是否为限功率网格的标识数组  
BBoxLimit = np.zeros((PNum, VNum), dtype=int)  
# 标记限功率网格  
for i in range(2, PNum - TopP):  
    if PowerLimit[i] == 1:
        BBoxLimit[i, int(PBoxMaxIndex[i] + CurveWidthR + 1):VNum] = 1

# 初始化数据异常需要剔除的网格标识数组  
BBoxRemove = np.zeros((PNum, VNum), dtype=int)  
# 标记需要剔除的网格  
for m in range(PNum - TopP): 
    for n in range(int(PBoxMaxIndex[m]) + int(CurveWidthR), VNum):  # 注意Python中的索引从0开始，因此需要减去1  
        BBoxRemove[m, n] = 1 
    # 功率主带左侧的超发网格，从最大索引向左直到第一个网格  
    
    for n in range(int(PBoxMaxIndex[m]) - int(CurveWidthL)+1, 0, -1):  # 使用range的步长参数来实现从右向左的迭代  
        BBoxRemove[m, n-1] = 2  # 注意Python中的索引从0开始，因此需要减去1

# 初始化变量  
CurveTop = np.zeros((2, 1), dtype=int)  
CurveTopValve = 1  # 网格的百分比阈值  
BTopFind = 0  
mm = 0  
#确定功率主带的左上拐点，即额定风速位置的网格索引
CurveTop = np.zeros((2, 1), dtype=int)  
CurveTopValve = 1  # 网格的百分比阈值  
BTopFind = 0  
mm = 0   
for m in range(PNum - TopP, 0, -1):  # 注意Python的range是左闭右开区间，所以这里从PNum-TopP开始到1（不包括0）  
    for n in range(int(np.floor(int(VCutIn) / intervalwindspeed)), VNum - 1):  # 使用floor函数来向下取整  
        if (VBoxPercent[m, n - 1] < VBoxPercent[m, n]) and (VBoxPercent[m, n] <= VBoxPercent[m, n + 1]) and (XBoxNumber[m, n] >= 3):   
            CurveTop[0] = m  
            CurveTop[1] = n  #[第80个,第40个]
            BTopFind = 1  
            mm = m  # mm是拐点所在功率仓，对应其index
            break  # 找到后退出内层循环  
    if BTopFind == 1:  
        break  # 找到后退出外层循环
        
IsolateValve = 3  #功率主带右侧孤立点占比功率仓阈值 3%
# 遍历功率仓和网格  
for m in range(PNum - TopP):    
    for n in range(int(PBoxMaxIndex[m]) + int(CurveWidthR), VNum):  
        # 检查PBoxPercent是否小于阈值，如果是，则标记BBoxRemove为1  
        if PBoxPercent[m, n] < IsolateValve:   
            BBoxRemove[m, n] = 1
#功率主带顶部宽度
CurveWidthT = np.floor((maxP - PRated) / intervalP) + 1  
# 标记额定功率以上的超发点(PNum-PTop之间)  
for m in range(PNum - TopP, PNum):   
    for n in range(VNum):  
        BBoxRemove[m, n] = 3
  
# 标记功率主带拐点左侧的欠发网格  
for m in range(mm-1, PNum - TopP): 
    for n in range(int(CurveTop[1]) - 2):
        BBoxRemove[m, n] = 2    # BBoxRemove数组现在包含了根据条件标记的超发点和欠发网格的信息

#以网格的标识，决定该网格内数据的标识。
# DzMarch809Sel数组现在包含了每个数据点的标识
DzMarch809Sel = np.zeros(nCounter1, dtype=int)  # 初始化标识数组   
nWhichP = 0  
nWhichV = 0  
nBadA = 0   
for i in range(nCounter1):  
    for m in range( PNum ):   
        if (DzMarch809[i, 1] > m * intervalP) and (DzMarch809[i, 1] <= (m+1) * intervalP):  
            nWhichP = m  #m记录的是index
            break  
    for n in range( VNum ):  # 注意Python的range是左闭右开区间，所以这里到VNum+1  
        if DzMarch809[i, 0] > ((n+1) * intervalwindspeed - intervalwindspeed/2) and DzMarch809[i, 0] <= ((n+1) * intervalwindspeed + intervalwindspeed / 2):  
            nWhichV = n  #index
            break  
    if nWhichP >= 0 and nWhichV >= 0:  
        if BBoxRemove[nWhichP, nWhichV] == 1:   
            DzMarch809Sel[i] = 1  
            nBadA += 1  
        elif BBoxRemove[nWhichP, nWhichV] == 2:  
            DzMarch809Sel[i] = 2  
        elif BBoxRemove[nWhichP , nWhichV] == 3:  
            DzMarch809Sel[i] = 0  # 额定风速以上的超发功率点认为是正常点，不再标识  
# DzMarch809Sel数组现在包含了每个数据点的标识

##############################滑动窗口方法
# 存储限负荷数据  
PVLimit = np.zeros((nCounter1, 3))  #存储限负荷数据  %第3列用于存储限电的点所在的行数
nLimitTotal = 0  
nWindowLength = 6   #滑动窗口长度设置为6
LimitWindow = np.zeros(nWindowLength)  #滑动窗口空列表
UpLimit = 0    #上限
LowLimit = 0   #下限
PowerStd = 30  # 功率波动方差  
nWindowNum = np.floor(nCounter1/nWindowLength) #6587
PowerLimitUp = PRated - 100  
PowerLimitLow = 100  

# 循环遍历每个窗口  
for i in range(int(nWindowNum)):  
    start_idx = i * nWindowLength  
    end_idx = start_idx + nWindowLength  
    LimitWindow = DzMarch809[start_idx:end_idx, 1]  # 提取当前窗口的数据  
      
    # 检查窗口内所有数据是否在功率范围内  
    bAllInAreas = np.all(LimitWindow >= PowerLimitLow) and np.all(LimitWindow <= PowerLimitUp)  
    if not bAllInAreas:  
        continue  
      
    # 计算方差上下限  
    UpLimit = LimitWindow[0] + PowerStd  
    LowLimit = LimitWindow[0] - PowerStd  
      
    # 检查窗口内数据是否在方差范围内  
    bAllInUpLow = np.all(LimitWindow >= LowLimit) and np.all(LimitWindow <= UpLimit)  
    if bAllInUpLow:  
        # 标识窗口内的数据为限负荷数据  
        DzMarch809Sel[start_idx:end_idx] = 4  
          
        # 存储限负荷数据  
        for j in range(nWindowLength):  
            PVLimit[nLimitTotal, :2] = DzMarch809[start_idx + j, :2]  
            PVLimit[nLimitTotal, 2] = Point_line[start_idx + j]  # 对数据进行标识  
            nLimitTotal += 1  
# PVLimit现在包含了限负荷数据，nLimitTotal是限负荷数据的总数


#将功率滑动窗口主带平滑化
# 初始化锯齿平滑的计数器  
nSmooth = 0  
# 遍历除了最后 TopP+1 个元素之外的所有 PBoxMaxIndex 元素  
for i in range(PNum - TopP - 1):  
    PVLeftDown = np.zeros(2)  
    PVRightUp = np.zeros(2)  
    # 检查当前与下一个 PBoxMaxIndex 之间的距离是否大于等于1  
    if PBoxMaxIndex[i + 1] - PBoxMaxIndex[i] >= 1:  
        # 计算左下和右上顶点的坐标  
        PVLeftDown[0] = (PBoxMaxIndex[i]+1 + CurveWidthR) * 0.25 - 0.125  
        PVLeftDown[1] = (i) * 25  
        PVRightUp[0] = (PBoxMaxIndex[i+1]+1 + CurveWidthR) * 0.25 - 0.125  
        PVRightUp[1] = (i+1) * 25  
          
        # 遍历 DzMarch809 数组  
        for m in range(nCounter1):  
            # 检查当前点是否在锯齿区域内  
            if (DzMarch809[m, 0] > PVLeftDown[0]) and (DzMarch809[m, 0] < PVRightUp[0]) and (DzMarch809[m, 1] > PVLeftDown[1]) and (DzMarch809[m, 1] < PVRightUp[1]):
                # 检查斜率是否大于对角连线  
                if (DzMarch809[m, 1] - PVLeftDown[1]) / (DzMarch809[m, 0] - PVLeftDown[0]) > (PVRightUp[1] - PVLeftDown[1]) / (PVRightUp[0] - PVLeftDown[0]):
                    # 如果在锯齿左上三角形中，则选中并增加锯齿平滑计数器  
                    DzMarch809Sel[m] = 0  
                    nSmooth += 1  
# DzMarch809Sel 数组现在＋包含了锯齿平滑的选择结果，nSmooth 是选中的点数


###################################存储数据
# 存储好点  
nCounterPV = 0  # 初始化计数器  
PVDot = np.zeros((nCounter1, 3))  # 初始化存储好点的数组  nCounter1是p>0的数
for i in range(nCounter1):  
    if DzMarch809Sel[i] == 0:  
        nCounterPV += 1  
        PVDot[nCounterPV-1, :2] = DzMarch809[i, :2]  
        PVDot[nCounterPV-1, 2] = Point_line[i]  # 好点 Point_line记录nCounter1在原始数据中的位置 
nCounterVP = nCounterPV  
 
# 对所有数据中的好点进行标注    
for i in range(nCounterVP):  
    Labeled_March809[int(PVDot[i, 2] - 1), (SM[1]-1)] = 1  # 注意Python的索引从0开始，并且需要转换为整数索引  
 
# 存储坏点  
nCounterBad = 0  # 初始化计数器  
PVBad = np.zeros((nCounter1, 3))  # 初始化存储坏点的数组  
for i in range(nCounter1):  
    if DzMarch809Sel[i] in [1, 2, 3]:  
        nCounterBad += 1  
        PVBad[nCounterBad-1, :2] = DzMarch809[i, :2]  
        PVBad[nCounterBad-1, 2] = Point_line[i]  
    
# 对所有数据中的坏点进行标注  
for i in range(nCounterBad):  
    Labeled_March809[int(PVBad[i, 2] - 1),(SM[1]-1)] = 5  # 坏点标识  

# 对所有数据中的限电点进行标注   
for i in range(nLimitTotal):  
    Labeled_March809[int(PVLimit[i, 2] - 1),(SM[1]-1)] = 4  # 限电点标识  
# 对所有的数据点进行标注  
# Labeled_March809是array，提取所第四列的值保存为dataframe
A = Labeled_March809[:,3]
A=pd.DataFrame(A,columns=['lab'])


mergedTable = pd.concat([scada_10min,A],axis=1)#合并dataframe
D = mergedTable[mergedTable['lab'] == 1]#选择为1的行

ws = D["风速"].values  #array
ap = D["变频器电网侧有功功率"]

# fig=plt.figure(figsize=(10,6),dpi=500)  #figsize是图形大小，dpi像素
fig=plt.figure()  #figsize是图形大小，dpi像素
plt.scatter(ws,ap,s=1,c='black',marker='.') #'.'比'o'要更小

# plt.scatter(x2,y2,s=10,c='b',marker='.',label='5.8-6.5建模噪声点')

x_major_locator=MultipleLocator(5)
y_major_locator=MultipleLocator(500)
ax=plt.gca()
ax.xaxis.set_major_locator(x_major_locator)
ax.yaxis.set_major_locator(y_major_locator)
plt.xlim((0,30))
plt.ylim((0,2200))
plt.tick_params(labelsize=8)

# plt.grid(c='dimgray',alpha=0.2)

plt.xlabel("V/(m$·$s$^{-1}$)",fontsize=8)
plt.ylabel("P/kW",fontsize=8)

# plt.savefig(r'D:\赵雅丽\研究生学习资料\学习资料\劣化度健康度\spyder\大论文\图\风速-功率.jpg',bbox_inches='tight')
plt.show()