一般来说，一个系统的LPV模型可以通过雅可比线性化、状态变换法以及函数替代法等理论推导的方法得到[3]。不过，在实际系统中就需要通过系统辨识的方法获得模型参数。现有LPV模型辨识方法分为两类：基于状态空间模型[4]和基于输入输出模型[5]。针对输入输出模型的辨识方法有一般有两类：全局非线性拟合的辨识算法以及局部滑动窗口辨识算法。第1类算法要求选定一组恰当的函数基底，再计算各基底的权值，如非线性最小二乘算法[6]。这类算法对参变函数的近似程度高，但是需要对系统的先验。第2类算法属于非参数化的方法，典型算法如基于支持向量机辨识算法[7]和基于贝叶斯的辨识算法[8]以及递推最小二乘算法等。此类算法在参变函数近似为分段常数函数的情况下，可以准确地辨识得到系统模型。但是，如果当参变函数是线性分段或者类似正弦函数等非线性的情况下，再利用常数近似就会引入较大的函数近似误差。另外，最小二乘算法求解的基础是要求辨识数据个数大于未知数个数。因而，当模型规模较大时，由于模型参数大量增加，这就需要辨识数据也大量增加。一方面，对于采样成本高或者系统响应时间长的情况，很难获得大量的辨识数据。另一方面，LPV模型的参数是一直在变化的，因而增加测量数据就意味着需要拟合更复杂的变化情况。如果近似模型维持原状，就会降低模型的近似程度，也就是模型的近似程度也限制了辨识数据的多少。在这种情况下，最小二乘算法很难同时提高计算精度，又减少因为数据量增多而带来的近似误差。但其实高阶模型中的参数是存在大量零值的，而非零值是十分有限的。因此，在解决高阶模型辨识问题上，利用NNG[9]和LASSO[10]等稀疏估计器是一个有效的办法。本文基于压缩测量辨识(Compress Measurement Identification，CMI)算法提出的动态压缩测量辨识(Dynamic CMI，DCMI)算法是从两个方面解决上述问题。其一，利用“匀速变化”和“非匀速变化”模型，提高对参变函数函数的近似精度，同时避免函数基底的选择。其二，利用压缩感知理论，通过减少所需的辨识数据个数，从而提高算法在辨识数据有限的情况下的辨识精度及参数规模。
1 压缩感知 压缩感知理论是在2006年由Donoho[11]、Candès和陶哲轩[12]提出的一个高效的采样理论。近年来，该理论在信号处理领域中，已经取得了大量令人瞩目的成绩[13-14]。信号的稀疏性是压缩感知的理论基础。对于一个离散信号x，其稀疏度可以定义为信号非零项的个数。所以，信号的稀疏度越小，该信号的稀疏程度就越高，那么理论上所需的采样数也就越少。而采样值yi是由测量向量ri分别与信号x进行内积运算得到的，即yi=〈ri, x〉。考虑多次测量，则测量值y可以写为矩阵形式，即y=Rx。其中，矩阵R被称为测量矩阵。如果R是一个单位矩阵，那么压缩感知采样就退化为了传统采样。在压缩感知理论中，测量值y的长度要远小于原信号x。为了保证这个采样过程是无损的，Candès和Romberg给出了有限等距性质(Restrict Isometric Property，RIP)[15]，即当测量矩阵R满足RIP时，采样过程就是无损的。S阶有限等距常数δS定义为使得测量矩阵R任选几列组成的子矩阵RT, T
 | (1) |
式中：x为任意稀疏信号。当有效等距常数(Restrict Isometry Constants, RIC)满足0 < δS < 1时，该测量矩阵认为是满足RIP的[16]。根据压缩感知理论可知，欠定方程y=Rx解空间中最稀疏的解就是原信号x。由于这是一个组合问题(NP-Hard)，因此该问题需要利用贪婪的思想进行求解。而文献[17]中提出的正交匹配追踪(Orthogonal Matching Pursuit, OMP)算法是其中的代表，该算法在每次迭代中寻找最匹配的向量，作为信号的支撑基，从而得到原信号的最佳近似。另外，Candès和陶哲轩[12]已经证明了，该问题可以转换为在解空间中寻找信号的最小1范数的解，如下：
 | (2) |
