CN101520642B - 火电厂单元机组协调控制系统被控对象动态参数整定方法 - Google Patents

Abstract

火电厂单元机组协调控制系统被控对象动态参数整定方法，属于工业控制技术领域。具体步骤为：一、根据现场扰动试验可复现性最好的原则，确定通过试验得到的参数为T，T₁，K₁，K₂；二、通过七个物理约束条件对其余十三个被控对象动态参数进行计算，得出十三个被控对象动态参数。本发明以最少和最容易测定的试验数据为前提，应用被控对象本身内在的物理规律和约束条件对其动态参数进行最优化求解，通过计算整定出被控对象全部动态参数。减少了确定协调控制系统被控对象动态参数过程中要求进行大量现场扰动试验的制约，为协调控制系统的仿真研究和现场整定提供了一种便捷的方法。

Description

火电厂单元机组协调控制系统被控对象动态参数整定方法

技术领域

本发明涉及一种动态参数的整定方法，特别是一种单元机组协调控制系统被控对象动态参数整定方法，属于工业控制技术领域。

背景技术

火电厂单元机组的被控对象锅炉-汽机是一个复杂的热力学过程，具有复杂的非线性、多变量和时变特征。一般只能采用现场扰动试验的方法来获取被控对象的动态参数。文献检索发现，沈邦盛、王世海等在青山电厂、焦作电厂分别进行了200MW单元机组被控对象动态特性试验(《热力发电》1982.1，1984.4)。由于被控对象本身内部多变量耦合及各种干扰因素，试验过程往往存在较大的偏差，获得的试验数据甚至无法表征出其正确的物理规律，获得有效数据很大程度依赖经验和偶然性。即便试验获取了一组较理想的参数，也很难满足不同机组、不同工况的仿真研究和现场调整需要。火电厂机组运行的安全性和经济性又不允许进行过多和反复的动态特性试验，因此，单纯由试验获取被控对象的动态参数受到很大的局限。进一步的文献检索发现，国内外对单元机组协调控制系统的分析和研究大多集中于协调控制回路和PID调节器参数的分析和整定。对被控对象的动态参数基本都局限于试验结果。迄今未能发现一种能够较好满足仿真研究或现场试验要求的整定方法。

发明内容

本发明针对现有技术的不足，提出了一种基于直接能量平衡原理的单元机组协调控制系统被控对象动态参数的整定方法。基于对一个2×4维的单元机组被控对象的多变量模型的分析，提出了在既有模型结构上，以最少和最容易测定的试验数据为前提，应用被控对象本身内在的物理规律和约束条件对其动态参数进行最优化求解，通过计算整定出其余的被控对象动态参数。

本发明是通过如下技术方案实现的，具体包括以下步骤：

步骤一，根据现场扰动试验可复现性最好的原则，确定试验参数。具体为采用汽包锅炉-汽机的单元机组作为协调控制系统的被控对象时，其动态模型一般可简化表达为2×4结构。其中，汽机调门开度μ_T的扰动可以视作被控对象的外部扰动，而锅炉燃烧率指令μ_B的扰动则作为单元机组本身的内部扰动。机组功率N_T的动态特性可以反映单元机组对电网负荷需求的响应能力。机前压力P_T则反映了单元机组被控对象之间的能量和质量平衡，表征了机组运行的稳定性。

被控对象的传递函数结构可分别由锅炉燃烧率和汽机调门开度扰动试验确定，并表达为：

(1)汽机调门开度μ_T扰动时各参数的传递函数：

N_T对μ_T扰动的传递函数： $W_{\mathrm{NT}}=\frac{N_T\left(s\right)}{{\mathrm{\μ}}_T\left(s\right)}=\frac{K_{5}s}{{(1+T_{4}s)}^2}---\mathrm{MW}/\%$

P_T对μ_T扰动的传递函数： $W_{\mathrm{PT}}=\frac{P_T\left(s\right)}{{\mathrm{\μ}}_T\left(s\right)}=-(K_2+\frac{K_1}{1+T_{1}s})---\mathrm{MPa}/\%$

