电熔镁炉是一种以菱镁矿为原料, 由电流控制器控制熔炼电流来生产电熔镁砂的重要设备.产品电熔镁砂是一种应用于冶金、化工、航天等领域的重要高级耐火材料.电熔镁群炉需量指当前时刻和当前时刻之前一定时间内群炉功率的平均值, 用于度量高耗能电熔镁群炉的用电量.在生产过程中需量不得超过规定的最大需量即需量峰值, 以限制电熔镁群炉的用电量.需量监控系统对群炉需量进行实时监控, 当需量超过需量峰值的限幅值时会切断某台炉供电, 以保证群炉需量不超过需量峰值; 当需量低于限幅值时再恢复该台炉供电, 使该炉继续生产.为了保证电熔镁砂质量, 需要电流控制器将电流控制在工艺设定值附近.在电流控制器的调节作用下, 熔炼过程中原料杂质成分含量增大和颗粒长度变大可能会导致需量尖峰, 即需量先增大超过限幅值而后下降低于限幅值.而需量尖峰会造成切断电熔镁炉的供电.然而, 切断供电会破坏炉内温度场吸热和放热平衡, 降低电熔镁砂质量, 因此对需量进行准确的预报对于避免尖峰时刻的错误拉闸显得十分重要.

近年来, 针对电力系统的功率预报问题相关学者开展了一系列研究, 多采用时间序列方法^[1]、支持向量回归机^[2]、神经网络^[3]和混合方法^[4]等.例如, 文献[1]为了制定某地区的发电量的日计划, 基于该地区过去14天用电总功率的数据(采样周期为15分钟), 采用相似形时间序列法先得到参考基准曲线, 然后预报未来一天用电总功率的曲线.文献[2]为了提高某地电力系统的可靠性, 使用该地过去7天用电总功率的数据(采样周期为1小时), 采用支持向量回归机的方法预报未来一天用电总功率曲线.文献[3]为了合理调度某地电力系统, 用该地过去7天风力发电总功率的数据(采样周期为0.5小时), 采用前馈神经网络和上下限估计的方法预报未来一天风力发电总功率的上下限.文献[4]为了降低某地电力市场的经济损失风险, 基于昨天和一周前的一天的用电功率数据和温度数据(采样周期为1小时), 采用遗传算法-径向基函数神经网络(Genetic algorithm --- radial basis function neural network, GA-RBFNN)方法预报未来一天的用电功率.文献[5]为了准确预报群炉需量变化趋势减少不必要的切断供电, 使用某电熔镁砂厂过去时刻的群炉功率数据(采样周期为7秒), 采用RBFNN方法预报下一时刻的群炉需量.

在电熔镁炉生产过程中, 控制系统通过调节电极位置, 改变电弧弧长, 进而控制熔炼电流稳定在电流设定值附近, 从而实现在满足产量约束的条件下尽量降低产品单吨能耗的控制目标^[6].电熔镁炉的熔炼电流通常在15 000 A左右, 电能消耗巨大.熔炼过程中, 当原料杂质成分含量增大和颗粒长度变大时, 工作电阻减小, 电流变大, 需量上升.此时电流控制器会调节弧长使工作电阻变大, 降低电流, 需量又随之下降, 这样就会出现需量先升高后下降的尖峰现象.因此根据功率变化特性建立需量的动态模型才能更准确地预报需量.而文献[1-4]预报的对象采样周期时间尺度较大, 与电熔镁群炉的功率变化特性不同, 而且只单纯依据过去功率的数据进行预报, 没有研究对象的动态特性分析.因此上述文献难以直接适用预报电熔镁群炉需量.文献[5]虽然对电熔镁群炉需量进行了预报, 但只是将历史功率数据和功率变化率整体视为非线性函数关系进行处理, 预报精度有待提高.

本文首先建立和分析需量动态模型, 提出由功率变化率的线性项、基于PACF输入变量决策的RBFNN未知非线性项估计和需量计算模型组成的需量预报方法, 通过某电熔镁砂厂实际数据的仿真实验和工业实验表明所提方法能够准确预报需量的变化趋势.

1.   (k + 1)时刻需量的动态模型

为了准确预报$(k+1)$时刻需量, 本文首先根据需量定义建立$(k+1)$时刻的需量模型, 然后建立功率变化率与电流控制系统输出电流之间的由线性项与未知非线性项组成的动态模型.

1.1   需量的定义

如图 1所示, 在熔炼电压$U$作用下, 电流控制器调节升降电机转速$u_i$, 使电流实际值$y_i$跟踪电流设定值$y^{*}$而产生用电功率.炉内原料吸收电弧释放的热量熔化形成逐渐上涨的MgO熔池, 熔炼结束后熔池经过冷却、结晶、破碎等工序形成产品电熔镁砂.熔炼过程中功率变送器测量电力变压器得到群炉功率数据$p(k)$, $p(k - 1)$, $\cdots$, $p(k - n + 1)$, 由需量计算装置基于定义1可得到当前时刻的群炉需量$\bar P(k)$.

图 1  电熔镁群炉需量监控原理图

Fig. 1  Schematic diagram of demand monitoring process for FMFG

下载: 全尺寸图片 幻灯片

定义1. 群炉需量$\bar P(k)$为$k$时刻和$(k-1)$, $\cdots$, $(k-n+1)$时刻群炉功率$p(k)$的平均值

其中, 群炉功率$p(k)$为

其中, $m$为电熔镁群炉台数, 各炉熔炼电压为$U$ (常量), ${y_i}(k)$为第$i$台电熔镁炉的熔炼电流, $\rm cos\varphi$为功率因数.

1.2   (k + 1)时刻需量模型

由定义1可递推得到$(k+1)$时刻的需量$\bar{P}(k+1)$.

其中, $\bar P(k)$, $p(k)$, $p(k + 1 - n)$都已知, 可见预报$\bar P(k + 1)$的关键在于估计未知的功率变化率$\Delta p(k)$.

