一、算法思想

用当前时间点和前面两点建立模型，并且这个模型可以用多个岭函数（平滑曲线）之和来表示。

二、建模过程

1.基本数据准备：将连续的实测功率值表示为如下形式，前两行为输入x，最后一行为输出y。
投影寻踪回归ppr - 图1
2.对输出y进行标准化处理，处理成y΄
3.找投影方向α
投影寻踪回归ppr - 图2
投影向量的行数为岭函数的个数，列数为输入x的维度。
找第一个投影方向(a,a)，这个投影方向并不是随机的，由输入x以及输出y的协方差（x和y的关系）和自方差（x和x的关系）之间的关系来确定确定权重，权重就是投影方向的系数，这样就找到第一个投影方向。（github上的α是随机的）。
4. 求岭函数f
有了第一个投影方向(a,a)之后，用所有的输入和这个投影方向作向量乘，将向量乘后的结果作为横坐标和y΄一一对应，这样就有了多对散点(αx, y)，拟合这些散点的过程就是求岭函数f的过程。求岭函数f首先是确定一个窗口，每个窗口里包含的点数是一致的，并且和光滑系数有关，每个窗口的点数=总点数0.5光滑系数+0.5（向下取整）。每次滑动的距离为一个窗口，每个窗口里包含多少点，就要平滑出多少点，最后的岭函数f其实就是这些平滑之后的散点。
5.求权重系数β
它等于协方差和除以自方差和，自方差和是岭函数的预测值的平方和，协方差和是岭函数和y΄的乘积和
6.找第二个岭函数的α，β，f
用第一轮每一个标准化后的实际值y΄和岭函数的预测值ŷ的偏差，用来作为下一轮的输出y
7. 迭代停止条件
（1）迭代次数>岭函数个数
（2）残差收敛：前一次的残差<后一次的残差

三、电场测试

投影寻踪回归ppr - 图3
江苏盐城大丰润龙
20天：0.579，30天：0.592，现在：0.564，邓老师（前）：0.704，邓老师：0.315，平滑：0.527，加入预测功率：0.654
5月调和平均准确率对比（20天数据建模）

电场名 ppr 现在
江苏大丰润龙 0.572 0.581
苜蓿台 0.606 0.688
望洋台 0.593 0.639
华电草湖二场 0.487 0.559

四、后续计划

处理方式 20天 30天 直接预测第8点 加入短期预测功率 平滑处理
调和平均准确率 0.579 0.592 0.515 0.653 0.514

复现邓老师的ppr算法基础上进行改进：
（1）参数测试，特别是建模时长和输入维数对超短期的准确率有一定的影响：
有提升，但是提升不太多。
（2）直接预测第8点：效果不好，没有直接对第一点进行建模效果好。
（3）和纯时序的超短期预测算法进行对比，看ppr算法是否是纯时序算法中效果最好的：由于现有的算法都加了短期预测功率，不是纯时序，从纯时序角度对比较困难。
（4）预测功率数据参与建模，方式是将这个预测功率数据加入到输入中：提升较多。
（5）数据预处理一：补断点，处理得比较简单，两点之间的均值。
（6）数据预处理二：实测数据进行平滑处理，避免波动较大的数据干扰模型:准确率比没有平滑更低，有可能和采用的平滑方式效果不佳有关，目前采用的是卷积平滑（np.convolve）。

对github上的ppr算法的改进
（7）测试github的ppr算法进行迭代训练、预测和加上短期预测功率的效果。
（8）对调和平均准确率异常的分析，华电草湖二场（652112）5月12日调和平均准确率0.088，现在的算法是0.272。

