0 引　言

舰船蒸汽动力装置运行时，由于自身环境和各种操作需要，负荷会出现大幅度且频繁的变化^[1]。常需主汽轮机在高转速区运行时突然降速，或从低转速区运行状态陡然加速到高转速，这种陡增陡降的运行方式会给主锅炉带来冲击负荷。舰用主锅炉因船舶机舱容量有限，尺寸和热容量较陆用锅炉小，需设计高效可靠的机炉协调系统来应对主汽轮机带来的负荷变化^[2]。机炉协调控制系统（Coordinated Control System，CCS）是将主锅炉和主汽轮机作为一个整体来控制，包含燃油控制、风量控制、给水控制和蒸汽压力控制等^[3-5]。系统热工动态过程非常复杂，具有强非线性、大时延、多变量耦合、约束多和干扰多的特点^[6-8]。目前，以机理建模为基础，采用PID策略的机炉协调系统无法得到令人满意的控制效果^[9-11]。为改善机炉协调控制质量，张铁军等^[12]以锅炉-汽机协调系统为研究对象，利用模糊集合理论建立了该系统的非线性离散全局模糊模型，设计了一种新型的扩展状态空间预测控制器，有效地改善了整个系统的控制品质。潘晖等^[13]提出一种支持向量机（SVM）的广义预测控制（GPC）建模方法和控制方法，提高了系统的稳定性和鲁棒性。

本文提出采用T-S模糊建模方法建立舰船协调系统模型，并基于该模型利用多模型广义预测控制方法设计合理的控制器^[14]，以实现对协调系统快速稳定的控制，提升蒸汽动力装置整体性能。

1 舰船机炉系统改进T-S模糊模型

以某训练用机组为试验对象，选择燃油量、给水量和主汽门开度作为输入变量，以螺旋桨转速、汽轮机输出功率和主蒸汽压力作为输出变量。每隔1 s采集一次现场DCS数据作为初始样本，共采集4 500组数据，为减小噪声等干扰因素对辨识效果产生的影响，在4500组数据中提取900组用于机组T-S模糊模型辨识，所提取的3个输出变量的900组数据分别如图1～图3所示。

为方便T-S模糊模型的搭建，按照输出变量的不同将实际机组分解成3个多输入单输出的子模型，分别为螺旋桨转速模型、输出功率模型和主蒸汽压力模型。

首先进行输入变量选择，根据启发性知识，选择燃油量、给水量和主汽门开度作为输入变量，并确定实际机组T-S模糊模型前件和后件参数中含有所有的输入变量和输出变量。由此初步确定的3个子模型结构分别为：

螺旋桨转速子模型

$ \begin{gathered} {y_1}(k) = {p_{10}} + {p_{11}}{u_1}(k - 1) + {p_{12}}{u_1}(k - 2) + {p_{13}}{u_2}(k - 1)+ \\ \begin{array}{*{20}{c}} {}&{}&{} \end{array}{p_{14}}{u_{\text{2}}}(k - {\text{2}}) + {p_{15}}{u_{\text{3}}}(k - 1) + {p_{16}}{y_1}(k - 1)。\\ \end{gathered} $

输出功率子模型

$ \begin{gathered} {y_2}(k) = {p_{20}} + {p_{21}}{u_1}(k - 1) + {p_{22}}{u_1}(k - 2) + {p_{23}}{u_2}(k - 1)+ \\ \begin{array}{*{20}{c}} {}&{}&{} \end{array}{p_{24}}{u_{\text{2}}}(k - {\text{2}}) + {p_{25}}{u_{\text{3}}}(k - 1) + {p_{26}}{y_1}(k - 1)。\\ \end{gathered} $

主蒸汽压力子模型

$ \begin{gathered} {y_3}(k) = {p_{30}} + {p_{31}}{u_1}(k - 1) + {p_{32}}{u_2}(k - 1) + {p_{33}}{u_3}(k - 1)+ \\ \begin{array}{*{20}{c}} {}&{}&{} \end{array}{p_{34}}{u_3}(k - 2) + {p_{35}}{y_3}(k - 1) + {p_{36}}{y_3}(k - 2)。\\ \end{gathered} $