1.3   功率变化率与电流之间的动态模型

式(3)中功率变化率$\Delta p(k)$的定义为

其中, $y_i(k+1)$和$y_i(k)$是电流闭环控制系统的输出.

首先建立以电极升降电机转速$u_i$为输入, 以电流$y_i$为输出的被控对象模型.电流$y_i$与工作电阻$R_i$之间关系有

其中, $R_i$^[7]为电弧电阻$R_{i, {\rm arc}}$和熔池电阻$R_{i, {\rm pool}}$之和, 即

其中, 电弧电阻$R_{i, {\rm arc}}$^[8]为

其中, 参数${G_0} = {( {{g_0}\pi r_{{{\rm arc}}}^2\exp ( {{{ - {T_0}} / {{T_1}}}} )} )^{-1}}$, $g_{0}$为电弧电导率常数, ${r_{{{\rm arc}}}}$为电弧弧柱半径, $T_{0}$为气体电离温度常数, $T_{1}$为电弧间隙温度, $L_{i, {\rm arc}}$为电弧长度, $h_{i, {\rm elec}}$为电极末端位置, $T$为运行时间, $u_i$为升降电机转速, $r$为升降机构等效齿轮半径^[9], $h_{i, {\rm pool}}(B_{i, 1}$, $B_{i, 2}, y_i)$为熔池高度, $B_{i, 1}$为原料杂质成分含量, $B_{i, 2}$为原料颗粒长度.

熔池电阻$R_{i, {\rm pool}}$^[10]为

其中, $\rho_{i, {{\rm pool}}}$为熔池电阻率, $D$为熔池直径.由式(7)和式(8), 得

对式(5)两边求导, 代入式(9)得到以升降电机转速$u_i$为输入, 电流$y_i$为输出的动态模型, 即

显然, 式(10)中${u_i}y_i^2$表明电流$y_i$与升降电机转速$u_i$之间为非线性关系, 且$y_i^2$系数中包含未知的熔池电阻率${\rho _{i, {{\rm pool}}}}({B_{i, 1}}, {B_{i, 2}})$和熔池高度非线性变化$\frac{{\rm d}{h_{i, {\rm pool}}}(B_{i, 1}, B_{i, 2}, y_i)}{{\rm d}t}$.由于工艺要求将电流控制在电流设定值附近, 电熔镁炉在工作点附近运行.将其在工作点附近线性化后, 被控对象(10)可由一阶线性模型与未建模动态的形式表示, 即

其中, $A(z^{-1})$, $B(z^{-1})$是关于$z^{-1}$的多项式, $A({z^{ - 1}})$ $=$ $1+{a_1}{z^{ - 1}}$, $B({z^{ - 1}}) = {b_{{\rm 0}}}$, ${\mathit{\boldsymbol{ h}}_{i, {\rm pool}}} \in {{\bf R} ^{{n_c}}}$, $n_c$未知. $v_i\left[ \cdot \right]$为高阶非线性项, 称未建模动态.

设计电流PID控制器为

其中, $H({z^{-1}}) = 1 - {z^{-1}}$, $G({z^{-1}}) = {g_0} + {g_1}{z^{-1}} + {g_2}{z^{-2}}$.将控制器(12)作用于被控对象(11), 得电流控制系统的闭环方程, 即

其中, $T({z^{-1}}) = (1-{z^{- 1}})A({z^{-1}}) + {z^{-1}}B({z^{-1}})$ $\times$ $G({z^{-1}})$, 由式(13)得到电流控制系统闭环方程的输出, 即

其中, ${T^*}({z^{ - 1}}) = T({z^{ - 1}}) - 1$.从式(14)可以看出, $y_i(k+1)$与过去时刻的电流$ - {T^*}({z^{ - 1}})y_i(k + 1)$、电流设定值$y^*$、未建模动态${v_i}[{y_i}(k), {u_i}(k), {\mathit{\boldsymbol{ h}}_{i, \rm {pool}}}]$相关.

将式(14)代入式(4)中可得功率变化率$\Delta p(k)$的动态模型.

其中, $C = \sqrt 3 U\cos \varphi$, $C\sum_{i = 1}^m (1 - {z^{ - 1}}){v_i}[{y_i}(k), $ ${u_i}(k), $ $\mathit{\boldsymbol{ h}}_{i, \rm {pool}}]$为高阶非线性项.式(15)表明了功率变化率$\Delta p(k)$与电流${y_i}(k)$之间的关系.

由式(5)可知影响电流${y_i}(k)$的直接因素是工作电阻${R_i}(k)$.将式(5)代入式(15)可得功率变化率$\Delta p(k)$与工作电阻${R_i}(k)$之间的关系.

再由式(7)和式(8)进一步得到功率变化率$\Delta p(k)$和电流影响因素(原料杂质成分含量$ B_{i, 1} $、颗粒长度$ B_{i, 2} $)之间的关系

其中, $ {f_{{R_i}}}\left( {{B_{i, 1}}(k), {B_{i, 2}}(k)} \right) $为关于原料杂质成分含量$ B_{i, 1} $、颗粒长度$ B_{i, 2} $的函数.

由式(3), (15)~(18)可以看出群炉需量与电流及电流因素之间的关系.当多台电熔镁炉出现原料杂质成分含量$ B_{i, 1}$增大和颗粒长度$ B_{i, 2} $变大的工况时, 会使熔池高度$ h_{i, {\rm pool}} $升高^[11], 电弧长度${L_{i, {{\rm arc}}}}$减小, 进而电阻$ R_i $减小, 电流$ y_i $变大, 功率变化率$ \Delta p$ $>$ $0 $, 需量$ \bar{P} $增大.由于电流控制系统调节电机转速$ u_i $使电弧长度$ L_{i, {\rm arc}} $增大, 又使电流$ y_i $减小到设定值$ y^* $附近, 需量$ \bar{P} $又减小, 导致需量$ \bar{P} $先增大然后减小, 出现需量尖峰.