上述问题可以进一步转化为典型的线性规划问题，那么也就有很多成熟的算法可以使用，比如内点法、单纯形法等。除此以外，Li和Zhang还提出了可以利用TV范数重建原信号[18]。
2 LPV模型 本文中采用LPV模型的输入输出的表示形式为
 | (3) |
式中：y(k)为系统的输出信号；u(k)为系统的输入信号；参数p为一个可测量的参量；ai(p)和bi(p)为关于参数p的函数；na和nb分别为模型的输入和输出阶数；当p=k时，LPV模型也就变化为常见的线性时变(Linear Time-Variant, LTV)模型。在实际系统中，还要考虑到系统的输入延迟以及噪声，在此用e(k)为一个零均值的随机噪声；d为系统时延。上述模型可以改写为
 | (4) |
式中：
 |
 |
3 动态压缩测量辨识算法 CMI算法是利用压缩感知理论对线性时不变系统进行辨识的一种算法。该算法在应用于LPV模型时，则需要加以改进。本文受到压缩感知运动成像理论的启发[19]，将参数变化函数的连续变化视作多个离散点之前的转换。假设系统状态函数θ(p)是满足唯一性的关于参数p的复杂函数，即对于每个参数p都有唯一的状态θ(p)与之对应。这个假设在大部分非线性系统中都是满足的，比如温度控制系统中的传热系数是与温度一一对应的。但是，当系统中包含继电器特性环节，该系统就不满足这个假设。对函数θ(p)进行采样就可得到一系列离散的状态点{θ(p1), θ(p2), …, θ(pn)}。为了保证采样的准确性，采样过程要满足Nyquist采样定理，即采样频率要大于θ(p)带宽的二倍以上。此时，连续函数θ(p)被n个离散点分割为n－1段函数，原辨识问题也就变为了对每段函数的辨识问题。如果可以准确地辨识得到其中一段的系统状态，那么也就可以通过滚动迭代的方法得到局部参变函数。以下将以辨识系统其中一段，即从离散状态θ(p1)到θ(p2)的过程为基础。
由于参变函数可能是非线性函数，所以在离散后的每一段函数中还是非线性的，也就很难直接对其进行辨识。为了简化这个问题，可以引入函数θ(p)是连续n阶可导的假设。此时，该函数可以利用泰勒级数在p1和p2的中点处展开，此时参变函数可以表示为
 | (5) |
式中：θ(n)(p)为θ(p)的n阶导数；Rn(p)为泰勒公式的余项。如果只保留一阶导数，得到原函数的一个近似线性函数为
 | (6) |
该线性函数和原函数的近似误差主要取决于高阶项的大小。如果采样频率远大于原信号频率或者原信号是缓慢变化的，那么利用式(6)来近似就可以得到比较好的结果。然而，式(6)等号右侧的各项都是未知的，所以可以将式(6)改写为
 | (7) |
改写为式(7)的好处在于
相应的，系统的输出也可以视为由各离散状态以及输入信号而产生的。此时，如果使用传统的辨识算法，辨识结果将无法收敛。为了保证近似模型的精度就要提高系统采样频率，但系统的参数一直在改变，所以对于每个状态下所获得的辨识数据M的个数可能小于待辨识的参数个数N。所以，该辨识问题相当于求解一个欠定的方程组。显然，用传统的辨识算法都很难处理这种情况。这就需要用到压缩感知理论，结合该理论的CMI算法正好适用于这种采样数不足的情况。不过，利用CMI算法对该模型进行辨识，仍然会存在很大的误差。这是因为CMI算法适用于模型不变的系统，而此处系统状态是一直在改变的。因此，还需要引入“运动”的思想，即系统状态一直在变化的，是从一个离散状态“运动”到另一个离散状态的过程。而系统每一时刻的输入是对这个“运动”状态的一次观测，系统的输出就是状态在“运动”过程中所产生的。也就是说，系统在任意时刻的输出，都可以表示为由相邻的两个状态与输入信号的卷积所产生的。从之前的叙述可知，系统某一时刻的采样值y(k)可以认为是由临近的两个离散状态所共同影响，并且两个状态对系统输出影响的比例是一直在变动的。此时，式(4)可以改写为
 | (8) |