当输入变量确定后，进行输入空间模糊划分，采用FCM聚类算法进行输入空间模糊划分，设定聚类数c为5，权值指数m为2，容许误差值 $ \varepsilon $ 取 $ {10^{ - 6}} $ ，所求出的3个子模型聚类中心矩阵如下：

螺旋桨转速子模型

$ {V_1} = \left[ \begin{gathered} {\text{92}}{\text{.53 92}}{\text{.49 448}}{\text{.38 448}}{\text{.36 0}}{\text{.79 296}}{\text{.71}} \\ {\text{67}}{\text{.64 67}}{\text{.58 424}}{\text{.01 423}}{\text{.97 0}}{\text{.74 257}}{\text{.93}} \\ {\text{59}}{\text{.15 59}}{\text{.07 415}}{\text{.99 415}}{\text{.90 0}}{\text{.74 235}}{\text{.83}} \\ {\text{84}}{\text{.35 84}}{\text{.32 441}}{\text{.35 441}}{\text{.30 0}}{\text{.75 287}}{\text{.14}} \\ {\text{76}}{\text{.40 76}}{\text{.36 432}}{\text{.95 432}}{\text{.92 0}}{\text{.75 275}}{\text{.09}} \\ \end{gathered} \right]。$

输出功率子模型

$ {V_2} = \left[ \begin{gathered} {\text{92}}{\text{.84 92}}{\text{.81 448}}{\text{.62 448}}{\text{.61 0}}{\text{.79 1264}}{\text{.62}} \\ {\text{69}}{\text{.50 69}}{\text{.44 425}}{\text{.79 425}}{\text{.74 0}}{\text{.74 834}}{\text{.02}} \\ {\text{77}}{\text{.83 77}}{\text{.79 434}}{\text{.40 434}}{\text{.35 0}}{\text{.75 990}}{\text{.73}} \\ {\text{85}}{\text{.19 85}}{\text{.16 442}}{\text{.26 442}}{\text{.22 0}}{\text{.75 1124}}{\text{.70}} \\ {\text{60}}{\text{.11 60}}{\text{.04 416}}{\text{.99 416}}{\text{.90 0}}{\text{.74 657}}{\text{.68}} \\ \end{gathered} \right] 。$

主蒸汽压力子模型

$ {V_3} = \left[ \begin{gathered} {\text{77}}{\text{.54 434}}{\text{.17 0}}{\text{.75 0}}{\text{.75 2}}{\text{.16 2}}{\text{.16}} \\ {\text{68}}{\text{.95 425}}{\text{.24 0}}{\text{.74 0}}{\text{.74 1}}{\text{.84 1}}{\text{.84}} \\ {\text{59}}{\text{.90 416}}{\text{.77 0}}{\text{.74 0}}{\text{.74 1}}{\text{.49 1}}{\text{.49}} \\ {\text{85}}{\text{.03 441}}{\text{.96 0}}{\text{.75 0}}{\text{.75 2}}{\text{.45 2}}{\text{.45}} \\ {\text{92}}{\text{.77 448}}{\text{.60 0}}{\text{.79 0}}{\text{.79 2}}{\text{.58 2}}{\text{.58}} \\ \end{gathered} \right]。$

然后进行参数辨识，先进行前件参数辨识，选择高斯型隶属度函数作为辨识函数，其表达式为：

$ {\mu _A}(x) = {e^{ - {{(\frac{{x - a}}{{\sqrt 2 b}})}^2}}}。$

(1)

式中： $ a $ 为数据集合中心； $ b $ 为数据标准差； $ {\mu _A}(x) $ 为样本 $ x $ 在子集A中的隶属度。

再进行后件参数辨识，采用递推最小二乘法计算出结论参数，设定规则数为5的情况下所计算出的3个子模型的结论参数矩阵如下：

寻优后螺旋桨转速子模型