$ u_i(k) $与过去时刻的电流$y_i(k)$相关, 可表示为

熔炼过程中熔池高度$ {\mathit{\boldsymbol{ h}}_{i, {\rm pool}}} $增大, 电流$y_i(k)$会随之增大, 因此$ {\mathit{\boldsymbol{ h}}_{i, {\rm pool}}} $与$y_i(k)$也相关, 即

所以$ {v_i}[{y_i}(k), {u_i}(k), {\mathit{\boldsymbol{ h}}_{i, {\rm pool}}}] $可简记为电流$ \mathit{\boldsymbol{ y}}_i $的非线性函数$ {V_i}({\mathit{\boldsymbol{ y}}_i}) $, $ {\mathit{\boldsymbol{ y}}_i} = {[{y_i}(k), {y_i}(k - 1), \cdots]^{{\rm T}}} $.将式(15)改写为

其中, $ C\sum_{i = 1}^m {(1 - {z^{ - 1}}){V_i}\left( {{\mathit{\boldsymbol{ y}}_i}} \right)} $为各炉电流的函数乘以熔炼电压$ U $和功率因数$ \cos\varphi $再累加, 与群炉功率有关.故可将线性组合$[ - z{T^*}({z^{ - 1}}) - 1]p(k) + B({z^{ - 1}})G({z^{ - 1}}){p^*} $之外的部分记为非线性函数$ \bar V(\mathit{\boldsymbol{ p}}) $, 则$ \Delta p(k) $的动态模型进一步改写为

其中, $ \bar V(\mathit{\boldsymbol{ p}}) $为关于$ \mathit{\boldsymbol{ p}} = {\left[ {p(k), p(k - 1), \cdots } \right]^{{\rm T}}} $和过去时刻的$ \bar V(\mathit{\boldsymbol{ p}}) $的未知非线性函数.至此, $ (k+1) $时刻的需量动态模型由式(3)和式(22)组成.

在需量的动态模型中, 影响需量的主要因素为原料杂质成分含量$ B_{i, 1} $和颗粒长度$ B_{i, 2} $.而熔池内的气泡、电弧闪变抖动等因素对需量的干扰本文视为均值为零, 方差有界的白噪声.在生产过程中电流控制器工作在稳定状态, $ B_{i, 1} $和$ B_{i, 2} $变化引起的电流变化是有界的, 即在一个闭集中.由于电流变化是有界的, 功率也有界, 可以采用神经网络的方法对功率变化率的非线性项进行估计.在此条件下, 本文提出一种数据与模型驱动的需量预报方法.

2.   需量预报方法

根据式(3)和式(22)组成的需量动态模型, 本文提出由线性模型、基于PACF输入变量决策的RBFNN未知非线性函数估计和需量计算模型组成的需量预报方法.

2.1   需量预报策略

需量动态模型中, 功率变化率$ \Delta p(k) $包含的线性部分$[ - z{T^*}({z^{ - 1}}) - 1]p(k) + B({z^{ - 1}})G({z^{ - 1}}){p^*} $已知, 将需量动态模型的式(22)改写为

其中, $ \Delta {p_1}(k) $和$ \bar V(k) $分别为$ \Delta p(k) $的线性部分和未建模动态, 即

其中, 参数$ a_1$, $b_0 $由式(11)结合实验确定, 控制器参数已知, 因此可计算求得线性部分$ \Delta p_1(k) $.因为$ \bar{V} (k) $的非线性函数$ f_v(\cdot) $的输入和输出都有界, 是一个闭集, 因此可基于文献[12-13]采用RBFNN估计出, 结合$\Delta p_1(k)$得到$ \Delta \hat{p}(k) $, 再代入式(3)中即可计算出$ (k+1) $时刻需量预报值$ \hat{\bar P}(k+1)$.

因此, 本文提出由线性模型、基于PACF输入变量决策的RBFNN未知非线性项估计和$(k+1)$时刻需量计算模型组成的预报模型结构, 如图 2所示.

图 2  电熔镁群炉需量预报方法结构框图

Fig. 2  The structure diagram of demand forecasting method for FMFG

下载: 全尺寸图片 幻灯片

1) 线性模型.线性部分$ \Delta p_1(k) $根据工艺实验确定被控对象的线性模型参数$ a_1$, $b_0 $和控制器参数$g_0$, $g_1$, $g_2 $计算.

2) 非线性项函数$ \bar V(k) $估计.根据采集的功率序列样本$ \{ p(i), i = 1, 2, \cdots , N\} $, 采用PACF方法分析, 决策$ \hat {\bar V}(k) $估计模型的功率历史数据输入个数$ n_f $和未建模动态历史数据输入个数$ n_v $.以过去的功率和非线性项为输入, 采用RBFNN估计出$\hat {\bar V}(k)$.

3) $ (k+1) $时刻需量计算模型.根据检测的$ p(k) $, $ p(k-n+1) $、非线性函数估计值$\hat {\bar V}(k)$和计算的线性部分$ \Delta p_1(k) $, 求出最终的需量预报值$\hat {\bar P}(k + 1)$.

2.2   需量预报算法

基于上述策略, 提出需量预报算法如下:

其中, 功率变化率的估计$ \Delta \hat p(k) $为线性部分$ \Delta p_1(k) $和非线性函数估计值$\hat {\bar V}(k)$之和, 即

2.2.1   线性模型

由式(22)可得$ \Delta p_1(k) $, 即

2.2.2   n_f和n_v决策算法

采用文献[3, 5]中的PACF决策$ n_f $和$ n_v $值.首先采集功率样本序列$ \{ p(i), i = 1, 2, \cdots , N\} $, 然后按文献[5]的计算步骤求得PACF函数$ \psi ({n_f}) $.当

则选择此时的$n_f$值.

决策$n_v$的方法类似, 具体决策条件如下:

2.2.3   V(k)的估计算法

