本发明公开了一种基于约束网格层次聚类的空分过程PWA模型辨识方法。该方法以精馏塔为研究对象，包括特征向量的估计、特征向量的聚类与区域划分的联合估计、以及PWA子模型的参数估计。本发明通过将回归域划分成网格，并且在层次聚类的过程中考虑实际回归域划分的约束条件，完成了数据聚类与区域划分的联立求解。相对于已有技术，本发明减轻了奇异点可能导致的回归域不可划分问题，且避免了复杂的切面方程参数求解工作，同时可应用于非线性划分的回归域，在较高噪信比时也能保持较高的辨识精度。
1.一种基于约束网格层次聚类的空分过程PWA模型辨识方法，其特征在于，该方法包括以下步骤：
步骤(1)：获取空分过程的辨识数据，包括：每个时刻精馏塔的空气进入量及空分结束后氧气、氮气和氩气各产品组分，根据PWA模型的要求，将空分过程的辨识数据的处理成数据点，再根据空分过程与PWA对象的特点和数据点的分布，设定局部数据集的大小，对每个数据点选择邻近的几个点组成该数据点的局部数据集，对局部数据集通过最小二乘方法估计局部模型作为特征向量，并且估计局部模型的协方差矩阵，将协方差矩阵作为该模型在后续计算中的权重；
步骤(2)：设定PWA模型回归域划分的精细度作为网格的数量，将步骤(1)得到的特征向量划分到各个网格中，将每个网格视为一个类，根据网格之间的相邻关系初始化距离矩阵，初始化距离矩阵后对每个网格进行聚类操作，聚类时需要满足实际回归域的完备和不重叠划分约束；聚类操作主要包括合并聚类和分裂聚类，执行合并聚类和分裂聚类的标准是比较生成随机数与给定值的大小，包括以下子步骤：(2.1)生成随机数，所述随机数在(0,1)之间服从均匀分布；
(2.2)若随机数的大小和聚类个数符合合并条件，执行聚类合并操作，更新聚类中包含的网格并更新距离矩阵；若随机数的大小和聚类个数符合分裂条件，执行聚类分裂操作，更新聚类并更新距离矩阵；
(2.3)判断聚类个数是否符合目标聚类个数，如不符合返回(2.1)，如符合执行(2.4)；
(2.4)不断执行分裂再合并操作调整网格划分的边界，直到达到收敛条件或达到循环次数上限，完成聚类；
步骤(3)：根据特征向量的聚类结果划分原始数据，对每个类中的数据用最小二乘方法估计ARX子模型；由子模型与区域的对应关系，得到空分过程的完整PWA模型。
2.如权利要求1所述的一种基于约束网格层次聚类的空分过程PWA模型辨识方法，其特征在于，所述步骤(1)具体包括以下子步骤：(1.1)采集空分过程的辨识数据，包括：每个时刻的精馏塔的空气进入量及空分结束后氧气、氮气和氩气各产品组分并预处理组成数据点(x(k),y(k)),k＝1,...,N，N为数据点数T
量，其中，调度变量x(k)＝[y(k‑1) u(k‑1)] ，u(k‑1)为系统在时刻k‑1的空气量，y(k‑1)为系统在时刻k‑1的空分结束获得的氧气量、氮气量和氩气量，根据PWARX模型表示为与x(k)的函数关系；
(1.2)选择周围的c‑1个数据点 作为数据点(x(k),y(k))的局部数据集Lk，局部数据集中的点满足：
L为局部数据集的集合，L\Lk表示数据集Lk以外其他数据集；
(1.3)对每个局部数据集中的数据通过最小二乘方法估计局部模型作为特征向量ξk：其中，回归向量矩阵
且x1,x2,...xc∈Lk， 为局部数据集Lk中包含的调度变量与输出向量，并计算协方差矩阵：
3.如权利要求1所述的一种基于约束网格层次聚类的空分过程PWA模型辨识方法，其特征在于，所述步骤(2)中，根据网格之间的相邻关系初始化距离矩阵的具体方法是：设置网格在y(k‑1)和u(k‑1)两个维度上的精度为m，将步骤(1)得到的每个特征向量按照原始数据的范围划分到每个网格中，将每个网格视作一个类，并初始化距离矩阵：不相邻的两个聚类之间的距离为无穷大，相邻聚类之间的距离用加权Ward距离衡量；首先计算每个聚类的平方和误差：
其中，Vk表示协方差矩阵，加权Ward距离为将两个相邻的聚类Ci、Cj合并成聚类Ct造成的信息损失值，即相邻类之间的距离通过下式计算：ΔESSi,j＝ESS(Ct)‑ESS(Ci)‑ESS(Cj)其中Ct为类Ci和Cj合并之后的类，且有ΔESSi,j＝ΔESSj,i，得到初始化的距离矩阵。
4.如权利要求1所述的一种基于约束网格层次聚类的空分过程PWA模型辨识方法，其特征在于，步骤(2.2)中，合并的条件为：生成的随机数在(0,p]区间内，分裂的条件为：生成的随机数在(p,1)区间内且聚类集合中包含具有多个网格的类，其中p为合并概率且0.5＜p＜
5.如权利要求1所述的一种基于约束网格层次聚类的空分过程PWA模型辨识方法，其特征在于，步骤(2.2)中每次合并聚类时，选择距离矩阵中加权Ward距离最小的两个类进行合并，合并时采用图的节点合并的方法更新合并后聚类的对外相邻关系，同时更新相对应的类和距离矩阵。
6.如权利要求1所述的一种基于约束网格层次聚类的空分过程PWA模型辨识方法，其特征在于，步骤(2.2)中分裂聚类时，首先筛选出包含多个网格的聚类，计算每个类的平方和误差，将平方和误差大小作为对该类执行分裂操作的概率，选择一个类分裂：在该聚类的边缘网格中，计算将该网格从类中分裂出来时引起的平方和误差的变化值，将该变化值作为选择该边缘网格执行分裂操作的概率，分裂出一个网格；执行分裂操作后，增加距离矩阵的维数、更新被分裂的类、新的类与相关类的对外连接关系和距离大小。
7.如权利要求1所述的一种基于约束网格层次聚类的空分过程PWA模型辨识方法，其特征在于步骤(3)中，根据约束网格层次聚类的聚类结果即得到的特征向量所在的类由特征向量和原始数据之间的双向映射关系，原始数据可按照如下方式同样被分到s个类中：
如果 则(x(k),y(k))∈Fi并且对每个原始数据类Fi,i＝1,...,s中的数据通过最小二乘方法估计PWA模型的子模型参数：
其中 和 代表属于数据类Fi的调度变量和输出向量；完成空分过程PWA模型辨识。
基于约束网格层次聚类的空分过程PWA模型辨识方法
技术领域
[0001] 本发明涉及空分过程的混杂系统辨识领域，特别是涉及一种基于约束网格层次聚类的分段仿射(PieceWise Affine，简称PWA)模型的辨识问题，其特点是在聚类过程中引入
区域划分的约束条件从而得到数据分类与区域划分的联立解，求解精度较高且求解速度较
快。
背景技术
[0002] 空分即空气分离，是通过分离空气以制取氧、氮和稀有气体等的过程，常用的方法包括化学法、电解法、吸附法、膜分离法、低温精馏分离法等，其中低温精馏分离法具有生产
成本低、技术成熟的优点，适合大规模的工业生产，因此称为工业上制取氧气等气体的重要
方法。如何对空分过程建立准确的模型，便于进行过程控制从而提高产量和质量，是一个具
有研究价值的问题。
[0003] 空分过程具有复杂的动态特性，难以用线性模型描述。混杂模型(Hybrid System，简称HS)是离散事件动态系统与连续变量动态系统统一而形成的，综合了由连续变量描述
的系统机理和以离散变量描述的逻辑约束，因此适用于空分过程的建模。
[0004] 分段仿射(PieceWise Affine，简称PWA)模型是混杂系统的一类重要子类，常用于非线性系统的建模，能以任意精度拟合非线性系统。该模型将系统的状态空间分割成若干
个多面体区域，不同的区域由不同的线性子模型描述其特性，并且系统状态在各区域的交
界处保持连续。经研究表明，其他混杂模型经过简单的变换也可以转换为PWA模型的形式，
因此PWA模型具有普遍适用性。相比于非线性模型，PWA模型能够简化控制器的设计和分析
问题，已经有学者实现了针对PWA模型的鲁棒控制、预测控制等控制方法；并且由于其在各
个子区域上的线性特性，可利用线性理论中的成熟理论应用于非线性系统的分析。基于上
述优点，PWA模型近年来引起了学者们的广泛关注，并成功应用于航空航天、化工生产、机械
制造等领域。
[0005] PWA模型的辨识主要包括两个方面的内容：对回归向量集的完备、不重叠分割和对子模型参数的估计。目前对PWA模型的辨识方法主要分为三大类，第一类是基于聚类的方
法，如基于K均值聚类的方法、基于单链接(Single Link)层次聚类的方法等；第二类是基于
优化的方法，如期望值最大法、l1优化方法等；第三类为一些其他辨识算法，如有界误差法
等。其中，聚类方法在在易用性(或易初始化)与建模精度之间达到了较好的平衡，因此聚类
方法在许多研究中都得到了应用。基于聚类的PWA模型辨识方法主要包括以下几个步骤：
[0006] (1)对每个数据点建立局部数据集，通过最小二乘方法估计局部模型作为特征向量，并估计局部模型的协方差矩阵；
[0007] (2)通过聚类算法对特征向量聚类；
[0008] (3)根据特征向量和原始数据之间的双射关系，原始数据也被分到几个类中，对每个类中的数据估计PWA模型的子模型
[0009] (4)根据原始数据的分布对回归空间做区域分割并估计切面参数，实现对状态空间的完整、不重叠划分。
[0010] 通过上述步骤完成传统聚类算法对PWA模型的辨识。但是由于在辨识算法中区域的划分与数据的聚类是分开进行的，因此聚类过程中产生的奇异点容易造成区域的不可划
分问题，从而严重影响最终的辨识结果；此外，回归域的划分通常是对回归向量构成的高维
空间的划分，切面参数的估计是比较复杂的，因此如何在聚类过程中的引入对回归域的完
整、不重叠划分，从而提高辨识的精度，是一个需要解决的问题。
发明内容
[0011] 基于上面所述的一些研究热点和问题，本发明提出一种基于约束网格层次聚类的空分过程PWA模型辨识方法，将回归域的完整、不重叠划分约束引入聚类优化过程，得到数
据聚类与区域划分的联立解，改进了基于聚类的传统PWA模型辨识方法中存在的一些问题，
使得空分装置的过程建模具有更高的精确性和快速性。
[0012] 本发明通过以下技术方案予以实现：
[0013] (1)获取空分过程的辨识数据，包括：每个时刻精馏塔的空气进入量及空分结束后氧气、氮气、氩气等各产品组分，根据PWA模型的要求，将空分过程的辨识数据的处理成数据
点，再根据PWA模型的特点和数据点的分布，设定局部数据集的大小，在每个数据点附近选
择距离最近的几个点组成局部数据集，对局部数据集中的数据通过最小二乘方法估计局部
模型作为特征向量，以更精确的表示该数据点的动态特性；并且估计局部模型的协方差矩
阵，将协方差矩阵作为该模型在后续计算中的权重，提高辨识的精度。上述参数的具体计算
方法将在具体实施方式中予以详细说明。
[0014] (2)设定PWA模型回归域划分的精细度即网格的数量，将步骤(1)得到的局部模型划分到各个网格中，将每个网格视为一个类并初始化距离矩阵，相邻聚类之间的距离使用
加权Ward方法衡量，不相邻的聚类之间距离为无穷大。初始化距离矩阵后对每个网格进行
聚类操作，聚类操作主要包括合并聚类、分裂聚类，执行合并聚类或分裂聚类的标准是比较
生成的随机数与给定值的大小，具体执行的方法为：
[0015] (2.1)生成随机数，随机数服在(0,1)之间服从均匀分布；
[0016] (2.2)若随机数的大小和聚类个数符合合并条件，执行聚类合并操作，更新聚类并更新距离矩阵；若随机数的大小和聚类个数符合分裂条件，执行聚类分裂操作，更新聚类并
更新距离矩阵；
[0017] (2.3)判断聚类个数是否符合目标聚类个数，如不符合返回(2.1)，如符合执行(2.4)；
[0018] (2.4)不断执行分裂再合并操作调整网格划分的边界，检测聚类是否达到收敛条件，或检测当前的迭代次数，若结果收敛或已达到最高迭代次数，则退出聚类过程。由于完
整的回归域被划分成不重叠的网格，完成对网格的聚类过程就自然而然地完成了对整个回
归域的划分。得到最终聚类结果
[0019] (3)由于特征向量和原始数据之间存在双向映射关系，根据特征向量的聚类结果将原始数据划分到对应的类中，对每个类中的数据用最小二乘方法估计ARX模型，由此得到
各个子区域上的模型，并得到完整回归域上的PWA模型。
[0020] 进一步地，所述步骤(1)具体包括以下子步骤：
[0021] (1.1)采集空分过程的辨识数据，包括：每个时刻精馏塔的空气进入量及空分结束后氧气、氮气、氩气等各产品组分，并预处理组成数据点(x(k),y(k)),k＝1,...,N，其中，调
T
度变量x(k)＝[y(k‑1) u(k‑1)]，u(k‑1)为系统在时刻k‑1的空气量，y(k‑1)为系统在时刻
k‑1的空分结束获得的氧气量、氮气量、氩气量等，根据PWARX模型表示为与x(k)的函数关
系。(1.2)选择周围的c‑1个数据点 作为数据点(x(k),y(k))的局部数据集Lk，局部数据
集中的点满足：
[0022]
[0023] L为局部数据集的集合，L\Lk表示数据集Lk以外其他数据集。
[0024] (1.3)对每个局部数据集中的数据通过最小二乘方法估计局部模型作为特征向量ξk：
[0025]
[0026] 其中，回归向量矩阵
[0027]
[0028] 且x1,x2,...xc∈Lk， 为局部数据集Lk中包含的调度变量与输出向量，并计算协方差矩阵：
[0029]
[0030] I为c×c维的单位矩阵。
[0031] 进一步地，所述步骤(2)中，根据网格之间的相邻关系初始化距离矩阵的具体方法是：
[0032] 设置网格在y(k‑1)和u(k‑1)两个维度上的精度为m，将步骤(1)得到的每个特征向量按照原始数据的范围划分到每个网格中，将每个网格视作一个类，并初始化距离矩阵：不
相邻的两个类之间的距离为无穷大，相邻聚类之间的距离用加权Ward距离衡量。首先计算
每个聚类的平方和误差：
[0033]
[0034] 其中，μi表示聚类Ci的均值，即
[0035]
[0036] 加权Ward距离为将两个相邻的聚类Ci、Cj合并成聚类Ct造成的信息损失值，即相邻类之间的距离通过下式计算：
[0037] ΔESSi,j＝ESS(Ct)‑ESS(Ci)‑ESS(Cj)
[0038] 其中Ct为类Ci和Cj合并之后的类，且有ΔESSi,j＝ΔESSj,i，得到初始化的距离矩阵。
[0039] 进一步地，步骤(2.2)中，合并的条件为：生成的随机数在(0,p]区间内，分裂的条件为：生成的随机数在(p,1)区间内且包含多个网格的聚类，其中p为概率且0.5＜p＜1。
[0040] 进一步地，步骤(2.2)中每次合并聚类时，选择距离矩阵中加权Ward距离最小的两个类进行合并，合并时采用图的节点合并的方法更新合并后聚类的对外相邻关系，同时更
新相对应的类和距离矩阵，保证类和类之间的相邻关系和距离关系的正确性。合并的概率
是p。
[0041] 进一步地，步骤(2.2)中分裂聚类时，首先筛选出包含多个网格的聚类，计算每个类的平方和误差，将平方和误差大小作为对该类执行分裂操作的概率，选择一个类分裂：在
该聚类的边缘网格中，计算将该网格从类中分裂出来时引起的平方和误差的变化值，将该
变化值作为选择该边缘网格执行分裂操作的概率，分裂出一个网格。执行分裂操作后，增加
距离矩阵的维数、更新新的类与相关类的对外连接关系和距离大小，保证类和类之间的相
邻关系和距离关系的正确性。分裂的概率为1‑p。
[0042] 进一步地，步骤(3)中，根据约束网格层次聚类的聚类结果即得到的特征向量所在的类 由特征向量和原始数据之间的双向映射关系，原始数据可按照如下方式
同样被分到s个类中：
[0043]
[0044] 并且对每个原始数据类Fi,i＝1,...,s中的数据通过最小二乘方法估计PWA模型的子模型参数：
[0045]
[0046] 其中 和 代表属于数据类Fi的调度变量和输出向量。完成空分过程PWA模型辨识。
[0047] 本发明的有益效果是：
[0048] 本发明以空分过程为研究对象，通过PWA模型描述空分过程的动态特性，并提出一种基于约束网格层次聚类的算法完成模型的估计。算法将调度变量的回归域划分成任意精
度的若干个网格，利用改进的层次聚类算法对数据进行分类。在聚类的过程中考虑回归域
上的完整、不重叠划分，从而完成数据聚类与区域划分的联合估计，减轻了奇异点可能引起
的区域不可划分问题从而对结果产生的严重影响，提高了辨识的精确度；同时避免了对区
域切面方程的进行复杂的多参数估计，简化了区域划分问题，提高了辨识的速度。该方法不
仅可以应用于线性划分的区域估计，也可应用于非线性划分的区域估计问题，在噪信比较
高时仍然能保持较好的辨识准确度，对于空分过程等工业生产过程能得到较好的建模结
果。
附图说明
[0049] 图1是基于约束网格层次聚类的二维调度变量的PWA模型辨识算法的流程框图；
[0050] 图2是一个二维调度变量的PWARX模型的示例；
[0051] 图3是一个3×3网格划分的无向图表示；
[0052] 图4是无向图的邻接矩阵表示形式；
[0053] 图5是聚类合并后的距离矩阵更新结果；
[0054] 图6是聚类分裂后的距离矩阵更新结果；
[0055] 图7为仿真PWA模型的区域划分
[0056] 图8为基于约束网格层次聚类的PWA模型辨识算法得到的区域划分结果；
[0057] 图9为传统基于聚类的PWA模型辨识算法得到的区域划分结果。
具体实施方式
[0058] 如图1所示，本发明采用一种基于约束网格层次聚类的空分过程PWA辨识算法对空分过程进行辨识，实施步骤如下：
[0059] (1)构建局部数据集和计算特征向量
[0060] 获取空分过程的辨识数据，包括：每个时刻精馏塔的空气进入量及空分结束获得的氧气、氮气、氩气等各产品组分。给定PWA模型完整回归域的上界和下界，给定局部数据集
的初始值大小为c，图2为一个二维调度变量的PWARX模型的示例，对于更高维的调度变量，
可以通过增加矩阵维度的方式进行辨识，辨识过程不失一般性，考虑如下PWARX模型的一般
表示形式：
[0061]
[0062] 其中x(k)＝[y(k‑1) u(k‑1)]T为调度变量，u(k‑1)为系统在时刻k‑1的输入，本实施例中对应于时刻k‑1的的空气进入量，y(k‑1)为系统在时刻k‑1的输出，对应于时刻k‑1空
分结束获得的氧气量、氮气量、氩气量等，ε(k)为均值为0的高斯噪声。子模型个数为s，θi为
子模型参数：
[0063]
[0064] n为模型的阶次，对于二维调度变量的PWA模型有n＝1，consti为PWA模型中第i个ARX子模型的常数项系数， 为第i个ARX子模型中y(k‑j)项的系数，
为第i个ARX子模型中u(k‑j)的系数，其中k＞j。通过工业过程采集到N个数据点，对每个数
据点(x(k),y(k)),k＝1,...,N选择其周围的c‑1个数据点 作为局部数据集Lk，选择的
数据点满足：
[0065]
[0066] L为局部数据集的集合，上式代表对数据点(x(k),y(k))，数据集Lk中的数据点到当前数据点的距离小于数据集Lk以外其他数据集L\Lk中的数据点 到当前点的
距离。对每个局部数据集中的数据通过最小二乘方法估计局部模型作为特征向量ξk，即：
[0067]
[0068] 其中回归向量矩阵
[0069]
[0070] 且x1,x2,...xc∈Lk， 为局部数据集Lk中包含的调度变量与输出向量，并按照下式计算协方差矩阵：
[0071]
[0072] I为c×c维的单位矩阵。
[0073] (2)执行约束网格层次聚类算法
[0074] 根据图2可知，对PWA模型的区域划分实际上是对x‑y平面的划分。设置网格在y(k‑1)和u(k‑1)两个维度上的精度为m，将步骤(1)得到的每个特征向量按照原始数据的范围划
分到每个网格中，将每个网格视作一个类，并初始化距离矩阵。不相邻的两个类之间的距离
为无穷大，相邻聚类之间的距离用加权Ward距离衡量。首先计算每个聚类的平方和误差：
[0075]
[0076] 其中，μi表示聚类Ci的均值，即
[0077]
[0078] 加权Ward距离为将两个相邻的聚类Ci、Cj合并成聚类Ct造成的信息损失值，即相邻类之间的距离通过下式计算：
[0079] ΔESSi,j＝ESS(Ct)‑ESS(Ci)‑ESS(Cj) (8)
[0080] 其中Ct为类Ci和Cj合并之后的类，且有ΔESSi,j＝ΔESSj,i，得到初始化的距离矩阵如下表所示：
[0081] 表1.初始化距离矩阵
[0082]
[0083] 其中，M＝m×m为网格总数量，inf代表对应的两个类之间的距离无穷大。
[0084] (3)生成随机数
[0085] 调用MATLAB中的rand函数生成一个随机数，生成的随机数在区间(0,1)上服从均匀分布。
[0086] (4)合并或分裂聚类并更新距离矩阵
[0087] 给定聚类合并的概率为p且0.5＜p＜1，若生成的随机数在(0,p]区间内，则执行聚类合并操作，每次合并时搜索距离矩阵中最小的距离对应的两个类Ci、Cj进行合并，将Cj中
的元素添加到Ci中，并且在类中删除类Cj；距离矩阵更新采用类似图的节点合并的方式进行
更新，以简单的3×3网格划分为例，可以用图3中所示的无向图表示类之间的邻接关系，无
向图的邻接关系可用图4中的邻接矩阵表示，在实际的运算过程中，图4矩阵中的1用加权
Ward距离替代。在节点合并时，找到合并的两个节点的对外邻接关系的并集作为新的类的
邻接关系，更新相应的加权Ward距离，并且删除类Cj在距离矩阵中对应的行和列，完成聚类
如图5所示的合并操作下的距离矩阵更新。
[0088] 若生成的随机数在(p,1)区间内，则执行聚类分裂操作。执行分裂操作时，先筛选出所有包含多个网格的类Cindex(i),i＝1,...,W，W为多网格类的数量，每个多网格类满足
size(Cindex(i))＞1，其中index为多网格类的下标集合，因为只有包含多个网格的类才能被
分裂。按照公式(6)计算所有多网格类的平方和误差，平方和误差越大则该类越分散，有较
大被分裂的可能性，对所有多网格类的平方和误差求和，即：
[0089]
[0090] 调用MATLAB中的rand函数，生成在区间(0,ESS)上服从均匀分布的一个随机数r，查找该随机数所在的范围参数q满足
[0091]
[0092] 此时选择类Cindex(q)作为被分裂的类，如图5中所示的分类和距离矩阵，由于只有一个多网格类因此选择第一个类分裂，在实际的运算过程中，图5矩阵中的1用加权Ward距离
替代。
[0093] 选择了被分裂的类后，需要在该类的边缘网格中选择一个被分裂出来的网格，令分裂后的类为C'index(q)，分裂出的网格为新的类C'，计算Ward值：
[0094] ΔESS＝ESS(Cindex(q))‑ESS(C'index(q))‑ESS(C') (11)
[0095] 该Ward值越大，则该网格有越大的可能性被分裂出来，被分裂网格选择的算法与被分裂类选择的方法相同。分裂出的网格更新到距离矩阵的最后一行和最后一列，更新类
C'index(q)和C'的对外邻接关系并计算Ward距离。图5中的聚类里，分裂网格1或网格4具有相
同的Ward值，因此随机选择一个网格分裂出去，将网格4作为新的类添加到距离矩阵的最后
一行和最后一列，如图6所示，更新第一行、第一列与最后一行、最后一列的对外连接关系。
[0096] (6)判断聚类个数与设定值是否相等
[0097] 若当前聚类的数量与设定值相等，则执行下一步，调整聚类的边界网格；若当前聚类的数量与设定值不相等，则生成随机数继续执行约束网格层次聚类算法，一直迭代直到
收敛到给定值。
[0098] (7)调整聚类的边缘网格
[0099] 完成大致的聚类过程后，执行新的循环体。在循环体中，先执行分裂操作，每个类都执行一次分裂操作，再通过合并操作将聚类的个数合并到设定值，若该循环执行了若干
次聚类结果都不发生变化，则认为聚类算法已经收敛，结束当前循环；若该循环的执行次数
已达到最大迭代上限，也结束当前循环，得到最终聚类结果
[0100] (8)估计子模型与搜索参数
[0101] 在上述聚类过程中，网格被划分后回归区域也完成了完整、不重叠划分。由于实际参与聚类的为特征向量，根据特征向量与辨识数据的双向映射关系，辨识数据也按照如下
的方式进行分类：
[0102]
[0103] 对每个类Fi,i＝1,...,s中数据通过最小二乘方法估计PWA模型的各个子参数模型：
[0104]
[0105] 其中 和 代表属于数据类Fi的调度变量和输出向量。由此得到如模型(1)所示的PWA模型。
[0106] 本发明进行仿真的硬件平台是装有Windows10系统的PC机，软件平台为MATLAB R2015b。
[0107] 本发明示例应用于一个二维调度变量的PWA模型，其模型参数为：
[0108]
[0109] x(k)＝[y(k‑1) u(k‑1)]T
[0110] 区域划分如图7所示。其中，噪声ε(k)～N(0,12)，w为工作点，噪信比为
[0111] N/S＝var(y(k)‑y0(k))/var(y0(k)‑w)＝3.26％
[0112] 其中y0(k)为无噪输出，var(·)代表信号方差。产生随机信号u(k)在区间[‑10,10]内均匀分布，并且经过滤波：
[0113]
[0114] 对上述模型产生的数据，模型辨识的局部数据集大小c＝10，网格划分精度为m×m＝20×20。分别采用本发明和基于聚类的传统PWA模型辨识方法对测试用例进行辨识，可得
本发明中对模型的区域划分结果如图8所示，各个子区域上的模型参数如下所示：
[0115]
[0116] x(k)＝[y(k‑1) u(k‑1)]T
[0117] 图9所示为传统基于聚类的PWA模型辨识方法得到的区域划分结果，模型参数的估计结果如下所示：
[0118]
[0119] x(k)＝[y(k‑1) u(k‑1)]T
[0120] 根据上述模型辨识结果，改进的方法具有更好的辨识结果。针对多调度变量的空分过程，辨识算法也能得到能充分描述其动态特性的PWA模型。
[0121] 上述实施例并非是对于本发明的限制，本发明并非仅限于上述实施例，只要符合本发明要求，均属于本发明的保护范围。