$ \begin{gathered} {y_{11}}(k) = -67.80 -0.03{u_1}(k - 1) + {\text{0.09}}{u_1}(k - 2) + {\text{0.02}}{u_2}(k - 1) + \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {\text{0}}{\text{.05}}{u_2}(k - 2)-{\text{ 27}}{\text{.73}}{u_3}(k - 1){\text{ + 1}}{\text{.17}}{y_1}(k - 1), \\ \end{gathered} $

$ \begin{gathered} {y_{12}}(k) = {{67}}{{.47 + 0}}{{.43}}{u_1}(k - 1){{ - 0}}{{.08}}{u_1}(k - 2){{ - 0}}{{.07}}{u_2}(k - 1)+ \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{ 0}}{{.15}}{u_2}(k - 2){{ - 147}}{{.02}}{u_3}(k - 1){{ + 0}}{{.94}}{y_1}(k - 1), \\ \end{gathered} $

$ \begin{gathered} {y_{13}}(k) = {{53}}{{.79 + 0}}{{.04}}{u_1}(k - 1){{ - 0}}{{.01}}{u_1}(k - 2){{ - 0}}{{.01}}{u_2}(k - 1)+ \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{ 0}}{{.04}}{u_2}(k - 2){{ - 84}}{{.64}}{u_3}(k - 1){{ + 0}}{{.96}}{y_1}(k - 1) , \\ \end{gathered} $

$ \begin{gathered} {y_{14}}(k) = {{44}}{{.85 - 0}}{{.03}}{u_1}(k - 1){{ + 0}}{{.05}}{u_1}(k - 2){{ - 0}}{{.04}}{u_2}(k - 1)+ \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{ 0}}{{.06}}{u_2}(k - 2){{ - 69}}{{.84}}{u_3}(k - 1){{ + 1}}{{.00}}{y_1}(k - 1) , \\ \end{gathered} $

$ \begin{gathered} {y_{15}}(k) = {{70}}{{.44 + 0}}{{.09}}{u_1}(k - 1){{ - 0}}{{.05}}{u_1}(k - 2){{ + 0}}{{.07}}{u_2}(k - 1)+ \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{ - 0}}{{.11}}{u_2}(k - 2){{ - 100}}{{.50}}{u_3}(k - 1){{ + 1}}{{.07}}{y_1}(k - 1)。\\ \end{gathered} $

寻优后输出功率子模型

$ \begin{gathered} {y_{21}}(k) = {{ - 517}}{{.55 - 0}}{{.52}}{u_1}(k - 1){{ + 1}}{{.63}}{u_1}(k - 2){{ + 0}}{{.32}}{u_2}(k - 1)+ \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{ 0}}{{.93}}{u_2}(k - 2){{ - 327}}{{.61}}{u_3}(k - 1){{ + 1}}{{.09}}{y_2}(k - 1), \\ \end{gathered} $

$ \begin{gathered} {y_{22}}(k) = {{3092}}{{.02 + 3}}{{.06}}{u_1}(k - 1){{ - 0}}{{.81}}{u_1}(k - 2){{ - 3}}{{.66}}{u_2}(k - 1)- \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} { {1}}{{.48}}{u_2}(k - 2){{ - 1537}}{{.18}}{u_3}(k - 1){{ + 1}}{{.10}}{y_2}(k - 1) , \\ \end{gathered} $

$ \begin{gathered} {y_{23}}(k) = {{784}}{{.81 + 0}}{{.43}}{u_1}(k - 1){{ - 0}}{{.40}}{u_1}(k - 2){{ + 0}}{{.79}}{u_2}(k - 1) - \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{ 0}}{{.97}}{u_2}(k - 2){{ - 997}}{{.39}}{u_3}(k - 1){{ + 1}}{{.04}}{y_2}(k - 1) ，\\ \end{gathered} $

$ \begin{gathered} {y_{24}}(k) = {{616}}{{.94 - 0}}{{.29}}{u_1}(k - 1){{ + 0}}{{.78}}{u_1}(k - 2){{ - 0}}{{.66}}{u_2}(k - 1)+ \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{ 0}}{{.83}}{u_2}(k - 2){{ - 956}}{{.95}}{u_3}(k - 1){{ + 0}}{{.99}}{y_2}(k - 1)，\\ \end{gathered} $