github的代码

  1. # -*- coding: utf-8 -*-
  2. import numpy
  3. from scipy.interpolate import UnivariateSpline
  4. from sklearn.base import BaseEstimator, RegressorMixin, TransformerMixin
  5. from sklearn.utils.validation import check_X_y, check_is_fitted
  6. from sklearn.utils import as_float_array, check_random_state, check_array
  7. from matplotlib import pyplot
  8. import numpy as np
  9. import MySQLdb
  10. import pandas as pd
  11. import datetime
  12. import time
  15. def mysql_conn(db='exportdb'):  # 数据库连接信息
  17.     user = 'dbadmin'
  18.     pwd = 'Goldwind@123'
  19.     ip = '10.12.17.97'
  20.     port = 3306
  21.     conn_modeling = MySQLdb.connect(host=ip, db=db, user=user, passwd=pwd, port=port, charset="utf8")
  23.     return conn_modeling
  26. def get_data(start_, end_):  # 获取实测数据
  28.     export_data = pd.read_sql("select rectime, power, predict_power from `%s` where rectime>='%s' and rectime<='%s'" %
  29.                               ('factory_' + wfid, start_, end_), mysql_conn())
  30.     return export_data
  33. def get_expand_data(data, dimension_):  # 将一维数据变为多维,最后一维为输出,前面为输入
  35.     refactor_sam = []
  36.     for i in range(dimension_ - 1, len(data)):
  37.         block = dimension_ * [0]
  38.         for j in range(dimension_):
  39.             block[j] = data[i]
  40.             i -= 1
  41.         block = list(reversed(block))
  42.         refactor_sam.append(block)
  43.     col = []
  44.     for i in range(dimension_):
  45.         if i != dimension_ - 1:
  46.             col.append(data.name + '_%s' % i)
  47.         else:
  48.             col.append('output')
  49.     expand_df = pd.DataFrame(refactor_sam, columns=col)
  50.     return expand_df
  53. def update_input(predict_result_, array_input_):  # 对输入值进行更新
  55.     array_input_ = np.delete(array_input_, 0)  # 删除掉numpy.ndarray的第一个值
  56.     array_input_ = np.append(array_input_, predict_result_)  # 插入预测值
  57.     return array_input_
  60. def rmse(real, pred, cap_=49500.0):  # 均方根误差准确率
  61.     s = 1 - np.sqrt(np.mean((real - pred) ** 2)) / cap_
  62.     return s * 100
  65. def harmonic_mean(sc, pre, cap_):  # 计算超短期准确率
  66.     index = (sc < 0.03 * cap_) & (pre < 0.03 * cap_)
  67.     sc = sc[~index]
  68.     pre = pre[~index]
  69.     if len(sc) > 0:
  70.         a = np.abs(sc / (sc + pre) - 0.5)
  71.         b = abs(sc - pre) / sum(abs(sc - pre))
  72.         th = 1 - 2 * sum(a * b)
  73.     else:
  74.         th = 1
  75.     return th
  78. def get_predict_data(ptime):
  79.     predict_power = pd.read_sql("select predict_power from `%s` where rectime='%s'" %('factory_' + wfid, ptime),
  80.                                 mysql_conn()).iloc[0, 0]
  81.     return predict_power
  84. def smooth_data(df):  # 数据平滑
  85.     x = df['power'].values
  86.     box = np.ones(10) / 10
  87.     y = np.convolve(x, box, mode='same')  # 两个数组长度可以不一样,参数为same时取更长的一个
  88.     y[:5], y[-5:] = x[:5], x[-5:]
  89.     df['new_power'] = y
  90.     return df
  93. def compensate_data(df, start, end):  # 补缺值
  94.     time_list = pd.date_range(start, end, freq='15min').tolist()
  95.     time_df = pd.DataFrame({'rectime': time_list})  # 首先构造一个时间完整的df
  96.     merge_df = pd.merge(time_df, df, on='rectime', how='outer')  # merge外链接
  97.     merge_df = merge_df.reset_index(drop=True)
  98.     for m in range(len(merge_df)):
  99.         if m == 0 and merge_df.iloc[0, 1] == np.nan:
  100.             merge_df.iloc[0, 1] = (merge_df.iloc[1, 1]+merge_df.iloc[2, 1])/2
  101.         elif m == len(merge_df) and merge_df.iloc[-1, 1] == np.nan:
  102.             merge_df.iloc[-1, 1] = (merge_df.iloc[-2, 1]+merge_df.iloc[-3, 1])/2
  103.         elif merge_df.iloc[m, 1] == np.nan:
  104.             merge_df.iloc[m, 1] = (merge_df.iloc[m+1, 1]+merge_df.iloc[m-1, 1])/2
  105.     # merge_df = merge_df.fillna(merge_df.mean())  # NAN值用平均来补
  106.     return merge_df
  109. class ProjectionPursuitRegressor(BaseEstimator, TransformerMixin, RegressorMixin):
  110.     def __init__(self, r=10, fit_type='polyfit', degree=3, opt_level='high',
  111.                  example_weights='uniform', out_dim_weights='inverse-variance',
  112.                  eps_stage=0.0001, eps_backfit=0.01, stage_maxiter=100,
  113.                  backfit_maxiter=10, random_state=None, show_plots=False,
  114.                  plot_epoch=50):
  116.         if not isinstance(r, int) or r < 1:  # 如果r不是int或者r<1,报错
  117.             raise ValueError('r must be an int >= 1.')
  118.         if fit_type not in ['polyfit', 'spline']:
  119.             raise NotImplementedError('fit_type ' + fit_type + ' not supported')
  120.         if not isinstance(degree, int) or degree < 1:
  121.             raise ValueError('degree must be >= 1.')
  122.         if opt_level not in ['low', 'medium', 'high']:
  123.             raise ValueError('opt_level must be either low, medium, or high.')
  124.         if not (isinstance(example_weights, str) and example_weights == 'uniform'):
  125.             try:
  126.                 example_weights = as_float_array(example_weights)
  127.             except (TypeError, ValueError) as error:
  128.                 raise ValueError('example_weights must be `uniform`, or array-like.')
  129.             if numpy.any(example_weights < 0):
  130.                 raise ValueError('example_weights can not contain negatives.')
  131.         if not (isinstance(out_dim_weights, str) and out_dim_weights in
  132.                 ['inverse-variance', 'uniform']):
  133.             try:
  134.                 out_dim_weights = as_float_array(out_dim_weights)
  135.             except (TypeError, ValueError) as error:
  136.                 raise ValueError('out_dim_weights must be either ' + \
  137.                                  'inverse-variance, uniform, or array-like.')
  138.             if numpy.any(out_dim_weights < 0):
  139.                 raise ValueError('out_dim_weights can not contain negatives.')
  140.         if eps_stage <= 0 or eps_backfit <= 0:
  141.             raise ValueError('Epsilons must be > 0.')
  142.         if not isinstance(stage_maxiter, int) or stage_maxiter <= 0 or \
  143.                 not isinstance(backfit_maxiter, int) or backfit_maxiter <= 0:
  144.             raise ValueError('Maximum iteration settings must be ints > 0.')
  146.         # Save parameters to the object
  147.         params = locals()
  148.         for k, v in params.items():
  149.             if k != 'self':
  150.                 setattr(self, k, v)
  152.     def transform(self, X):
  154.         check_is_fitted(self, '_alpha')
  155.         X = check_array(X)
  156.         return numpy.dot(X, self._alpha)
  158.     def predict(self, X):
  160.         P = self.transform(X)
  161.         Y = sum([numpy.outer(self._f[j](P[:, j]), self._beta[:, j]) for j in range(self.r)])
  162.         return Y if Y.shape[1] != 1 else Y[:, 0]
  164.     def fit(self, X, Y):
  166.         X, Y = check_X_y(X, Y, multi_output=True)  # 标准估算器的输入验证,检查X和y的长度是否一致,强制X为2D和y为1D
  167.         if Y.ndim == 1:  # standardize Y as 2D so the below always works 判断Y的维度
  168.             Y = Y.reshape((-1, 1))  # reshape returns a view to existing data
  170.         self._random = check_random_state(self.random_state)  # 随机产生随机数, 没有限制
  172.         if isinstance(self.example_weights, str) and self.example_weights == 'uniform':
  173.             self._example_weights = numpy.ones(X.shape[0])
  174.         elif isinstance(self.example_weights, numpy.ndarray):
  175.             if X.shape[0] != self.example_weights.shape[0]:
  176.                 raise ValueError('example_weights provided to the constructor' +
  177.                                  ' have dimension ' + str(self.example_weights.shape[0]) +
  178.                                  ', which disagrees with the size of X: ' + str(X.shape[0]))
  179.             else:
  180.                 self._example_weights = self.example_weights
  182.         if isinstance(self.out_dim_weights, str) and self.out_dim_weights == 'inverse-variance':
  183.             variances = Y.var(axis=0)
  184.             # There is a problem if a variance for any column is zero, because
  185.             # its relative scale will appear infinite. 
  186.             if max(variances) == 0:  # if all zeros, don't readjust weights
  187.                 variances = numpy.ones(Y.shape[1])
  188.             else:
  189.                 # Fill zeros with the max of variances, so the corresponding
  190.                 # columns have small weight and are not major determiners of loss.
  191.                 variances[variances == 0] = max(variances)
  192.             self._out_dim_weights = 1. / variances
  193.         elif isinstance(self.out_dim_weights, str) and \
  194.                 self.out_dim_weights == 'uniform':
  195.             self._out_dim_weights = numpy.ones(Y.shape[1])
  196.         elif isinstance(self.out_dim_weights, numpy.ndarray):
  197.             if Y.shape[1] != self.out_dim_weights.shape[0]:
  198.                 raise ValueError('out_dim_weights provided to the constructor' +
  199.                                  ' have dimension ' + str(self.out_dim_weights.shape[0]) +
  200.                                  ', which disagrees with the width of Y: ' + str(Y.shape[1]))
  201.             else:
  202.                 self._out_dim_weights = self.out_dim_weights
  204.         self._alpha = self._random.randn(X.shape[1], self.r)  # p x r
  205.         self._beta = self._random.randn(Y.shape[1], self.r)  # d x r
  206.         self._f = [lambda x: x * 0 for j in range(self.r)]  # zero functions
  207.         self._df = [None for j in range(self.r)]  # no derivatives yet
  209.         for j in range(self.r):  # for each term in the additive model
  210.             self._fit_stage(X, Y, j, True)
  212.             if self.opt_level == 'high':
  213.                 self._backfit(X, Y, j, True)
  214.             elif self.opt_level == 'medium':
  215.                 self._backfit(X, Y, j, False)
  216.         return self
  218.     def _fit_stage(self, X, Y, j, fit_weights):
  219.         P = self.transform(X)  # the n x r projections matrix
  220.         R_j = Y - sum([numpy.outer(self._f[t](P[:, t]), self._beta[:, t].T) for t
  221.                        in range(self.r) if t is not j])  # the n x d residuals matrix
  223.         # main alternating optimization loop
  224.         itr = 0  # iteration counter
  225.         prev_loss = -numpy.inf
  226.         loss = numpy.inf
  227.         p_j = P[:, j]  # n x 1, the jth column of the projections matrix
  228.         # 总的投影次数不能大于stage_maxiter,并且前后两次的损失差值不能小于eps_stage
  229.         while abs(prev_loss - loss) > self.eps_stage and itr < self.stage_maxiter:
  230.             # To understand how to optimize each set of parameters assuming the
  231.             # others remain constant, see math.pdf section 3.
  233.             # find the f_j
  234.             beta_j_w = self._out_dim_weights * self._beta[:, j]  # weighted beta
  235.             targets = numpy.dot(R_j, beta_j_w) / (
  236.                     numpy.inner(self._beta[:, j], beta_j_w) + 1e-9)
  237.             # Fit the targets against projections.
  238.             self._f[j], self._df[j] = self._fit_2d(p_j, targets, j, itr)
  240.             # find beta_j
  241.             f = self._f[j](p_j)  # Find the n x 1 vector of function outputs.
  242.             f_w = self._example_weights * f  # f weighted by examples
  243.             self._beta[:, j] = numpy.dot(R_j.T, f_w) / (numpy.inner(f, f_w) + 1e-9)
  245.             # find alpha_j
  246.             if fit_weights:
  247.                 # Find the part of the Jacobians that is common to all
  248.                 J = -(self._df[j](p_j) * numpy.sqrt(self._example_weights) * X.T).T
  249.                 JTJ = numpy.dot(J.T, J)
  250.                 A = sum([self._out_dim_weights[k] * (self._beta[k, j] ** 2) * JTJ
  251.                          for k in range(Y.shape[1])])
  252.                 # Collect all g_jk vectors in to a convenient matrix G_j
  253.                 G_j = R_j - numpy.outer(self._f[j](p_j), self._beta[:, j].T)
  254.                 b = -sum([self._out_dim_weights[k] * self._beta[k, j] *
  255.                           numpy.dot(J.T, G_j[:, k]) for k in range(Y.shape[1])])
  257.                 delta = numpy.linalg.lstsq(A.astype(numpy.float64, copy=False),
  258.                                            b.astype(numpy.float64, copy=False), rcond=-1)[0]
  259.                 # TODO implement halving step if the loss doesn't decrease with
  260.                 # this update.
  261.                 alpha = self._alpha[:, j] + delta
  262.                 # normalize to avoid numerical drift
  263.                 self._alpha[:, j] = alpha / numpy.linalg.norm(alpha)
  265.             # Recalculate the jth projection with new f_j and alpha_j
  266.             p_j = numpy.dot(X, self._alpha[:, j])
  268.             # Calculate mean squared error for this iteration
  269.             prev_loss = loss
  270.             # Subtract updated contribution of the jth term to get the
  271.             # difference between Y and the predictions. 
  272.             diff = R_j - numpy.outer(self._f[j](p_j), self._beta[:, j].T)
  273.             # square the difference, multiply rows by example weights, multiply
  274.             # columns by their weights, and sum to get the final loss
  275.             loss = numpy.dot(numpy.dot(self._example_weights, diff ** 2),
  276.                              self._out_dim_weights)
  277.             itr += 1
  279.     def _backfit(self, X, Y, j, fit_weights):
  280.         itr = 0
  281.         prev_loss = -numpy.inf
  282.         loss = numpy.inf
  283.         while abs(prev_loss - loss) > self.eps_backfit and itr < self.backfit_maxiter:
  284.             for t in range(j):
  285.                 self._fit_stage(X, Y, t, fit_weights)
  287.             prev_loss = loss
  288.             diff = Y - self.predict(X).reshape(Y.shape)
  289.             loss = numpy.dot(numpy.dot(self._example_weights, diff ** 2), self._out_dim_weights)
  290.             itr += 1
  292.     def _fit_2d(self, x, y, j, itr):
  293.         if self.fit_type == 'polyfit':
  294.             coeffs = numpy.polyfit(x.astype(numpy.float64, copy=False),
  295.                                    y.astype(numpy.float64, copy=False), deg=self.degree,
  296.                                    w=self._example_weights)
  297.             fit = numpy.poly1d(coeffs)
  298.             deriv = fit.deriv(m=1)
  299.         elif self.fit_type == 'spline':
  300.             order = numpy.argsort(x)
  301.             fit = UnivariateSpline(x[order], y[order], w=self._example_weights, k=self.degree, s=len(y) * y.var(),
  302.                                    ext=0)
  303.             deriv = fit.derivative(1)
  305.         if self.show_plots and (itr % self.plot_epoch == 0):
  306.             pyplot.scatter(x, y)
  307.             pyplot.title('stage ' + str(j) + ' iteration ' + str(itr))
  308.             pyplot.xlabel('projections')
  309.             pyplot.ylabel('residuals')
  310.             xx = numpy.linspace(min(x), max(x), 100)
  311.             yy = fit(xx)
  312.             pyplot.plot(xx, yy, 'g', linewidth=1)
  313.             pyplot.show()
  315.         return fit, deriv
  318. if __name__ == '__main__':
  320.         # 需要配置如下6个变量
  321.         wfid = '320907'
  322.         ridge_number = 20  # 岭函数个数
  323.         dimension = 5  # 需要的数据总维度,如果是二阶输入的话,就是2+8=10,后面8个为输出
  324.         train_data_sets = 96 * 30  # 训练的时候用到的输入和输出对数
  325.         predict_start = '2019-05-01 00:00:00'  # 要预测的开始时间
  326.         predict_points = 96*2  # 要预测的点数(一天96点)
  328.         start_time = time.time()
  329.         cap = pd.read_sql("select powercap from `wind_farm_info` where wfid='%s'"
  330.                           % wfid, mysql_conn('WindFarmdb')).iloc[0, 0]
  332.         train_end = datetime.datetime.strptime(str(predict_start), '%Y-%m-%d %H:%M:%S') + \
  333.                     datetime.timedelta(hours=-2)  # 超短期提前两小时预测当前的
  334.         # train_start由train_end和train_data_sets以及dimension决定,保证有192对输入输出
  335.         train_start = datetime.datetime.strptime(str(train_end), '%Y-%m-%d %H:%M:%S') + \
  336.                       datetime.timedelta(minutes=-15 * (train_data_sets + dimension - 2))
  338.         predict_result_dict = {}
  339.         for i in range(predict_points):  # 每次训练和预测1个点
  341.             new_train_start = datetime.datetime.strptime(str(train_start), '%Y-%m-%d %H:%M:%S') + \
  342.                               datetime.timedelta(minutes=15 * i)
  343.             new_train_end = datetime.datetime.strptime(str(train_end), '%Y-%m-%d %H:%M:%S') + datetime.timedelta(
  344.                 minutes=15*i)
  346.             predict_result_list = []
  348.             for j in range(1, 9):  # 每一个时间点进行8次迭代训练和预测(之前16次,节约时间)
  349.                 predict_power_time_t = datetime.datetime.strptime(str(new_train_end), '%Y-%m-%d %H:%M:%S') + \
  350.                                        datetime.timedelta(minutes=15 * (j - 1))  # 取预测功率的时间,用于训练
  351.                 predict_power_time_p = datetime.datetime.strptime(str(new_train_end), '%Y-%m-%d %H:%M:%S') + \
  352.                                        datetime.timedelta(minutes=15 * j)  # 取预测功率的时间,用于预测
  353.                 if j == 1:  # 第1个点进行训练和预测
  355.                     train_data = get_data(new_train_start, new_train_end)
  356.                     expanded_real_data = get_expand_data(train_data['power'], dimension)  # 首先取实测功率数据并扩展为dimension-1维
  357.                     predict_data = train_data['predict_power'].iloc[dimension - 1:].reset_index(
  358.                         drop=True)  # 预测功率的时间和输出的时间一致
  360.                     expanded_data = pd.merge(expanded_real_data, predict_data.to_frame(), left_index=True,
  361.                                              right_index=True)
  362.                     expanded_data.insert(dimension - 1, 'p_power', expanded_data['predict_power'],
  363.                                          allow_duplicates=False)
  364.                     expanded_data = expanded_data.drop(['predict_power'], axis=1)
  366.                     train_input = expanded_data.iloc[:, 0: -1]  # 前几列为训练的输入,最后一列为训练的输出
  367.                     train_output = expanded_data.iloc[:, -1]
  369.                     estimator = ProjectionPursuitRegressor(ridge_number)
  370.                     estimator.fit(train_input.values, train_output.values)  # 建立模型
  372.                     init_predict_input = expanded_data.iloc[-1, np.r_[1:dimension - 1, dimension]].values.tolist()
  373.                     predict_power = get_predict_data(predict_power_time_p)  # 预测数据的预测功率部分
  374.                     init_predict_input.append(predict_power)
  376.                     predict_input_array = np.array(init_predict_input).reshape(1, -1)  # 转换为数组并reshape
  377.                     predict_result = estimator.predict(predict_input_array)
  378.                     predict_result_list.append(predict_result[0])
  380.                 else:  # 第2~16个点进行训练和预测
  381.                     if j == 2:  # 第一次取j=1时expanded_data最后一行
  382.                         last_col = {}
  383.                         for k in range(dimension):  # 先构建dict
  384.                             if 0 <= k < dimension-2:
  385.                                 last_col['power_%s' % k] = expanded_data.iloc[-1, k + 1]
  386.                             elif k == dimension-2:
  387.                                 last_col['power_%s' % (dimension - 2)] = expanded_data.iloc[-1, -1]
  388.                             else:
  389.                                 last_col['output'] = predict_result
  390.                         predict_power = get_predict_data(predict_power_time_t)
  391.                         last_col['p_power'] = predict_power
  393.                         columns = ['power_%s' % i for i in range(dimension - 1)]
  394.                         columns.append('p_power')  # 最后一行的columns
  395.                         columns.append('output')
  396.                         last_df = pd.DataFrame(last_col, columns=columns, index=[0])
  397.                         new_expanded_data = pd.concat([expanded_data, last_df])  # 新构建的和前面的expanded_data连接
  398.                         new_expanded_data = new_expanded_data.iloc[1:, :].reset_index(drop=True)  # 删掉第一行并更新索引
  400.                     else:
  401.                         last_col = {}
  402.                         for k in range(dimension):  # 先构建dict
  403.                             if 0 <= k < dimension - 2:
  404.                                 last_col['power_%s' % k] = new_expanded_data.iloc[-1, k + 1]
  405.                             elif k == dimension - 2:
  406.                                 last_col['power_%s' % (dimension - 2)] = new_expanded_data.iloc[-1, -1]
  407.                             else:
  408.                                 last_col['output'] = predict_result
  409.                         predict_power = get_predict_data(predict_power_time_t)
  410.                         last_col['p_power'] = predict_power
  412.                         columns = ['power_%s' % i for i in range(dimension - 1)]
  413.                         columns.append('p_power')  # 最后一行的columns
  414.                         columns.append('output')
  415.                         last_df = pd.DataFrame(last_col, columns=columns, index=[0])
  417.                         new_expanded_data = pd.concat([new_expanded_data, last_df])  # 新构建的和前面的expanded_data连接
  418.                         new_expanded_data = new_expanded_data.iloc[1:, :].reset_index(drop=True)  # 删掉第一行并更新索引
  420.                     train_input = new_expanded_data.iloc[:, 0: -1]  # 前几列为训练的输入,最后一列为训练的输出
  421.                     train_output = new_expanded_data.iloc[:, -1]
  423.                     estimator = ProjectionPursuitRegressor(ridge_number)
  424.                     estimator.fit(train_input.values, train_output.values)  # 建立模型
  426.                     predict_input = new_expanded_data.iloc[-1, np.r_[1:dimension - 1, dimension]].values.tolist()
  427.                     predict_power = get_predict_data(predict_power_time_p)  # 预测数据的预测功率部分
  428.                     predict_input.append(predict_power)
  429.                     predict_input_array = np.array(predict_input).reshape(1, -1)
  431.                     predict_result = estimator.predict(predict_input_array)
  432.                     predict_result_list.append(predict_result[0])
  434.             for p in range(len(predict_result_list)):  # 最后对小于0和大于cap的情况进行处理
  435.                 if predict_result_list[p] > cap:
  436.                     predict_result_list[p] = cap
  437.                 elif predict_result_list[p] < 0:
  438.                     predict_result_list[p] = 0
  439.             predict_result_dict[new_train_end + datetime.timedelta(hours=2)] = predict_result_list
  440.         print predict_result_dict
  441.         key_list = []
  442.         predict_list_1, predict_list_2, predict_list_3, predict_list_4, predict_list_5, predict_list_6, predict_list_7, \
  443.         predict_list_8 = [], [], [], [], [], [], [], []
  444.         for k, v in predict_result_dict.items():
  445.             key_list.append(str(k))
  446.             predict_list_1.append(v[0])
  447.             predict_list_2.append(v[1])
  448.             predict_list_3.append(v[2])
  449.             predict_list_4.append(v[3])
  450.             predict_list_5.append(v[4])
  451.             predict_list_6.append(v[5])
  452.             predict_list_7.append(v[6])
  453.             predict_list_8.append(v[7])
  454.         df = pd.DataFrame({'time': key_list, 'predict_1': predict_list_1, 'predict_2': predict_list_2,
  455.                            'predict_3': predict_list_3, 'predict_4': predict_list_4,
  456.                            'predict_5': predict_list_5, 'predict_6': predict_list_6,
  457.                            'predict_7': predict_list_7, 'predict_8': predict_list_8})
  459.         df = df[['time', 'predict_1', 'predict_2', 'predict_3', 'predict_4', 'predict_5',
  460.                  'predict_6', 'predict_7', 'predict_8']].sort_values('time').reset_index(drop=True)
  461.         df.to_csv(r"C:\Users\35760\PycharmProjects\LocalProject\PPR\datas\320907\predict_df_%s.csv" %
  462.                   datetime.datetime.now().strftime('%m_%d_%H_%M'))
  463.         # df.to_csv(r"predict_df_%s.csv" % wfid)  # 服务器部署
  464.         df['time'] = pd.to_datetime(df['time'], format="%Y-%m-%d %H:%M:%S")
  465.         df['day'] = df['time'].apply(lambda x: '%s-%02d-%02d' % (x.year, x.month, x.day))
  466.         df = df.rename(columns={'predict_8': 'predict_power'})
  467.         harmonic_acc_list = []
  468.         day_list = []
  469.         rmse_list = []
  470.         data = pd.read_sql("select rectime, power from `%s` where rectime>='2019-05-01' and rectime<='2019-05-31'" %
  471.                            ('factory_' + wfid), mysql_conn('exportdb'))
  472.         data = data.rename(columns={'rectime': 'time'})
  473.         merge_data = pd.merge(data, df, on='time', how='inner')
  474.         print merge_data
  475.         for day, df_ in merge_data.groupby('day'):  # 计算调和平均准确率
  476.             harmonic_acc = harmonic_mean(df_['power'], df_['predict_power'], cap)
  477.             accuracy = rmse(df_['power'], df_['predict_power'], cap)
  478.             harmonic_acc_list.append(harmonic_acc)
  479.             rmse_list.append(accuracy)
  480.             day_list.append(day)
  481.         df_harv = pd.DataFrame({'day': day_list, 'harv': harmonic_acc_list[0], 'rmse': rmse_list[0]})
  482.         df_harv.to_csv(r"C:\Users\35760\PycharmProjects\LocalProject\PPR\datas\320907\ppr_result_%s.csv" %
  483.                        datetime.datetime.now().strftime('%m_%d_%H_%M'))
  484.         # df_harv.to_csv(r"ppr_result_%s.csv" % wfid)
  485.         end_time = time.time()
  486.         print 'time:%s' % (end_time - start_time)  # 计算一下程序执行的时间,如果过长,要考虑优化代码

