东华理工大学理学院，江西 南昌

收稿日期：2021年9月13日；录用日期：2021年10月6日；发布日期：2021年10月15日

摘要

灰色预测模型是一类处理小样本数据预测的有效方法，其中单变量灰色预测模型 $GM\left(1,1,t^{\alpha }\right)$ 是一个重要的研究对象。 $FAGM\left(1,1,t^{\alpha }\right)$ 模型是在基于一阶累加的 $GM\left(1,1,t^{\alpha }\right)$ 模型的基础上引入分数阶累加形式的预测模型。然而，该模型的精度不够高，且容易存在过拟合现象。本文结合Lasso回归中的L1正则化思想，对分数阶累加的灰色时间幂模型 $FAGM\left(1,1,t^{\alpha }\right)$ 进行正则化，提出正则化分数阶灰色时间幂预测模型 $LFAGM\left(1,1,t^{\alpha }\right)$，使用坐标下降算法替代最小二乘估计来求解模型的参数。同时，使用灰狼优化算法(GWO)搜索 $LFAGM\left(1,1,t^{\alpha }\right)$ 模型的最优非线性参数。并基于中国农业耕地灌溉面积情况(2008~2019年)进行算例分析，结果表明， $LFAGM\left(1,1,t^{\alpha }\right)$ 模型具有更高的预测精度。

关键词

分数阶灰色预测模型，L1正则化，灰狼优化算法，坐标下降

L1 Regularized Fractional Grey Time Power Forecasting Model Based on GWO Optimization

Bangqiang Dong, Yanfei Zhang^*, Muhua Ding, Ping Chen

School of Science, East China University of Technology, Nanchang Jiangxi

Received: Sep. 13^th, 2021; accepted: Oct. 6^th, 2021; published: Oct. 15^th, 2021

ABSTRACT

Grey prediction model is an effective method for dealing with small sample data prediction, in which univariate grey prediction model $GM\left(1,1,t^{\alpha }\right)$ is an important research object. The $FAGM\left(1,1,t^{\alpha }\right)$ model is a prediction model based on the first-order cumulative $GM\left(1,1,t^{\alpha }\right)$ model by introducing the fractional-order cumulative form. However, the accuracy of the model is not high enough and it is prone to over-fitting. This paper combines the L1 regularization idea in Lasso regression to regularize the fractional cumulative gray time power model $FAGM\left(1,1,t^{\alpha }\right)$, and proposes a regularized fractional gray time power prediction model $LFAGM\left(1,1,t^{\alpha }\right)$, using coordinate descent algorithm instead of least squares estimation to solve the model parameter. At the same time, Grey Wolf Optimization (GWO) is used to search the optimal nonlinear parameters of $LFAGM\left(1,1,t^{\alpha }\right)$ model. Based on the agricultural farmland irrigation area in China (2008~2019), the results show that the $LFAGM\left(1,1,t^{\alpha }\right)$ model has higher prediction accuracy.

Keywords:Fractional Gray Prediction Model, L1 Regularized, Grey Wolf Optimization Algorithm, Coordinate Descent

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

灰色预测是对既含有已知信息又含有不确定信息的系统进行预测，对在一定范围内变化的、与时间有关的灰色过程进行预测控制的理论。灰色预测模型由邓聚龙教授 [1] 首次提出，在专家学者的不懈努力和完善下，弥补了由于数据少所带来的预测问题，尤其在实际应用中如何提高预测模型的精度一直是广大学者研究的热点问题。在初始值优化角度方面，Ding [2] 等人对初始值进行了非线性优化，使得所构建的模型的初始值各分量加权系数都具备了可调节性的特点；王义闹 [3] 等人同样以误差平方和最小值为目标函数，优化了指数模型的一个参数，并且得到了最小的平方误差和；Wang [4] 等人提出了一种基于矩阵的算法，建立了初始条件和发展系数间的关系。在模型参数优化的角度方面，通过怎么样的方式去优化灰色预测模型是当下灰色预测模型较为重要的内容，例如Liu [5] 等人利用Weibull累积分布函数构建了一个双形参数来替代原来较为传统的常数参数值，并且通过拟合实例验证了所构建的模型具有较高的拟合度；Zhao [6] 等人使用蚁狮优化算法优化了灰色预测模型参数，很大程度地提升了模型的预测精度，同时证明了该算法可行性；Yi-Shian-Lee [7] 设置误差的最小绝对值为目标函数，利用遗传算法来提高预测模型的精度。从灰色预测模型的累加方式上看，其中Wen Kunli [8] 提出了局部的序列灰色累加方式；刘解放 [9] 等人提出并使用了分数阶反向累加的离散灰色模型，通过实例证明了其可行性和优越性；谢波 [10] 等人通过在 $GM\left(1,1,t\right)$ 预测模型基础上选取合适的数据维度来增加数据长度，然后取多期的预测值的平均值作为最终的预测值，进一步提高了模型的精度，减少了数据波动；陈英超 [11] 等人则是通过将累加生成变为卷积变换，从而建立带有时间项的 $GM\left(1,1,t\right)$ 灰色预测模型。但随着灰色系统越发复杂，传统的灰色预测模型已经不能单单只应用指数特征去进行描述，所以应该在这些较为复杂的系统中进行内在的演化与分析，钱吴永 [12] 等人建立了带有时间幂次项的 $GM\left(1,1,t^{\alpha }\right)$ 模型并研究讨论了参数 $\alpha$ 在不同取值下模型的性质，并对沿海高速路软土地地基沉降进行了模拟分析；而Wu Lifeng [13] 等人将灰色预测模型累加次数从整数过度到分数，构建了分数阶灰色 $FGM\left(1,1\right)$ 预测模型；Wu [14] 等人引入了分数阶的累加方式，建立了分数阶含有时间幂次项的灰色 $FAGM\left(1,1,t^{\alpha }\right)$ 预测模型，并对中国卫生支出进行了研究分析。