$ \begin{gathered} {y_{25}}(k) = {{359}}{{.20 + 0}}{{.31}}{u_1}(k - 1){{ - 0}}{{.13}}{u_1}(k - 2){{ - 0}}{{.14}}{u_2}(k - 1)+ \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{ 0}}{{.29}}{u_2}(k - 2){{ - 564}}{{.07}}{u_3}(k - 1){{ + 0}}{{.98}}{y_2}(k - 1) 。\\ \end{gathered} $

寻优后主蒸汽压力子模型

$ \begin{gathered} {y_{31}}(k) = {{ - 0}}{{.69 + 0}}{{.01}}{u_1}(k - 1) + {{0}}{{.01}}{u_2}(k - 1){{ - 0}}{{.76}}{u_3}(k - 1)+ \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{0}}{{.71}}{u_3}(k - 2) + {{1}}{{.37}}{y_3}(k - 1){{ - 0}}{{.41}}{y_3}(k - 2) ，\\ \end{gathered} $

$ \begin{gathered} {y_{32}}(k) = {{ - 0}}{{.10 + 0}}{{.01}}{u_1}(k - 1){{ + 0}}{{.01}}{u_2}(k - 1){{ - 0}}{{.53}}{u_3}(k - 1)+ \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{0}}{{.57}}{u_3}(k - 2) + {{1}}{{.56}}{y_3}(k - 1){{ - 0}}{{.58}}{y_3}(k - 2) ，\\ \end{gathered} $

$ \begin{gathered} {y_{33}}(k) = {{0}}{{.01 + 0}}{{.01}}{u_1}(k - 1) + {{0}}{{.01}}{u_2}(k - 1){{ - 0}}{{.74}}{u_3}(k - 1) + \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{0}}{{.70}}{u_3}(k - 2) + {{1}}{{.40}}{y_3}(k - 1){{ - 0}}{{.41}}{y_3}(k - 2)，\\ \end{gathered} $

$ \begin{gathered} {y_{34}}(k) = {{ - 0}}{{.73 - 0}}{{.01}}{u_1}(k - 1){{ + 0}}{{.01}}{u_2}(k - 1){{ - 0}}{{.45}}{u_3}(k - 1) + \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{1}}{{.04}}{u_3}(k - 2) + {{1}}{{.89}}{y_3}(k - 1){{ - 0}}{{.92}}{y_3}(k - 2)，\\ \end{gathered} $

$ \begin{gathered} {y_{35}}(k) = {{0}}{{.21 + 0}}{{.01}}{u_1}(k - 1){{ + 0}}{{.01}}{u_2}(k - 1){{ - 0}}{{.77}}{u_3}(k - 1) + \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{} \end{array} {{0}}{{.61}}{u_3}(k - 2) + {{1}}{{.36}}{y_3}(k - 1){{ - 0}}{{.45}}{y_3}(k - 2)。\\ \end{gathered} $

2 多变量广义预测控制 2.1 多变量广义预测控制基本算法

多变量广义预测控制算法具有与单输入单输出预测控制算法一致的控制思想和设计过程，主要应用于多输入多输出系统中。以 $ n $ 输入 $ n $ 输出系统为例，所采用的CARIMA数学模型可用离散差分方程表示为：

$ A({z^{ - 1}})\Delta y(t) = B({z^{ - 1}})\Delta u(t - 1) + C({z^{ - 1}})\xi (t)。$

(2)

其中： $ A({z^{ - 1}}) $ 、 $ B({z^{ - 1}}) $ 和 $ C({z^{ - 1}}) $ 是 $ n \times n $ 的矩阵多项式。具体表示如下：

$ A({z^{ - 1}}) = I + {A_1}{z^{ - 1}} + \cdots + {A_{{n_A}}}{z^{ - {n_A}}} ，$

(3)