P₁对μ_T扰动的传递函数： $W_{1T}=\frac{P_1\left(s\right)}{{\mathrm{\μ}}_T\left(s\right)}=\frac{K_{3}s}{{(1+T_{2}s)}^2}---\mathrm{MPa}/\%$

P_D对μ_T扰动的传递函数： $W_{\mathrm{DT}}=\frac{P_D\left(s\right)}{{\mathrm{\μ}}_T\left(s\right)}=-\frac{K_4}{{(1+T_{3}s)}^2}---\mathrm{MPa}/\%$

(2)锅炉燃烧率指令μ_B扰动时各参数的传递函数：

N_T对μ_B扰动的传递函数： $W_{\mathrm{NB}}=\frac{N_T\left(s\right)}{{\mathrm{\μ}}_B\left(s\right)}=\frac{a}{(1+{\mathrm{\η}}_3\mathrm{Ts})(1+\mathrm{Ts})}---\mathrm{MW}/t/h$

P_T对μ_B扰动的传递函数： $W_{\mathrm{PB}}=\frac{P_T\left(s\right)}{{\mathrm{\μ}}_B\left(s\right)}=\frac{b}{{(1+\mathrm{Ts})}^2}---\mathrm{MPa}/t/h$

P₁对μ_B扰动的传递函数： $W_{1B}=\frac{P_1\left(s\right)}{{\mathrm{\μ}}_B\left(s\right)}=\frac{c}{{(1+{\mathrm{\η}}_1\mathrm{Ts})}^2}---\mathrm{MPa}/t/h$

P_D对μ_B扰动的传递函数： $W_{\mathrm{DB}}=\frac{P_D\left(s\right)}{{\mathrm{\μ}}_B\left(s\right)}=\frac{d}{{(1+{\mathrm{\η}}_2\mathrm{Ts})}^2}---\mathrm{MPa}/t/h$

式中：N_T      ______机组功率

P_T        ______汽机机前压力

P₁        ______汽机一级压力

P_D，P_d    ______汽包压力

对燃煤机组W_NB的二个极点基本相等，而燃油机组两极点相差较大。因此，对燃煤机组，上述传递函数中可以整理成具有二重极点的二阶惯性环节，即T′＝η₃T。

对上述模型的阶跃响应求解析式：

(1)在汽机调门μ_T扰动下：

${\mathrm{\ΔP}}_T\left(t\right)=L^{-1}[-(\frac{K_1}{1+T_{1}s}+K_2)\×\frac{1}{s}]=-K_2-K_1(1-e^{-\frac{t}{T_1}})$

${\mathrm{\ΔP}}_1\left(t\right)=L^{-1}[\frac{K_{3}s}{{(1+T_{2}s)}^2}\×\frac{1}{s}]=\frac{K_3}{{T_2}^2}\×t\×e^{-\frac{t}{T_2}}$

${\mathrm{\ΔP}}_d\left(t\right)=L^{-1}[-\frac{K_4}{1+T_{3}s}\×\frac{1}{s}]=-K_4(1-e^{-\frac{t}{T_3}})$

${\mathrm{\ΔN}}_T\left(t\right)=L^{-1}[\frac{K_{5}s}{{(1+T_{4}s)}^2}\×\frac{1}{s}]=\frac{K_5}{{T_4}^2}\×t\×e^{-\frac{t}{T_5}}$

(2)在锅炉燃烧率指令μ_B扰动下：

从模型结构中可以看出，需确定的参数共有十七个，其中：九个比例增益系数K₁～K₅、a、b、c、d，五个时间常数T、T₁～T₄，三个时间常数的系数η₁～η₃。这些模型参数都包含了一定的热力学规律，因此，首先从这些众多的参数中选择出数目最少，同时又最容易测量的参数作为试验测量量。

根据实践经验，在μ_B和μ_T扰动下，ΔP_T的复现度较高。这样，确定四个参数T，T₁，K₁，K₂作为试验测量量。

步骤二，通过下述七个物理约束条件对十三个被控对象动态参数进行整定计算。

条件一：μ_B阶跃扰动下，P₁/P_T理论上代表了汽机调门的开度，从动态扰动的角度看，它应不随时间变化。将上述条件改写成数学表达式：即：

(P_T0和P₁₀指基准工况下的运行参数)