以上文献都是基于传统的灰色预测模型进行优化建模，并于不同优化角度都保留了最小二乘法，对于现实数据处理往往存在过拟合问题。本研究则通过在分数阶累加灰色时间幂预测模型 $FAGM\left(1,1,t^{\alpha }\right)$ 建模基础上使用L1正则化，即在传统的最小二乘法中的损失函数中引入惩罚项，构建正则化预测模型 $LFAGM\left(1,1,t^{\alpha }\right)$，使用坐标下降算法求解模型的参数。并对于超参数问题，通过灰狼优化算法(GWO)优化调整L1正则化参数 $\lambda$，$LFAGM\left(1,1,t^{\alpha }\right)$ 模型的非线性参数 $\alpha$ 以及分数阶累加阶数r，从而避免由最小二乘法带来的过拟合问题。最后，以中国农业耕地灌溉面积情况(2008~2019年)进行了仿真模拟，验证了方法的可行性。

2. 分数阶灰色预测模型

2.1. $GM\left(1,1,t^{\alpha }\right)$ 模型

$GM\left(1,1,t^{\alpha }\right)$ 模型是由钱吴永等人首次提出。该模型将灰色作用量由常数b变为了一个关于时间呈现非线性变化的幂函数形式，扩展了单变量灰色预测模型的适用范围，很大程度上提高了模型的拟合和预测精度。设 $x^{\left(0\right)}={\left(x^{\left(0\right)}\left(1\right),x^{\left(0\right)}\left(2\right),\cdots ,x^{\left(0\right)}\left(n\right)\right)}^{\text{T}},n\ge 4$ 为非负原始数据列。对 $x^{\left(0\right)}={\left(x^{\left(0\right)}\left(1\right),x^{\left(0\right)}\left(2\right),\cdots ,x^{\left(0\right)}\left(n\right)\right)}^{\text{T}}$ 进行一次累加得到新的生成数据列 $x^{\left(1\right)}$ ( $x^{\left(0\right)}$ 的1-AGO序列)为

$x^{\left(1\right)}={\left(x^{\left(1\right)}\left(1\right),x^{\left(1\right)}\left(2\right),\cdots ,x^{\left(1\right)}\left(n\right)\right)}^{\text{T}}$ (1.1)

其中 $x^{\left(1\right)}\left(k\right)={\displaystyle {\sum }_{i=1}^k{x}^{\left(0\right)}}\left(i\right),k=1,2,\cdots ,n$。令 $z^{\left(1\right)}$ 为一阶累加序列 $x^{\left(1\right)}$ 的紧邻均值生成序列，即

$z^{\left(1\right)}={\left(z^{\left(1\right)}\left(2\right),z^{\left(1\right)}\left(3\right),\cdots ,z^{\left(1\right)}\left(n\right)\right)}^{\text{T}}$ (1.2)

其中 $z^{\left(1\right)}\left(m\right)=\delta x^{\left(1\right)}\left(m\right)+\left(1-\delta \right)x^{\left(1\right)}\left(m-1\right),m=2,3,\cdots ,n$ 且 $\delta =0.5$。

设灰色 $GM\left(1,1,t^{\alpha }\right)$ 模型的白化方程为

$\frac{\text{d}{x}^{\left(1\right)}\left(t\right)}{\text{d}t}+a{x}^{\left(1\right)}\left(t\right)=b{t}^{\alpha }+c,\alpha >0$ (1.3)