$ B({z^{ - 1}}) = {B_0} + {B_1}{z^{ - 1}} + \cdots + {B_{{n_B}}}{z^{ - {n_B}}}，$

(4)

$ C({z^{ - 1}}) = I + {C_1}{z^{ - 1}} + \cdots + {C_{{n_c}}}{z^{ - {n_c}}}。$

(5)

式中： $ u(t) $ 为 $ n \times {\text{1}} $ 维的输入向量； $ y(t) $ 表示 $ n \times {\text{1}} $ 维的输出向量； $ \Delta $ 为 $ n \times n $ 对角差分矩阵，有 $ \Delta = diag\{ 1 - {z^{ - 1}}\} $ ， $ 1 - {z^{ - 1}} $ 为差分算子； $ \xi (t) $ 为均值为0的白噪声向量。

选取以下函数作为性能指标：

$ J = \sum\limits_{j = {N_0}}^N {{{\left\| {y(t + j) - {y_r}(t + j)} \right\|}_Q}^2} + \sum\limits_{j = 1}^{{N_u}} {\left\| {\Delta u(t + j - 1)} \right\|_\lambda ^2}。$

(6)

式中： $ {N_0} $ 为最小预测时域，通常情况下取 $ {N_0} = 1 $ ； $ N $ 为最大预测时域； $ {N_u} $ 为控制时域。 $ \Delta u(t + j) = 0 $ ， $j= {N}_{u}，\cdots N$ ，为在 $ {N_u} $ 步后控制量不再发生变化。 $\left\| x \right\|_{{\lambda}} ^2 = {x^{\rm{T}}}{\boldsymbol{\lambda}} x$ ， $\left\| x \right\|_Q^2 = {x^{\rm{T}}}{\boldsymbol{Q}}x$ ，其中 $ \lambda $ 和 ${\boldsymbol{Q }}$ 是 $ n \times n $ 对角矩阵，都具有控制加权的作用，分别为 ${\boldsymbol{ \lambda }}={\rm{diag}} ({\lambda }_{1}，\cdots ，{\lambda }_{n})$ 和 ${\boldsymbol{Q}}={\rm{diag}}({Q}_{1}，\cdots ，{Q}_{n})$ 。 $ {y_r}(t) $ 是n×1维有界的设定值向量。

在 $ n \times 1 $ 维 $ {y_r}(t) $ 中含有 $ n $ 个设定值，每一个设定值对应一个输出变量，对于绝大多数的工业生产，采用的是恒值控制，所以一般将这 $ n $ 个值设定为常值。令第一个输出变量的设定值表达式为 $ {y_{1r}}(t + j) $ ， $ j=1，2，\cdots $ ，所对应的设定常值为 $ {y_{1r}} $ ，当前时刻的输出为 $ {y_{\text{1}}}(t) $ ，则有如下一阶滤波方程：

$ \left\{ \begin{gathered} {y_{1r}}(t) = {y_{\text{1}}}(t)，\\ {y_{1r}}(t + j) = \alpha {y_{1r}}(t + j - 1) + (1 - \alpha ){y_{1r}} 。\\ \end{gathered} \right. $

(7)

式中： $ \alpha $ 为柔化系数，主要是用来平衡控制中的快速性和鲁棒性。较小 $ \alpha $ 值的控制系统具有更快的跟踪速度，但鲁棒性较差。反之， $ \alpha $ 值较大会使系统具有更好的鲁棒性，但跟踪速度会下降。因此，在选择 $ \alpha $ 值时应当综合考虑。

为方便建立预测模型，引入以下矩阵多项式方程，该矩阵多项式方式是由多个Diophantine方程组成：

$ I = {E_j}({z^{ - 1}})A({z^{ - 1}})\Delta + {z^{ - j}}{F_j}({z^{ - 1}})，$

(8)

$ {E_j}({z^{ - 1}})B({z^{ - 1}}) = {G_j}({z^{ - 1}}) + {z^{ - j}}{H_j}({z^{ - 1}})。$

(9)

式中： $ j=\text{1}，\cdots N $ ，并且有

$ {E_j}({z^{ - 1}}) = {E_0} + {E_1}({z^{ - 1}}) + \cdots + {E_{j - 1}}{z^{ - j + 1}} ，$

(10)

$ {F_j}({z^{ - 1}}) = F_0^j + F_1^j({z^{ - 1}}) + \cdots + F_{{n_A}}^j{z^{ - {n_A}}} ，$

(11)

$ {G_j}({z^{ - 1}}) = {G_0} + {G_1}({z^{ - 1}}) + \cdots + {G_{j - 1}}{z^{ - j + 1}}，$

(12)

$ {H_j}({z^{ - 1}}) = H_0^j + H_1^j({z^{ - 1}}) + \cdots + H_{{n_B} - 1}^j{z^{ - {n_B} + 1}}。$

(13)

联立式(2)、式(7)和式(8)可得：

$ \begin{gathered} y(t + j) = {G_j}({z^{ - 1}})\Delta u(t + j - 1) + {F_j}({z^{ - 1}})y(t)+ \\ \begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{} \end{array}}&{}&{} \end{array} {H_j}({z^{ - 1}})\Delta u(t - 1) 。\\ \end{gathered} $