考虑M次采样，系统输出可以表示为
 | (9) |
式中：pi和qi分别为在k时刻状态θ(p1)和θ(p2)在系统输出中所占的权重；Φ代表测量矩阵；P、Q为权重矩阵；e为噪声向量。其中，权重矩阵P、Q与Φ的“*”运算为Hadamard内积，即对应项相乘。而且，从本节描述可知矩阵P、Q与变化速度v相关。对于“匀速”变化，pi和qi的大小由变化速度v所决定；对于非匀速变化，只需要按照“非匀速”速度计算得出各项值即可。在此，本文将[P*??Φ, Q*Φ]称为比例测量矩阵，用B表示。而[θ(p1), θ(p2)]整体变为了一个新的待辨识参量，其长度为原信号的两倍，用x表示。那么，式(9)就变为了常见的压缩感知的形式，y=Bx+e。通过上述方法，就可以在采集到所有辨识数据后获得系统模型，即离线的对系统进行辨识。而此方法是可以很容易推广到在线辨识的情况的。因为，本文对假设的各离散状态的位置是任意的。所以，只要在采集到了最新的M个辨识数据后，就可以辨识得到两个离散状态。而且，通过滚动辨识的方法就可以实时的获取最新的系统状态，利用该状态得到的系统输出预测值也是最准确的。在此给出，系统的一步预测器：
 | (10) |
式中：y*(k+1)为一步预测器在k时刻的预测值；
为了保证采样的准确性，就需要测量矩阵满足RIP。在压缩感知理论中，Candès和陶哲轩[12]证明了独立同分布的随机测量矩阵是满足RIP的。不过，由于线性系统的输出是用卷积的形式表示的，所以测量矩阵是一个有结构的Toeplitz矩阵。Bajwa等证明了当Toeplitz矩阵各项服从高斯分布或者伯努利分布时，Toeplitz矩阵有很大可能满足RIP，只要测量数M>O(S2lb(N/S))[20]，N为信号总长度。随后Rauhut等利用Dudley不等式证明了，当测量数M>O(S1.5·lb(N1.5))[21]，测量矩阵也是满足RIP的。不过，随机测量矩阵对测量数的要求只要M>O(Slb(N/S))[22]，这是因为Toeplitz矩阵各列之间的相关度要大于随机测量矩阵。对于本文中的比例测量矩阵，该矩阵是由两个Toeplitz矩阵分别乘以不同的系数所组成，因而该矩阵相比于测量矩阵，其各列之间的相关性被进一步放大。因此，要想获得同等准确的采样值，所需的测量数肯定要更多。
4 仿真试验 假设待辨识系统模型参数是稀疏的，即在{ai(p)}和{bi(p)}中一共只有S个位置是非零值，其他位置均为0。鉴于实际系统中S个非零项的大小和位置都是未知的，因此在每一次仿真中非零项的位置和大小在都会随机改变。由于DCMI方法并不要求提前设定准确的模型阶数，因而在仿真中选择的模型阶数足够大即可。仿真中设置稀疏度S=6，输入信号阶数na与输出信号阶数nb均为20，即信号总长度N=40。
系统的输入信号u(k)各项是从标准的正态分布中随机选取。为了接近真实情况，观测值是包含噪声的。该噪声的大小由信噪比(SNR)所衡量。本文比较了DCMI算法、CMI算法以及最小二乘算法对θ(p1)和θ(p2)的变化过程的辨识结果。此外，由于最小二乘算法在欠定条件下无法获得唯一解，因而本文在目标函数中增加了正则项来解决。其中，DCMI和CMI的恢复算法均选用OMP算法来重建模型参数信号。为了保证仿真结果的普适性，本文进行了200次蒙特卡罗仿真，并且在每次仿真中都会重新随机生成一个新的LPV系统模型。考虑到真实系统中的模型真值是未知的，因而本文除了比较模型本身外，还比较了利用模型预测的系统输出以及实测系统输出的差距来评价模型的准确性。具体指标为均方误差(Mean Square Error，MSE)和平均接近程度(Average Fit Rate，AFR)，该指标的具体形式如下：
 | (11) |