由于$\bar V (k)$和过去时刻的$ \bar V(k - 1), \cdots, \bar V(k - {n_v}) $相关, 也和$ p(k), \cdots, p(k - {n_f} + 1) $相关.因此, $ \bar V(k) $表示以$ n_f $个过去时刻功率和$ n_v $个过去时刻非线性项为输入的未知非线性函数${f_v}( \cdot ) $, 即

根据文献[5], RBFNN输出当前的$ \hat{\bar V}（(k)$为

其中, $H$为RBFNN的隐含层节点的个数, ${\mathit{\boldsymbol{ d}}_j} =$ ${\rm diag}\{\omega ^a_{1j}, \omega ^a_{2j}, \cdots , \omega ^a_{({n_f} + {n_v})j}\}$, $\omega^a _{ij}$表示输入层神经元$i$与隐含层神经元$j$之间的连接权, 全部元素$\omega ^a_{ij}$组成连接权矩阵${{{\omega}}^a}$, $n_f+n_v$行$H$列即${{{\omega}}^a}$ $\in$ ${{\bf R} ^{({n_f} + {n_v}) \times H}}$, $\mathit{\boldsymbol{ x}} = [p(k), \cdots, p(k - {n_f} + 1), \bar V(k - 1)$, $\cdots$, $\bar V(k - {n_v})]^{\rm T}$向量的维数为${n_f} + {n_v}$, $\mathit{\boldsymbol{ x}}\in{{\bf R} ^{{n_f} + {n_v}}}$, $\sigma$为高斯函数的宽度, $\left\| \cdot \right\|$为欧几里得范数, $\mathit{\boldsymbol{ C}}_j$是第$j$个高斯函数中心点, ${{\mathit{\boldsymbol{\omega}}}^{b}}$为输出层权值向量, 维数为$H$, ${\mathit{\boldsymbol{\omega}} ^b} \in {{\bf R} ^H}$, $\beta$为偏置.

1) 输入层到隐含层的连接权矩阵${{{\omega}}^a}$的设定

由文献[14]知${{{\omega}}^a}$对RBFNN拟合精度的作用很小, 所以将${{{\omega}}^a}$所有元素设为1, 即

则${\phi _j}{{\rm (}}{\mathit{\boldsymbol{ d}}_j}\mathit{\boldsymbol{ x}}{{\rm ) = }}{\phi _j}{{\rm (}}\mathit{\boldsymbol{ x}}{{\rm )}}$.

2) 隐含层节点数$H$及对应中心点$\mathit{\boldsymbol{ C}}_j$和高斯函数宽度$\sigma$的选择

对训练样本$\{( {\bar V(k), {\mathit{\boldsymbol{ x}}_k}} ), k = 1, \cdots , N \}$进行5折交叉验证实验, 将5折后每4份数据样本数记为$M$, 即$M = \frac{4}{5}N$, 并作为隐含层中心点的候选集, 剩下$\frac{1}{5}N$的训练样本作为验证集. RBFNN对训练样本集合$\{ ( {\bar V(k), {\mathit{\boldsymbol{ x}}_k}} ), k = 1, \cdots , N\}$的回归映射表示如下^[15]:

其中,

其中, ${\mathit{\boldsymbol{\bar V}}} \in {{\bf R} ^N}$为导师信号, 即期望输出向量, ${\Phi}$ $\in$ ${{\bf R} ^{N \times {H}}}$为RBFNN隐含层输出矩阵, ${{\bf \Xi }} \in {{\bf R} ^N}$为RBFNN对训练样本建模的误差向量. $ {\phi _i}({\mathit{\boldsymbol{ x}}_j}) $为

${\mathit{\boldsymbol{ x}}_j} = [ p(j), \cdots , p(j - {n_f} + 1), \bar V(j - 1), \cdots, \bar V(j-$ ${n_v})]^{\rm T} $, ${\mathit{\boldsymbol{ C}}_i} = [ p(i), p(i - 1), \cdots , p(i - {n_f} + 1)$, $\bar V(i-$ $1)$, $\cdots$, $\bar V(i - {n_v})]^{\rm T}$.对矩阵${\Phi}$用Gram-Schmidt方法^[16]进行正交分解来逐一确定中心点, 步骤如下:

步骤1. 对于所有$1\le i\le M$, 计算如下公式:

得到$M$个候选样本作为中心点的贡献程度指标$[{ err}]_1^{(i)}$, 最大的贡献度指标为

选择

以上标$i_1$对应的样本${\mathit{\boldsymbol{ x}}_{{i_1}}}$作为第1个中心点

步骤H (2 ≤ H ≤ H_max). 对于所有的$1\le i$ $\le$ $M$, $i \ne {i_1}, \cdots, i \ne {i_{H - 1}}$, 其中, ${i_1}, \cdots, {i_{H - 1}}$为前面$H - 1$步已经选择作为中心点的样本下标.计算如下公式:

$M$个候选样本中的贡献度指标第$H$大的值为

选择

以上标${i_H}$对应的样本${\mathit{\boldsymbol{ x}}_{{i_H}}}$作为第$H$个中心点

当$H = {H_{\max }}$时, ${H_{\max }}$个中心点选择完毕.

节点数$H$逐一增加, 对候选集拟合的误差越来越小, 但当$H$增至某临界值$H_c$值时, 验证集拟合的误差开始增大, 出现过拟合现象.因此权衡预报精度和避免过拟合的要求, 选择节点数$H$.

当节点数$H$确定之后, 对完整训练样本进行正交分解得到贡献度较大的前$H$个中心点.

隐层节点的高斯函数采用相同的$\sigma$, 采用文献[17]的方法, 在一个合理区间$[{\sigma _{\min }}, {\sigma _{\max }}]$内进行5折交叉验证实验, 选择验证集预报误差收敛速度变化不明显的临界值${\sigma _c}$.

3) 输出层权重向量${{\mathit{\boldsymbol{\omega}}}^b}$和偏置$\beta$的求解

采用最小二乘法求如下方程:

