水轮机稳定工况识别与相关性分析

library(tidyverse)
library(lubridate)

预处理

读入初始数据：

inp <- read_delim('wqx.csv', delim = ',')

Parsed with column specification:
cols(
  time = col_datetime(format = ""),
  hydrahead = col_double(),
  activepwr = col_double(),
  vibz = col_double()
)

设置参数，剔除无效数据：

# 有效功率下限，单位：mW
min_pwr_thr <- 0

# 数据采样周期
# sample_intv <- '5 mins'
sample_intv <- 'minute'

# 有效功率波动最大范围
# 超过此范围则认为工况发生了变化，单位：mW
step_thr <- 10
# 最小有效窗口宽度：20分钟
min_win_width <- 20

inpf <- inp %>% filter(activepwr > min_pwr_thr)

按照指定时间间隔生成采样数据：

inpf$intv <- floor_date(inpf$time, sample_intv)
spd <- inpf %>%
  group_by(intv) %>%
  summarise(hyd = median(hydrahead), atp = median(activepwr), vib = median(vibz))

数据分析和可视化

有效功率的时序分布和密度图：

ggplot(spd, aes(x = intv, y = atp)) +
  geom_point()

ggplot(spd, aes(atp)) + geom_density()

整合所有稳定窗口，计算每个特征在稳定运行期内的统计特征：

source('../filter_stable_wins.R')
stables <- filter_stable(spd, 'atp', step_thr, min_win_width)
summary(stables)

      intv                          hyd             atp                 vib          
 Min.   :2019-12-10 17:45:00   Min.   :40.16   Min.   :  0.05765   Min.   :  0.2510  
 1st Qu.:2019-12-17 17:51:30   1st Qu.:44.41   1st Qu.:  0.10207   1st Qu.:  0.3231  
 Median :2019-12-25 13:52:00   Median :45.43   Median :  0.12152   Median :  0.3851  
 Mean   :2019-12-25 13:41:43   Mean   :45.94   Mean   : 40.17701   Mean   :  9.1524  
 3rd Qu.:2020-01-02 03:30:30   3rd Qu.:48.13   3rd Qu.:  0.17324   3rd Qu.:  7.0675  
 Max.   :2020-01-09 16:13:00   Max.   :52.98   Max.   :252.73363   Max.   :206.7480  

三个特征的分布和相关性：

stables %>% select(-intv) %>% PerformanceAnalytics::chart.Correlation()

统一时序图：

scaled <- stables %>% select(-intv) %>% scale()
longtbl <- reshape2::melt(scaled)
ggplot(longtbl, aes(x = Var1, y = value, color = Var2)) + geom_line()

最长稳定窗口分析

stables %>% select(-intv) %>% cor()  # 3个特征的相关性

           hyd        atp        vib
hyd  1.0000000 -0.5111608 -0.3216112
atp -0.5111608  1.0000000  0.5224529
vib -0.3216112  0.5224529  1.0000000

windows <- max_windows(spd$atp, 10, 4)
longest <- spd[which.max(windows) : (which.max(windows) + max(windows) - 1), ]
ggplot(longest, aes(x = intv, y = hyd)) + geom_line()

ggplot(longest, aes(x = intv, y = atp)) + geom_line()

ggplot(longest, aes(x = intv, y = vib)) + geom_line()

两两特征散点图：

ggplot(stables, aes(x = hyd, y = atp)) + geom_point()

ggplot(stables, aes(x = atp, y = vib)) + geom_point()

ggplot(stables, aes(x = hyd, y = vib)) + geom_point()

输出稳定工况数据：

write_csv(stables, "wqx_stables.csv")