Trajectory Control of Quadrotor With Cable-Suspended Load via Dynamic Feedback Linearization
-
摘要: 三维空间下的四旋翼吊挂运输系统是一种欠驱动、强耦合、多变量的非线性系统. 根据系统的动力学特点, 将系统分解为双质点系绳连接子系统和四旋翼姿态控制子系统. 选择与系统自由度维数相同的广义坐标并基于虚位移原理计算对应的广义力, 从而建立系统的拉格朗日动力学方程. 利用微分平滑特性证明了运输系统存在平凡零动态, 因此可通过动态反馈转化为线性和能控系统. 经过2次动态扩展和变量代换, 原系统扩展为总相对阶等于系统状态维度的线性能控系统. 基于赫尔维茨稳定性判据, 设计了跟踪误差指数收敛的动态反馈控制律. 该方法可作为一类非线性系统控制器设计的标准方法. 最后以三维空间的螺旋曲线及水平面内频率变化的圆周曲线为参考轨迹进行仿真, 仿真结果验证了控制系统的有效性.Abstract: A quadrotor with cable-suspended load in 3-D space is considered, which is underactuated, strongly coupling and nonlinear. According to the dynamic feature, the system is decoupled into quadrotor attitude control subsystem and double points link subsystem. The generalized coordinate with the same dimension of the system freedom is selected and the generalized force is calculated based on the principle of virtual displacement, then the Lagrange dynamic equation of the system is established. The quadrotor-load system is proved to have ordinary zero dynamics based on differentially flat property, so it can be transformed into linear and controllable system by dynamic feedback. After 2 dynamic expansions and variable substitutions, the original system is extended to a linear controllable system whose total relative orders are equal to the system state dimensions. Based on Hurwitz stability criterion, a dynamic feedback controller with exponential convergence of position error is designed. This method can be used as a standard method for a class of nonlinear systems. Finally, the spiral curve in 3-D space and the circular curve in the horizontal plane with varying frequency are taken as the reference trajectory for simulation. Simulation results demonstrate the effectiveness of the proposed method.
-
Key words:
- Quadrotor /
- transport /
- zero dynamics /
- dynamic feedback
-
电熔镁砂是一种生产航天、航空和工业所需耐火材料的原料. 电熔镁砂生产过程采用多台电熔镁炉, 熔炼菱镁矿石, 生产电熔镁砂. 电熔镁炉是一种重大耗能设备, 耗电成本占总生产成本的60%以上[1-2]. 需量是电熔镁砂生产过程的耗电指标, 供电部门规定需量为当前时刻和过去29个采样时刻的电熔镁砂生产用电功率的均值, 其采样周期为7秒[3]. 电熔镁砂生产企业采用需量监控系统监控需量值, 当需量达到其限幅值(22 300 kW)时会切断其中一台电熔镁炉的供电, 使需量不超过限幅值, 然而切断供电会破坏炉内温度场吸热和放热平衡, 降低产品质量[4].
在电熔镁砂生产过程中, 当原料杂质成分含量增大和颗粒大小变化等造成阻抗减小, 导致熔化电流增大、功率增大, 进而需量上升. 由于熔化电流设定值不变, 因此熔化电流控制系统使熔化电流下降, 需量下降, 出现需量先升后降的尖峰现象. 当需量尖峰值达到其限幅值会造成拉闸断电, 为了避免需量尖峰造成的不必要的断电, 需要对需量尖峰进行识别, 因此需要对需量进行多步预报.
需量多步预报的研究集中在城市用电领域. 文献[5]利用采样周期为1小时的用电需量数据, 提出一种自回归模型与三次样条曲线相结合的混合负荷多步预报方法, 预报大西洋城未来一天每小时的用电需量. 文献[6]利用采样周期为0.5小时的用电需量数据, 提出基于信号滤波和季节调整的多输出前馈神经网络与经验模态分解相结合的多步预报模型, 预报澳大利亚新南威尔士州未来一天每半小时的用电需量. 文献[7]利用采样周期为1分钟的需量数据, 提出基于非线性自回归神经网络的多步预报方法, 预报未来1分钟、10分钟、1小时和2小时的住宅用电需量. 文献[8]利用采样周期为0.5小时的用电需量数据, 提出基于多目标优化算法的最小二乘支持向量机的多步预报方法, 预报电力市场未来0.5小时、1小时、2小时和3小时的用电需量. 文献[9]利用采样周期为1天的用电需量数据, 提出基于支持向量机的多步预报模型, 预报某配电公司未来30天的每天最大用电需量. 文献[10]利用采样周期为0.5小时的用电需量数据, 提出基于支持向量机、极限学习机和多层循环神经网络的用电需量多步预报方法, 预报了未来1.5小时的每半小时的城市用电需量. 文献[11]利用采样周期为1小时的用电需量数据, 提出基于自回归滑动平均模型的某地区未来一天的用电需量预报方法. 文献[12]利用采样周期为1个月的用电需量数据, 提出基于多层神经网络的未来12个月的地区用电需量预报模型. 文献[13]利用采样周期为1个月的用电需量数据, 提出基于神经网络的未来12个月的地区用电需量预报方法. 文献[14]利用采样周期为15分钟的用电需量数据, 提出基于变分模态分解的某地区未来1分钟、3分钟和5分钟的用电需量预报方法. 文献[15]利用采样周期为30秒的用电需量数据, 提出基于集成经验模态分解和卷积神经网络的某地区未来30秒、1分钟、1.5分钟和2分钟的用电需量预报方法. 电熔镁砂生产用电需量的变化是快变化的动态系统, 需量的采样周期为7秒, 发生需量尖峰的整个时间在70秒之内. 由于电熔镁砂生产过程是由多台电熔镁炉运行组成, 其用电需量变化过程是模型结构与系统阶次未知的非线性动态系统, 而城市用电需量变化是一个慢变化的系统, 可以采用上述静态建模方法[5-15]. 电熔镁砂生产用电需量的多步预报难以采用文献[5-15]方法进行需量尖峰的准确预报.
文献[16-17]利用采样周期为7秒的用电需量数据, 提出数据与模型驱动的电熔镁砂生产用电需量单步预报方法. 为了提高需量预报精度, 文献[18]采用用电需量大数据提出系统辨识与自适应深度学习相结合的电熔镁砂生产用电需量一步预报方法. 由于需量尖峰是需量先升高后下降, 需要对需量进行多步预报.
本文利用电熔镁砂生产过程熔化电流闭环控制系统方程, 建立了由线性模型和未知非线性动态系统组成的群炉需量多步预报模型, 将系统辨识与自适应深度学习相结合, 采用端边云协同结构, 提出了电熔镁砂生产用电需量多步智能预报方法. 采用电熔镁砂生产过程的工业大数据验证了所提方法可以准确预报用电需量的变化趋势.
1. 电熔镁砂需量监控过程描述
如图1所示, 电熔镁砂的生产过程由多台电熔镁炉控制系统, 供电系统和需量监控系统组成, 每台电熔镁炉控制系统由电熔镁炉、拉闸系统、加料系统、熔化电流控制系统组成[2]. 需量监控系统由功率采集装置和需量监控计算机组成. 电熔镁炉采用埋弧方式, 通过加料系统将原矿送入电熔镁炉内, 电流控制系统控制三相电极与原料之间的距离使产生的电弧电流达到生产工艺规定的熔化电流, 形成熔池. 边熔化边加料, 使熔池增高至炉口, 熔炼过程结束. 熔炼过程持续约10小时, 每台电熔镁炉每炉次耗电约40 MWh, 因此电熔镁炉是高耗能设备.
为了节能减排, 电力部门规定当前时刻和过去29个时刻的功率的平均值作为当前时刻的需量值并规定了需量的限幅值, 当需量实际值超过需量限幅值对用电企业罚款, 因此生产企业设立需量监控系统, 操作人员监视需量实际值, 当超过限幅值, 通过拉闸系统切断其中一台电熔镁炉的电源. 由于电熔镁砂生产过程常出现尖峰现象, 导致不必要的拉闸. 为避免尖峰导致的不必要的拉闸断电, 需要对尖峰进行预报, 由于尖峰现象是需量先升高后下降, 因此需要对需量进行多步预报.
2. 电熔镁砂生产用电需量多步预报模型
需量$ \bar{p}(k) $为$ k $时刻群炉功率$ p(k) $与过去29个采样周期的群炉功率的均值, 即
$$ \begin{split} \bar{p}(k)=\;& \frac{1}{30}\sum_{j = 0}^{29} p(k-j) =\\ &\bar{p}(k-1)+\frac{p(k)-p(k-30)}{30} \end{split} $$ (1) 其中, $ k = 1 $表示采样周期7秒, $ (k+n) $时刻的需量$\bar{p}(k+n)$为
$$ \begin{split} \bar{p}(k+n) = &\frac{1}{30}\sum_{j = 0}^{29} p(k+n-j)=\frac{1}{30}\sum_{i = 0}^n{p}(k+n-i)\;+\\&\frac{1}{30}\sum_{j = 0}^{29-n} p(k-j)=\bar{p}(k)\;-\\ &\frac{1}{30}\sum_{i = 0}^n{p}(k-30+i)+\frac{1}{30}\sum_{i = 1}^n p(k+i)\ \end{split} $$ (2) 其中, 群炉功率$ p(k+1), \cdots, p(k+n) $未知, 因此需要建立群炉功率多步预报模型.
2.1 群炉功率多步预报模型
首先需建立群炉功率$ p(k) $的动态模型. 电熔镁炉群炉功率$ p(k) $为 $M$个电熔镁炉功率之和$ \sum_{s = 1}^M p_s(k) $.
采用文献[18]的方法建立第$ s $个电熔镁炉的功率$ p_s(k) $的动态模型. 第$ s $个电熔镁炉的熔化电流$ y_s(k) $为电流闭环控制系统的输出, 其动态模型为
$$ \begin{align} T_s(z^{-1})y_s(k) = G_s(z^{-1})y^*+v_s(k) \end{align} $$ (3) 其中, $ y^* $是已知的熔化电流设定值, $ v_s(k) $为未知非线性项, $ T_s(z^{-1}) $ 和 $ G_s(z^{-1}) $为关于$ z^{-1} $的多项式, 即$ T_s(z^{-1}) = 1+t_{s1}(z^{-1})+t_{s2}(z^{-2})+t_{s3}(z^{-3}) $, $G_s(z^{-1}) = g_{s0}+g_{s1}(z^{-1})+g_{s2}(z^{-2})$.
第$ s $个电熔镁炉的功率动态模型为
$$ \begin{split} p_s(k) =\; &\sqrt{3}Uy_s(k)\cos \phi =-\;t_{s1}p_s(k-1)\;-\\ &t_{s2}p_s(k-2)-t_{s3}p_s(k-3)\;+\\ &d_{s0}p^*+v_s(k) \end{split} $$ (4) 其中, $ p^* $是熔化电流设定值 $ y^* $对应的功率; $p_s(k- 1) $, $p_s(k-2) $, $ p_s(k-3)$为$ (k-1) $, $ (k-2) $, $(k- 3) $时刻的第$s$个电熔镁炉功率; $d_{s0} =g_{s0}+ g_{s1}+ g_{s2} $.
采用未知常数$ t_{1}, t_{2}, t_{3}, d_{0} $代替式(4)中的$ t_{s1}, t_{s2}, t_{s3}, d_{s0} $, 群炉功率$ p(k) $为
$$ \begin{split} p(k) =\; &\sum_{s = 1}^{M} p_s(k)= -\;t_{1}p(k-1)-t_{2}p(k-2)\;-\\ &t_{3}p(k-3)+d_{0}p^*+v(k) \end{split} $$ (5) 其中, $ v(k) $为$v(k) = \sum_{s = 1}^{M}\left((t_1-t_{s1})p(k-1)+\right. (t_2-$ $t_{s2})p(k-2) $ $+\;(t_3\;-\;t_{s3})p(k-3)\;-\;(d_0\;-\;d_{s0})p^*\;+$ $v_s(k) ). $
由式(5)可得$ (k+1) $时刻的群炉功率$ p(k+1) $为
$$ \begin{split} {p}(k+1) =\; &{\boldsymbol{\varphi}}(k){{\boldsymbol{\theta}}}^1+{v}(k+1)=\\ &{\boldsymbol{\varphi}}(k){{\boldsymbol{\theta}}}^1+{r}(k+1) \end{split} $$ (6) 其中, ${\boldsymbol{\varphi}}(k) \,= \,( p(k),\, p(k-1),\, p(k-2), p(k-3) )$, $ {{\boldsymbol{\theta}}}^1 = ( {\theta}_1^1, {\theta}_2^1, {\theta}_3^1, {\theta}_4^1)^{\rm{T}} = (-t_{1}, -t_{2}, -t_{3}, d_{0})^{\rm{T}} $, ${r}(k\;+ 1) = {v}(k+1).$ $ (k+2) $时刻的群炉功率 $ p(k+2) $为
$$ \begin{split} {p}(k+2) =\; &{\boldsymbol{\varphi}}(k+1){\boldsymbol{\theta}}^1+v(k+2)=\\ &({\boldsymbol{\varphi}}(k){{\boldsymbol{\theta}}}^1+{v}(k+1), p(k), p(k-1),\\ & p^* ){\boldsymbol{\theta}}^1+v(k+2)=\\ &p(k)((\theta_1^1)^2+\theta_2^1)+p(k-1)(\theta_1^1\theta_2^1+\theta_3^1)\;+\\ &p(k-2)\theta_3^1\theta_1^1+p^*(\theta_1^1\theta_4^1+\theta_4^1)\;+\\ &r(k+1)\theta_1^1+v(k+2)=\\ &{\boldsymbol{\varphi}}(k){{\boldsymbol{\theta}}}^2(k)+r(k+2) \\[-10pt]\end{split} $$ (7) 其中, 线性模型参数${{\boldsymbol{\theta}}}^2 = ({\theta}_1^2, {\theta}^2_{2}, {\theta}^2_{3}, {\theta}^2_{4})^{\rm{T}},{\theta}_1^2 = (\theta_1^1)^2\, +$ $\theta_2^1, {\theta}^2_{2} = \theta_1^1\theta_2^1+\theta_3^1, {\theta}^2_{3} = \theta_3^1\theta_1^1, {\theta}^2_{4} = \theta_1^1\theta_4^1+\theta_4^1$, 未知非线性项$ r(k+2) $ 为 ${r}(k+2) \;=\; {r}(k+1)\theta_1^1(k+1)\;+ v(k+2)$.
采用数学归纳法可以证明假设$ (k+n-1) $时刻的群炉功率$ p(k+n-1) $为式(8)成立, 则$ (k+n) $时刻的群炉功率$ p(k+n) $表达式(9)成立
$$ {p}(k+n-1) = {\boldsymbol{\varphi}}(k){{\boldsymbol{\theta}}}^{n-1}+ {r}(k+n-1) $$ (8) $$ {p}(k+n) = {\boldsymbol{\varphi}}(k){{\boldsymbol{\theta}}}^n+ {r}(k+n) $$ (9) 其中, ${{\boldsymbol{\theta}}}^{n-1} = ({\theta}_1^{n-1}, {\theta}^{n-1}_{2}, {\theta}^{n-1}_{3}, {\theta}^{n-1}_{4})^{\rm{T}}$, ${{\boldsymbol{\theta}}}^n = ({\theta}_1^n,$${\theta}_2^n, {\theta}^n_{3}, {\theta}^n_{4})^{\rm{T}}, {\theta}_1^n = {\theta}_1^{n-1}\theta_1^1+{\theta}_1^{n-2}\theta_2^1+{\theta}_1^{n-3}\theta_3^1$, ${\theta}_2^n = $ ${\theta}_2^{n-1}\theta_1^1+{\theta}_2^{n-2}\theta_2^1+{\theta}_2^{n-3}\theta_3^1$, $ {\theta}_3^n = {\theta}_3^{n-1}\theta_1^1+{\theta}_3^{n-2}\theta_2^1\;+ $ ${\theta}_3^{n-3}\theta_3^1 $, $ {\theta}_4^n = {\theta}_4^{n-1}\theta_1^1+{\theta}_4^{n-2}\theta_4^1+{\theta}_4^{n-3}\theta_3^1+\theta_4^1$, 未知非线性项$ {r}(k+n) $为${r}(k+n) = {r}(k+n- 1)\theta_1^1+{r}(k\;+ n-2)\theta_2^1+{r}(k+n-3)\theta_3^1+{v}(k+n)$. 因此群炉功率多步预报模型为
$$ \begin{align} {P}(k+n) = &\Theta^{\rm{T}}\varphi^{\rm{T}}(k)+ {R}(k+n) \end{align} $$ (10) 其中, $ {P}(k+n) = (p(k+1), \cdots, p(k+n))^{\rm{T}} $, 模型参数$ \Theta = ({\boldsymbol{\theta}}^1, \cdots, {\boldsymbol{\theta}}^n) $, $R(k+n) = (r_1(k+1), \cdots, r_i(k\,+ i), \cdots, r_n(k+n))^{\rm{T}}$.
由式(10)可知模型参数$ \Theta $的辨识方程为
$$ \begin{align} P_s(k) = \begin{pmatrix} {\boldsymbol{\varphi}}(k-1){\boldsymbol{\theta}}^1\\ \vdots\\ {\boldsymbol{\varphi}}(k-i){\boldsymbol{\theta}}^i\\ \vdots\\ {\boldsymbol{\varphi}}(k-n){\boldsymbol{\theta}}^n \end{pmatrix} + {R_s}(k) \end{align} $$ (11) 其中, $ n\times 1 $的向量$ P_s(k+n) $为
$$ P_s(k+n) = (p(k), \cdots, p(k))^{\rm{T}} $$ $$ {R_s}(k) = (r_1(k), \cdots, r_i(k), \cdots, r_n(k))^{\rm{T}} $$ 采用最小二乘算法[19]离线辨识模型参数$ \Theta $, 由式(10)和式(11)可得群炉功率预报模型为
$$ \begin{align} {P}(k+n) = \hat{\Theta}^{\rm{T}}{\boldsymbol{\varphi}}^{\rm{T}}(k)+ \bar{R}(k+n) \end{align} $$ (12) 其中, $\hat{\Theta}^{\rm{T}} = ({\hat{\boldsymbol{\theta}}}^1, \cdots, {\hat{\boldsymbol{\theta}}}^n), \bar{R}(k+n) = (\Theta-\hat{\Theta}^{\rm{T}})\varphi^{\rm{T}}(k)\,+$ $R(k+n) = (\bar{r}_1(k+1), \cdots, \bar{r}_i(k+i), \cdots, \bar{r}_n(k + n))^{\rm{T}}$,
$$ \begin{split} \bar{r}_i(k+i) =\;& f_i(p(k), p(k-1), \cdots, \\ &\bar{r}_i(k), \bar{r}_i(k-1), \cdots) \end{split} $$ (13) 其中, $ f(\cdot) $是模型结构与阶次未知的未知非线性函数, $i = 1, \cdots, n$. 由式(11)知
$$ \begin{align} {\bar{R}_s}(k) = P_s(k)- \begin{pmatrix} {\boldsymbol{\varphi}}(k-1){\hat{\boldsymbol{\theta}}}^1\\ \vdots\\ {\boldsymbol{\varphi}}(k-i){\hat{\boldsymbol{\theta}}}^i\\ \vdots\\ {\boldsymbol{\varphi}}(k-n){\hat{\boldsymbol{\theta}}}^n \end{pmatrix} \end{align} $$ (14) 其中, $ \bar{R}_s(k) = \left(\bar{r}_{1}(k), \cdots, \bar{r}_{i}(k), \cdots, \bar{r}_{n}(k)\right) $.
2.2 需量多步预报模型
由式(2)和式(12)可以得到需量多步预报模型为
$$ \begin{split} \bar{P}\left( k+n \right) =\; &\left( \begin{matrix} \bar{p}(k) \\ \bar{p}(k) \\ \vdots \\ \bar{p}(k) \\ \end{matrix} \right)-\frac{1}{30}\left( \begin{matrix} p(k-29) \\ \displaystyle\sum\limits_{i = 1}^{2}{p(k-30+i)} \\ \vdots \\ \displaystyle\sum\limits_{i = 1}^{n}{p(k-30+i)} \\ \end{matrix} \right)+\\ & \frac{1}{30}\left( \begin{matrix} p(k\text+1) \\ \displaystyle\sum\limits_{i = 1}^{2}{p(k+i)} \\ \vdots \\ \displaystyle\sum\limits_{i = 1}^{n}{p(k+i)} \end{matrix} \right) =\\ &\left( \begin{matrix} \bar{p}(k) \\ \bar{p}(k) \\ \vdots \\ \bar{p}(k) \end{matrix} \right)-\frac{1}{30}\left( \begin{matrix} p(k-29) \\ \displaystyle\sum\limits_{i = 1}^{2}{p(k-30+i)} \\ \vdots \\ \displaystyle\sum\limits_{i = 1}^{n}{p(k-30+i)} \end{matrix} \right)+\\ &\frac{1}{30}\left( \begin{matrix} {\boldsymbol{\varphi}} (k){{{\hat{\theta }}}^{1}} \\ \displaystyle\sum\limits_{i = 1}^{2}{{\boldsymbol{\varphi}} (k){{{\hat{\theta }}}^{i}}} \\ \vdots \\ \displaystyle\sum\limits_{i = 1}^{n}{{\boldsymbol{\varphi}} (k){{{\hat{\theta }}}^{i}}} \\ \end{matrix} \right)+\frac{1}{30}\left( \begin{matrix} \bar{r}(k) \\ \displaystyle\sum\limits_{i = 1}^{2}{\bar{r}(k+i)} \\ \vdots \\ \displaystyle\sum\limits_{i = 1}^{n}{\bar{r}(k+i)} \\ \end{matrix} \right) =\\ &{\bar{P}_S(k)}+\frac{1}{30}AP(k+n-30)\;+\\ &\frac{1}{30}A\hat{\Theta}^{\rm{T}}{\boldsymbol{\varphi}}^{\rm{T}}(k)+ \frac{1}{30}A\bar{R}(k+n) \\[-10pt]\end{split} $$ (5) 其中, $\bar{P}(k+n) = (\bar{p}(k+1), \cdots, \bar{p}(k+n))^{\rm{T}}$, ${\bar{P}_S(k)} = (\bar{p}(k), \cdots, \bar{p}(k))^{\rm{T}}$, $P(k+n-30) = (p(k-29), \cdots, p(k-n+30))^{\rm{T}}$, $ \bar{R}(k+n) $ 为 $(\bar{r}_1(k), \cdots, \bar{r}_i(k+ i), \cdots, \bar{r}_n(k+n))^{\rm{T}}$, $ A $为非0元素为1的$ n\times n $维下三角阵. 式(15)中$ {\bar{P}_S(k)}, P(k+n-30), \hat{\Theta}^{\rm{T}}{\boldsymbol{\varphi}}^{\rm{T}}(k) $已知, 向量$ \bar{R}(k+n) $是模型结构阶次与参数未知的非线性动态系统.
3. 需量多步智能预报方法
3.1 需量多步预报策略
采用文献[18]方法建立式(15)的未知非线性动态系统的由在线深度学习多步预报模型、自校正深度学习多步预报模型和自校正机制组成的自适应深度学习多步预报模型.
采用图2所示的端边云协同结构实现电熔镁砂生产用电需量多步智能预报算法, 其中端−需量监控系统实时采集电熔镁砂生产过程中的群炉功率$ p(k) $与需量$ \bar{p}(k) $, 边−需量预报计算机执行数据处理、在线计算$ \hat{\Theta}^{\rm{T}}\varphi^{\rm{T}}(k) $, 采用窗口长度为$ N $的数据和$ \bar{R}(k\;+ n) $的在线深度学习多步预报模型获得$ \hat{\bar{R}}(k+n) $, 采用式(15)求取需量多步预报值$ \hat{\bar{P}}(k+n) $, 云−数据服务器和人工智能计算平台采用$ k $时刻以及以前所有时刻的输入输出数据和$ \bar{R}(k+n) $的自校正深度学习多步预报模型更新自校正深度学习多步预报模型的全部权值参数和偏置参数. 自校正机制在线监控$ \bar{R}(k+n) $的在线深度学习多步预报模型的预报精度, 当不满足精度要求时, 采用自校正深度学习多步预报模型的参数校正在线深度学习多步预报模型的参数, 从而保证需量的预报精度.
3.2 $ \bar{R}(k+n) $的自适应深度学习多步预报算法
首先建立$ \bar{R}(k+n) $的自适应深度学习多步预报算法. 采用长短周期记忆(Long short-term memory, LSTM)[20-22]的网络架构提出如图3所示的需量第$ i $步预报$ \bar{r}(k+i) $的深度学习预报模型结构. 结合$ \bar{r}(k+i) $动态特性, 将其输入变量作为单个神经元的输入, 未知阶次用神经元个数$ t $来表示, 由式(13)和式(14)可知$ \bar{r}(k+i) $的深度学习网络的第$ j $个神经元的输入为 ${\boldsymbol{x}}_i(k-j+1) = (p(k- j\;+ 1), \bar{r}_i(k-j+1))^{\rm{T}}$, 其中 $\bar{r}_i(k-j+1) = p(k- j\;+ 1)-{\boldsymbol{\varphi}}(k- i-j+1){\hat{\boldsymbol{\theta}}}^i$, $ j = 1, \cdots, t $. 采用文献[18]的训练方法, 利用30个炉次(150000组)的群炉功率数据离线训练$ \bar{r}(k+i) $的深度学习模型结构, 确定神经元个数$ t = 25 $、神经元的节点个数$ \bar{h} = 200 $、隐藏层数$ L = 3 $、在线训练数据窗口长度$ N = 2\,000 $. 采用该深度学习预报模型结构, 提出由在线深度学习预报模型、自校正深度学习预报模型和自校正机制组成的$ \bar{r}(k+i) $的自适应深度学习预报算法.
采用窗口长度$ N $的实时数据在线校正$ \bar{r}(k+i) $的在线深度学习预报模型的输出层权值与偏置参数. 校正算法为
$$ \begin{align} \hat{\bar{r}}_i(k+i) = {\hat{\boldsymbol{\beta}}}_i(k){\boldsymbol{h}}_i^3(k)+\hat{b}_i(k) \end{align} $$ (16) 其中,
$$ {\hat{\boldsymbol{\beta}}}_i(k) = {\hat{\boldsymbol{\beta}}}_i(k-1)-\alpha\frac{\partial L_i(k)}{\partial {{\boldsymbol{\beta}}}_i(k-1)} $$ (17) $$ \hat{b}_i(k) = \hat{b}_i(k-1)-\alpha\frac{\partial L_i(k)}{\partial{b}_i(k-1)} $$ (18) $$ L_i(k) =\frac{1}{N}\sum_{m = 0}^{N-1}||\Delta \bar{r}_i(k-m)||_2 $$ (19) 其中, $ \Delta\bar{r}_i(k-m) = \bar{r}_i(k-m)-\hat{\bar{r}}_i(k-m) $.
采用当前时刻$ k $以及以前所有时刻的实时数据校正自校正深度学习预报模型的全部权值和偏置参数. 其输出层的权值与偏置参数采用式(16) ~ (19)校正. 第$ l $层第$ j $个神经元的输出
$$ \begin{split} {\boldsymbol{h}}_i^l(k-j+1) =\;& {\boldsymbol{o}}_i^l(k-j+1) \;\odot\\ &\tanh( C_i^l(k-j+1)) \end{split} $$ (20) 其中, $ \odot $为哈达玛积[22], $ \tanh(\cdot) $为双曲正切函数, $ ( C_i^l(k-j+1)) $和$ {\boldsymbol{o}}_i^l(k-j+1) $表示第$ l $层第$ j $个神经元的状态和输出门的输出, 采用文献[12]方法计算, 其中, 第$ l $层第$i $个神经元的输出门的权值与偏置参数的校正算法为
$$ W_i^l(k-j+1) = W_i^l(k-j)-\alpha\frac{\partial L_i(k)}{\partial W_i^l(k-j)} $$ (21) $$ b_i^l(k-j+1) = b_i^l(k-j)-\alpha\frac{\partial L_i(k)}{\partial b_i^l(k-j)} $$ (22) $$L_i(k) =\frac{1}{k}\sum_{m = 0}^{k-1}||\Delta \bar{r}_i(k-m)||_2 $$ (23) 其中, $ \Delta\bar{r}_i(k-m) = \bar{r}_i(k-m)-\tilde{r}_i(k-m) $.
为了准确预报需量尖峰, 需要保证需量预报误差精度和需量变化趋势预报精度, 由式(15)知需量预报精度取决于$ \bar{r}(k+i) $的预报精度. 采用自校正机制监控$ \bar{r}(k+i) $的在线深度学习预报模型的预报误差和变化趋势的预报精度, 当不满足精度要求时, 采用自校正深度学习预报模型的各层权值参数与偏置参数校正在线深度学习预报模型的权值参数与偏置参数.
自校正机制采用$ \bar{r}(k+i) $的预报误差$ \Delta \bar{r}_i(k) $, 未知非线性系统第$ i $步上升趋势预报准确率$ { TPR_i(k)} $和第$ i $步下降趋势预报准确率$ { TNR_i(k)} $三项指标, 即
$$ |\Delta\bar{r}_i(k)| = |\bar{r}_i(k)-\hat{\bar{r}}_i(k)| $$ (24) $$ { TPR_i(k)} = \frac{\sum\limits_{m = 0}\limits^{N-1}{ TP}_i(k-m)}{\sum\limits_{m = 0}\limits^{N-1}{ TP}_i(k-m)+\sum\limits_{m = 0}\limits^{N-1}{ FP}_i(k-m)} $$ (25) $$ { TNR_i(k)} =\frac{\sum\limits_{m = 0}\limits^{N-1}{ TN}_i(k-m)}{\sum\limits_{m = 0}\limits^{N-1}{ TN}_i(k-m)+\sum\limits_{m = 0}\limits^{N-1}{ FN}_i(k-m)} $$ (26) 其中, $ {TP}_i(k), { FP}_i(k), { TN}_i(k), { FN}_i(k) $的计算方式如表1所示.
表 1 $ { TP}_i(k), { FP}_i(k), { TN}_i(k), {FN}_i(k)$的计算方式Table 1 Formula mode of ${ TP}_i(k), { FP}_i(k), $ ${ TN}_i(k), {FN}_i(k) $$\hat{\bar{r} }_i(k)-\hat{\bar{r} }_i(k-1)\geq 0$ $\hat{\bar{r} }_i(k)-\hat{\bar{r} }_i(k-1) < 0$ $\bar{r}_i(k)-\bar{r}_i(k-1)\geq 0$ ${TP}_i(k)=1$ ${FP}_i(k)=1$ $\bar{r}_i(k)-\bar{r}_i(k-1)< 0$ ${FN}_i(k)=1$ ${TN}_i(k)=1$ $ k $时刻在线深度学习预报模型的预报误差$ |\Delta\bar{r}_i(k)|\geq \delta_i $ ($ \delta_i $为预报误差上界), 自校正深度学习预报模型的预报误差$ |\Delta\bar{r}_i(k)|< \delta_i $, 且在线深度学习预报模型的上升趋势和下降趋势预报准确率均小于自校正深度学习预报模型的上升趋势和下降趋势预报准确率, 采用自校正深度学习预报模型的全部权值参数和偏置参数校正在线深度学习预报模型的全部权值参数和偏置参数.
3.3 需量多步智能预报算法
由式(15)可得需量多步预报模型为
$$ \begin{split} \hat{\bar{P}}(k+n) =\; &{\bar{P}_S(k)}-\frac{1}{30}AP(k+n-30)\;+\\ &\frac{1}{30}A\hat{\Theta}^{\rm{T}}\varphi^{\rm{T}}(k)+ \frac{1}{30}A\hat{\bar{R}}(k+n) \end{split} $$ (27) 其中,
$$ \begin{split} \hat{\bar{P}}(k+n) =\;& (\hat{\bar{p}}_1(k+1), \cdots, \hat{\bar{p}}_i(k+i), \cdots,\\ &\hat{\bar{p}}_n(k+n))^{\rm{T}}\\ \hat{\bar{R}}(k+n) =\;&(\hat{\bar{r}}_1(k+1), \cdots, \hat{\bar{r}}_i(k+i), \cdots,\\ &\hat{\bar{r}}_n(k+n))^{\rm{T}} \end{split} $$ 端边云协同的需量多步智能预报算法:
1)端−需量监控系统实时采集电熔镁砂生产过程中的群炉功率$ p(k) $与需量$ \bar{p}(k) $;
2) 边−需量预报计算机执行数据处理和需量在线多步预报模型. 在线计算$ \hat{\Theta}^{\rm{T}}\varphi^{\rm{T}}(k) $, 由式(14)计算$ \hat{\bar{R}}(k+n) $. 采用窗口长度为$ N $的输入输出数据由$ \bar{R}(k+n) $的在线深度学习多步预报模型得其预报值$ \hat{\bar{R}}(k+n) $, 由需量多步预报模型式(27)得需量多步预报值$ \hat{\bar{P}}(k+n) $;
3)云−数据服务器和人工智能计算平台采用$ k $时刻以及以前所有时刻的输入输出数据和$ \bar{R}(k+n) $的自校正深度学习多步预报模型得其预报值$ \hat{\bar{R}}(k+ n) $. 采用自校正机制的三项指标式(24) ~ (26), 当$ \bar{R}(k+n) $的在线深度学习多步预报模型的预报误差超过预报精度上界时, 采用自校正深度学习多步预报模型的权值参数和偏置参数校正在线深度学习多步预报模型的权值参数和偏置参数.
4. 实验结果
采用某电熔镁砂生产企业的实际功率和需量数据进行了本文提出的需量多步智能预报方法的实验, 并与文献[9]提出的基于支持向量机、极限学习机和循环神经网络的用电需量多步预报方法进行了对比实验.
采用30个炉次的150000组功率数据离线确定$ \bar{r}(k+i) $的深度学习预报模型结构, 神经元个数$ t = 25 $、神经元的节点个数$ \bar{h} = 200 $、隐藏层数$ L = 3 $. 校正算法中的式(17)、式(18)、式(21)和式(22)中的学习率为$ \alpha = 0.1 $, 在线训练数据窗口长度$ N = 2\,000 $. 自校正机制中的预报误差上界$ \delta_i = 100 $.
由于发生需量尖峰的整个时间小于70秒, 而需量的采样周期为7秒, 因此选择需量预报步数$ i = 1, \cdots, 10 $. 采用线性模型参数辨识方程式(11)得模型参数$ \hat{\Theta} $为
$$ \begin{split} &{{{\hat{\boldsymbol{\theta}} }}^{1}} = {{\left( {0}{.140, 0}{.075, -0}{.157, 0}{.869} \right)}^{\rm T}}\\ & {{{\hat{\boldsymbol{\theta}} }}^{2}} = {{\left( {0}{.201, 0}{.135, -0}{.042, 0}{.587} \right)}^{\rm T}} \\ & {{{\hat{\boldsymbol{\theta}} }}^{3}} = {{\left( {0}{.204, 0}{.146, 0}{.061, 0}{.431} \right)}^{\rm T}}\\ &{{{\hat{\boldsymbol{\theta}} }}^{4}} = {{\left( {0}{.201, 0}{.109, 0}{.103, 0}{.411} \right)}^{\rm T}} \\ & {{{\hat{\boldsymbol{\theta}} }}^{5}} = {{\left( {0}{.203, 0}{.103, 0}{.068, 0}{.427} \right)}^{\rm T}}\\ &{{{\hat{\boldsymbol{\theta}} }}^{6}} = {{\left( {0}{.207, 0}{.103, 0}{.063, 0}{.409} \right)}^{\rm T}} \\ & {{{\hat{\boldsymbol{\theta}} }}^{7}} = {{\left( {0}{.208 , 0}{.103 , 0}{.066, 0}{.389} \right)}^{\rm T}} \\ &{{{\hat{\boldsymbol{\theta}} }}^{8}} = {{\left( {0}{.204, 0}{.103, 0}{.069, 0}{.374} \right)}^{\rm T}}\\ &{{{\hat{\boldsymbol{\theta}} }}^{9}} = {{\left( {0}{.203, 0}{.097, 0}{.071, 0}{.364} \right)}^{\rm T}}\\ &{{{\hat{\boldsymbol{\theta}} }}^{10}} = {{\left( {0}{.202, 0}{.095, 0}{.067, 0}{.358} \right)}^{\rm T}}\end{split} $$ 采用上述30个炉次的150000组需量数据离线建立文献[9]的用电需量多步预报模型:
$$ \begin{split} \hat{\bar{p}}(k+i) =\; &{{W}_{1}}{{\hat{\bar{p}}}_{svm}}(k+i)+{{W}_{2}}{{\hat{\bar{p}}}_{elm}}(k+i)\;+\\ &{{W}_{3}}{{\hat{\bar{p}}}_{rnn}}(k+i) \end{split} $$ 其中, $ \hat{\bar{p}}_{svm}, \hat{\bar{p}}_{elm}, \hat{\bar{p}}_{rnn} $分别为支持向量机、极限学习机和循环神经网络的输出, 上述模型的阶次$ t = 25 $, 支持向量机中径向基函数的期望为150、方差为0.08, 循环神经网络的每层节点数为$ \bar{h} = 150 $, 隐藏层数为$ L = 2 $, 权参数$W_1 = 0.16, W_2 = 0.34, W_3 = 0.71$.
采用实时采集的70个炉次的350000组的需量与功率数据对本文所提需量多步预报算法与文献[9]多步预报算法进行了预报步数$ i = 1, \cdots, 10 $的对比实验. 采用文献[23]的均方根误差(Root mean square error, RMSE)式(28), 文献[24]的拟合优度$ R^2 $式(29), 文献[17]的平均绝对百分比误差(Mean absolute percentage error, MAPE)式(30), 需量第$ i $步预报的上升趋势预报准确率$ {TPR}_i $式(25)和下降趋势预报准确率$ { TNR}_i $式(26)指标对本文所提需量多步预报算法和文献[9]多步预报算法的实验结果进行评估. 实验结果见表2.
表 2 需量预报精度对比Table 2 Precision comparison of demand forecast预报步数$i$ 1 2 3 4 5 6 7 8 9 10 $R^2_i\;(\%)$ 本文 99.96 99.62 99.59 99.47 99.39 99.31 98.99 98.51 98.03 97.95 文献[9] 90.34 90.05 89.77 89.54 88.73 88.48 88.01 87.76 87.33 86.94 ${{RMSE} }_i$ 本文 9.93 11.06 11.99 13.03 13.89 14.73 16.05 16.83 17.93 18.78 文献[9] 24.92 30.01 34.49 39.99 44.79 50.23 54.93 60.05 65.32 70.64 ${{MAPE} }_i\;(\%)$ 本文 0.04 0.05 0.05 0.06 0.06 0.07 0.07 0.08 0.08 0.08 文献[9] 0.11 0.13 0.15 0.18 0.20 0.22 0.24 0.27 0.29 0.32 $TPR_i\;(\%)$ 本文 94.88 93.21 92.19 91.42 90.17 89.77 88.21 90.05 91.55 89.66 文献[9] 86.12 82.11 80.05 80.11 78.94 79.33 79.11 77.06 80.15 79.02 $TNR_i\;(\%)$ 本文 93.22 94.67 92.19 92.01 94.21 93.18 90.96 89.99 88.12 90.01 文献[9] 81.12 80.04 80.67 83.72 79.99 80.15 77.56 86.77 80.15 76.91 $$ {{RMSE}}_i = \sqrt{\frac{1}{N_t}\sum\limits_{k = 1}^{N_t} \left( \bar{p}(k)-\hat{\bar{p}}_i(k)\right)^2} $$ (28) $$ { R}_i^2 = \frac{\sum\limits_{k = 1}\limits^{N_t}\left( \hat{\bar{p}}_i(k)-{\tilde{p}}_i\right)^2 }{\sum\limits_{k = 1}\limits^{N_t}\left(\bar{p}_i(k)-\tilde{p}_i \right)^2 }\times 100{\text{%}} $$ (29) 其中, $ \tilde{p}_i = \frac{1}{N_t}\sum_{k = 1}^{N_t} \bar{p}(k+i-1) $.
$$ {{MAPE}}_i = \frac{1}{N_t}\sum\limits_{k = 1}^{N_t}\left|\frac{\bar{p}(k)-\hat{\bar{p}}_i(k)}{\bar{p}(k)}\right|\times 100{\text{%}} $$ (30) 为了能清楚地对比实验结果, 采用图4 ~ 图6表示从在线的350000组数据的实验结果中随机抽取100组$ i = 1, 5, 10 $步需量预报结果, 虚线为本文需量预报结果, 实线为文献[9]需量预报结果, 点线为需量真实值. 可以看出本文需量预报方法与文献[9]方法相比, 预报精度明显提高.
从表2可以看出, 本文的方法与文献[9]方法的$ i = 1, 5, 10 $步预报结果相比, $ R^2 $提高$ 11.01{\text{%}} $, $ {{RMSE}}_i $降低51.86, $ {{MAPE}}_i $降低$ 0.24{\text{%}} $, 上升趋势$ { TPR}_i $提高$ 10.23{\text{%}} $, 下降趋势$ { TNR}_i $提高$ 14.22{\text{%}} $. 对比实验结果表明本文方法比文献[9]方法的需量多步预报精度和需量变化趋势的预报精度明显提高.
5. 结论
本文通过电熔镁砂生产过程熔化电流闭环控制系统方程建立需量动态模型, 在此基础上建立了由线性模型和未知非线性动态系统组成的需量多步预报模型, 采用文献[18]方法建立了未知非线性动态系统的自适应深度学习多步预报模型, 在此基础上提出了端边云协同的电熔镁砂生产过程需量多步智能预报方法. 采用70个炉次的电熔镁砂生产过程的实际数据的实验结果表明本文的方法与文献[9]方法的$ i = 1, 5, 10 $ 步预报结果相比, $ R^2 $提高$ 11.01{\text{%}} $, $ {{RMSE}}_i $降低51.86, $ {{MAPE}}_i $降低$ 0.24{\text{%}} $, 上升趋势$ { TPR}_i $提高$ 10.23{\text{%}} $, 下降趋势$ { TNR}_i $提高$ 14.22{\text{%}} $, 验证了本文所提的预报方法可以准确预报需量的变化趋势, 为实现需量尖峰的准确预报和控制创造了条件.
-
表 1 仿真中使用的模型参数
Table 1 Model parameters in the simulations
变量 参数 单位 mq 0.4 kg ml 0.1 kg l1 0.8 m l2 0.2 m g −9.8 m·s−2 -
[1] 梁潇, 方勇纯, 孙宁. 平面四旋翼无人飞行器运送系统的轨迹规划与跟踪控制器设计. 控制理论与应用, 2015, 32(11): 1430−1438Liang Xiao, Fang Yong-Chun, Sun Ning. Trajectory planning and tracking controller design for a planar quadrotor unmanned aerial vehicle transportation system. Control Theory and Applications, 2015, 32(11): 1430−1438 [2] Alothman Y, Gu D B. Quadrotor transporting cablesuspended load using iterative Linear Quadratic regulator (iLQR) optimal control. In: Proceedings of the 8th IEEE Conference on Computer Science and Electronic Engineering. Colchester, UK: IEEE, 2017. 168−173 [3] Qian L H, Liu H H T. Dynamics and control of a quadrotor with a cable suspended payload. In: Proceedings of the 30th IEEE Canadian Conference on Electrical and Computer Engineering. Windsor, Canada: IEEE, 2017. 1−4 [4] 王诗章, 鲜斌, 杨森. 无人机吊挂飞行系统的减摆控制设计. 自动化学报, 2018, 44(10): 45−54Wang Shi-Zhang, Xian Bin, Yang Sen. Anti-swing controller design for an unmanned aerial vehicle with a slung-load. Acta Automatica Sinica, 2018, 44(10): 45−54 [5] Cruz P J, Oishi M, Fierro R. Lift of a cable-suspended load by a quadrotor: A hybrid system approach. In: Proceedings of the 2015 IEEE American Control Conference. Chicago, IL, USA: IEEE, 2015. 1887−1892 [6] Cruz P J, Fierro R. Cable-suspended load lifting by a quadrotor UAV: Hybrid model, trajectory generation, and control. Autonomous Robots, 2017, 41(8): 1629−1643 doi: 10.1007/s10514-017-9632-2 [7] de Angelis E L, Giulietti F, Pipeleers G. Two-time-scale control of a multirotor aircraft for suspended load transportation. Aerospace Science and Technology, 2019, 84(2019): 193−203 [8] Isidori A, Wang B, Zhuang S X. Nonlinear Control System I. Beijing: Electronics Industry Press, 2012. [9] 梁晓, 胡欲立. 基于微分平滑的四旋翼运输系统轨迹跟踪控制. 控制理论与应用, 2019, 36(4): 525−532Liang Xiao, Hu Yu-Li. Trajectory control of a quadrotor with a cable-suspended load based on differential Flatness. Control Theory and Applications, 2019, 36(4): 525−532 [10] Lee T, Leok M, McClamroch N H. Geometric tracking control of a quadrotor UAV on SE(3). In: Proceedings of the 49th IEEE Conference on Decision and Control (CDC). Atlanta, GA, USA: IEEE, 2010. 5420−5425 [11] 王宁, 王永. 基于模糊不确定观测器的四旋翼飞行器自适应动态面轨 迹跟踪控制. 自动化学报, 2018, 44(4): 685−695Wang Ning, Wang Yong. Fuzzy uncertainty observer based adaptive dynamic surface control for trajectory tracking of a quadrotor. Acta Automatica Sinica, 2018, 44(4): 685−695 [12] 陈乐生, 王以伦. 多刚体动力学基础. 哈尔滨: 哈尔滨工程大学出版社, 1995.Chen Le-Sheng, Wang Yi-Lun. Multi-body System Dynamics. Harbin: Harbin Engineering University Press, 1995. [13] 赵杰梅, 胡忠辉. 基于动态反馈的AUV水平面路径跟踪控制. 浙江大学学报(工学版), 2018, 52(18): 1467−1481Zhao Jie-Mei, Hu Zhong-Hui. Path following control of AUV in horizontal plane based on dynamic feedback control. Journal of Zhejiang University (Engineering Science), 2018, 52(18): 1467−1481 [14] 苏善伟, 朱波, 向锦武, 林岩. 非线性非最小相位系统的控制研究综 述. 自动化学报, 2015, 41(1): 9−21Su Shan-Wei, Zhu Bo, Xiang Jin-Wu, Lin Yan. A survey on the control of nonlinear non-minimum phase systems. Acta Automatica Sinica, 2015, 41(1): 9−21 [15] 易国, 毛建旭, 王耀南, 郭斯羽, 缪志强. 非完整移动机器人目标环绕动态反馈 线性化控制. 控制理论与应用, 2017, 34(7): 895−902 doi: 10.7641/CTA.2017.60962Yi Guo, Mao Jian-Xu, Wang Yao-Nan, Guo Si-Yu, Miao Zhi-Qiang. Circumnavigation of a target with nonholonomic mobile robots via dynamic feedback linearization. Control Theory and Applications, 2017, 34(7): 895−902 doi: 10.7641/CTA.2017.60962 [16] Fliess M. Flatness and defect of non-linear systems: Introductory theory and examples. International Journal of Control, 1995, 61(6): 1327−1361 doi: 10.1080/00207179508921959 [17] Taniguchi T, Eciolaza L, Sugeno M. Quadrotor control using dynamic feedback linearization based on piecewise bilinear models. In: Proceedings of the 2014 IEEE Symposium on Computational Intelligence in Control and Automation (CICA). Orlando, USA: IEEE, 2014. 1−7 [18] Choi I H, Bang H C. Quadrotor-tracking controller design using adaptive dynamic feedback-linearization method. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2013, 228(12): 2329−2342 [19] Rego B S, Raffo G V. Suspended load path tracking control based on zonotopic state estimation using a tilt-rotor UAV. In: Proceedings of the 19th IEEE International Conference on Intelligent Transportation Systems (ITSC). Rio de Janeiro, Brazil: IEEE, 2016. 1445−1451 期刊类型引用(7)
1. 张喜铭,徐欢,杨秋勇,高伟,张睿喆. 考虑电力系统数据治理智能化的数据库生成方法研究. 制造业自动化. 2024(02): 160-165+171 . 百度学术
2. 祝健杨,辛明勇,代奇迹. 基于深度神经网络的数字电网边缘侧数据迁移. 电子设计工程. 2024(08): 55-58+63 . 百度学术
3. 张福胜,张扬,薛志胜,彭驰. 基于物联网的煤矿井下机电设备状态智能监测. 机械与电子. 2024(04): 45-49 . 百度学术
4. 于秋玲,梁锦照,陈康平. 基于云边协同的交直流混联电网线路过负荷协调控制方法. 沈阳工业大学学报. 2024(03): 241-247 . 百度学术
5. 陈琳,赵冬,尹建兵,王明昶,徐航. QoS数据不确定性下电力业务服务自适应选取模型. 电子设计工程. 2024(13): 131-134+139 . 百度学术
6. 李远征,任潇,葛磊蛟,彭靖轩,徐秋实,李曦. 基于可逆固体氧化物电池的电氢耦合微电网全生命周期规划-运营研究. 中国电机工程学报. 2024(13): 5169-5185 . 百度学术
7. 金义,宋永春,马路遥,曹雨,宋树,贺琳. 高维时序空间大数据关联特征精准挖掘算法. 电子设计工程. 2024(16): 161-165 . 百度学术
其他类型引用(3)
-