其中, $\mathit{\boldsymbol{ s}} = {[1, \cdots , 1]^{\rm T}}$, $\mathit{\boldsymbol{ s}} \in {{\bf R} ^N}$, 求得${{\mathit{\boldsymbol{\omega}}}^b}$和$\beta$.

其中, ${\left[ \cdot \right]^† }$为广义逆矩阵, 由奇异值分解方法解得.至此, $ \hat{\bar V}(k) $由RBFNN的估计算法如下:

其中, $n_f$和$n_v$由式(28)和式(29)给出, $H$, $\mathit{\boldsymbol{ C}}_j$和$ \sigma $由式(44)~(46)给出, ${{\mathit{\boldsymbol{\omega}}}^b}$和$\beta$由式(48)给出.

4) $\bar V(k)$预报模型的在线更新

由于群炉生产状态会发生迁移, 需根据当前预报误差对预报模型的参数${\mathit{\boldsymbol{\omega}}^b}$进行更新.设${{\mathit{\boldsymbol{\omega}}}^b}$在$k$时刻为${\mathit{\boldsymbol{\omega}}}^b_k$, 求$(k+1)$时刻的权值向量$\mathit{\boldsymbol{\omega }}^b_{k + 1}$.

采用文献[18]递推正交最小二乘法的更新步骤如下:

步骤1. $k$时刻${\mathit{\boldsymbol{\omega}}^b}$的值为$\mathit{\boldsymbol{\omega}}^b_k$, 隐含层输出矩阵为${{\Phi} _k}$, 对${{\Phi} _k}$正交分解

其中, ${{Q}_k} \in {{\bf R} ^{N \times N}}$是正交矩阵, ${{R}_k} \in {{\bf R} ^{{H} \times {H}}}$是上三角矩阵, ${O} \in {{\bf R} ^{(N - H) \times {H}}}$是零矩阵.

步骤2. $k+1$时刻递推正交分解求得${{R}_{k + 1}}$.

其中, $\mathit{\boldsymbol{\vartheta}} _{\;k + 1}^{\rm T} = [{\phi _1}({\mathit{\boldsymbol{ x}}_{k{{\rm + }}1}}), \cdots, {\phi _H}({\mathit{\boldsymbol{ x}}_{k{{\rm + }}1}})]$.再由${Q}_{k + 1}$计算${ \Omega} _{k + 1}$

其中, $\varepsilon (k{{\rm + }}1) =\bar V(k + 1{{\rm ) - }}\mathit{\boldsymbol{\vartheta}} _{\;k + 1}^{\rm T}\mathit{\boldsymbol{\omega }}^b_k$.

步骤3. 由最小二乘法求得$\Delta \mathit{\boldsymbol{\omega }}^b_k$.

代入式(50)得到$\mathit{\boldsymbol{ \omega }}^b_{k + 1}$, 更新$k=k+1$, 返回步骤2.

综上所述, 本文预报方法由式(25)给出, 其中功率变化率预报值$ \Delta \hat p (k) $由式(26)给出.式(26)中的线性部分$ \Delta p_1 (k)$由式(27)给出. $\hat {\bar V} (k)$为$\Delta p (k)$的未建模动态的神经网络估计, 由式(49)给出.

3.   仿真实验与工业实验

3.1   仿真实验

将本文所提预报方法对某电熔镁砂厂1号电力变压器的1个炉次的群炉需量数据进行仿真实验, 以验证方法的有效性.实验数据中功率的采样周期为7秒, 需量参数$ n=30 $, 群炉生产台数为4台, 当天设定的需量限幅值为22 100 kVA.从该炉次4 620组数据中选取4 000组数据$p(i)$, $i = 1, \cdots, 4 000 $进行仿真实验, 其中$p(i)$, $i = 1, \cdots, 2 000$为训练集, $p(i)$, $i = 2 001, \cdots, 4 000$为测试集.根据下式求得建模所用的$ \Delta p(i)$

3.1.1   线性模型∆p₁(i)

根据现场工艺实验在工作点$ y^*=15 000 $ A, $ u^*$ $=$ $6.1918\times10^{-4} $ rad $\cdot $ s$^{-1}$处线性化被控对象模型, 计算得到参数$ a_1$, $b_0 $如下:

控制器参数$ g_0 = 6.300 035 $, $g_1 = -11.9$, $g_2 = 5.6$, 则功率变化率线性模型$ \Delta p_1(i) $计算公式为

其中, ${p^*} = \sqrt 3 mU\cos \varphi {y^*}$, $m = 4$, $U = 190$, $\cos \varphi$ $=$ $0.92$.

3.1.2   n_f和n_v的决策

对2 000组训练数据组成的功率时间序列求PACF序列值为

PACF的95 %的置信区间为$[ - {{\rm 0}}{{\rm .0448}}, {{\rm 0}}{{\rm .0448}}]$, 当$n_f$ $=4$时, PACF值落入95 %置信区间内, 因此设定

类似地, 对2 000组训练数据组成的未建模动态求PACF序列值为

PACF的95 %的置信区间为$[ - {{\rm 0}}{{\rm .0448}}, {{\rm 0}}{{\rm .0448}}]$, 设定

3.1.3   基于RBFNN的V(i)估计

根据下式整理RBFNN训练所需的输入和输出数据:

其中, $ {\bf x} = [p(i), \cdots , p(i - {n_f} + 1), \bar V(i -1), \cdots$, $\bar V(i$ - ${n_v})]^{{\rm T}}$.先将$ \sigma $由经验选为0.8, 将节点数$ H $由1逐一增加至$ H_{\max} = 200 $.

如图 3所示, 进行5折交叉验证实验, 可以看出一开始随着节点个数增加, 候选集和验证集的预报误差指标(MAPE, RMSE)都在减小, 但是当节点数$ H >60$时, 虽然候选集的预报误差指标依然在减小, 但是验证集的预报误差指标开始增大.因此根据实验结果选择隐含层节点数

图 3  隐层节点数的交叉验证

Fig. 3  Cross-validation of the number of hidden nodes

