伺服压力机与传统机械压力机相比，滑块运动轨迹可以进行较大范围调整，能够实现对冲压成形速度和成形压力的精确控制^[1]，从而保证冲压过程的高效率和冲压件的高质量，具有巨大的应用潜力。伺服压力机，特别是大型伺服压力机，工作过程中主驱电机加减速频繁，负载变化大，功率和电压波动大，易对电网造成冲击，同时也会造成电能的损失^[2]。以1台公称压力为2 000 t的大型伺服压力机为例，其峰值功率可达3 000 kW，而平均功率只有300 kW左右。如果缺少有效的能量管理系统，电网侧需要配置容量大于2 000 kW的变压器以满足功率的剧烈波动，造成很大的电能浪费。因此，能量管理系统是大型伺服压力机和冲压线的重要功能模块，通过合理管控能量主回路各单元，“削峰填谷”平衡峰值功率，降低供电(进线)功率要求，避免对电网冲击和电能损失。

目前，伺服冲压线能量管理主要通过添加储能模块实现能量的有序回收与释放。常用的储能介质有电容^[3]和飞轮^[4]。电容储能采用大容量电容作为能量储存介质，实现电能的快速吸收与释放。飞轮储能采用储能电机驱动大惯量飞轮，通过控制储能电机减速或加速，使飞轮的动能发生变化。储能电机交替处于发电和电动状态，实现机械动能和电能的转换。电容储能的动态特性好，但是体积较大，造价较高；飞轮储能的能量密度大，更适合大型伺服冲压线的能量回收与释放^[5]。在基于飞轮储能的能量管理系统中，其关键技术是飞轮储能电机和主驱动电机的高精度动态同步，使两者速度反向变化，完成峰值能量的吸收与释放^[6]。

目前，飞轮储能电机的控制方法是：根据伺服压力机的理论运动规划，给定飞轮电机的参考转速，然后采用传统的比例积分微分(PID)反馈控制^[7]，根据直流母线的电压波动经过反馈环节对参考转速作出修正，从而控制飞轮的实际转速，稳定直流母线电压，达到能量管理的目标。上述反馈控制方案中，PID参数的调整难度大；同时，无法避免转速动态误差的产生；此外，由于飞轮惯量较大，储能电机动态特性较差，难以快速响应大范围的调整^[8]。针对这些问题，本文提出了基于飞轮转速动态规划的控制算法，根据压力机的实际功率需求，提前完成飞轮运动的精确规划，求解出储能电机的转速控制参数。

规划控制方案对伺服压力机功率模型的准确性要求较高。大型伺服压力机工作于周期性工况，各周期内的功率及负载情况变动较小。为提高模型的准确性，研究者们引入了自适应算法，通过迭代建立精确的压力机功率模型。在自适应算法中，最小方差法^[9]、零极点相消法、极点配置法^[10]等都可以较好地解决线性系统的自适应问题，但对非线性系统的兼容性差^[11]。针对非线性系统的自适应问题，研究学者们提出了最小二乘算法^[12]、神经元模型^[13]和模糊算法^[14]等，并进行了实验研究；在此基础上，侯忠生等提出了梯度投影算法，克服了非线性系统参数估计中可能出现的奇异问题，具有良好的鲁棒性^[15]。

本文采用梯度投影算法，迭代修正伺服压力机的实际功率模型，以建立精确的功率模型。在此基础上，将冲压线的进线功率波动最小作为优化目标，进行飞轮转速的规划求解，并采用B-样条曲线进行转速拟合，从而保证飞轮转速的光滑和连续。最后，使用Simulink建立了冲压线的理论仿真模型，对该能量管理方法的效果进行仿真分析和验证。

1 飞轮转速规划原理 1.1 基于梯度投影算法的功率模型

本小节采用梯度投影算法建立修正因数的计算准则，实现对进线功率模型的修正，从而建立精确的伺服压力机功率模型。

在不考虑电容充放电功率的假设下，伺服压力机理想的进线功率模型为

$ {P_{{\rm{ In }}}} = {P_{{\rm{ Flywheel }}}} + {P_{{\rm{ Prsdrv }}}}. $

(1)

即进线功率P_In等于飞轮功率P_Flywheel与压力机主驱动功率P_Prsdrv之和。实际进线功率模型为