可见只要选取η₁≈1即可满足条件1。同样道理，因机组功率响应几乎与P₁响应同步，故应选取η₃≥η₁，且η₃≈η₁。

条件二：因为μ_T＝P₁/P_T，所以在μ_T阶跃扰动下，P₁/P_T应同比阶跃，即：

$\frac{P_{10}+{\mathrm{\ΔP}}_1\left(t\right)}{P_{T0}+{\mathrm{\ΔP}}_T\left(t\right)}=\frac{P_{10}+{\mathrm{\ΔP}}_1(\∞)}{P_{T0}+{\mathrm{\ΔP}}_T(\∞)}$

$\⇔\frac{P_{10}+{\mathrm{\ΔP}}_1(\∞)}{P_{T0}+{\mathrm{\ΔP}}_T(\∞)}=\frac{P_1\left(t\right)-{\mathrm{\ΔP}}_1(\∞)}{P_T\left(t\right)-{\mathrm{\ΔP}}_T(\∞)}$

$\⇔\frac{P_{10}}{P_{T0}-(K_1+K_2)}=\frac{\frac{K_3}{T_2}\×t\×e^{-\frac{t}{T_2}}}{K_1{e}^{-\frac{t}{T_1}}}$

化简可得： $\frac{P_{10}}{P_{T0}-(K_1+K_2)}\×\frac{K_1}{K_3}=\frac{1}{{T_2}^2}\×t\×e^{(\frac{t}{T_1}-\frac{t}{T_2})}$ ①

上式左端是一常数，故应该满足： $\frac{d}{\mathrm{dt}}[\frac{1}{{T_2}^2}\×t\×e^{(\frac{1}{T_1}-\frac{1}{T_2})t}]\≈0$

$\⇒\frac{1}{{T_2}^2}\×e^{(\frac{1}{T_1}-\frac{1}{T_2})t}[1+t\×(\frac{1}{T_1}-\frac{1}{T_2})]\≈0$

令： $m=\frac{T_2}{T_1},$ 得函数： $F_1\left(t\right)=\frac{1}{{T_1}^2}\×\frac{1}{m^2}\×e^{(1-\frac{1}{m})\frac{t}{T_1}}[1+\frac{t}{T_1}(1-\frac{1}{m})]\≈0$

选用平方误差准则，即要求出一最佳值，使得

为最小，即F₁(t)最接近零。显然，m应满足0＜m＜1，否则，F₁(t)不收敛。

${\∫}_0^{\∞}{F_1}^2\left(t\right)\mathrm{dt}={\∫}_0^{\∞}\frac{1}{{T_1}^4}\×\frac{1}{m^4}\×e^{2(1-\frac{1}{m})\frac{t}{T_1}}\×[1+\frac{2T}{T_1}(1-\frac{1}{m})+\frac{t^2}{{T_1}^2}{(1-\frac{1}{m})}^2]\mathrm{dt}$

上式中再令 $A=(1-\frac{1}{m})\frac{1}{T_1},$ 则有：

${\∫}_0^{\∞}{F}_1^2\left(t\right)\mathrm{dt}=\frac{1}{{T_1}^4}\frac{1}{m^4}{\∫}_0^{\∞}{e}^{2\mathrm{At}}[1+2\mathrm{At}+{\left(\mathrm{At}\right)}^2]\mathrm{dt}=\frac{1}{{T_1}^4}\×\frac{1}{m^4}[-\frac{1}{2A}+2A\frac{1}{{\left(2A\right)}^2}-A^2\frac{2}{{\left(2A\right)}^2}]$

$=-\frac{1}{4{\mathrm{Am}}^4{T_1}^4}=-\frac{1}{{T_1}^3}\×\frac{1}{{4m}^3(m-1)}=\frac{1}{{T_1}^3}\×\frac{1}{{4m}^3(1-m)}$

要使

最小，必须满足： $\frac{d}{\mathrm{dt}}\left[\frac{1}{m^3(1-m)}\right]=0$

即：-3m^-4(1-m)^-1-m^-3(1-m)^-2＝0

$m=\frac{3}{4}=0.75$

∴T₂＝0.75T₁