下载: 全尺寸图片 幻灯片

如图 4所示, 在区间范围$[{\sigma _{\min}}, {\sigma _{\max }}] =$ $[0.01$, $2.00]$内, 通过5折交叉验证实验, 随着$ \sigma $增大候选集和验证集的预报误差指标都相应减小, 当$ \sigma >1.2$时预报误差指标减小的幅度很小.因此根据实验结果选择高斯函数宽度

图 4  高斯函数宽度的交叉验证

Fig. 4  Cross-validation of the width of Gaussian function

下载: 全尺寸图片 幻灯片

确定$H$和$\sigma$之后, 对整个训练集数据进行正交分解, 得到$H$个中心点$ \mathit{\boldsymbol{ C}}_j $的初始值${\mathit{\boldsymbol{ C}}_j}(0)$, $j =1$, $2$, $\cdots$, $60 $, 即

由式(48)求出权值向量${{\mathit{\boldsymbol{\omega}}}^b}$的初始值${{\mathit{\boldsymbol{\omega }}}^b(0)}$和偏置$\beta$初始值$\beta(0)$为

3.1.4   ∆p(k)预报算法验证

将离线训练后的预报初模型对剩余的2 000组数据进行验证, 如图 5所示.可以看出功率变化率的预报值基本包括功率变化率中有规律的信息.功率变化率预报误差为$\Delta {e_p}(k) = \Delta p(k) - \Delta \hat p(k)$, 可以看出误差$ \Delta {e_p}(k) $序列明显呈现一个白噪声的特性.统计的$\Delta p(k)$预报误差指标如表 1所示, 预报误差序列的方差为1.1481E+6, 均方根误差为1 071.3.

图 5  ∆p(k)预报验证曲线

Fig. 5  Forecast validation curves of ∆p(k)

下载: 全尺寸图片 幻灯片

表 1  ∆p(k)预报误差指标

Table 1  Forecast error indicators of ∆p(k)

方差	RMSE
1.1481E+6	1 071.3

下载: 导出CSV 

| 显示表格

3.1.5   P(k + 1)预报算法验证

验证实验中需量有三次较明显的先上升后下降的趋势, 分别是: $k=171$~$220$时段, 如图 6所示; $k$ $=571$~$620$时段, 如图 7所示; $k=1 356$~$1 405$时段, 如图 8所示.

图 6  时段1需量曲线

Fig. 6  Demand curve for the 1st time period

下载: 全尺寸图片 幻灯片

图 7  时段2需量曲线

Fig. 7  Demand curve for the 2nd time period

下载: 全尺寸图片 幻灯片

图 8  时段3需量曲线

Fig. 8  Demand curve for the 3rd time period

下载: 全尺寸图片 幻灯片

图 6中, $k=171$时, 需量实际值为$\bar P(k)=21 427$ kW, 对应需量预报值$\hat{\bar P}(k)=21 421$ kW; 当$k$ $=187$时, 需量上升至该时段内最大值$\bar P(k)=21 833$ kW, 对应需量预报值$\hat{\bar P}(k)=21 842$ kW; 之后由于电流控制系统的调节作用, $\bar P(k)$开始下降, 当$k=194$时, 需量降低至$\bar P(k)=21 598$ kW, 对应需量预报值为$\hat{\bar P}(k)=21 620$ kW.

图 7中, $k=571$时, 需量实际值为$\bar P(k)=21 678$ kW, 对应需量预报值$\bar P(k)=21 687$ kW; $k$ $=$ $592$时, 需量为该时段最大值$\bar P(k)=22 067$ kW, 需量预报值$\hat{\bar P}(k)=22 080$ kW; 之后$\bar P(k)$开始下降, $k=620$时, 需量$\bar P(k)=21 730$ kW, 对应需量预报值为$\hat{\bar P}(k)=21 704$ kW.

图 8中, $k=1 356$时, 需量实际值为$\bar P(k)=21 450$ kW, 对应需量预报值$\hat{\bar P}(k)=21 445$ kW; $k$ $=1 385$时, 需量为该时段最大值$\bar P(k)=22 011$ kW, 需量预报值$\hat{\bar P}(k)=22 006$ kW; 之后$\bar P(k)$开始下降, $k=1 399$时, 需量$\bar P(k)=21 667$ kW, 对应需量预报值为$\hat{\bar P}(k)=21 633$ kW.

从图 6~8可以看出, 在三个时段中需量实际值和需量预报值的上升和下降趋势基本类似.为定量分析预报方法的性能, 本文采用误差方差、预报精度百分比(Percent better, PB)^[19]、均方根误差(Root mean square error, RMSE)^[20-23]、平均绝对百分误差(Mean absolute percentage error, MAPE)^[3]作为评估预报性能的指标, 计算公式如下:

其中, $ {\bar P_{0.5}}(k) = (1 + 0.5 \% ) \times \bar P(k) $, $N$为误差序列的样本数量.对包含图 6~8时间段测试集2 000个采样点的预报实验误差序列由式(68)进行计算, 对应预报性能指标见表 2, 误差序列的方差为1 275.7, 预报精度百分比为97.55 %, RMSE为35.7104, MAPE为0.1054 %, 与文献[5]中预报方法对比, 本文的方法明显减小了预报误差的方差和RMSE, 提高了预报精度百分比.由式(25)可知, 需量预报算法的复杂度和精度取决于$ \Delta p(k) $的估计.由于利用了所建立的$\Delta p(k)$动态模型式(22), 将$\Delta p(k)$的估计分为可精确计算的线性部分$ \Delta p_1(k) $和未建模动态$ \bar V(k) $.其中, $ \Delta p_1(k) $由式(27)计算, 未知的$ \bar V(k) $由式(49)采用神经网络估计.因此, 与采用神经网络整体估计$\Delta p(k)$的文献[5]相比, 本文所提方法对线性部分利用已知信息进行精确四则运算, 没有增加神经网络估计算法的复杂度, 同时使预报精度明显提高, 满足了现场预报的要求.