$ P_{{\rm{In}}}^\prime = P_{{\rm{Flywheel}}}^\prime + P_{{\rm{Prsdrv}}}^\prime + {P_{{\rm{ others }}}}. $

(2)

其中：P′_In为实际进线功率，P′_Flywheel为飞轮实际功率，P′_Prsdrv为压力机实际功率，P_others为其他功率损耗。

实际工作中难以测得P′_Flywheel、P′_Prsdrv与P_others的真实值，因此考虑对实际进线功率模型进行拟合求解。

引进函数β(t)和θ(t)，使得：

$ {\beta (t){P_{{\rm{FIywheel}}}}(t) \to P_{{\rm{Flywheel}}}^\prime (t).} $

(3)

$ {\theta (t){P_{{\rm{Prsdrv}}}}(t) \to P_{{\rm{Prsdrv}}}^\prime (t) + {P_{{\rm{others}}}}(t).} $

(4)

压力机功率曲线P_Prsdrv(t)通常连续且存在零点，在零点附近P_Prsdrv(t)的值将远小于其他工作状态时的值，若此时恰好出现冲击载荷，即功率存在阶跃性波动，则θ(t)取值需非常大，才能满足式(4)逼近，这将极大减慢迭代速度。据此本文提出改进式(5)，

$ \theta (t)[{P_{{\rm{Prsdrv}}}}(t) + {\rm{lft}}] - {\rm{lft}} \to P_{{\rm{Prsdrv}}}^\prime (t) + {P_{{\rm{ others }}}}(t). $

(5)

其中lft=－min(P_Prsdrv(t))+1 000，确保与θ(t)的被影响量处于同一数量级，以保证迭代的收敛速度在各部分基本保持一致。由于飞轮功率为求解所得，函数β(t)和θ(t)的值可通过计算进行实时更新，无须进行额外调整。

因此，伺服压力机的进线功率模型可修正为：

$ \begin{array}{l} {P_{{\rm{ In }}}}(t) = \beta (t){P_{{\rm{ Flywheel }}}}(t) + \theta (t)[{P_{{\rm{ Prsdrv }}}}(t) + {\rm{lft}}] - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{Ift}} = \beta (t){P_{{\rm{Flywheel}}}}(t) + \theta (t){P_{{\rm{PrsdrvCor}}}}(t) - {\rm{lft}}.\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {P_{{\rm{ PrsdrvCor }}}}(t) = {P_{{\rm{ Prsdrv }}}}(t) + {\rm{lft}}. \end{array} $

(6)

进一步定义系统的进线功率模型误差为

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{error}}(t + 1) = {P_{{\rm{In}}}}(t) - P_{{\rm{In}}}^\prime (t) = \\ \beta (t){P_{{\rm{ Flywheel }}}}(t) + \theta (t){P_{{\rm{ Prsdrv }}}}(t) - [P_{{\rm{ In }}}^\prime (t) + {\rm{lft}}] = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \beta (t)[{P_{{\rm{ Flywheel }}}}(t) - {\bf{Fct}} (t){\mathit{\boldsymbol{P}}_{{\rm{Ccl}}}}(t)]. \end{array} $

(7)

其中：

$ \begin{array}{*{20}{c}} { {\bf{Fct}} (t) = [{\beta ^{ - 1}}(t), - {\beta ^{ - 1}}(t)\theta (t)],}\\ {{\mathit{\boldsymbol{P}}_{{\rm{Ccl}}}}(t) = {{[P_{{\rm{In}}}^\prime (t) + {\rm{lft}} ,{P_{{\rm{PrsdrvCor}}}}(t)]}^{\rm{T}}}.} \end{array} $

(8)

称Fct为修正因数。

基于梯度投影算法对修正因数Fct提出估计准则，

$ \begin{array}{*{20}{l}} {J({\bf{Fct}}) = {{[{P_{{\rm{Flywheel}}}}(t) - {\bf{Fct}}(t + T){\mathit{\boldsymbol{P}}_{{\rm{Ccl}}}}(t)]}^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} c{{[{\bf{Fct}}(t + T) - {\bf{Fct}}(t)]}^2}.} \end{array} $

(9)

式中c[Fct(t+T)－Fct(t)]²项的引入是为了惩罚修正因数的过度变化。

由于