(14)

将上式用向量形式表示如下：

$ {\boldsymbol{y}} = G{\boldsymbol{u}} + {\boldsymbol{F}}({z^{ - 1}})y(t) + {\boldsymbol{H}}({z^{ - 1}})\Delta u(t - 1) 。$

(15)

其中：

$ {\boldsymbol{y}} = {\left[ {y{{(t + 1)}^{\rm{T}}}, \cdots ,y{{(t + N)}^{\rm{T}}}} \right]^{\rm{T}}}，$

(16)

$ {\boldsymbol{u}} = {\left[ {\Delta u{{(t)}^{\rm{T}}}, \cdots ,\Delta u{{(t + {N_u} - 1)}^{\rm{T}}}} \right]^{\rm{T}}} ，$

(17)

$ {\boldsymbol{F}}({z^{ - 1}}) = {\left[ {{F_1}{{({z^{ - 1}})}^{\rm{T}}}, \cdots ,{F_N}{{({z^{ - 1}})}^{\rm{T}}}} \right]^{\rm{T}}}，$

(18)

$ {\boldsymbol{H}}({z^{ - 1}}) = {\left[ {{H_1}{{({z^{ - 1}})}^{\rm{T}}}, \cdots ,{H_N}{{({z^{ - 1}})}^{\rm{T}}}} \right]^{\rm{T}}} ，$

(19)

$ G = {\left[ {\begin{array}{*{20}{c}} {{G_0}}&0&0&0 \\ {{G_1}}&{{G_0}}&0&0 \\ \vdots & \vdots & \vdots & \vdots \\ {{G_{{N_u} - 1}}}&{{G_{{N_u} - 2}}}& \cdots &{{G_0}} \\ {{G_{N - 1}}}&{{G_{N - 2}}}& \cdots &{{G_{N - {N_u}}}} \end{array}} \right]_{{N_{\text{1}}}n \times {N_u}n}}。$

(20)

将设定参考轨迹定义如下：

$ {{\boldsymbol{y}}_r} = {\left[ {{y_r}{{(t + 1)}^{\rm{T}}}, \cdots ,{y_r}{{(t + N)}^{\rm{T}}}} \right]^{\rm{T}}}。$

(21)

于是可将式(5)改写为：

$ J = {\left\| {{\boldsymbol{y}} - {{\boldsymbol{y}}_r}} \right\|^2} + \left\| {\boldsymbol{u}} \right\|_\Lambda ^2 ，$

(22)

式中， $\Lambda = {\rm{diag}}\{ \lambda \}$ 。

把式(14)代入到式(21)，并使 $ J $ 对 $ u $ 求导，当 $ J $ 取得最小值时有：

$ \begin{split} {\boldsymbol{u}} =& {({G^{\rm{T}}}QG + \Lambda )^{ - 1}}{G^{\rm{T}}}Q[{{\boldsymbol{y}}_r} - {\boldsymbol{F}}({z^{ - 1}})y(t) - \\ & {\boldsymbol{H}}({z^{ - 1}})\Delta u(t - 1)] 。\end{split} $