公式(1.3)中的a为 $GM\left(1,1,t^{\alpha }\right)$ 模型的发展系数， $b{t}^{\alpha }+c$ 为灰色作用量。因此， $GM\left(1,1,t^{\alpha }\right)$ 模型的基本形式为

$x^{\left(0\right)}\left(k\right)+a{z}^{\left(1\right)}\left(k\right)=b{k}^{\alpha }+c,k=2,3,\cdots ,n$ (1.4)

方程(1.4)同时也被称为灰色微分方程。注意到，线性方程(1.4)为一个超定方程组，根据最小二乘法，我们可得参数向量 $\beta ={\left(a,b,c\right)}^{\text{T}}$ 的估计式为

${\left(a,b,c\right)}^{\text{T}}={\left({\Delta }^{\text{T}}\Delta \right)}^{-1}\left({\Delta }^{\text{T}}\Psi \right)$ (1.5)

其中

(1.6)

其中n表示用于建立 $GM\left(1,1,t^{\alpha }\right)$ 模型的样本个数。

2.2. 分数阶 $FAGM\left(1,1,t^{\alpha }\right)$ 预测模型

一般地，一阶累加方式为分数阶累加方式的一个特例，并且分数阶累加能够更加体现数据序列内在的规律，更不易破坏数据间潜在的关联。另外，分数阶累加相比一阶累加更加灵活，能适应更多的数据。因此，在 $GM\left(1,1,t^{\alpha }\right)$ 模型的基础上，将分数阶累加考虑进来，能够显著提高 $GM\left(1,1,t^{\alpha }\right)$ 模型的预测精度和扩大使用范围。由于一阶累加的 $GM\left(1,1,t^{\alpha }\right)$ 模型与分数阶累加的 $FAGM\left(1,1,t^{\alpha }\right)$ 模型在建模过程中只是累加方式不同，由此我们先引入分数阶累加和累减的定义。

定义1：设 $X^{\left(0\right)}={\left(x^{\left(0\right)}\left(1\right),x^{\left(0\right)}\left(2\right),\cdots ,x^{\left(0\right)}\left(n\right)\right)}^{\text{T}}$ 为n个观测到的非负原始数据序列，T为转置符号。 $X^{\left(r\right)}$ 是 $X^{\left(0\right)}$ 的r阶累加生成序列(r-AGO)序列。其中 $x^{\left(r\right)}\left(k\right)={\displaystyle {\sum }_{i=1}^k{x}^{\left(r-1\right)}\left(i\right)}$，$k=1,2,\cdots ,n$。注意到，为了形式更加紧凑，该累加可以改写为更加方便的矩阵格式

$X^{\left(r\right)}=A^r{X}^{\left(0\right)}$ (1.7)

其中

$A^r={\left(\begin{array}{ccccc}\left[\begin{array}{c}r\\ 0\end{array}\right]& 0& 0& \cdots & 0\\ \left[\begin{array}{l}r\hfill \\ 1\hfill \end{array}\right]& \left[\begin{array}{l}r\hfill \\ 0\hfill \end{array}\right]& 0& \cdots & 0\\ \left[\begin{array}{l}r\hfill \\ 2\hfill \end{array}\right]& \left[\begin{array}{l}r\hfill \\ 1\hfill \end{array}\right]& \left[\begin{array}{l}r\hfill \\ 0\hfill \end{array}\right]& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ \left[\begin{array}{c}r\\ n-1\end{array}\right]& \left[\begin{array}{c}r\\ n-2\end{array}\right]& \left[\begin{array}{c}r\\ n-3\end{array}\right]& \cdots & \left[\begin{array}{l}r\hfill \\ 0\hfill \end{array}\right]\end{array}\right)}_{n\times n}$ (1.8)

并且 $\left[\begin{array}{l}r\hfill \\ i\hfill \end{array}\right]=\frac{r\left(r+1\right)\cdots \left(r+i-1\right)}{i!}=\left(\begin{array}{c}r+i-1\\ i\end{array}\right)=\frac{(r+i-1)!}{i!\left(r-1\right)!},\left[\begin{array}{l}0\hfill \\ i\hfill \end{array}\right]=0,\left[\begin{array}{l}0\hfill \\ 0\hfill \end{array}\right]=\left(\begin{array}{l}0\hfill \\ 0\hfill \end{array}\right)=1$。

定义2：设第r阶逆累加生成(累减)定义为 $x^{\left(r-1\right)}\left(k\right)=x^{\left(r\right)}\left(k\right)-x^{\left(r\right)}\left(k-1\right),k=2,3,\cdots ,n$ 且 $x^{\left(r\right)}\left(1\right)=x^{\left(0\right)}\left(1\right),r>0$。用矩阵 $D^r$ 表示满足 $X^{\left(0\right)}=D^r{X}^{\left(r\right)}$ 的r-IAGO矩阵，其中