$ \begin{array}{*{20}{c}} { {\bf{Fct}} (t + T){\mathit{\boldsymbol{P}}_{{\rm{Ccl}}}}(t) = {\bf{Fct}} (t){\mathit{\boldsymbol{P}}_{{\rm{Ccl}}}}(t) + }\\ {[{\bf{Fct}}(t + T) - {\bf{Fct}}(t)]{\mathit{\boldsymbol{P}}_{{\rm{Ccl}}}}(t),} \end{array} $

(10)

代入估计准则(9)中，得到

$ \begin{array}{*{20}{c}} {J({\bf{Fct}}) = {{[{P_{{\rm{ Flywheel }}}}(t) - {\bf{Fct}}(t + T){\mathit{\boldsymbol{P}}_{{\rm{Ccl}}}}(t)]}^2} + }\\ {[c + {{\left\| {{\mathit{\boldsymbol{P}}_{{\rm{Ccl}}}}(t)} \right\|}^2}]{{[{\bf{Fct}}(t + T) - {\bf{Fct}}(t)]}^2}.} \end{array} $

(11)

令一阶导数为零，求解得到Fct(t)的周期更新律为

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\bf{Fct}}(t) = \\ \left\{ \begin{array}{l} {\bf{Fct}}(t - T) + \frac{{\alpha \cdot \mathit{\boldsymbol{P}}_{{\rm{Ccl}}}^{\rm{T}}(t - T) {\rm{error}} (t - T + 1)}}{{c + \mathit{\boldsymbol{P}}_{{\rm{Ccl}}}^{\rm{T}}(t - T){\mathit{\boldsymbol{P}}_{{\rm{Ccl}}}}(t - T)}},\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} t \ge T;\\ {\bf{Fc}}{{\bf{t}}^0}(t),\;\:0 < t < T. \end{array} \right. \end{array} $

(12)

其中：Fct⁰(t)为有界函数；引入α作为学习率，且c>0，0 < α·β_max < 2，α和β的符号相同，β_max为β(t)的上界。

1.2 飞轮转速求解和描述

基于修正后的进线功率模型，考虑到进线功率的平均值在具体负载确定的情况下是一个未知常量，为了便于建立优化模型，将进线功率的平方和最小作为优化目标，以反映进线功率的波动性，

$ \begin{array}{*{20}{c}} {J = \int_0^{{T_{{\rm{ cycle }}}}} {P_{{\rm{In}}}^2} {\rm{d}}t = }\\ {\int_0^{{T_{{\rm{ cycle }}}}} ( \beta (t){J_{\rm{F}}}{{\dot \omega }_{\rm{F}}}(t){\omega _{\rm{F}}}(t) + }\\ {\theta (t){P_{{\rm{PrsdrvCor}}}}(t) - {{({\rm{lft}})}^2}){\rm{d}}t.} \end{array} $

(13)

其中：T_cycle为周期运动时间，J_F为飞轮转动惯量，ω_F为飞轮转速。${P_{{\rm{Flywheel }}}}(t) = {J_{\rm{F}}}{\dot \omega _{\rm{F}}}(t){\omega _{\rm{F}}}(t)$，满足

$ {\rm{min}}J = {\rm{min}}\left( {\int_0^{{T_{{\rm{ cycle }}}}} {P_{{\rm{ In }}}^2} {\rm{d}}t} \right). $

(14)

求解量为飞轮转速曲线。约束条件如下：

$ \begin{array}{l} \begin{array}{*{20}{c}} {{\omega _{{\rm{min}}}} < {\omega _{\rm{F}}}(t) < {\omega _{{\rm{max}}}},}\\ {|{J_{\rm{F}}}{{\dot \omega }_{\rm{F}}}(t)| < {T_{{\rm{max}}}},}\\ {|{{\ddot \omega }_{\rm{F}}}(t)| < {{\ddot \omega }_{{\rm{max}}}},} \end{array}\\ |{J_{\rm{F}}}{{\dot \omega }_{\rm{F}}}(t){\omega _{\rm{F}}}(t)| < {P_{{\rm{max}}}}. \end{array} $

(15)

即飞轮的转速、飞轮的角加速度、飞轮的角加加速度、飞轮的功率存在范围约束。

规划算法的具体流程(如图 1所示)如下：