(23)

把 ${({G^{\rm{T}}}G + \Lambda )^{ - 1}}{G^{\rm{T}}}$ 的前 $ n $ 行写作：

$ {P^{\rm{T}}} = \left[ {{P_1}, \cdots ,{P_N}} \right]。$

(24)

定义：

$ P({z^{ - 1}}) = {P_N} + {P_{N - 1}}{z^{ - 1}} + \cdots + {P_1}{z^{ - N + 1}} ，$

(25)

式中， $ {P_i}(i = 1,2, \cdots ,N) $ 为 $ n \times n $ 矩阵。

根据式(22)～式(24)，可得广义预测控制律为：

$ \begin{split} \Delta u(t) = & {P^{\rm{T}}}[{{\boldsymbol{y}}_r} - {\boldsymbol{F}}({z^{ - 1}})y(t) - {\boldsymbol{H}}({z^{ - 1}})\Delta u(t - 1)]=\\ & P({z^{ - 1}}){{\boldsymbol{y}}_r}(t + N) - \alpha ({z^{ - 1}})y(t)- \\ &\beta ({z^{ - 1}})\Delta u(t - 1)。\end{split} $

(26)

其中 $ \alpha ({z^{ - 1}}) $ 和 $ \beta ({z^{ - 1}}) $ 为参数矩阵多项式，分别为：

$ \alpha ({z^{ - 1}}) = \sum\limits_{j = 1}^N {{P_j}{F_j}({z^{ - 1}})} ，$

(27)

$ \beta ({z^{ - 1}}) = \sum\limits_{j = 1}^N {{P_j}{H_j}({z^{ - 1}})}。$

(28)

2.2 多模型预测控制加权策略

由于所搭建机组模型包含多个局部线性子模型，因此，需将多个子模型进行合理的加权以实现更优的控制。考虑到本文研究对象需在不同工况下切换，为避免出现模型加权得到的整体模型与实际偏差过大的问题，采用基于控制器加权的预测控制方法设计2个系统控制器。

基于控制器加权方法的基本原理是分别对所建立的多个子模型设计对应的控制器，再根据不同时刻各子模型输出与实际输出的接近程度，计算出该模型对应控制器输出占总控制器输出的权重系数。随着时间的变化，权重系数也会不断发生改变，从而使总控制器输出合理，以满足控制的快速性和稳定性。基于控制器加权策略的原理如图4所示。

假定某系统具有 $ n $ 个子模型，且每一个子模型含有 $ m $ 个输出变量，在 $ k $ 时刻的综合加权系数计算步骤如下：

先计算每一个子模型在 $ k $ 时刻的输出误差平方和：

$ {\varepsilon _{ki}} = \sum\limits_{j = 1}^m {{\eta _j}{{({y_{kj}} - {y_{rkj}})}^2}}。$

(29)

由式（28）得到该时刻各子模型的高斯隶属度为：

$ {\mu _{ki}} = \exp ( - \frac{1}{2}{\varepsilon _{ki}}/{\beta ^2})。$

(30)

由此可得到第 $ i $ 个子模型控制器的权重系数为：

$ {v_{ki}} = {\mu _{ki}}/\sum\limits_{i = 1}^n {{\mu _{ki}}}。$

(31)

则 $ k $ 时刻所需要的总控制器输出：

$ {u_k} = \sum\limits_{i = 1}^n {{v_{ki}}} {u_{ki}}。$

(32)