$D^r={\left(\begin{array}{ccccc}\left[\begin{array}{c}-r\\ 0\end{array}\right]& 0& 0& \cdots & 0\\ \left[\begin{array}{c}-r\\ 1\end{array}\right]& \left[\begin{array}{c}-r\\ 0\end{array}\right]& 0& \cdots & 0\\ \left[\begin{array}{c}-r\\ 2\end{array}\right]& \left[\begin{array}{c}-r\\ 1\end{array}\right]& \left[\begin{array}{c}-r\\ 0\end{array}\right]& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ \left[\begin{array}{c}-r\\ n-1\end{array}\right]& \left[\begin{array}{c}-r\\ n-2\end{array}\right]& \left[\begin{array}{c}-r\\ n-3\end{array}\right]& \cdots & \left[\begin{array}{c}-r\\ 0\end{array}\right]\end{array}\right)}_{n\times n}$ (1.9)

并且 $\left[\begin{array}{c}-r\\ i\end{array}\right]=\frac{-r\left(-r+1\right)\cdots \left(-r+i-1\right)}{i!}={\left(-1\right)}^{r\left(r-1\right)\cdots \left(r-i+1\right)}i!,\left[\begin{array}{c}-r\\ i\end{array}\right]=0,i>r$。

由公式(1.8)和(1.9)可以明显看出一阶累加是分数阶累加的一个特殊例子，即当 $r=1$ 时为一阶累加。接下来介绍一下分数阶累加和累减算子的一个基本定理。

定理1：矩阵 $A^r$ 和矩阵 $D^r$ 满足 $A^r{D}^r=I_n$

证明1：由定义1可知 $det\left(A^r\right)=1$，这意味着 $A^r$ 为可逆矩阵。从 $X^{\left(r\right)}=A^r{X}^{\left(0\right)}$ 可知 ${\left(A^r\right)}^{-1}{X}^{\left(r\right)}={\left(A^r\right)}^{-1}{A}^r{X}^{\left(0\right)}=X^{\left(0\right)}$。根据定义2，有 $X^{\left(0\right)}=D^r{X}^{\left(r\right)}$。因此，显然有 $A^r{D}^r=I_n$，证毕。

由前所述， $GM\left(1,1,t^{\alpha }\right)$ 模型与 $FAGM\left(1,1,t^{\alpha }\right)$ 模型的建模过程只是累加方式发生了些许改变，因此 $FAGM\left(1,1,t^{\alpha }\right)$ 模型的白化微分方程为

$\frac{\text{d}{x}^{\left(r\right)}\left(t\right)}{\text{d}t}+a{x}^{\left(r\right)}\left(t\right)=b{t}^{\alpha }+c,r>0,\alpha >0$ (1.10)

将等式(1.10)在区间 $\left[k-1,k\right]$ 进行积分可得

${\displaystyle {\int }_{k-1}^k\text{d}}{x}^{\left(r\right)}\left(t\right)+a{\displaystyle {\int }_{k-1}^k{x}^{\left(r\right)}}\left(t\right)\text{d}t=b{\displaystyle {\int }_{k-1}^k{t}^{\alpha }}\text{d}t+c{\displaystyle {\int }_{k-1}^k\text{d}}t$ (1.11)

根据梯形公式以及 $z^{\left(r\right)}\left(k\right)=0.5x^{\left(r\right)}\left(k-1\right)+0.5x^{\left(r\right)}\left(k-1\right)+0.5x^{\left(r\right)},k=2,3,\cdots ,n$，上述公式(1.11)可进一步简化为

$x^{\left(r\right)}\left(k\right)-x^{\left(r\right)}\left(k-1\right)+\alpha z^{\left(r\right)}\left(k\right)=b\frac{k^{1+\alpha }-{\left(k-1\right)}^{1+\alpha }}{1+\alpha }+c$ (1.12)

同样地，根据最小二乘法可得 $FAGM\left(1,1,t^{\alpha }\right)$ 模型的参数向量 $\beta ={\left(a,b,c\right)}^{\text{T}}$ 的估计式为

${\left(a,b,c\right)}^{\text{T}}={\left(B{}^{\text{T}}B\right)}^{-1}\left(B^{\text{T}}Y\right)$ (1.13)

其中

(1.14)

定理2：假设观测序列、累加和累减方式满足前面的定义1和定义2，则 $FAGM\left(1,1,t^{\alpha }\right)$ 模型的时间响应序列为

