系统辨识MATLAB仿真

设一非线性系统如下所示:()()()()()111y k y k e u k k βαε--=-+-+0.75α= 0.35β=0.25γ=,()k ε是零均值，方差为0.01的噪声序列（均匀白噪声）。

（1）试设计一种激励信号能持续激励系统的各工作点（平衡点）（2）用适当的方法辨识出系统的等价模型（用另一组数据来检验模型的泛化性） 说明：下面讨论的都是离散系统，所以时间坐标均采用离散时间节点k 。

解：(1) 线性化处理寻找系统的合理输入信号 可以求得系统的平衡点为：()0.75y k α== (1.1)按题意要求最后系统必须收敛于平衡点附近，即：()lim 0.75k y k →∞= (1.2)为了找出系统的合理输入信号，使得系统最终工作在平衡点附近，这里可以将系统线性化处理，将上述非线性系统进行泰勒展开得：()()()()()()()2323111111112!3!!n n y k y k y k y k y k u k k n αβαβαβαβγε--+---++-=-+因为 ，后面()1y k -的高阶项都可以扔掉（只作为寻找输入信号使用）， 所以系统可以化为下式：此时不妨设系统输出()y k 的最后的极限为A , 从式1.2得0.75A =。

那么应该满足()()lim lim 1k k y k y k A →∞→∞=-= (1.5)从而有 ()()()11A u k k αβγε-=-+ (1.6)同时为了抵消系统的部分噪声，这里采用MA TLAB 软件编程产生另一服从同一分布的均匀噪声()1k ε，将式1.6变形得：()()()01111u k u k k αβεγγ--=-(1.7)式1.7中()0u k 是一个最后收敛于系统平衡点0.75的基本信号，这里可以采用一阶线性系统的阶跃响应曲线作为基本信号()0u k ，同时考虑系统的平衡点，即设计为：()/00.75k Tu k e-=- (1.8)T 是一阶线性系统对应的时间常数，反应到输入基本信号()0u k 上就是过零点作()0u k 对应()()()()11y k y k u k k αβγε--=-+01αβ<（1.3）(1.4)曲线的切线的斜率1KT =。

图1 基本信号产生框图很显然有()lim0.75tu k→∞=(1.9)从而最后的输入信号设计为：()()()011111u k u k kαβεγγ-=+-+(1.10)式1.10中,,αβγ题中已知，()1kε满足：()1E kε=⎡⎤⎣⎦，()10.01D kε=⎡⎤⎣⎦，()u k 的产生参照图1。

代入式1.8得：()()(1)/111(0.75)1k Tu k e kαβεγγ-+-=--+(1.11)产生()u k过程中的时间常数T由仿真效果决定，如果从系统的性能指标而言，T尽可能小，仿真如下：图2 噪声随机抵消情况下的原系统响应曲线系统输入基本信号()u k阶跃信号11Ts+图3噪声完全抵消情况下的原系统响应曲线仿真分析：图2与图3中都是原非线性系统的响应曲线图，两者的输入信号有差别，分别为()()12,u k u k ，如下式表示：()()()1011111u k u k k αβεγγ-=+-+ (1.12)()()()201111u k u k k αβεγγ-=+-+ (1.13)()0u k 参照式1.8。

式1.12中()11k ε+是MA TLAB 中产生的符合同一均匀分布的给定随机噪声，从而式1.12中输入信号()1u k 实现噪声随机抵消，结果系统的响应()y k 呈现出很大的噪声干扰，这是符合对未知系统认识的逻辑过程的，如果要得出理想的响应曲线，也就是要求系统最终在平衡点附近趋于稳定，那么需要设计相应的滤波器，这里不再讨论，后续会尝试设计滤波器滤波。

式1.13中()1k ε+是系统给定的噪声，那么设计的输入信号()u k 实现了对噪声的完全抵消，得到的理想响应曲线如图3所示，从图3中可以看出，()y k 很快的收敛于0.75附近，且上升时间以及超调量等性能指标均符合理想响应要求，这里就没有具体计算，显然符合控制系统的快稳准的三要素。

(2) 最小二乘法辨识模型通过MATLAB 编程从原非线性系统中采集一些输入输出数据，利用最小二乘法辨识出系统近似的线性模型，实现框图如图4。

时不变SISO 系统动态过程的数学模型为()()()()()11A z y k B z u k k ε--=+ （2.1）其中()u k 是系统的输入，()y k 是系统的输出，()()11,B A z z --的展开式分别如下：()112121a a n n A z a z a z a z ----=+++ （2.2） ()11212b b n n B z b z b z b z ----=++ （2.3）参数1212,,,,,,,a b n n a a a b b b 待辨识。