表 2  需量预报误差指标

Table 2  Forecast error indicators of demand

方法	方差	PB (%)	RMSE	MAPE (%)
文献[5]	1 533.0	97.05	39.1921	0.0979
本文	1 275.7	97.55	35.7104	0.1054

下载: 导出CSV 

| 显示表格

3.2   工业实验

将提出的预报方法应用于某电熔镁砂厂1号电力变压器负载的电熔镁群炉需量监控过程, 如图 1所示.该过程主要由1号电力变压器(型号: SF9-22500/66, 额定容量: 22 500 kVA), 高压断路器(10 kV, DJS-10机械闭锁), 5台电炉变压器(型号: HKS-4500/10, 额定容量: 4 500 kVA), 5台电熔镁炉(直径2.7 m, 高3 m), 升降电机(型号: YVP160M-4, 标称功率11 kW), PID电流控制器(CPU型号: 313-6CF03-0AB0).当天设定需量限幅值为21 800 kW, 群炉生产台数为4台.该厂对需量的管控动作分为切断和恢复两种.当需量实际值超过需量限幅值且超过时间大于4个采样周期(28秒)时, 进行切断操作.当需量实际值低于需量限幅值时, 对该台电熔镁炉进行恢复供电的操作.

3.2.1   需量预报系统简介

采用本文所提的预报方法研制了需量预报系统, 硬件平台为图 1中的研华IPC-7120需量监控计算机以及Siemens CP5621通讯板卡.

需量预报软件平台包括: STEP7-Micro/WIN编程软件, PC Access OPC服务器软件, SIMATIC WinCC Explorer过程监视软件, jdk1.7版本的JAVA软件开发工具包, JAVA集成开发环境Eclipse各1套以及Windows 7操作系统.基于上述软件平台和本文所提算法研制的需量预报软件的界面如图 9所示, 其中实线表示需量实际值, 虚线表示需量预报值.

图 9  需量预报软件界面

Fig. 9  The interface of demand forecasting software

下载: 全尺寸图片 幻灯片

3.2.2   预报模型的参数选择

预报模型的参数为仿真实验结束时的预报模型参数, 即模型输入变量个数${n_f} = 4$, ${n_v} = 7$, 隐含层节点数$H=60$, 高斯函数宽度$\sigma = 1.2$, 其他参数值具体如下, 中心点${ {\mathit{\boldsymbol{ C}}}_j}$, $j = 1, 2, \cdots, 60$为

输出层的权值向量${{\mathit{\boldsymbol{\omega}}}^b}$为

偏置$\beta$为

3.2.3   实验结果

在时段21 : 00 : 00~05 : 22 : 29进行的工业实验中, 正常情况下需量实际值和预报值曲线如图 10~12所示.

图 10  时段1需量曲线

Fig. 10  Demand curve for the 1st time period

下载: 全尺寸图片 幻灯片

图 11  时段2需量曲线

Fig. 11  Demand curve for the 2nd time period

下载: 全尺寸图片 幻灯片

图 12  时段3需量曲线

Fig. 12  Demand curve for the 3rd time period

下载: 全尺寸图片 幻灯片

在时段1 (22 : 48 : 23~22 : 54 : 13)中, $\bar P(k)$由22 : 48 : 23的20 801 kW上升至22 : 51 : 11的21 689 kW, $\hat{\bar P}(k)$由22 : 48 : 23的20 810 kW上升至22 : 51 : 11的21 690 kW; 之后$\bar P(k)$下降至22 : 54 : 13的21 104 kW, $\hat{\bar P}(k)$降至21 123 kW, 如图 10所示.

在时段2 (23 : 18 : 43~23 : 24 : 33)中, $\bar P(k)$由23 : 18 : 43的21 253 kW上升至23 : 21 : 10的21 783 kW, $\hat{\bar P}(k)$由23 : 18 : 43的21 247 kW上升至23 : 21 : 10的21 786 kW; 之后$\bar P(k)$下降至23 : 24 : 33的21 041 kW, $\hat{\bar P}(k)$降至21 049 kW, 如图 11所示.

在时段3 (00 : 46 : 13~00 : 52 : 03)中, $\bar P(k)$由00 : 46 : 13的21 122 kW上升至00 : 49 : 08的21 780 kW, $\hat{\bar P}(k)$由00 : 46 : 13的21 127 kW上升至00 : 49 : 08的21 813 kW; 之后$\bar P(k)$下降至00 : 52 : 03的21 375 kW, $\hat{\bar P}(k)$降至21 369 kW, 如图 12所示.

超限拉闸情况下, 需量实际值和预报值曲线如图 13所示, 在22 : 38 : 00为21 826 kW时出现了拉闸操作.拉闸前后一段时间(10个采样周期)需量实际值与需量预报值及误差见表 3.

图 13  超限拉闸时段需量曲线

Fig. 13  Demand curve for cut off time period

下载: 全尺寸图片 幻灯片

表 3  超限拉闸时段需量预报误差

Table 3  Demand forecast errors during cut off time period

时间	需量实际值(kW)	需量预报值(kW)	误差(kW)
22 : 36 : 50	21 626	21 635	-9
22 : 36 : 57	21 654	21 650	4
22 : 37 : 04	21 692	21 685	7
22 : 37 : 11	21 728	21 725	3
22 : 37 : 18	21 754	21 751	4
22 : 37 : 25	21 788	21 787	1
22 : 37 : 32	21 812	21 829	-17
22 : 37 : 39	21 834	21 849	-15
22 : 37 : 46	21 835	21 839	-4
22 : 37 : 53	21 833	21 829	4
22 : 38 : 00	21 826	21 822	4
22 : 38 : 07	21 691	21 819	-128
22 : 38 : 14	21 510	21 463	47
22 : 38 : 21	21 335	21 318	17
22 : 38 : 28	21 191	21 176	15
22 : 38 : 35	21 049	21 056	-7
22 : 38 : 42	20 908	20 907	1
22 : 38 : 49	20 751	20 769	-18
22 : 38 : 56	20 581	20 577	4
22 : 39 : 03	20 407	20 385	22