式中：ti为第i次仿真的系统观测向量；ti*为一步预测器获得的系统输出值；D为总的试验次数。由于系统的输出值是有噪声的，所以需要用MSE来评价的预测值的偏差和方差。而AFR描述的是预测值和观测值的近似程度，因而该指标可以更直观地反映系统预测的准确度。
4.1 匀速变化模型 假设两个系统状态θ(p1)和θ(p2)之间的变化是线性的，LPV模型真值和不同算法辨识结果的对比如图 1所示。虽然已知状态变化是线性，但是状态中的各参数的变化速度是不同的，而且其大小也是是未知的。
 |
| 图 1 LPV模型真值和不同算法辨识结果对比 Fig. 1 Comparison between truth value of LPV model and identification result of different algorithms |
| 图选项 |
此时，权重矩阵P和Q由系统测量值个数M给出：
 | (12) |
图 1中：Idx表示θ(p)向量中的位置，f(p)表示θ(p)中下标为Idx处的参数值大小。该图中展现了200次试验中的两次结果，其中图 1(a1)、(b1)为LPV模型θ(p)的真值，其他各分图分别表示各辨识算法计算得到的模型参数值。可以看出，最小二乘算法和CMI算法的结果是用常数来近似表示线性函数，只有DCMI算法辨识结果与模型真值是一致的。另外，最小二乘算法的结果在很多原本是零的位置，得到了大量的非零值，即导致了模型的过拟合，引入了大量的误差。
鉴于在高阶LPV模型的辨识问题中，系统输入输出的测量值个数是主要的影响因素，因而随后重点分析测量数对各算法的影响。由于在实际系统中，模型真值一般是未知的，因而也就无法直接对模型进行比较。另一个评价的办法就是给定一个输入信号，对比系统输出的实测值和利用模型得到的预测值，该预测值由第3节中提到的一步预测器给出。图 2显示了3种算法在不同测量数条件下的预测结果。其中，最小二乘算法效果最差，而且随着测量数的增加，预测精度反而下降。此外，该算法在测量数M=40时，由于过拟合导致极大的预测误差。相比较而言，CMI算法要强于最小二乘算法，但同样随着测量数的增加，算法的预测精度反而逐渐下降。这同样是由于算法没有考虑到参数是动态变化的所导致的。从图中可以看到，DCMI算法的预测精度显然是最高的，但该算法在测量数不足时表现较差。因为，该算法需要同时重建两个时刻的模型参数，其长度是其他两个算法的2倍。由于观测本身存在误差，而且一步预测器利用的是上一时刻的模型参数，所以DCMI的预测值与观测值还是存在差距的。同时，由于算法计算时间同测量数多少直接相关。因而，在仿真中记录了各算法所花费的时间，并发现算法计算花费时间均在在1 ms以内，但会随着模型规模的增加而增加。
 |
| 图 2 测量数对不同算法预测精度的影响 Fig. 2 Influence of measurement number on prediction accuracy of different algorithms |
| 图选项 |
接下来，本文对3种算法的鲁棒性进行了测试。为了避免测量数在结果中的影响，在此设置测量数M=50，仿真结果如图 3所示。从图中可以看出，CMI和DCMI两种算法对白噪声都有一定的抗干扰能力，而最小二乘算法的对噪声比较敏感，其结果偏差较大。当噪声分贝降低到40 dB左右时，噪声基本不会影响DCMI算法的预测的准确性。这是因为CMI算法考虑到了信号本身是稀疏的，所以大部分噪声在辨识结果中被抑制了。而且重建算法OMP采用了贪婪的思想，即只考虑了测量矩阵中影响最大的几列，从而降低了噪声对恢复结果的干扰。
 |