1) 在冲压线启动前，通过主驱动电机转速与扭矩之间的关系，构建初始的进线功率模型，储能电机与主驱动电机的能量流动方向相反，得到初始飞轮转速曲线。

2) 在冲压线启动后，运行初始飞轮转速曲线，每半个周期更新并修正一次主驱动电机的负载模型，同时相应地求解新的飞轮储能电机转速曲线。半个周期的设定方式可以保证飞轮运动的连续性。

采用4次B-样条曲线进行飞轮速度拟合。实际生产中，一个冲压周期时长约3~5 s。在Intel Core i5-8250U CPU(1.60 GHz)测试平台下，进行一次完整周期的飞轮转速曲线规划需要约2 s的CPU时间。由于计算需要花费一定时间，因此无法在完整周期执行结束后进行更新律的计算，这会大大减慢修正因数的更新速率，并且不能保证飞轮运动的连续性。为此，本文提出将每个冲压周期分割为两部分，每半个冲压周期进行一次更新律的计算，在下一个周期的对应阶段完成修正因数的更新。由于计算规模减小，每次计算花费的CPU时间仅约为1 s，通过选用高性能工控机或工作站，可进一步缩短计算时间。此外，由于B-样条曲线的局部性，每半个周期计算得到新的转速后，另半个周期内的飞轮转速不发生变化，可实现新-旧飞轮曲线的衔接。

图 2为飞轮转速更新过程示意图。其中：C_k(t)表示当前飞轮转速曲线，C_k+1(t)表示下一次更新后的飞轮转速曲线。对于每半个冲压周期，计算该半个冲压周期内的误差以进行修正因数的更新，并求解得到新的飞轮转速曲线，另半个冲压周期仍使用旧的飞轮转速曲线，因此新曲线与旧曲线总是能存在重合部分，从而能够实现飞轮运动的连续，并提高模型迭代拟合的效率。

2 约数周期正弦扰动下的进线功率仿真

由于冲压线能量系统内的各分量(例如送料系统)的周期与冲压周期相关性很强，也就是说，与冲压周期相关的扰动占多数，为此在模拟工况中引入约数周期正弦扰动，即扰动信号的周期T与冲压周期T_cycle满足nT=T_cycle，其中n=1, 2, …。下面进一步观察在约数周期正弦扰动下模型的性能。

采用半周期自适应飞轮转速规划算法分析冲压线在该扰动下的相应控制结果。大型伺服线中包含4台压力机，后3台压力机与第1台压力机的虚轴相位差及功率比率分别为(60°, 0.6)、(120°, 0.5)、(180°, 0.3)。

利用MATLAB程序仿真，仿真参数设置如表 1所示。

仿真使用的是整线功率因数为1.1、飞轮功率因数为0.9的实际负载模型，引入幅值为600 kW、周期为2.5 s(即冲压周期的一半，n=2)的正弦波扰动，即

$ \begin{array}{*{20}{c}} {P_{{\rm{In}}}^\prime = 0.9{P_{{\rm{Flywheel}}}}(t) + 1.1{P_{{\rm{PrsLn}}}}(t) + }\\ {600{\kern 1pt} {\kern 1pt} {\rm{sin}}\left( {\frac{{2\pi t}}{{2.5}}} \right).} \end{array} $

(16)

图 3为第1台压力机(P_Prsdrv)与整线(P_PrsLn)的功率曲线。由于第1台压力机的功率占比最大，整线功率曲线与第1台压力机接近，故可将第一台压力机的功率作为初始负载模型提供给飞轮转速规划程序。规划程序从仿真结果中获得当前时刻的进线功率值，并向仿真环境输出当前时刻的飞轮转速值。

由图 4的仿真结果可以看出，在几个周期后，拟合的进线功率与实际进线功率高度吻合，且进线功率的波动趋于平缓。同时，由图 5可以看出，飞轮转速也因负载模型的修正而逐渐产生变化，且转速曲线连续光滑。

定义整个周期的误差(周期误差)为误差函数平方在整个周期时间内的积分，即

$ {\rm{erro}}{{\rm{r}}_{{\rm{ cycle }}}} = \int_0^{{T_{{\rm{ cycle }}}}} { {\rm{erro}}{{\rm{r}}^2}} (t){\rm{d}}t. $

(17)