${\hat{x}}^{\left(r\right)}\left(k\right)=\left(x^{\left(0\right)}\left(1\right)-\frac{c}{a}\right){\text{e}}^{-a\left(k-1\right)}+\frac{c}{a}+\frac{b}{2}{\text{e}}^{-a\left(k-1\right)}{\displaystyle \underset{i=1}{\overset{k-1}{\sum }}\left(i^{\alpha }{\text{e}}^{a\left(i-1\right)}+{\left(i+1\right)}^{\alpha }{\text{e}}^{ai}\right)},k=2,3,\cdots ,n$ (1.15)

且 $FAGM\left(1,1,t^{\alpha }\right)$ 模型的预测值 ${\hat{x}}^{\left(0\right)}\left(k\right)$ 为

${\hat{x}}^{\left(0\right)}=D^r{\hat{x}}^{\left(r\right)}$ (1.16)

证明2：已知公式(1.10)的解为齐次方程的通解及其特解的和。因此，由常数变异法可得

$\frac{\text{d}{x}^{\left(r\right)}\left(t\right)}{\text{d}t}+a{x}^{\left(r\right)}\left(t\right)=0$ (1.17)

令 $x^{\left(r\right)}\left(t\right)=U{\text{e}}^{-at}$，其中U为常数，并将 $x^{\left(r\right)}\left(t\right)=U{\text{e}}^{-at}$ 代入公式(1.10)得到

$\frac{\text{d}U\left(t\right)}{\text{d}t}={\text{e}}^{at}\left(b{t}^{\alpha }+c\right)$ (1.18)

对等式(1.18)在区间 $\left[1,t\right]$ 上进行积分可以得到

${\displaystyle {\int }_1^t\text{d}U\left(s\right)}=b{\displaystyle {\int }_1^t{s}^{\alpha }}{\text{e}}^{as}\text{d}s+c{\displaystyle {\int }_1^t{\text{e}}^{as}}\text{d}s$ (1.19)

进一步可以有

$U\left(t\right)=U\left(1\right)+b{\displaystyle {\int }_1^t{s}^{\alpha }}{\text{e}}^{as}\text{d}s+\frac{c}{a}\left({\text{e}}^{at}-{\text{e}}^a\right)$

又因为 $U\left(1\right)=x^{\left(r\right)}\left(1\right){\text{e}}^a=x^{\left(0\right)}\left(1\right){\text{e}}^a$，根据数值积分公式中的两点梯形近似公式可以得到

$\begin{array}{c}{x}^{\left(r\right)}\left(t\right)={\text{e}}^{-at}\left\{x^{\left(0\right)}\left(1\right){\text{e}}^a+b{\displaystyle {\int }_1^t{s}^{\alpha }}{\text{e}}^{as}\text{d}s+\frac{c}{a}\left({\text{e}}^{at}-{\text{e}}^a\right)\right\}\\ =\left(x^{\left(0\right)}\left(1\right)-\frac{c}{a}\right){\text{e}}^{-a\left(t-1\right)}+\frac{c}{a}+b{\text{e}}^{-a\left(t-1\right)}{\displaystyle {\int }_1^t{s}^{\alpha }}{\text{e}}^{a\left(s-1\right)}\text{d}s\\ =\left(x^{\left(0\right)}\left(1\right)-\frac{c}{a}\right){\text{e}}^{-a\left(t-1\right)}+\frac{c}{a}+\frac{b}{2}{\text{e}}^{-a\left(t-1\right)}{\displaystyle \underset{i=1}{\overset{t-1}{\sum }}\left(i^{\alpha }{\text{e}}^{a\left(i-1\right)}+{\left(i+1\right)}^{\alpha }{\text{e}}^{ai}\right)}\end{array}$ (1.20)

所以 $FAGM\left(1,1,t^{\alpha }\right)$ 的时间响应函数为：

${\hat{x}}^{\left(r\right)}\left(k\right)=\left(x^{\left(0\right)}\left(1\right)-\frac{c}{a}\right){\text{e}}^{-a\left(k-1\right)}+\frac{c}{a}+\frac{b}{2}{\text{e}}^{-a\left(k-1\right)}{\displaystyle \underset{i=1}{\overset{k-1}{\sum }}\left(i^{\alpha }{\text{e}}^{a\left(i-1\right)}+{\left(i+1\right)}^{\alpha }{\text{e}}^{ai}\right)}$ (1.21)

3. $FAGM\left(1,1,t^{\alpha }\right)$ 模型的正则化

3.1. 线性回归及L1正则化

在灰色预测模型中，通过求解(2.1)来得到模型的估计参数

${\hat{\beta }}_{linear}={\mathrm{argmin}}_{\beta }\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(B_i\beta -y_i\right)}^2}=\frac{1}{2}{\Vert B\beta -y\Vert }^2$ (2.1)