将T₂带入式①，取该式右端为最大值，可解出 $K_3=\frac{3e}{16}\×\frac{P_0}{P_{T_0}-(K_1+K_2)}\×K_1{T}_1$ 条件三：在μ_T扰动下，达到稳态时，因蒸汽流量不变，故有ΔP_d(∞)＝ΔP_T(∞)。

$\left.\begin{array}{c}{\mathrm{\ΔP}}_d(\∞)=K_4\\ {\mathrm{\ΔP}}_T(\∞)=K_1+K_2\end{array}\right\}\⇒K_4=K_1+K_2$

K₁，K₂值可由试验得出，这样同时可求得K₄之值。

条件四：在μ_T扰动下，热量 $Q\left(t\right)=P_1+C_k\×\frac{{\mathrm{dP}}_d}{\mathrm{dt}}$ 应保持不变，即：

$\mathrm{\ΔQ}\left(t\right)={\mathrm{\ΔP}}_1\left(t\right)+C_k\×\frac{\mathrm{d\Δ}{P}_d}{\mathrm{dt}}\≈0$

$\⇒\mathrm{\ΔQ}\left(t\right)=\frac{K_3}{T_2^2}t\×e^{-\frac{t}{T_2}}-C_k\frac{K_4}{T_3}{e}^{-\frac{t}{T_3}}\≈0$

仍用平方误差准则，即寻求合适C_k和T₃值，使得为最小，即接近于零。

同样，令： $\frac{T_3}{T_2}=\mathrm{\ρ},$ 则有：

${\∫}_0^{\∞}\mathrm{\Δ}{Q}^2\left(t\right)\mathrm{dt}={\∫}_0^{\∞}{(-\frac{C_k{K}_4}{\mathrm{\ρ}{T}_2}\×e^{-\frac{t}{\mathrm{\ρ}{T}_2}}+\frac{K_3}{T_2^2}\×t\×e^{-\frac{t}{T_2}})}^2\mathrm{dt}$

上式中再令： $\frac{C_k{K}_4}{T_2}=M,$ $\frac{K_3}{T_2^2}=N.$ 则：

${\∫}_0^{\∞}\mathrm{\Δ}{Q}^2\left(t\right)\mathrm{dt}={\∫}_0^{\∞}[N^2\×t^2\×e^{-\frac{2t}{T_2}}+\frac{M^2}{{\mathrm{\ρ}}^2}\×e^{-\frac{2t}{\mathrm{\ρ}{T}_2}}-\frac{2\mathrm{MN}}{\mathrm{\ρ}}\×t\×e^{-(1+\frac{1}{\mathrm{\ρ}})\frac{t}{T_2}}]\mathrm{dt}$

$=N^2\×\frac{2}{{\left(\frac{2}{T_2}\right)}^3}+\frac{M^2}{{\mathrm{\ρ}}^2}\×\frac{1}{\frac{2}{{\mathrm{\ρT}}_2}}-\frac{2\mathrm{MN}}{\mathrm{\ρ}}\frac{1}{{\left[(1+\frac{1}{\mathrm{\ρ}})\frac{1}{T_2}\right]}^2}$

$=\frac{{C_k}^2{K_4}^2}{2\mathrm{\ρ}{T}_2}-\frac{{2K}_3{K}_4}{T_2}{C}_k\frac{\mathrm{\ρ}}{{(1+\mathrm{\ρ})}^2}+\frac{{K_3}^2}{{4T}_2}$

要使

为最小，则必须满足：