式中： $ {\varepsilon _{ki}} $ 为 $ k $ 时刻第 $ i $ 个子模型的输出误差平方和，其中， $ k = 0,1,2, \cdots $ ， $ i{\text{ = 1,2,}} \cdots ,n $ ； $ {y_{kj}}(j = 1,2, \cdots ,m) $ 为第 $ i $ 个子模型 $ k $ 时刻的第 $ j $ 个变量输出值； $ {y_{rkj}}(j = 1, 2, \cdots ,m) $ 为子模型 $ k $ 时刻的第 $ j $ 个变量设定值； $ {\eta _j}(j = 1,2, \cdots ,m) $ 为第 $ j $ 个输出变量的加权系数； $ {\mu _{ki}} $ 为 $ k $ 时刻第 $ i $ 个子模型的高斯隶属度； $ \beta $ 为常数，用来控制隶属度函数的分辨率，本文取0.3； $ {u_{ki}} $ 为 $ k $ 时刻第 $ i $ 个子模型控制器的输出。

3 仿真结果分析

为了验证基于多模型广义预测控制的舰船机炉协调控制有效性，本文设计单模型广义预测控制进行仿真对比，选择组成的5个CARIMA子模型中的第1个模型进行对照试验，并模拟阶跃扰动过程和变工况过程。表1为预测控制算法设定参数值。

首先进行输出功率阶跃扰动试验，以45%工况点为初始工况点，待机组平稳运行在45%工况点后，在2800 s将功率设定值从590 kW阶跃至992.25 kW。在该试验中，预测控制的输出、输入响应曲线分别如图5和图6所示。

可知，在输出功率发生阶跃扰动时，螺旋桨转速和主蒸汽压力会在短时间内小幅度下降，然后在燃油阀门、给水阀门以及主汽阀门开度不断调节的过程中先升高再下降，最后回到设定值附近。通过计算各条曲线可得，在多模型预测控制下，输出功率从开始发生变化到最终稳定的时长为93 s，而单模型预测控制下的对应时长为115 s，比前者长22 s。同时，前者作用下的转速和主蒸汽压力在调节过程中的峰值分别为232.5 r/min和2.474 MPa，比后者的相应变量峰值低0.8 r/min和0.197 MPa。在3个阀门开度的变化上，前者比后者的开度变化范围分别低0.0018、0.0049和0.0084。

再进行变工况试验，以45%工况点为初始工况点，从45%工况缓慢增加至75%工况，再缓慢下降到45%工况。该试验中，当t=2280 s时，螺旋桨转速设定值先以0.3352 （r/min）/s从229 r/min缓慢升至250.9 r/min，输出功率设定值先以0.01825 kW/s从590 kW升至992.25 kW，主蒸汽压力设定值则以 ${\text{6.92}} \times {\text{1}}{{\text{0}}^{-{\text{4}}}}$ MPa/s从1.34 MPa升至2.17 MPa。待工况稳定，在t=4280 s时，螺旋桨转速设定值以0.3352 （r/min）/s从250.9 r/min缓慢降至229 r/min，输出功率设定值以0.01825 kW/s从992.25 kW降至590 kW，主蒸汽压力设定值以 ${\text{6}}{\text{.92}} \times {\text{1}}{{\text{0}}^{-{\text{4}}} }$ MPa/s从2.17 MPa降至1.34 MPa。图7和图8分别为实际机组在变工况时的预测控制输出、输入响应曲线。

可以看出，多模型预测控制方法能使机组的3个输出变量平稳且快速的跟踪设定值的变化，而单模型预测控制方式在螺旋桨转速和主蒸汽压力的调节中会产生较大的波动。经计算可得，前者作用下的3个输出变量均无超调，而后者调节过程中，转速最高值比设定值高4.6 r/min，最低值比设定值低3.7 r/min，压力最高值比设定值高0.234 MPa，最低值比设定值低0.135 MPa。同时，单模型预测控制方式中的3个阀门开度变化幅度也比多模型控制方法中的对应值更大。

4 结　语

针对舰船机炉协调系统锅炉侧和汽机侧运行特性差异显著、综合调节困难的问题，本文根据某型训练用机组输入输出数据建立了舰船机炉系统T-S模糊模型，并采用多模型广义预测控制算法设计了协调控制器。经验证，与单模型预测控制方法相比，该控制策略提高了响应速度，缩短了调节时间，减少了调节过程中的波动量，使输出更加平稳，有效提升蒸汽动力装置的整体性能。