| 图 3 观测噪声对不同算法预测精度的影响 Fig. 3 Influence of measurement noise on prediction accuracy of different algorithms |
| 图选项 |
4.2 非匀速变化模型 在4.1节中考虑了参数“匀速”变化的情况，本节将以一个参数以“匀加速”变化的系统作为例子来说明该方法在参数非线性变化的LPV系统中的使用。先假设系统两个相邻的离散状态θ(p1)和θ(p2)，参数变化速度v以加速度a均匀增加。由于参数变化是非匀速的，所以权重矩阵也就需要根据测量数的不同而重新计算。已知，两个离散状态之间的任意状态可以由式(7)表示。那么，对于测量初始时刻ts的状态θs和最后时刻te的状态θe有：
 | (13) |
式中：λe和λs分别为初始时刻和最后时刻两个状态的影响比例。
将式(13)整理可得
 | (14) |
再将式(14)代入式(7)，可得θs和θe之间任意状态θ的表示为
 | (15) |
根据式(15)就可以计算得出DCMI算法所需的权重矩阵P和Q。由于4.1节已经分析了噪声的影响，所以为了减少噪声对结果的影响，设置噪声大小为40 dB。将图 4与图 2对比可以发现，在线性变化的系统中最小二乘算法和CMI算法的AFR与DCMI算法差距在10%左右，而在非线性变化的系统中的该差距被拉大到20%以上。这是因为最小二乘和CMI算法的本质都是静态算法，因而算法辨识出的模型参数是各状态的均值。而该参数均值与当前值的差距，在非线性变化过程中被拉大了，所以利用该均值得到的预测值与观测值的差距也就被拉大了。这同样影响到了DCMI算法的准确性，相比于线性变化情况，此处最佳的AFR要降低了2%左右。而且，可以看出此时算法所需的测量数要高于线性变化的情况。其原因可能是由于非线性的权重矩阵使得测量矩阵性质变差。
 |