式2.1写成最小二乘格式()()()Ty k hk k θε=+ （2.4）式2.4中：()()()()()[y 1,,y ,u 1,,u ]T a b h k k k n k k n =------ （2.5）1212,,,,,,,a b Tn n a a a b b b θ⎡⎤=⎣⎦ （2.6）对于1,2,,,k L =方程2.4构成一个线性方程组，可以把它写成式2.7，L L L Y H E θ=+ （2.7）其中()()()1,y 2,,y TL Y y L =⎡⎤⎣⎦ （2.8） ()()()1,2,,TL E L εεε=⎡⎤⎣⎦ （2.9）()()()()()()()()()()()()()()()101011212211T a b T a b L T a b h y y n u u n y y n u u n h H y L y L n u L u L n h L ⎡⎤----⎡⎤⎢⎥⎢⎥----⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥------⎣⎦⎣⎦ （2.10）由题意知()k ε为均匀分布的噪声所以{}(){}(){}(){}120L E E E E E L εεε⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦（2.11）{}{}(){}()(){}()(){}()(){}(){}()(){}()(){}()(){}(){}22221121212212T L L L E E E L E E E L Cov E E E E I E L E L E L εεεεεεεεεεεσεεεεε⎡⎤⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦（2.12）式中20.01εσ=，I 是单位矩阵。

这里认为噪声()k ε与输入()u k 是不相关的，从而()(){}0,,E k u k l k l ε-=∀ （2.13）这里需要考虑的是记忆长度的选择，上述矩阵方程包含()a b n n +个未知数，一般要求a b L n n >+，这样才可能确定一个“最优”的模型参数θ，而且为了保证辨识的精度，L 必须充分大，由AIC 定阶准则知L 同样决定着所辨识系统的阶次，下面会阐述相应理论。

定义准则函数：()()()()21LTk J k y k h k θθ=⎡⎤=Λ-⎣⎦∑ （2.14）其中()k Λ是加权因子，它衡量的是采样数据的可信度，这里对数据的可信度难以肯定，就简单的将取()1k Λ=。

将准则函数写成二次型的形式如式2.15。

()()()TL L L L J Y H Y H θθθ=-- （2.15）设wls θ∧使得()|min wls J θθ=()()()||0T T L L L L wls wls J Y H Y H θθθθθθθ∧∧∂∂=--=∂∂ （2.16） 利用下面的两个向量微分公式展开得()()2T TT T x xx Ax x A x αα∂⎧=⎪⎪∂⎨∂⎪=⎪∂⎩（2.17） 其中A 是对称阵，得到正则方程()TTwls LL L L HH H Y θ∧= （2.18）当TL L H H 为正定矩阵时，有：()1T Twls L L L L H H H Y θ∧-= （2.19）从而()22|20wlsTL L J H H θθθ∧∂=>∂ （2.20） 结合式2.16到2.20知，wls θ∧使()|min wls J θθ=，并且wls θ∧是唯一的。

从而求得最小二乘估计值 ()1T T LS LLL L H H H Y θ∧-= （2.21）定阶：系统的等价线性模型的阶次确定可以基于AIC 准则来确定，首先将系统写成如式2.7的形式。

输出变量在n θ条件下的似然函数为：()()()22212exp 2L T n n n n n n n L Y H Y H θπσθθσ-∧∧⎧⎫⎛⎫⎪⎪=---⎨⎬ ⎪⎝⎭⎪⎪⎩⎭（2.22） 对应的对数似然函数为：()()()()221ln ln 2ln 222Tn n n n n n n n L L L l Y H Y H θθπσθθσ∧∧==----- （2.23） 其中L 是数据长度，根据极大似然原理得到模型参数θ及白噪声()k ε方差的极大似然估计值如下：()121T TML n n n n T ML ML n n n n H H H Y Y H Y H L θσθθ∧-∧∧∧⎧=⎪⎪⎨⎛⎫⎛⎫⎪=-- ⎪ ⎪⎪⎝⎭⎝⎭⎩（2.24） 将ML θ∧与2σ∧回代式2.23可得：2ln 2ML L l const θσ∧∧⎛⎫=- ⎪⎝⎭（2.25）AIC 准则的标准式为：2ln 2ML AIC N L N θ∧∧∧⎛⎫⎛⎫=-+ ⎪ ⎪⎝⎭⎝⎭（2.26）将式2.25代入AIC 准则的标准式2.26中得到：2ln 4AIC n L n σ∧∧∧⎛⎫=+ ⎪⎝⎭（2.27）上面各式中n ∧表示的就是模型要估计的阶次，最后找到使min AIC n ∧⎛⎫= ⎪⎝⎭的n ∧作为模型的阶次。

本次辨识采样的数据长度1000L =, 很显然2σ∧足够小，当数据长度取定时，显然n ∧可以尽可能的小，暂且取n ∧=4。

仿真分析：图5 定4阶情况下系统的辨识效果图6 定4阶情况下更换输入为原来的14情况下系统的辨识效果图7 定4阶情况下更换输入为原来的4倍情况下系统的辨识效果图5与图6及图7中分别给出了定4阶情况下的系统利用最小二乘法辨识的结果，图5中的辨识效果在幅值上出现减小，有一定的误差，但基本可以与原系统同步，为了检验模型的泛化性，在MATLAB编程过程中，将输入信号减小为原来的14以及增大为原来的4倍，图6中得到了较好的辨识效果，但图7却进一步增大幅度误差，从这里看出所辨识模型的泛化性并不好，即随着输入信号的增大，模型与系统的幅度差值增大。