邓老师的ppr算法改写

  1. # -*- coding: utf-8 -*-
  2. """
  3. @author: Haojie Shu
  4. @time: 
  5. @description:ppr算法加入风速特征
  6. """
  7. import numpy as np
  8. import MySQLdb
  9. import pandas as pd
  10. import datetime
  11. import math
  12. import time
  13. from sklearn.metrics import make_scorer
  14. from sklearn.metrics import r2_score
  15. from sklearn import svm
  16. from sklearn.model_selection import GridSearchCV
  19. def mysql_conn(db='exportdb'):  # 数据库连接信息
  21.     user = 'dbadmin'
  22.     pwd = 'Goldwind@123'
  23.     ip = '10.12.17.97'
  24.     port = 3306
  25.     conn_modeling = MySQLdb.connect(host=ip, db=db, user=user, passwd=pwd, port=port, charset="utf8")
  27.     return conn_modeling
  30. def get_data(start_, end_):  # 获取实测数据
  32.     export_data = pd.read_sql("select rectime, power, predict_power from `%s` where rectime>='%s' and rectime<='%s'" %
  33.                               ('factory_' + wfid, start_, end_), mysql_conn())
  34.     return export_data
  37. def get_expand_data(data, dimension_):  # 将一维数据变为多维,最后一维为输出,前面为输入
  39.     refactor_sam = []
  40.     for i in range(dimension_ - 1, len(data)):
  41.         block = dimension_ * [0]
  42.         for j in range(dimension_):
  43.             block[j] = data[i]
  44.             i -= 1
  45.         block = list(reversed(block))
  46.         refactor_sam.append(block)
  47.     col = []
  48.     for i in range(dimension_):
  49.         if i != dimension_ - 1:
  50.             col.append(data.name + '_%s' % i)
  51.         else:
  52.             col.append('output')
  53.     expand_df = pd.DataFrame(refactor_sam, columns=col)
  54.     return expand_df
  57. def update_input(predict_result_, array_input_):  # 对输入值进行更新
  59.     array_input_ = np.delete(array_input_, 0)  # 删除掉numpy.ndarray的第一个值
  60.     array_input_ = np.append(array_input_, predict_result_)  # 插入预测值
  61.     return array_input_
  64. def rmse(real, pred, cap_=49500.0):  # 均方根误差准确率
  65.     s = 1 - np.sqrt(np.mean((real - pred) ** 2)) / cap_
  66.     return s * 100
  69. def harmonic_mean(sc, pre, cap_):  # 计算超短期准确率
  70.     index = (sc < 0.03 * cap_) & (pre < 0.03 * cap_)
  71.     sc = sc[~index]
  72.     pre = pre[~index]
  73.     if len(sc) > 0:
  74.         a = np.abs(sc / (sc + pre) - 0.5)
  75.         b = abs(sc - pre) / sum(abs(sc - pre))
  76.         th = 1 - 2 * sum(a * b)
  77.     else:
  78.         th = 1
  79.     return th
  82. def get_predict_data(ptime):
  83.     predict_power = pd.read_sql("select predict_power from `%s` where rectime='%s'" %('factory_' + wfid, ptime),
  84.                                 mysql_conn()).iloc[0, 0]
  85.     return predict_power
  88. def smooth_data(df):  # 数据平滑
  89.     x = df['power'].values
  90.     box = np.ones(10) / 10
  91.     y = np.convolve(x, box, mode='same')  # 两个数组长度可以不一样,参数为same时取更长的一个
  92.     y[:5], y[-5:] = x[:5], x[-5:]
  93.     df['new_power'] = y
  94.     return df
  97. def compensate_data(df, start, end):  # 补缺值
  98.     time_list = pd.date_range(start, end, freq='15min').tolist()
  99.     time_df = pd.DataFrame({'rectime': time_list})
  100.     merge_df = pd.merge(time_df, df, on='rectime', how='outer')
  101.     merge_df = merge_df.reset_index(drop=True)
  102.     merge_df = merge_df.fillna(merge_df.mean())  # NAN值用平均来补
  103.     return merge_df
  106. class PPR:
  108.     def __init__(self, r):
  109.         """
  110.         r:岭函数个数, int
  111.         """
  112.         self.r = r
  113.         if not isinstance(self.r, int) or self.r < 1:
  114.             raise ValueError('r>1')
  115.         self.smooth_ratio = 0.1  # 拟合岭函数的光滑系数
  117.     def fit(self, x, y):  # 训练
  118.         """
  119.         1. x:输入数组, numpy.ndarray
  120.         2. y:输出数组, numpy.ndarray
  121.         """
  122.         ave = np.average(y)
  123.         std = np.std(y)
  124.         y_list = []
  125.         for i in range(y.shape[0]):  # 进行标准化处理
  126.             new_y = (y[i] - ave) / std
  127.             y_list.append(new_y)
  128.         y = np.array(y_list)
  129.         # 投影, 岭函数模型, 权重, 损失, 拟合值, x带索引的list
  130.         projection_list, ridge_weight_list, loss_list, smo_list, enu_x_list, x_list_list = [], [], [], [], [], []
  132.         projection = self.find_projection(x, y)  # 第一个投影方向
  133.         projection_list.append(projection)
  135.         x_list, y_list, smo, loss, enu_x = self.create_one_f(x, y)  # smo,拟合值列表(list), enu_x带顺序的列表
  137.         loss_list.append(loss)
  138.         smo_list.append(smo)  # smo_list是list的list
  139.         enu_x_list.append(enu_x)  # enu_x_list也是list的list
  140.         x_list_list.append(x_list)
  142.         ridge_weight = self.get_ridge_weight(x, y)  # 第一个权重
  143.         ridge_weight_list.append(ridge_weight)
  145.         for i in range(1, self.r):  # 从第二个岭函数开始循环(如果是4个岭函数,不能这么循环)
  146.             new_y_list = []
  147.             for j in range(y.shape[0]):
  148.                 for k in range(len(enu_x_list[i-1])):
  149.                     if enu_x_list[i-1][k][0] == j:  # k用来找索引
  150.                         new_y = y[j] - ridge_weight_list[i-1] * smo_list[i-1][k]  # 偏差 = y- 权重*smo
  151.                 new_y_list.append(new_y)  # 每一行得到一个新的y
  153.             new_y_array = np.array(new_y_list)  # 转为一维数组
  154.             projection = self.find_projection(x, new_y_array)  # 将偏差作为下一轮的y
  155.             projection_list.append(projection)
  157.             x_list, y_list, smo, loss, enu_x = self.create_one_f(x, new_y_array)
  158.             loss_list.append(loss)
  159.             smo_list.append(smo)
  160.             enu_x_list.append(enu_x)
  161.             x_list_list.append(x_list)
  163.             ridge_weight = self.get_ridge_weight(x, new_y_array)
  164.             ridge_weight_list.append(ridge_weight)
  166.             if loss_list[i] > loss_list[i-1]:  # 后一次的残差大于前一次的残差,迭代停止
  167.                 break
  168.         return projection_list, ridge_weight_list, ave, std, smo_list, x_list_list
  170.     def predict(self, x, y, input_data_):  # ppr函数里的predict函数
  172.         projection_list, ridge_weight_list, ave, std, new_smo_list, x_list_list = self.fit(x, y)
  174.         # 判断是哪种方式迭代停止的
  175.         if len(x_list_list) == self.r:  # 岭函数个数=r
  176.             output_list = []
  177.             for j in range(len(x_list_list)):  # 遍历每一个岭函数
  178.                 input_sum = 0
  179.                 for k in range(input_data_.shape[1]):
  180.                     input_sum += input_data_[0][k] * projection_list[j][k]  # 投影*输入
  181.                 # print 'input_sum:%s' % input_sum
  182.                 min_x, max_x = min(x_list_list[j]), max(x_list_list[j])
  183.                 # 都是等比例输出
  184.                 if input_sum > max_x:
  185.                     if x_list_list[j][-1] == x_list_list[j][-2]:  # 最大的两个值刚好相等(防止分母为0)
  186.                         output = new_smo_list[j][-1]
  187.                         new_output = output * ridge_weight_list[j]
  188.                         output_list.append(new_output)
  189.                     else:
  190.                         ratio = (input_sum - x_list_list[j][-2]) / (x_list_list[j][-1] - x_list_list[j][-2])
  191.                         output = ratio * (new_smo_list[j][-1] - new_smo_list[j][-2]) + new_smo_list[j][-2]
  192.                         # print 'max_output:%s' % output
  193.                         new_output = output * ridge_weight_list[j]
  194.                         output_list.append(new_output)
  195.                 elif input_sum < min_x:
  196.                     if x_list_list[j][1] == x_list_list[j][0]:
  197.                         output = new_smo_list[j][0]
  198.                         new_output = output * ridge_weight_list[j]
  199.                         output_list.append(new_output)
  200.                     else:
  201.                         ratio = (x_list_list[j][1]-input_sum)/(x_list_list[j][1]-x_list_list[j][0])
  202.                         output = new_smo_list[j][1]-(new_smo_list[j][1]-new_smo_list[j][0])*ratio
  203.                         # print 'min_output:%s' % output
  204.                         new_output = output * ridge_weight_list[j]
  205.                         output_list.append(new_output)
  206.                 else:  # 最大最小之间
  207.                     for h in range(len(x_list_list[j])):
  208.                         if x_list_list[j][h] == input_sum:  # input_sum和x_list_list某个值刚好相等
  209.                             output_list.append(new_smo_list[j][h] * ridge_weight_list[j])
  210.                             # print 'mid_output:%s' % ridge_weight_list[j]
  211.                             break  # 有可能相同的不止一个,取第一个就行
  212.                         elif x_list_list[j][h] < input_sum < x_list_list[j][h + 1]:  # 如果不相等,必定在两个之间
  213.                             ratio = (input_sum - x_list_list[j][h]) / (x_list_list[j][h+1] - x_list_list[j][h])
  214.                             output = (new_smo_list[j][h+1] - new_smo_list[j][h]) * ratio + new_smo_list[j][h]
  215.                             # print 'mid_output:%s' % output
  216.                             new_output = output * ridge_weight_list[j]
  217.                             output_list.append(new_output)
  218.                         else:
  219.                             pass
  221.         else:
  222.             output_list = []
  223.             for j in range(len(x_list_list)):  # 遍历每一个岭函数
  224.                 input_sum = 0
  225.                 for k in range(input_data_.shape[1]):
  226.                     input_sum += input_data_[0][k] * projection_list[j][k]  # 投影*输入
  227.                 min_x, max_x = min(x_list_list[j]), max(x_list_list[j])
  228.                 # 都是等比例输出
  229.                 if input_sum > max_x:
  230.                     if x_list_list[j][-1] == x_list_list[j][-2]:
  231.                         output = new_smo_list[j][-1]
  232.                         new_output = output * ridge_weight_list[j]
  233.                         output_list.append(new_output)
  234.                     else:
  235.                         ratio = (input_sum - x_list_list[j][-2]) / (x_list_list[j][-1] - x_list_list[j][-2])
  236.                         output = ratio * (new_smo_list[j][-1] - new_smo_list[j][-2]) + new_smo_list[j][-2]
  237.                         new_output = output * ridge_weight_list[j]
  238.                         output_list.append(new_output)
  239.                 elif input_sum < min_x:
  240.                     ratio = (x_list_list[j][0] - input_sum) / (x_list_list[j][1] - input_sum)
  241.                     output = (ratio * new_smo_list[j][1] - new_smo_list[j][0]) / ratio - 1
  242.                     new_output = output * ridge_weight_list[j]
  243.                     output_list.append(new_output)
  244.                 else:  # 最大最小之间
  245.                     for h in range(len(x_list_list[j])):
  246.                         if x_list_list[j][h] == input_sum:  # input_sum和x_list_list某个值刚好相等
  247.                             output_list.append(new_smo_list[j][h] * ridge_weight_list[j])
  248.                         elif x_list_list[j][h] < input_sum < x_list_list[j][h + 1]:  # 如果不相等,必定在两个之间
  249.                             ratio = (input_sum - x_list_list[j][h]) / (x_list_list[j][h+1] - x_list_list[j][h])
  250.                             output = (new_smo_list[j][h+1] - new_smo_list[j][h]) * ratio + new_smo_list[j][h]
  251.                             new_output = output * ridge_weight_list[j]
  252.                             output_list.append(new_output)
  253.                         else:
  254.                             pass
  255.         ult_output = sum(output_list)*std+ave
  256.         return ult_output
  258.     @staticmethod
  259.     def cal_var(x, y):
  261.         avg_x = np.mean(x, axis=0)  # 求数组x每一列的均值
  263.         cov_x_y_list = []
  264.         for p in range(x.shape[1]):
  265.             s = 0.0
  266.             for q in range(x.shape[0]):
  267.                 s = s + y[q] * (x[q, p] - avg_x[p])
  268.             cov_x_y_list.append([p, s / x.shape[0]])  # 计算x/y协方差
  270.         var_x_x_list = []
  271.         for i in range(x.shape[1]):
  272.             for j in range(i, x.shape[1]):
  273.                 s = 0.0
  274.                 for k in range(x.shape[0]):
  275.                     tt = (x[k, j] - avg_x[j]) * (x[k, j] - avg_x[i])
  276.                     s = s + tt
  277.                 var_x_x_list.append([i, j, s / x.shape[0]])
  278.         return cov_x_y_list, var_x_x_list
  280.     def find_projection(self, x, y):  # 由自方差和协方差来寻找投影方向(权重系数)
  282.         cov_list, var_list = self.cal_var(x, y)
  283.         project_dict = {}  # 投影方向的dict
  284.         for i in range(x.shape[1]):
  285.             project_dict[i] = 1
  287.         for i in range(x.shape[1]):
  288.             k1 = i + (i * (i + 1) / 2)  # 取自方差的索引,按照C#代码写的
  289.             sum_var = project_dict[i] * var_list[k1][2]
  290.             for j in range(i + 1, x.shape[1]):
  291.                 k2 = i + (j * (j + 1) / 2)
  292.                 sum_var += project_dict[i] * var_list[k2][2]
  293.             for k in range(i):
  294.                 k3 = i + (k * (k + 1) / 2)
  295.                 sum_var += project_dict[i] * var_list[k3][2]
  296.             project_dict[i] = cov_list[i][1] / sum_var
  297.         '''
  298.         用自方差和协方差求完投影方向后,还要进行一个简单的处理
  299.         '''
  300.         sum_project = 0
  301.         new_projection_dict = {}
  302.         for v in project_dict.values():
  303.             sum_project += v * v
  304.         sum_project = 1.0 / math.sqrt(sum_project)
  305.         for k, v in project_dict.items():
  306.             new_projection = v * sum_project
  307.             new_projection_dict[k] = new_projection
  308.         return new_projection_dict
  310.     def get_project_label(self, x, y):  # 获取投影后的坐标list, x和y是对应起来的
  312.         project_dict = self.find_projection(x, y)
  313.         x_list, y_list = [], []
  314.         for i in range(x.shape[0]):
  315.             x_label = 0
  316.             for j in range(x.shape[1]):
  317.                 x_label += x[i][j] * project_dict[j]  # 投影方向和输入作乘法,得到的值作为横坐标
  318.             x_list.append(x_label)
  319.             y_list.append(y[i])
  320.         return x_list, y_list
  322.     @staticmethod
  323.     def performance_metric(y_true, y_predict):
  324.         score = r2_score(y_true, y_predict)
  325.         return score
  327.     def fit_curve(self, x, y):  # 拟合曲线
  329.         x_list, y_list, smo_list, loss, index_list = self.create_one_f(x, y)
  330.         x_array = np.array(x_list)
  331.         smo_array = np.array(smo_list)
  332.         reshape_x_array = x_array.reshape(x_array.shape[0], -1)
  333.         scorig_fnc = make_scorer(self.performance_metric)
  334.         params = {'C': [10], 'gamma': [0.1]}
  335.         svc = svm.SVR()
  336.         grid = GridSearchCV(svc, param_grid=params, scoring=scorig_fnc, cv=5)  # 网格搜索
  337.         model = grid.fit(reshape_x_array, smo_array)
  338.         return model, loss
  340.     def create_one_f_back(self, x, y):  # 创建一个岭函数
  342.         x_list, y_list = self.get_project_label(x, y)
  343.         enu_x = list(enumerate(x_list))
  344.         df = pd.DataFrame({'x': x_list, 'y': y_list, 'enu_x': enu_x})
  345.         sorted_df = df.sort_values(by='x')  # 按照x进行升序排序
  347.         new_x_list = sorted_df['x'].tolist()
  348.         new_y_list = sorted_df['y'].tolist()
  349.         new_enu_x = sorted_df['enu_x'].tolist()
  351.         smo_list = self.super_filter(new_x_list, new_y_list)  # 岭函数拟合后的纵坐标
  353.         loss_sum = 0
  354.         for j in range(len(smo_list)):
  355.             loss = math.pow((new_y_list[j] - smo_list[j]), 2)
  356.             loss_sum += loss  # 总的误差
  358.         return new_x_list, new_y_list, smo_list, loss_sum, new_enu_x
  360.     def create_one_f(self, x, y):  # 创建一个岭函数
  362.         x_list, y_list = self.get_project_label(x, y)
  363.         enu_x = list(enumerate(x_list))
  364.         df = pd.DataFrame({'x': x_list, 'y': y_list, 'enu_x': enu_x})
  365.         sorted_df = df.sort_values(by='x')  # 按照x进行升序排序
  367.         new_x_list = sorted_df['x'].tolist()
  368.         new_y_list = sorted_df['y'].tolist()
  369.         new_enu_x = sorted_df['enu_x'].tolist()
  371.         smo_list = self.super_filter(new_x_list, new_y_list)  # 岭函数拟合后的纵坐标
  373.         loss_sum = 0
  374.         for j in range(len(smo_list)):
  375.             loss = math.pow((new_y_list[j] - smo_list[j]), 2)
  376.             loss_sum += loss  # 总的误差
  377.         avg_smo_list = sum(smo_list)/len(smo_list)
  378.         v = 0.0
  379.         for k in range(len(smo_list)):
  380.             smo_list[k] = smo_list[k]-avg_smo_list
  381.             v += smo_list[k]*smo_list[k]
  382.         if v > 0:
  383.             v = 1/math.sqrt(v/len(smo_list))
  384.             for s in range(len(smo_list)):
  385.                 smo_list[s] = smo_list[s]*v
  386.         return new_x_list, new_y_list, smo_list, loss_sum, new_enu_x
  388.     def super_filter(self, x_list, y_list):  # 超级滤波器
  389.         """
  390.         :param x_list: 经过排序的list,投影方向和输入x的乘积
  391.         :param y_list: y_list和x_list是对应的
  392.         在c#中
  393.         :return:
  394.         """
  395.         len_list = len(x_list)
  396.         scale = x_list[3 * len_list / 4] - x_list[len_list / 4]
  397.         vsm = scale * scale
  399.         smo_list = self.smooth_window(x_list, y_list, vsm)
  401.         new_y_list = []
  402.         for i in range(len_list):
  403.             new_y_list.append(smo_list[i])
  405.         new_smo_list = self.smooth_window(x_list, new_y_list, vsm)  # 调用两次,增加平滑度
  407.         return new_smo_list
  409.     def smooth_window(self, x_list, y_list, vsm):  # 对窗口内的点作平滑处理
  411.         """
  412.         :param x_list:经过排序的list,投影方向和输入x的乘积
  413.         :param y_list:y_list和x_list是对应的
  414.         :param vsm:滑动距离参数
  415.         需要注意的是,在C#中smooth_window函数虽然有smo_list,acvr_list这两个参数,但是它们并不起作用
  416.         :return:
  417.         """
  419.         win_x, win_y, fbw, xvar, cvar = 0.0, 0.0, 0.0, 0.0, 0.0  # 岭函数滑动窗口的一些初始值
  420.         n_points = int(0.5 * self.smooth_ratio * len(x_list) + 0.5)  # 每次窗口滑动要包含的点数
  422.         if n_points < 2:  # 每个窗口内的点数不能少于2
  423.             n_points = 2
  425.         n_iter = int(2 * n_points + 1)  # 循环次数
  426.         if n_iter > len(x_list):  # 循环次数不能大于横坐标的总个数
  427.             n_iter = len(x_list)
  429.         for i in range(n_iter):  # 第一次循环滤波处理, 计算波型边界参数
  431.             list_index = i  # 取x_list, y_list的一个索引
  432.             xti = x_list[list_index]
  434.             if list_index < 1:
  435.                 list_index = len(x_list) - 1 + list_index
  436.                 xti = x_list[list_index] - 1.0
  438.             fbo = fbw
  439.             fbw = fbw + 1
  441.             if fbw > 0:
  442.                 win_x = (fbo * win_x + xti) / fbw
  443.                 win_y = (fbo * win_y + y_list[list_index]) / fbw
  444.             tmp = 0.0
  446.             if fbo > 0:
  447.                 tmp = fbw * (xti - win_x) / fbo
  448.             xvar = xvar + tmp * (xti - win_x)
  449.             cvar = cvar + tmp * (y_list[list_index] - win_y)
  450.         smo_list, acvr_list = [], []
  451.         for p in range(len(x_list)):  # 第二次循环滤波处理,计算SMO[0..N],ACVR[0..N]
  452.             out = p - n_points - 1
  453.             inp = p + n_points
  454.             if out >= 1 and inp <= len(x_list):
  455.                 if out < 1:
  456.                     out = len(x_list) + out
  457.                     xto = x_list[out] - 1.0
  458.                     xti = x_list[inp]
  459.                 else:
  460.                     if inp > len(x_list) - 1:
  461.                         inp = inp - len(x_list)
  462.                         xti = x_list[inp] + 1.0
  463.                         xto = x_list[out]
  464.                     else:
  465.                         xto = x_list[out]
  466.                         xti = x_list[inp]
  467.                 fbo = fbw
  468.                 fbw = fbw - 1
  469.                 tmp = 0.0
  470.                 if fbw > 0.0:
  471.                     tmp = fbo * (xto - win_x) / fbw
  472.                 xvar = xvar - tmp * (xto - win_x)
  473.                 cvar = cvar - tmp * (y_list[out] - win_y)
  474.                 if fbw > 0.0:
  475.                     win_x = (fbo * win_x - xto) / fbw
  476.                     win_y = (fbo * win_y - y_list[out]) / fbw
  477.                 fbo = fbw
  478.                 fbw = fbw + 1
  479.                 if fbw > 0.0:
  480.                     win_x = (fbo * win_x + xti) / fbw
  481.                     win_y = (fbo * win_y + y_list[inp]) / fbw
  482.                 tmp = 0.0
  483.                 if fbo > 0.0:
  484.                     tmp = fbw * (xti - win_x) / fbo
  485.                 xvar = xvar + tmp * (xti - win_x)
  486.                 cvar = cvar + tmp * (y_list[inp] - win_y)
  488.             a = 0.0
  489.             if xvar > vsm:
  490.                 a = cvar / xvar
  491.             smo = a * (x_list[p] - win_x) + win_y
  492.             smo_list.append(smo)
  494.             h = 0.0
  495.             if fbw > 0.0:
  496.                 h = 1.0 / fbw
  497.             if xvar > vsm:
  498.                 h = h + (x_list[p] - win_x) * (x_list[p] - win_x) / xvar
  500.             acvr = 0
  501.             a = 1.0 - h
  502.             if a > 0.0:
  503.                 acvr = abs(y_list[p] - smo_list[p]) / a
  504.             else:
  505.                 if a > 1:
  506.                     acvr = acvr_list[p - 1]
  507.             acvr_list.append(acvr)
  509.         jj = 0
  510.         while jj <= len(x_list) - 1:
  511.             jo = jj
  512.             sy = smo_list[jj]
  513.             j = jj
  514.             fbw = 1
  515.             while True:
  516.                 j = j + 1
  517.                 if j >= len(x_list) - 1:
  518.                     break
  519.                 if x_list[j - 1] <= x_list[j]:
  520.                     break
  521.                 sy = sy + smo_list[j - 1]
  522.                 fbw = fbw + 1
  523.                 if j > jo:
  524.                     a = 0.0
  525.                     if fbw > 0.0:
  526.                         a = sy / fbw
  527.                     for q in range(j):
  528.                         smo_list[q] = a
  530.             jj = jj + 1
  531.         return smo_list
  533.     def get_ridge_weight(self, x, y):  # 岭函数的权重,也即β参数
  535.         x_list, y_list, smo_list, loss, index_list = self.create_one_f(x, y)
  536.         var_sum, cov_sum = 0.0, 0.0
  537.         for i in range(len(smo_list)):
  538.             cov_sum += y_list[i] * smo_list[i]  # 协方差和
  539.             var_sum += smo_list[i] * smo_list[i]  # 自方差和
  540.         ridge_weight = cov_sum / var_sum
  541.         return ridge_weight
  544. if __name__ == '__main__':
  546.     # 需要配置如下6个变量
  547.     wfid = '320907'
  548.     ridge_number = 4  # 岭函数个数
  549.     dimension = 5  # 需要的数据总维度,如果是二阶输入的话,就是2+8=10,后面8个为输出
  550.     train_data_sets = 96*30  # 训练的时候用到的输入和输出对数
  551.     predict_start = '2019-05-01 00:00:00'  # 要预测的开始时间
  552.     predict_points = 96*1  # 要预测的点数(一天96点)
  554.     start_time = time.time()
  555.     cap = pd.read_sql("select powercap from `wind_farm_info` where wfid='%s'"
  556.                       % wfid, mysql_conn('WindFarmdb')).iloc[0, 0]
  558.     train_end = datetime.datetime.strptime(str(predict_start), '%Y-%m-%d %H:%M:%S') + \
  559.                 datetime.timedelta(hours=-2)  # 超短期提前两小时预测当前的
  560.     # train_start由train_end和train_data_sets以及dimension决定,保证有192对输入输出
  561.     train_start = datetime.datetime.strptime(str(train_end), '%Y-%m-%d %H:%M:%S') + \
  562.                   datetime.timedelta(minutes=-15*(train_data_sets+dimension-2))
  564.     predict_result_dict = {}
  565.     for i in range(predict_points):  # 每次训练和预测1个点
  567.         new_train_start = datetime.datetime.strptime(str(train_start), '%Y-%m-%d %H:%M:%S') + \
  568.                           datetime.timedelta(minutes=15*i)
  569.         new_train_end = datetime.datetime.strptime(str(train_end), '%Y-%m-%d %H:%M:%S') + datetime.timedelta(
  570.             minutes=15*i)
  572.         predict_result_list = []
  574.         for j in range(1, 9):  # 每一个时间点进行8次迭代训练和预测(之前16次,节约时间)
  575.             predict_power_time_t = datetime.datetime.strptime(str(new_train_end), '%Y-%m-%d %H:%M:%S') + \
  576.                                  datetime.timedelta(minutes=15*(j-1))  # 取预测功率的时间,用于训练
  577.             predict_power_time_p = datetime.datetime.strptime(str(new_train_end), '%Y-%m-%d %H:%M:%S') + \
  578.                                  datetime.timedelta(minutes=15*j)  # 取预测功率的时间,用于预测
  579.             if j == 1:  # 第1个点进行训练和预测
  581.                 train_data = get_data(new_train_start, new_train_end)
  582.                 expanded_real_data = get_expand_data(train_data['power'], dimension)  # 首先取实测功率数据并扩展为dimension-1维
  583.                 predict_data = train_data['predict_power'].iloc[dimension-1:].reset_index(drop=True)  # 预测功率的时间和输出的时间一致
  585.                 expanded_data = pd.merge(expanded_real_data, predict_data.to_frame(), left_index=True, right_index=True)
  586.                 expanded_data.insert(dimension-1, 'p_power', expanded_data['predict_power'], allow_duplicates=False)
  587.                 expanded_data = expanded_data.drop(['predict_power'], axis=1)
  589.                 train_input = expanded_data.iloc[:, 0: -1]  # 前几列为训练的输入,最后一列为训练的输出
  590.                 train_output = expanded_data.iloc[:, -1]
  592.                 estimator = PPR(ridge_number)
  593.                 estimator.fit(train_input.values, train_output.values)  # 建立模型
  595.                 init_predict_input = expanded_data.iloc[-1, np.r_[1:dimension-1, dimension]].values.tolist()
  596.                 predict_power = get_predict_data(predict_power_time_p)  # 预测数据的预测功率部分
  597.                 init_predict_input.append(predict_power)
  599.                 predict_input_array = np.array(init_predict_input).reshape(1, -1)  # 转换为数组并reshape
  600.                 predict_result = estimator.predict(train_input.values, train_output.values, predict_input_array)
  601.                 predict_result_list.append(predict_result)
  603.             else:  # 第2~16个点进行训练和预测
  604.                 if j == 2:  # 第一次取j=1时expanded_data最后一行
  605.                     last_col = {}
  606.                     for k in range(dimension):  # 先构建dict
  607.                         if 0 <= k < dimension-2:
  608.                             last_col['power_%s' % k] = expanded_data.iloc[-1, k+1]
  609.                         elif k == dimension-2:
  610.                             last_col['power_%s' % (dimension-2)] = expanded_data.iloc[-1, -1]
  611.                         else:
  612.                             last_col['output'] = predict_result
  613.                     predict_power = get_predict_data(predict_power_time_t)
  614.                     last_col['p_power'] = predict_power
  616.                     columns = ['power_%s' % i for i in range(dimension-1)]
  617.                     columns.append('p_power')  # 最后一行的columns
  618.                     columns.append('output')
  619.                     last_df = pd.DataFrame(last_col, columns=columns, index=[0])
  621.                     new_expanded_data = pd.concat([expanded_data, last_df])  # 新构建的和前面的expanded_data连接
  622.                     new_expanded_data = new_expanded_data.iloc[1:, :].reset_index(drop=True)  # 删掉第一行并更新索引
  624.                 else:
  625.                     last_col = {}
  626.                     for k in range(dimension):  # 先构建dict
  627.                         if 0 <= k < dimension-2:
  628.                             last_col['power_%s' % k] = new_expanded_data.iloc[-1, k+1]
  629.                         elif k == dimension-2:
  630.                             last_col['power_%s' % (dimension-2)] = new_expanded_data.iloc[-1, -1]
  631.                         else:
  632.                             last_col['output'] = predict_result
  634.                     predict_power = get_predict_data(predict_power_time_t)
  635.                     last_col['p_power'] = predict_power
  637.                     columns = ['power_%s' % i for i in range(dimension-1)]
  638.                     columns.append('p_power')  # 最后一行的columns
  639.                     columns.append('output')
  640.                     last_df = pd.DataFrame(last_col, columns=columns, index=[0])
  642.                     new_expanded_data = pd.concat([new_expanded_data, last_df])  # 新构建的和前面的expanded_data连接
  643.                     new_expanded_data = new_expanded_data.iloc[1:, :].reset_index(drop=True)  # 删掉第一行并更新索引
  645.                 train_input = new_expanded_data.iloc[:, 0: -1]  # 前几列为训练的输入,最后一列为训练的输出
  646.                 train_output = new_expanded_data.iloc[:, -1]
  648.                 estimator = PPR(ridge_number)
  649.                 estimator.fit(train_input.values, train_output.values)  # 建立模型
  651.                 predict_input = new_expanded_data.iloc[-1, np.r_[1:dimension-1, dimension]].values.tolist()
  652.                 predict_power = get_predict_data(predict_power_time_p)  # 预测数据的预测功率部分
  653.                 predict_input.append(predict_power)
  654.                 predict_input_array = np.array(predict_input).reshape(1, -1)
  656.                 predict_result = estimator.predict(train_input.values, train_output.values, predict_input_array)
  657.                 predict_result_list.append(predict_result)
  659.         for p in range(len(predict_result_list)):  # 最后对小于0和大于cap的情况进行处理
  660.             if predict_result_list[p] > cap:
  661.                 predict_result_list[p] = cap
  662.             elif predict_result_list[p] < 0:
  663.                 predict_result_list[p] = 0
  664.         predict_result_dict[new_train_end + datetime.timedelta(hours=2)] = predict_result_list
  665.     """
  666.     计算调和平均准确率
  667.     """
  668.     key_list = []
  669.     predict_list_1, predict_list_2, predict_list_3, predict_list_4, predict_list_5, predict_list_6, predict_list_7, \
  670.     predict_list_8 = [], [], [], [], [], [], [], []
  671.     for k, v in predict_result_dict.items():
  672.         key_list.append(str(k))
  673.         predict_list_1.append(v[0])
  674.         predict_list_2.append(v[1])
  675.         predict_list_3.append(v[2])
  676.         predict_list_4.append(v[3])
  677.         predict_list_5.append(v[4])
  678.         predict_list_6.append(v[5])
  679.         predict_list_7.append(v[6])
  680.         predict_list_8.append(v[7])
  681.     df = pd.DataFrame({'time': key_list, 'predict_1': predict_list_1, 'predict_2': predict_list_2,
  682.                        'predict_3': predict_list_3, 'predict_4': predict_list_4,
  683.                        'predict_5': predict_list_5, 'predict_6': predict_list_6,
  684.                        'predict_7': predict_list_7, 'predict_8': predict_list_8})
  686.     df = df[['time', 'predict_1', 'predict_2', 'predict_3', 'predict_4', 'predict_5',
  687.              'predict_6', 'predict_7', 'predict_8']].sort_values('time').reset_index(drop=True)  # 调整df列的顺序并按照time进行排序
  688.     # df.to_csv(r"C:\Users\35760\PycharmProjects\LocalProject\PPR\datas\320907\predict_df_%s.csv" %
  689.     #           datetime.datetime.now().strftime('%m_%d_%H_%M'))
  690.     df.to_csv(r"predict_df_%s.csv" % wfid)  # 服务器部署
  691.     df['time'] = pd.to_datetime(df['time'], format="%Y-%m-%d %H:%M:%S")
  692.     df['day'] = df['time'].apply(lambda x: '%s-%02d-%02d' % (x.year, x.month, x.day))
  693.     df = df.rename(columns={'predict_8': 'predict_power'})
  694.     harmonic_acc_list = []
  695.     day_list = []
  696.     rmse_list = []
  697.     data = pd.read_sql("select rectime, power from `%s` where rectime>='2019-05-01' and rectime<='2019-05-31'" %
  698.                        ('factory_' + wfid), mysql_conn('exportdb'))
  699.     data = data.rename(columns={'rectime': 'time'})
  700.     merge_data = pd.merge(data, df, on='time', how='inner')
  702.     for day, df_ in merge_data.groupby('day'):  # 计算调和平均准确率
  703.         harmonic_acc = harmonic_mean(df_['power'], df_['predict_power'], cap)
  704.         accuracy = rmse(df_['power'], df_['predict_power'], cap)
  705.         harmonic_acc_list.append(harmonic_acc)
  706.         rmse_list.append(accuracy)
  707.         day_list.append(day)
  708.     df_harv = pd.DataFrame({'day': day_list, 'harv': harmonic_acc_list, 'rmse': rmse_list})
  709.     # df_harv.to_csv(r"C:\Users\35760\PycharmProjects\LocalProject\PPR\datas\320907\ppr_result_%s.csv" %
  710.     #                datetime.datetime.now().strftime('%m_%d_%H_%M'))
  711.     df_harv.to_csv(r"ppr_result_%s.csv" % wfid)
  712.     end_time = time.time()
  713.     print 'time:%s' % (end_time - start_time)  # 计算一下程序执行的时间,如果过长,要考虑优化代码