其中，y为真实值， $B\beta$ 为拟合值。由前所述，为了得到 $GM\left(1,1,t^{\alpha }\right)$ 模型和 $FAGM\left(1,1,t^{\alpha }\right)$ 模型中的参数，必须求解一个超定线性方程组(1.5)或者(1.13)。显然(1.5)和(1.13)是一个标准的线性回归。然而，线性回归可能存在过拟合现象，且不能进行变量选择，从而导致模型也较为复杂。其次，矩阵 $\Delta$ 和B还可能会存在多重共线性，以及当变量间的相关系数为−1或者1时，模型将会存在完全多重共线性的问题。

正则化方法是解决过拟合的一个有效方法。常见的正则化方法有Ridge回归、Elastic Net回归以及Lasso回归等。其中，Ridge回归也称作L2正则化，Elastic Net回归也被称作弹性网络，Lasso回归也称作L1正则化。而本文引入L1正则化法来优化改进 $FAGM\left(1,1,t^{\alpha }\right)$，构建新的 $LFAGM\left(1,1,t^{\alpha }\right)$ 模型。该模型与模型 $FAGM\left(1,1,t^{\alpha }\right)$ 相比，是在原来的损失函数基础上添加了一个参数的L1正则化项，即增加了一个惩罚函数以弥补线性回归所带来的不足。新的损失函数如(2.2)所示。

${\hat{\beta }}_{lasso}=\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(B_i\beta -y_i\right)}^2}+\lambda {\Vert \beta \Vert }_1=\frac{1}{2}{\Vert B\beta -y\Vert }^2+\lambda {\Vert \beta \Vert }_1$ (2.2)

其中 $p\left(\lambda \right)=\lambda {\Vert \beta \Vert }_1\left(\lambda >0\right)$ 为惩罚函数， $\lambda$ 为正则参数。

3.2. $LFAGM\left(1,1,t^{\alpha }\right)$ 模型的线性参数估计

求解最优化问题(2.2)已经有很多有效的方法，其中坐标下降法为一个简单有效的求解方法。该方法为一种非梯度优化的方法，目前已被广泛应用于机器学习中求解大规模数据优化问题。在坐标下降法中，每一步的迭代都是沿一个坐标的方向进行搜索，通过循环使用不同的坐标方法来达到目标函数的极小值。对于一个凸函数 $f\left(x\right),x\in R^n$，该算法的计算步骤总结如下

$\begin{array}{l}{x}_1^{\left(k\right)}=\underset{x_1}{\arg \min }f\left(x_1,x_2^{\left(k-1\right)},x_3^{\left(k-1\right)},\cdots ,x_n^{\left(k-1\right)}\right)\\ x_2^{\left(k\right)}=\underset{x_2}{\arg \min }f\left(x_1^{\left(k\right)},x_2,x_3^{\left(k-1\right)},\cdots ,x_n^{\left(k-1\right)}\right)\\ x_3^{\left(k\right)}=\underset{x_3}{\arg \min }f\left(x_1^{\left(k\right)},x_2^{\left(k\right)},x_3,\cdots ,x_n^{\left(k-1\right)}\right)\\ x_n^{\left(k\right)}=\underset{x_n}{\arg \min }f\left(x_1^{\left(k\right)},x_2^{\left(k\right)},x_3^{\left(k\right)},\cdots ,x_n\right)\end{array}$ (2.3)

该算法每一次迭代都会更新参数x的一个维度。其原理是通过不停的迭代来构造序列 $x^0,x^1,x^2,\cdots$ 来求解问题。最终收敛到期望的极值点。因此，问题(2.2)的坐标下降迭代格式为

${\hat{\beta }}_i^{lasso}=S_{\lambda /{\Vert B_i\Vert }^2}\left(\frac{B_i^{\text{T}}\left(Y-B_{-i}{\beta }_{-i}\right)}{B_i^{\text{T}}{B}_i}\right)$ (2.4)

其中S为软阈值函数。

3.2.1 $LFAGM\left(1,1,t^{\alpha }\right)$ 模型的超参数优化

注意到，上述所有的计算都是在分数阶r，幂指数 $\alpha$ 和正则参数 $\lambda$ 已知的情况下进行计算的。但是在实际计算时，模型 $LFAGM\left(1,1,t^{\alpha }\right)$ 对分数阶r，幂指数 $\alpha$ 和正则参数 $\lambda$ 这三个超参数非常敏感。因此，本文设计了一个优化模型(2.5)来搜索这三个参数的最优值。

$\begin{array}{l}{\min }_{r,\lambda ,\alpha }\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left|\left({\hat{x}}^{\left(r\right)}\left(k\right)-x^{\left(0\right)}\left(k\right)\right)/x^{\left(0\right)}\left(k\right)\right|\times 100\%}\\ \text{s}\text{.t}\text{.}r>0,\lambda >0,\alpha >0\end{array}$ (2.5)