下载: 导出CSV 

| 显示表格

由表 3看出, 在22 : 36 : 50~22 : 37 : 53拉闸前半段时间内需量实际值为上升趋势, 对应的需量预报值同样为上升趋势, 趋势相同.在22 : 37 : 32时$\bar P(k)$为21 812 kW, 开始超出需量限幅值, 对应$\hat {\bar P}(k)$为21 829 kW; 22 : 37 : 32之后连续4次采样时刻需量为21 834 kW, 21 835 kW, 21 833 kW, 21 826 kW, 对应的预报值为21 849 kW, 21 839 kW, 21 829 kW, 21 822 kW, 拉闸前4次采样时刻的需量实际值和需量预报值的趋势相同, 工厂在22 : 38 : 00进行拉闸操作, 拉闸时的$\bar P(k)$为21 826 kW, $\hat {\bar P}(k)$为21 822 kW; 从22 : 38 : 07~22 : 39 : 03拉闸后半段时间$\bar P(k)$下降迅速, 对应的$\hat {\bar P}(k)$也开始下降.

上述需量变化过程中, 虽然需量实际值超过了需量限幅值, 但没有超过需量峰值.从图 13和表 3可以看出, 需量预报值在拉闸时刻附近趋于缓慢下降趋势, 如果不拉闸需量实际值有可能会下降到限幅值以下, 那么这次拉闸就可能是不必要的拉闸.因此, 将本文所提预报方法与现有需量监控系统结合将有助于减少尖峰引起的不必要拉闸, 对提高需量限幅值的设定值, 从而提高生产过程的电能利用率有一定的指导作用.

超限拉闸时段的需量预报误差变化曲线如图 14所示, 由于预报方法是基于拉闸时和之前的功率数据来预报, 预报的是不拉闸情况的下一时刻需量, 所以在拉闸后下一时刻需量预报误差与之前预报误差相比发生了较大的波动变化; 随着时间推移, 预报模型的输入变量开始包含拉闸之后的功率数据, 在经过2个采样时刻之后, 需量预报误差又逐渐减小, 说明预报模型可以及时地随着需量动态变化而在线调整.

图 14  超限拉闸时段需量预报误差变化曲线

Fig. 14  Demand forecast error curve during cut off time period

下载: 全尺寸图片 幻灯片

根据现场实验的数据, 对恢复动作时间段的需量预报进行了仿真实验, 如图 15所示.图中$a$处22 : 39 : 17时刻对断电的电熔镁炉进行恢复供电的动作; 图中$b$处22 : 40 : 48时刻该台电熔镁炉由于刚恢复供电炉内工况不稳, 电弧闪灭出现了跳闸.图中$c$处22 : 43 : 22该台电熔镁炉再次恢复供电.

图 15  恢复供电动作下的需量预报

Fig. 15  Demand forecast curve during restore operations

下载: 全尺寸图片 幻灯片

在$a$, $c$两处恢复供电的动作之后的需量预报误差变化如图 16所示. $a$处首次恢复供电, 由于该炉停止供电一段时间造成炉内熔池温度下降, 所以恢复供电后的熔炼电流值要低于切断之前的电流值.所以需量会有上升但不会达到断电之前的需量值.当3个采样周期之后需量预报误差减小至正常范围内.在$b$处跳闸之后需量又开始突然下降, 需量预报误差经过2个采样周期之后减小至正常范围.在$c$处再次恢复供电, 需量开始减缓下降并逐渐开始上升, 需量预报误差经过3个采样周期减小至正常范围.说明对于恢复供电的动作, 本文提出的预报方法的预报误差也会产生较大波动, 但随着时间推移预报误差又减小到正常范围内.

图 16  恢复供电动作下的需量预报误差

Fig. 16  Demand forecast error curve during restore operations

下载: 全尺寸图片 幻灯片

对上述21 : 00 : 00~05 : 22 : 29时间段的实验结果进行性能分析, 实验的需量预报误差的概率分布和白度分析如图 17所示, 可以看出需量预报误差大部分在$ [-70, 70] $之内, 误差序列的自相关系数大部分在白度测试^[24]的95 %置信区间$ [-0.0299$, $0.0299] $之内, 因此可将预报误差视为白噪声序列, 说明了需量预报值的可靠性.

图 17  工业实验需量预报误差白度分析

Fig. 17  The whiteness analysis of demand forecast error in industrial experiment

下载: 全尺寸图片 幻灯片

根据式(68)计算得到工业实验的需量预报方法性能指标见表 4, 方差为1 049.8, 预报精度百分比为97.68 %, RMSE为32.3981, MAPE为0.0996 %.从表 4可以看出, 采用本文提出的需量预报方法能够较为准确地预报出下一时刻的需量值, 对提高需量限幅值的设定值具有指导意义.

表 4  工业实验需量预报误差指标

Table 4  Demand forecast error indicators of industrial experiment

方差	PB (%)	RMSE	MAPE
1 049.8	97.68	32.3981	0.0996

下载: 导出CSV 

| 显示表格

4.   结论

本文建立了$(k+1)$时刻电熔镁群炉需量模型, 利用$(k+1)$时刻需量$ \bar P(k) $取决于功率变化率$ \Delta p(k) $, 以及$\Delta p(k)$取决于电流控制系统输出电流的特点, 提出了电熔镁群炉需量预报方法.该方法由线性模型、基于PACF输入变量个数决策的RBFNN未知非线性函数估计和$(k+1)$时刻需量$\bar P(k+1)$计算模型组成.通过某电熔镁砂厂实际数据的仿真实验和工业实验表明, 该方法可准确预报需量变化趋势, 不仅对预报需量尖峰防止不必要拉闸有实际意义, 而且对于工业过程控制系统的运行指标的预报具有一定参考价值.