$\left\{\begin{array}{c}\frac{\∂}{{\∂C}_k}{\∫}_0^{\∞}\mathrm{\Δ}{Q}^2\left(t\right)\mathrm{dt}=0\\ \frac{\∂}{\∂\mathrm{\ρ}}{\∫}_0^{\∞}\mathrm{\Δ}{Q}^2\left(t\right)\mathrm{dt}=0\end{array}\right.$

即： $\left\{\begin{array}{c}\frac{{2C}_k{K_4}^2}{{2\mathrm{\ρT}}_2}-\frac{{2K}_3{K}_4}{T_2}\×\frac{\mathrm{\ρ}}{{(1+\mathrm{\ρ})}^2}=0\\ -\frac{{C_k}^2{K_4}^2}{{2\mathrm{\ρ}}^2{T}_2}-\frac{{2K}_3{K}_4}{T_2}\×C_k\×\frac{1-\mathrm{\ρ}}{{(1+\mathrm{\ρ})}^3}=0\end{array}\right.$

解此方程组得：

$\mathrm{\ρ}=3\⇒T_3={3T}_2=\frac{9}{4}{T}_1$

$C_k=\frac{{2K}_3}{K_4}{\left(\frac{\mathrm{\ρ}}{1+\mathrm{\ρ}}\right)}^2=\frac{9}{16}\×\frac{{2K}_3}{K_4}=\frac{9}{8}\frac{K_3}{K_4}$

条件五：在μ_B扰动下，热量信号近似为阶跃响应。

假设

的阶跃值为b。

可寻求一个最佳的η₂值，使得

为最小。

要使

为最小，必须满足：

假设η₁≈1，η₂≈1， ${\mathrm{\η}}_2=\frac{1}{\sqrt[3]{\frac{{24\mathrm{b\η}}_1^{2}T}{C_{k}c}}-1}$

条件六：作为一种无因次线性化模型，当μ_B和μ_T分别从0增至100％时，即：

$\left\{\begin{array}{c}a-(K_1+K_2)=P_{T0}\\ b+0=P_{10}\\ c-K_4=P_{d0}\\ d+0=N_{T0}\end{array}\right.=\left\{\begin{array}{c}a=P_{T0}+(K_1+K_2)\\ b=P_{10}\\ c=P_{d0}+K_4\\ d=N_{T0}\end{array}\right.$

条件七、因μ_T扰动时，P₁与N_T的响应特性几乎相同，故有：

$\frac{\mathrm{\Δ}{P}_{1\max }}{P_{1_0}}=\frac{{\mathrm{\ΔN}}_{T\max }}{N_{T_0}}$

$\stackrel{\·}{\·\·}\frac{K_3}{K_5}\frac{T_4}{T_2}=\frac{P_{1_0}}{N_{T_0}}\⇒K_5=\frac{T_4}{T_2}{K}_3\frac{N_{T_0}}{P_{1_0}}$

至此，根据步骤一，由μ_B和μ_T试验测定动态参数T，T₁，K₁，K₂，其余待定参数均由本方法确定，各参数整定结果可以归结如下：

(1)η₃≥η₁≥1，η₃≈1，η₁≈1

(2)T₂＝0.75×T₁

(3) $K_3=\frac{3e}{16}\×\frac{P_{1_0}}{P_{T_0}-(K_1+K_2)}\×K_1\×T_1$

(4) $T_3=\frac{9}{4}\×T_1$

(5)K₄＝K₁+K₂

(6) $C_K=\frac{9}{8}\×\frac{K_3}{K_4}$

(7)T₄≥T₂，T₄≈T₂

(8) $K_5=\frac{T_4}{T_4}\×K_3\×\frac{N_0}{T_{1_0}}$

(9)a＝P_T0+(K₁+K₂)

(10)b＝P₁₀

(11)c＝P_d0+K₄

(12)d＝N_T0

(13) ${\mathrm{\η}}_2=\frac{1}{\sqrt[3]{\frac{24b{\mathrm{\η}}_1^{2}T}{C_{k}c}}-1}$

本发明的有益效果：

本发明的方法以最少和最容易测定的试验数据为前提，应用被控对象本身内在的物理规律和约束条件对其动态参数进行最优化求解，通过计算整定出被控对象其余动态参数。减少了确定协调控制系统被控对象动态参数过程中要求进行大量现场扰动试验的制约。为协调控制系统的仿真研究和现场整定提供了一种便捷的方法。

附图说明

图1是单元机组协调控制系统被控对象2×4维动态模型示意图；

图2是通过本发明整定出的协调控制系统被控对象动态特性示意图。

具体实施方式

下面结合附图对本发明的具体实施作进一步描述，如图1、图2所示。本发明的具体实施步骤如下：

步骤一，分别由μ_B和μ_T试验获取如图一所示的被控对象动态参数T，T₁，K₁，K₂。

步骤二，通过计算公式可以求得其余的被控对象动态参数，

(1)η₃≥η₁≥1，η₃≈1，η₁≈1

(2)T₂＝0.75×T₁

(3) $K_3=\frac{3e}{16}\×\frac{P_{1_0}}{P_{T_0}-(K_1+K_2)}\×K_1\×T_1$

(4) $T_3=\frac{9}{4}\×T_1$

(5)K₄＝K₁+K₂

(6) $C_K=\frac{9}{8}\×\frac{K_3}{K_4}$

(7)T₄≥T₂，T₄≈T₂

(8) $K_5=\frac{T_4}{T_4}\×K_3\×\frac{N_0}{T_{1_0}}$

(9)a＝P_T0+(K₁+K₂)

(10)b＝P₁₀

(11)c＝P_d0+K₄

(12)d＝N_T0

(13) ${\mathrm{\η}}_2=\frac{1}{\sqrt[3]{\frac{{24\mathrm{b\η}}_1^{2}T}{C_{k}c}}-1}$

可以得到图2所示的满足被控对象仿真研究和现场整定所需的动态特性。

Claims (1)

1.一种火电厂单元机组协调控制系统被控对象动态参数整定方法，其特征在本发明是基于直接能量平衡原则，从单元机组被控对象物理机理的约束条件入手，通过最优化方法求得协调控制系统被控对象的十七个动态参数，其步骤为：

步骤一，根据现场扰动试验可复现性最好的原则，确定通过试验得到的参数为T，T₁，K₁，K₂；

步骤二，通过七个物理约束条件对其余十三个被控对象动态参数进行计算，得出十三个被控对象动态参数为：

(1)    η₃＞η₁＞1，η₃≈1，η₁≈1

(2)    T₂＝0.75×T₁

$\left(3\right)K_3=\frac{3e}{16}\×\frac{P_{10}}{P_{T0}-(K_1+K_2)}\×K_1\×T_1$

$\left(4\right)T_4=\frac{9}{4}\×T_1$

(5)    K₄＝K₁+K₂

$\left(6\right)C_K=\frac{9}{8}\×\frac{K_3}{K_4}$

(7)    T₄＞T₂；T₄≈T₂

$\left(8\right)K_5=\frac{T_4}{T_2}\×K_3\×\frac{N_{T0}}{P_{10}}$

(9)    a＝PT₀+(K₁+K₂)

(10)   b＝P₁₀

(11)   c＝P_d0+K₄

(12)   d＝NT₀

$\left(13\right){\mathrm{\η}}_2=\frac{1}{\sqrt[3]{\frac{24{\mathrm{\η}}_1^{2}T}{C_{k}c}}-1}$

其中：μ_T为汽机调门开度，μ_B为锅炉燃烧率指令；

所述的七个物理约束条件是指：

条件一：μ_B阶跃扰动下，P₁/P_T理论上代表了汽机调门的开度，从动态扰动的角度看，它应不随时间变化；

条件二：因为μ_T＝P₁/P_T，所以在μ_T阶跃扰动下，P₁/P_T同比阶跃；

条件三：在μT扰动下，达到稳态时，因蒸汽流量不变，故有ΔP_d(∞)＝ΔP_T(∞)；

条件四：在μ_T扰动下，热量

保持不变；

条件五：在μ_b扰动下，热量信号近似为阶跃响应；

条件六：当μ_B和μ_T分别从0增至100％时，即：

条件七、因μ_T扰动时，P₁与N_T的响应特性几乎相同，故有：

$\frac{{\mathrm{\ΔP}}_{1\max }}{P_{10}}=\frac{{\mathrm{\ΔN}}_{T\max }}{N_{T0}}$

$\stackrel{\·}{\·\·}\frac{K_3}{K_5}\frac{T_4}{T_2}=\frac{P_{10}}{N_{T0}}\⇒K_5=\frac{T_4}{T_2}{K}_3\frac{N_{T0}}{P_{10}}$

其中：K₁-K₅、a、b、c和d为九个比例增益系数，T、T₁-T₄为五个时间常数，η₁-η₃为三个时间常数的系数，P_T0和P₁₀为基准工况下的运行参数，N_T为机组功率，P_T为汽机机前压力，P₁为汽机的一级压力，P_d为汽包压力。

Images