由于模型(2.5)的具有高度非线性，求解复杂，传统的数值优化方法求解困难甚至难以求解。因此本文使用启发式群智能优化算法—灰狼优化算法(GWO)来求解优化模型(2.5)中的分数阶r，幂指数 $\alpha$ 和正则参数 $\lambda$。

灰狼算法 [15]，模拟了灰狼在也野外捕食场景。首先将灰狼分为四个阶级 $\theta$，$\xi$，$\eta$，w。其关系从上到下依次递减如图1。

图1. 灰狼等级制度

并且灰狼捕食过程可分为三个部分；

1) 跟踪，接近猎物；

2) 追逐，包围猎物逼停猎物；

3) 攻击猎物。

3.2.2. GWO算法模型

将灰狼种群中最好的三匹狼(最优解)分别设为 $\theta$，$\xi$ 和 $\eta$，它们可以引导全狼群对目标进行搜索。而其他的狼群定义为w，并且跟随 $\theta$，$\xi$ 和 $\eta$ 来更新位置。

第1步：包围猎物

在狩猎过程中，将灰狼的狩猎过程定义为：

$\stackrel{\to }{D}=\left|\stackrel{\to }{C}\cdot {\stackrel{\to }{X}}_p\left(t\right)-\stackrel{\to }{X}\left(t\right)\right|$ (2.6)

$\stackrel{\to }{X}\left(t+1\right)={\stackrel{\to }{X}}_p\left(t\right)-\stackrel{\to }{A}\cdot \stackrel{\to }{D}$ (2.7)

其中(2.6)表示为狼群个体与猎物之间的距离，(2.7)为狼群更行位置。t为当前迭代数。 $\stackrel{\to }{A}$ 和 $\stackrel{\to }{C}$ 为系数向量， ${\stackrel{\to }{X}}_p$ 和 $\stackrel{\to }{X}$ 为猎物的位置以及灰狼的位置。 $\stackrel{\to }{A}$ 和 $\stackrel{\to }{C}$ 计算公式为：

$\stackrel{\to }{A}=2\stackrel{\to }{\epsilon }\cdot \stackrel{\to }{r_1}-\stackrel{\to }{\epsilon }$ (2.8)

$\stackrel{\to }{C}=2\cdot \stackrel{\to }{r_2}$ (2.9)

其中， $\stackrel{\to }{\epsilon }$ 为收敛因子，并且会随着迭代次数从2减为0， $\stackrel{\to }{r_1}$ 和 $\stackrel{\to }{r_2}$ 为[0, 1]上的随机数。

第2步：狩猎

当灰狼识别猎物位置之后 $\xi$ 和 $\eta$ 会在 $\theta$ 的领导下包围猎物，灰狼跟踪猎物的位置可以表示为：

$\begin{array}{l}\stackrel{\to }{D_{\theta }}=\left|\stackrel{\to }{C_1}\cdot \stackrel{\to }{X_{\theta }}-\stackrel{\to }{X}\right|\\ \stackrel{\to }{D_{\xi }}=\left|\stackrel{\to }{C_2}\cdot \stackrel{\to }{X_{\xi }}-\stackrel{\to }{X}\right|\\ \stackrel{\to }{D_{\eta }}=\left|\stackrel{\to }{C_3}\cdot \stackrel{\to }{X_{\eta }}-\stackrel{\to }{X}\right|\end{array}$ (2.10)

其中， $\stackrel{\to }{D_{\theta }}$，$\stackrel{\to }{D_{\xi }}$ 和 $\stackrel{\to }{D_{\eta }}$ 分别为 $\theta$，$\xi$ 和 $\eta$ 与其他个体间的距离。 $\stackrel{\to }{X_{\theta }}$，$\stackrel{\to }{X_{\xi }}$ 和 $\stackrel{\to }{X_{\delta }}$ 为 $\theta$，$\xi$ 和 $\eta$ 的当前位置，且 $\stackrel{\to }{C_1}$，$\stackrel{\to }{C_2}$ 和 $\stackrel{\to }{C_3}$ 为随机向量， $\stackrel{\to }{X}$ 为灰狼当前位置。

$\begin{array}{l}\stackrel{\to }{X_1}=\stackrel{\to }{X_{\theta }}-A_1\cdot \left(\stackrel{\to }{D_{\theta }}\right)\\ \stackrel{\to }{X_2}=\stackrel{\to }{X_{\xi }}-A_2\cdot \left(\stackrel{\to }{D_{\xi }}\right)\\ \stackrel{\to }{X_3}=\stackrel{\to }{X_{\eta }}-A_3\cdot \left(\stackrel{\to }{D_{\eta }}\right)\end{array}$ (2.11)