| 图 4 线性变化情况下测量数对不同算法预测精度的影响 Fig. 4 Influence of measurement number on output prediction accuracy of different algorithms in case of linear variation |
| 图选项 |
5 结论 1) 根据仿真试验结果表明，DCMI算法中使用的匀速变化模型可以无损地对线性参变函数近似，而“非匀速变化”模型对非线性参变函数近似程度明显好于传统算法。
2) 本文提出的DCMI算法在测量数据量不足的情况下，仍然能够得到准确的LPV模型。另外，DCMI算法可以很好地避免参数过拟合问题，而无需根据先验知识选定恰当的模型阶数。
3) DCMI算法可以有效地抵抗白噪声对辨识结果的干扰。
参考文献
| [1] | HOFFMANN C, WERNER H. A survey of linear parameter-varying control applications validated by experiments or high-fidelity simulations[J]. IEEE Transactions on Control Systems Technology, 2015, 23(2): 416-433. DOI:10.1109/TCST.2014.2327584 |
| [2] | TURK D, GILLIS J, PIPELEERS G, et al. Identification of linear parameter-varying systems:A reweighted l2, 1-norm regularization approach[J]. Mechanical Systems and Signal Processing, 2018, 100: 729-742. DOI:10.1016/j.ymssp.2017.07.054 |
| [3] | 王明昊, 刘刚, 赵鹏涛, 等. 高超声速飞行器的LPV变增益状态反馈H∞控制[J]. 宇航学报, 2013, 34(4): 488-495. WANG M H, LIU G, ZHAO P T, et al. LPV H∞ control for hypersonic vehicle[J]. Journal of Astronautics, 2013, 34(4): 488-495. DOI:10.3873/j.issn.1000-1328.2013.04.006 (in Chinese) |
| [4] | VAN J W, VERHAEGEN M. Subspace identification of bilinear and LPV systems for open-and closed-loop data[J]. Automatica, 2009, 45(2): 372-381. DOI:10.1016/j.automatica.2008.08.015 |
| [5] | BAMIEB B, GIARRE L. Identification of linear parameter varying models[J]. International Journal of Robust and Nonlinear Control, 2002, 12(9): 841-853. |
| [6] | DE CAIGNY J, CAMINO J F, SWEVERS J. Interpolating model identification for SISO linear parameter-varying systems[J]. Mechanical Systems and Signal Processing, 2009, 23(8): 2395-2417. DOI:10.1016/j.ymssp.2009.04.007 |
| [7] | TOTH R, LYZELL C, ENQVIST M, et al.Order and structural dependence selection of LPV-ARX models using a nonnegative garrote approach[C]//Proceedings of the 48h IEEE Conference on Decision and Control(CDC)Held Jointly with 200928th Chinese Control Conference.Piscataway, NJ: IEEE Press, 2009: 7406-7411. |
| [8] | ROJAS C R, HJALMARSSON H.Sparse estimation based on a validation criterion[C]//Proceedings of 201150th IEEE Conference on Decision and Control and European Control Conference.Piscataway, NJ: IEEE Press, 2011: 2825-2830. |
| [9] | LAURAIN V, TOTH R, ZHENG W X, et al. Nonparametric identification of LPV models under general noise conditions:An LS-SVM based approach[J]. IFAC Proceedings Volumes, 2012, 45(16): 1761-1766. DOI:10.3182/20120711-3-BE-2027.00230 |
| [10] | GOLABI A, MESKIN N, TOTH R, et al. A Bayesian approach for LPV model identification and its application to complex processes[J]. IEEE Transactions on Control Systems Technology, 2017, 25(6): 2160-2167. DOI:10.1109/TCST.2016.2642159 |
| [11] | DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582 |
| [12] | CANDèS E J, TAO T. Decoding by linear programming[J]. IEEE Transactions on Information Theory, 2005, 51(12): 4203-4215. DOI:10.1109/TIT.2005.858979 |
| [13] | ROSSI M, HAIMOVICH A M, ELDAR Y C. Spatial compressive sensing for MIMO radar[J]. IEEE Transactions on Signal Processing, 2014, 62(2): 419-430. DOI:10.1109/TSP.2013.2289875 |
| [14] | WILLETT R M, DUARTE M F, DAVENPORT M A, et al. Sparsity and structure in hyperspectral imaging[J]. IEEE Signal Processing Magazine, 2014, 31(1): 116-126. DOI:10.1109/MSP.2013.2279507 |
| [15] | CANDèS E J, ROMBERG J. Sparsity and incoherence in compressive sampling[J]. Inverse Problems, 2006, 23(3): 969-985. |
| [16] | CANDèS E J. The restricted isometry property and its implications for compressed sensing[J]. Comptes Rendus Mathematique, 2008, 346(9-10): 589-592. DOI:10.1016/j.crma.2008.03.014 |
| [17] | PATI Y C, REZAⅡFAR R, KRISHNAPRASAD P S.Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[C]//Proceedings of Conference Record of the 27th Asilomar Conference on Signals, Systems and Computers.Piscataway, NJ: IEEE Press, 2002: 40-44. |
| [18] | LI C, ZHANG F. AVA inversion based on the L1-norm-based likelihood function and the total variation regularization constraint[J]. Geophysics, 2017, 82(3): 1-54. DOI:10.1190/geo2017-0321-TIOgeo.1 |
| [19] | 王盼盼, 姚旭日, 刘雪峰, 等. 基于行扫描测量的运动目标压缩成像[J]. 物理学报, 2017, 66(1): 76-83. WANG P P, YAO X R, LIU X F, et al. Motion target imaging using compress sensing[J]. Acta Physic Sinica, 2017, 66(1): 76-83. (in Chinese) |
| [20] | HAUPT J, BAJWA W U, RAZ G, et al. Toeplitz compressed sensing matrices with applications to sparse channel estimation[J]. IEEE Transactions on Information Theory, 2010, 56(11): 5862-5875. DOI:10.1109/TIT.2010.2070191 |
| [21] | RAUHUT H, ROMBERG J, TROPP J A. Restricted isometries for partial random circulant matrices[J]. Applied and Computational Harmonic Analysis, 2012, 32(2): 242-254. DOI:10.1016/j.acha.2011.05.001 |
| [22] | BARANIUK R, DAVENPORT M, DEVORE R, et al. A simple proof of the restricted isometry property for random matrices[J]. Constructive Approximation, 2008, 28(3): 253-263. DOI:10.1007/s00365-007-9003-x |