周期误差随系统运行的变化曲线如图 6所示。可以看出，从迭代开始，周期误差指标呈现指数级快速衰减，直至稳定在10^-4下，说明对于约数周期正弦扰动下的进线功率模型，程序拟合效果较好，程序拟合的进线功率曲线非常接近实际的进线功率曲线。

对于其他约数周期，进行数值仿真可以得到相似的结论。

3 飞轮转速自适应规划Simulink仿真 3.1 Simulink建模

本节将飞轮储能的完整过程在Simulink中进行仿真，检验算法在复杂状况下的工作性能，观察以进线功率波动最小为优化目标的算法引起的母线电压波动和母线电流波动状况。仿真模型主要包括飞轮转速计算模块、电机模拟模块、信号输入模块、运动规划与反馈模块等。

运动规划主要由飞轮速度曲线计算模块完成，每半个周期对曲线进行一次更新；反馈模块包括母线电压波动反馈和母线电流波动反馈。给母线电压和母线电流一个设定值，母线电压围绕设定值的波动通过控制器输入到飞轮伺服电机的速度环，母线电流的波动值通过控制器输入到飞轮伺服电机的电流环。加入反馈环节可以使系统形成闭环控制，实现飞轮电机动态抵消母线电压和电流的波动。

3.2 仿真结果分析 3.2.1 周期载荷下进线功率对比

对不同控制情况下的功率流进行观测，对比功率流波动情况。如图 7所示，根据伺服压力机理想轨迹进行飞轮转速规划得到的功率流曲线为蓝色，采用本文提出的自适应转速规划算法的功率流曲线为橙色，红色曲线为电压和电流反馈控制下的功率流。

可以发现，蓝色曲线功率流的波动为0~630 kW(周期运动阶段的数据，下同)，且存在大量的功率突变，在功率为零时的停滞时间较长。在进行储能电机自适应转速规划后，功率流的波动为0~ 580 kW，功率分布更为集中，功率为零的停滞时间明显变短，且有随周期增加迭代缩短的趋势，说明在算法的迭代下，飞轮转速不断调整，功率流的波动得到改善。与自适应规划相比，添加电压反馈和电流反馈控制后，功率曲线变动较小，具有减小的趋势但不明显。

3.2.2 周期载荷下母线电压波动对比

母线电压波动在自适应规划前后的变化更加显著，如图 8所示。蓝色曲线表示根据伺服压力机理想轨迹进行飞轮转速规划的电压波动，红色曲线表示采用自适应转速规划后的电压波动。未进行自适应规划和飞轮储能时，电压从600 V波动到900 V (周期运动阶段的数据，下同)，而在采用了迭代的自适应规划算法后，电压波动稳定后为600~ 650 V。这是由于采用运动规划后，功率波动减小，功率为零的停滞时间更短，充电时间更短，因此母线电压的升幅会减小。

从对典型仿真结果的对比和分析中可以发现，运动规划迭代算法在面对周期性载荷时，可以在开始的几个周期内对飞轮转速作出调整，收敛速度快，在实际仿真模型中算法调整效果显著，功率流波动在2~3个周期内就已经有了明显收敛，并在后续阶段进行微调。

在母线电压波动调整方面，由于功率流较快得到调整，而且最终调整目标是功率波动为一恒定正值，因此每个周期中功率近似为零的时长会随迭代而减小，即充电时间随着周期迭代快速减少，因而母线电压波动得到很大改善。

4 结论

本文将大型压力机能量管理系统中的自动控制问题转化为飞轮转速规划问题，利用自适应算法对伺服压力机的功率模型进行迭代求解，建立精确的功率模型，以进线功率波动最小为优化目标，求解飞轮转速，并采用B-样条曲线拟合，保证转速的连续和光滑，建立了基于飞轮转速自适应规划的大型伺服冲压线能量管理方法。随后对算法的合理性进行了验证，利用MATLAB程序仿真，评估了算法在引入约数周期正弦扰动下的性能，结果表明算法具有较好的优化效果，以半个运行周期为迭代周期，优化速度更快，也避免了飞轮预期转速突变的问题。利用Simulink建立了冲压线的理论仿真模型，仿真结果表明飞轮转速自适应规划算法对于进线功率、母线电压波动都有显著的改善。本研究为伺服压力机的能量管理提供了一种有效的思路，有利于减小系统峰谷差值，降低能耗。