$\stackrel{\to }{X}\left(t+1\right)=\frac{\stackrel{\to }{X_1}+\stackrel{\to }{X_2}+\stackrel{\to }{X_3}}{3}$ (2.12)

其中，(2.8)分别设定了狼群中w的朝向，以及θ，ξ和η前进步长和方向。公式(2.12)为w的最终位置。

第3步：攻击猎物

猎物，停止移动之后，灰狼便会开始攻击， $\stackrel{\to }{\epsilon }$ 收敛因子的值会逐渐减小，所以A的波动范围也会随之减小。在迭代的过程中当 $\stackrel{\to }{\epsilon }$ 的值从2下降到0时，而 $\stackrel{\to }{A}$ 的值会在 $\left[-\epsilon ,\epsilon \right]$ 之间变化，并且当 $\left|\stackrel{\to }{A}\right|<1$ 时狼群会发起攻击(局部最优)，当 $\left|\stackrel{\to }{A}\right|>1$ 时，狼群会脱离猎物，寻找更合适的猎物(全局最优)。

其中由公式(2.9)可知， $\stackrel{\to }{C}$ 为[0, 2]之间的随机数，表示为狼群所在位置对猎物的影响的随机权重，当 $C>1$ 时表示权重大， $C<1$ 时表示权重小。这样可以避免算法更加随机探索，同时避免陷入局部最优。

GWO算法流程图如图2所示。

图2. 灰狼优化算法流程图

4. 实例分析

为验证 $LFAGM\left(1,1,t^{\alpha }\right)$ 模型的有效性和实用性，本文分别使用 $FAGM\left(1,1,t^{\alpha }\right)$ 模型和 $LFAGM\left(1,1,t^{\alpha }\right)$ 模型对中国农业耕地灌溉面积情况进行模拟预测，数据来源于《中国统计年鉴》(2009~2020年)。其中，前7组为训练集，后5组为测试集，相关计算结果见表1。对于 $FAGM\left(1,1,t^{\alpha }\right)$ 模型，灰狼优化算法得到的最优非线性参数 $\alpha =6.364$，累加阶数 $r=0.026$，参数向量 $\beta ={\left(0.167,-0.03,1.229\times {10}^4\right)}^{\text{T}}$，对于 $LFAGM\left(1,1,t^{\alpha }\right)$ 模型，灰狼优化算法得到的最优非线性参数 $\alpha =0.728$，累加阶数 $r=0.010$，正则化参数 $\lambda =74.575$，参数向量 $\beta ={\left(-0.028,-186.322,0\right)}^{\text{T}}$。

表1. $FAGM\left(1,1,t^{\alpha }\right)$ 模型以及 $LFAGM\left(1,1,t^{\alpha }\right)$ 模型精度比较

由表1可见，经过L1正则化后的 $LFAGM\left(1,1,t^{\alpha }\right)$ 模型在测试组预测精度上会更高一些，可见正则化可使灰色预测模型具有更高的预测精度，能更好的避免过拟合问题，提高模型的泛化性能。同时，该数值实例也验证了 $LFAGM\left(1,1,t^{\alpha }\right)$ 模型有更好的实用性。

5. 小结

本文将借助Lasso的思想，将L1正则化方法考虑到灰色时间幂预测模型。然而，基于一阶累加的 $GM\left(1,1,t^{\alpha }\right)$ 模型的拟合精度和预测精度都没有分数阶累加高。因此，在 $GM\left(1,1,t^{\alpha }\right)$ 模型的基础上，本文将分数阶累加添加到 $GM\left(1,1,t^{\alpha }\right)$ 模型的建模过程中，并基于L1正则化方法，构建新的灰色 $LFAGM\left(1,1,t^{\alpha }\right)$ 模型。该模型在实际预测中具有更高的预测精度。可见，L1正则化和分数阶累加进一步提高了基于一阶累加的 $GM\left(1,1,t^{\alpha }\right)$ 模型的预测精度。然而，该方法是基于短期数据的预测处理，对于大量规模的数据处理还有待于进一步验证有效性和实用性。

基金项目

国家自然科学基金项目“数据驱动的高可解释迁移模糊系统预测建模新方法及应用研究” (71961001)。

文章引用

董帮强,张延飞,丁木华,陈 萍. 基于GWO优化的L1正则化分数阶灰色时间幂预测模型
L1 Regularized Fractional Grey Time Power Forecasting Model Based on GWO Optimization[J]. 应用数学进展, 2021, 10(10): 3277-3287. https://doi.org/10.12677/AAM.2021.1010344

参考文献