动态相关分析 系统建模理论与方法 教学课件

第7章动态相关分析系统建模理论与方法教学课件7.1.2根本定理

设所研究的系统为zn＝f(Zn-1,Ｕn,θn)+ξn，n∈Ｎ(7-3)式中，zi表示时刻i的系统观测值，且Zn-1＝［zn-1，zn-2，…，zn-s］Ｔ；Un＝［un,un-1，…，un-m］Ｔ，ui表示时刻i的系统输入值；θn是r维参数向量，且θn＝［α１n,α２n,…，αrn］Ｔ，αjn为n时刻的第j个参数，j＝1，2，…，r；ξn是零均值弱平稳随机序列，它表示系统所受的干扰。7.1.3数据处理方法1.增长记忆法

2.渐消记忆法

3.限定记忆法1.增长记忆法

当样本长度为n时，所估计的参数值记为θn，那么θn中应该包含n组样本的有用信息。假设样本容量增加为n+1，那么相应的参数估计值记为θn+1，它包含着n+1组样本的信息，对于时变系统模型，E(θn)≠E(θn+1)。这种差异是由于新增加了一组样本而引起的，也表达了新息对参数的修正作用。由于每次增加样本容量时，所有过去时刻的信息都保存在参数估计值中，因此称这种信息取用的方式为增长记忆法。用增长记忆法处理数据时，“老〞的信息总是包含在所估计的参数中。这对于描述时变系统完全没有必要。因为老的样本信息实际上在抵抗新息对参数的修正作用，阻碍了参数对现实系统的跟踪作用。特别是当参数估计中对老的信息记忆过多时，甚至于会“淹没〞新息对参数的影响。新息不能修改参数估计值的现象称为“数据饱和〞，如果参数估计时出现了数据饱和，那么新息不再对参数估计值有修正作用。所以，动态相关分析不宜采用增长记忆法来处理数据。2.渐消记忆法

由于系统是时变的，很容易想象：样本越老，那么它偏离系统的现实情况越远。因此，为了充分反映系统当前的情况，应该重视新的样本而将老样本“遗忘〞。遗忘的方法就是对样本进行加权处理，让新样本对新参数估计值作较大的奉献。也就是通过权重来强调当前样本的作用，逐渐消去老样本数据对新参数估计值的影响。这种信息的取用方式称为渐消记忆法。设μ为一个小于1而大于0的实数，用μ组成加权矩阵W＝μn-1μn-2⋱1使用加权最小二乘法，其准那么函数为Jw＝eＴWe＝∑ni=1μn-ie２(i)(7-17)由式(7-17)得到的参数估计值就可以表达渐消记忆的效果。通常称μ为遗忘因子。3.限定记忆法在渐消记忆法中，遗忘因子对模型参数有一定的影响，特别是对于残留在参数序列中的随机噪声，影响更为明显。一般来讲，用增长记忆法或渐消记忆法所得到的参数序列，其中的随机噪声不会是弱平稳的。这种缺陷常常增加了建立参数模型的复杂性，或者使第二代参数模型得不到好的估计值。采用限定记忆法在一定程度上能弥补这一缺陷。限定记忆的概念主要指事先规定每次进行参数估计的样本长度m(即记忆区间)，在参数估计过程中，每增加一组新样本就去掉一组老样本。这样一来在任何时候估计参数时所使用的样本都是最新的m个样本。本书所讨论的动态相关分析主要采用这种数据处理方法。7.2参数跟踪7.2.1线性时变系统

7.2.2非线性时变系统

7.2.3多层递阶模型7.2.1线性时变系统

当式(7-18)中的函数f(·)具有线性形式时，它可以写成如下形式的线性时变系统：z(k)＝ΦＴkθk+ξ(k)(7-19)式中，ΦＴk＝［Zk-1，Uk］。对于式(7-19)描述的系统，根据第2章讲述的最小二乘法参数估计迭代计算的原理，可以构成迭代计算公式P－１k＝P－１k-1＋ΦkΦＴkθ^k＝θ^k－１＋PkΦk［z(k)-ΦＴkθ^k-1］(7-20)为了讨论方便，把修正矩阵Pk记为Pk＝δ‖Φk‖２(7-21)式中，δ为由新息引起的修正系数，它可以由式(7-20)和式(7-21)计算出来。这样，就得到时变参数θk的参数估计值跟踪公式θ^k＝θ^k-1+δ‖Φk‖２Φk［z(k)-ΦＴkθk-1］(7-22)对于式(7-22)的跟踪公式存在如下定理。定理7-2当线性时变系统式(7-19)的参数估计值跟踪公式(7-20)中取δ＝１时，那么由它所确定的参数估计值必满足e［k,θ^k］＝z(k)－ФＴkθ^k＝０∀k≥1(7-23)称e［k,θ^k］为后验残差。7.2.2非线性时变系统

当系统具有式(7-18)表示的一般形式时，其模型参数也存在与式(7-22)类似的跟踪公式。设式(7-18)中的函数f(·)关于θk的具有一阶连续偏导数，其梯度记为θf(k,θ)，那么参数跟踪公式为θ^k＝θ^k-1+δ‖θf(k,θ^k-1)‖２θf(k,θ^k-1)｛z(k)-f［Zk-1,Uk，θ^k-1,k］｝(7-24)为了评价式(7-18)的跟踪能力，也可以建立与定理7-2类似的残差定理。但这个定理的证明比较复杂，本书不作更进一步的讨论，仅给出有关的结论。有兴趣了解其证明过程的读者可以查阅本书所列的参考文献。定理7-3设非线性系统式(7-18)中的观测值z(k)几乎处处有界，在跟踪公式(7-24)中，对δ＞0一致的有limk→∞θf(k,θ^∗k)Ｔθf(k,θ^k-1)‖θf(k,θ^k-1)‖２＝μ＞0式中，θ^∗k为θ^k和θ^k-1的加权和，即对0＜αk＜1，有θ^∗k＝αkθ^k+(1-αk)θ^k-1。那么对于任何ε＞0，必有δ＞0和Ｎ＞0，使当k≥N时，｛θ^k｝满足｜e(k,θ^k)｜＝｜z(k)-f［Zk-1,Uk，θ^k，k］｜＜ε(a.s.)(7-25)7.2.3多层递阶模型

从上面的两个定理可以看出，在一定的条件下，由参数跟踪公式所计算的参数序列｛θ１，θ2，…，θk｝对时变系统具有较好的跟踪作用。这些成果已经被一些研究人员所证实。动态相关分析在此根底上还要求进一步研究参数序列的变化规律。根据样本空间和参数空间同步性的观点，可以建立参数模型来反映时变系统样本空间的变化。如果想知道未来样本空间的状态，可以先利用参数模型对相应时刻的参数空间进行预测。如果第二代参数仍然是时变的，那么可以继续对第二代参数跟踪，引出的新一代参数称为第三代参数。这是一种运用多层模型描述时变系统的方法，有些文献称之为多层递阶预报方法。一般来讲，可以对参数序列｛θk｝建立模型θk＝g［Θk-1，Uk，βk,k］+ηk(7-26)式中，Θk-1＝［θ０，θ１，…，θk-1］Ｔ；Uk表示由影响参数θk的因素组成的向量；βk为第二代参数；ηk为序列｛θk｝中的随机干扰。比较式(7-26)和式(7-18)，可以看出它们数学形式根本相同，因此可以仿照式(7-22)或者式(7-24)建立第二代参数βk的跟踪公式。值得注意的是，在式(7-18)中，假定随机干扰序列｛ξk｝是零均值白噪声是允许的，而在式(7-26)中，不能再对随机干扰｛ηk｝做任何假设。正如在7.1节曾强调的那样，｛ηk｝是否为白噪声由参数估计方法所决定。在实际应用中，保证｛ηk｝仍然为白噪声是相当困难的事情，所以使用多层递阶预报方法时，随机干扰不能太大，递阶层次也不宜过多。7.3变参数回归方程7.3.1参数子模型

7.3.2差分法

7.3.3预报—校正法

7.3.4机理分析法7.3.1参数子模型

对于模型式(7-27)，采用小节介绍的数据处理方法，或者采用7.2节介绍的参数跟踪方法，都可以得到参数向量序列｛θ1，θ2，…，θk｝。对参数估计向量中的每个元素都可以建立一个数学模型，其模型结构可以从该元素对时间的散点图上分析出来。例如，对其中的第i个元素αik,k∈N建立模型αik＝gi(βi，k)+ηik(7-28)式中，βi＝［γ1i,γ2i,…，γqi］Ｔ表示第二代参数向量，q表示它的维数；ηik为随机序列｛ηk｝的第i个分量。称式(7-28)为第二代参数向量的第i个参数子模型，或者简称为参数子模型。理解参数子模型的概念并不是一件很困难的事情，实施这一方法也很容易。但是，为了获得成功必须注意两个问题：一个问题是要谨慎地选用系统模型式(7-27)的参数估计方法，使参数序列｛θk｝具有良好的统计特性，以保证第二代参数估计能顺利地进行；另一个问题是由于模型只是现实世界的近似，子模型的误差会成倍地反响到系统模型之中，所以使用子模型法时，系统模型和参数模型都不应该太复杂，子模型的层次也不要太多。参数子模型法的主要缺陷是计算工作量较大。系统模型的每个参数都对应一个参数子模型，需要实施一次参数估计。如果系统模型的参数比较多，计算工作量会很大。差分法可以大大降低计算时变参数的工作量。7.3.2差分法由于微分算子D＝d/dt具有核c(λ)＝iλ(i为虚数单位)，因此D∈Lc。根据定理7-1和推论7-2，D把弱平稳过程映射为弱平稳过程。作者在1985年提出的时变参数差分修正法，在某些情况下能够提高模型的精度。差分法的根本原理如下：设时变系统为zt＝f［Xt，θt，t］+ξt(7-29)对式(7-29)两边微分，得到dzt＝ƏfƏXＴtdXt+ƏfƏθＴtdθt+ƏfƏt+dξt(7-30)这是一个以dθt为未知量的线性方程，由于当ξt为零均值白噪声时，dξt也是零均值白噪声，所以可以用最小二乘法来求解dθt，利用dθt来修正已经得到的系统模型参数估计值θ^t，即可到达动态回归的目的。实际使用这种方法时，往往用差分来代替微分，因此把式(7-29)改写为Δzt＝ƏfƏXＴtΔXt+ƏfƏθＴtΔθtƏfƏt＋Δξt(7-31)令y＝Δzt－ƏfƏXＴtΔＸ－ƏfƏt，ηt＝Δξt，把它们代入式(7-31)，得到y＝ƏfƏθＴtΔθt+ηt(7-32)7.3.3预报—校正法上面介绍的两种方法都是考虑建立参数空间的子模型，通过对子模型规律的把握增加系统模型的可靠性和准确度，从而到达动态建模的效果。实际上也可以不用建立参数子模型的方法，而直接使用修正样本空间的方法来实现动态建模。这种方法的根本思想是：先按传统的建模方法对系统建模，然后分析模型使用时所发生的误差及其原因，建立预报—校正模型，求解初始模型的校正值。所以，也称这种方法为预报—校正法。7.3.4机理分析法

机理分析是很重要建模方法，但是机理分析需要具备丰富的专业知识，是相应的专业课程讨论的重要内容。系统建模理论研究的重点是通过数据(或信息)建立数学模型，它不需要掌握很多的学科专业知识。在动态相关的研究中机理分析有时显得更重要一些，因为机理分析能帮助建模人员理解参数的含义，从而确定参数方程的结构。在动态相关分析中，寻找样本空间或者参数空间的变化规律是重要但有时又是很困难的工作。由于实际情况限制，理论上要求的条件总是只能局部或者近似地得到满足。这往往成为引起模型误差的主要原因。有时，不追求建模理论上的严密性，通过简单的机理分析或者凭经验确定参数的变化规律，也能得到较好的模型。特别是当参数的物理含义比较清楚又有相关学科的理论支持时，用机理分析确定该参数的变化规律往往是很有效的方法。7.4二次回归分析7.4.1二次回归的概念

7.4.2二次回归模型的建立过程7.4.3对二次回归分析的说明

7.4.1二次回归的概念

当讨论时变线性方程对复杂系统的跟踪时，曾经讲到可以用时变系数的直线来逼近一个非线性系统。下面进一步讨论这个观点的数学根底。设有非线性系统zt＝f(Ｘt，θ)+ξt(7-38)式中，Ｘt表示由p个解释变量组成的向量。假定函数f(·)对于变量Ｘt中各元素的偏导数都存在，且在Ｘ０连续，那么可以把函数在Ｘ０附近展开为泰勒级数f(Ｘ,θ)＝f(Ｘ０,θ)+ƏfƏＸＴ(Ｘ-Ｘ０)+(Ｘ-Ｘ０)ＴƏ２fƏＸƏＸＴ(Ｘ-Ｘ０)+…(7-39)在式(7-39)中，所有的X都是时间的函数，为了书写简单，略去了所有的时间下标；此外，所有的偏导数都表示在X0点的计算值。7.4.2二次回归模型的建立过程1.构造参数序列2.构造参数子模型的样本序列

3.建立参数子模型并进行变量筛选

4.建立二次回归方程1.构造参数序列

设系统由动态回归方程z＝a0＋∑pi=1aixi+ξ(7-45)描述，其中p为解释变量的个数，并且z,a,x，ξ都是时间的函数，为了表达方式比较简单，这里省略了对时间的标记。把式(7-45)写成向量的表达式z＝ΦＴθ+ξ式中，ΦＴ＝［１,x1,…，xp］；θ＝［a0,a1，…，ap］Ｔ。2.构造参数子模型的样本序列

根据参数估计的方法，可以看到参数估计值序列｛θ^１，θ^2，…，θ^l｝中的每一个向量都包含了m个原始样本数据的信息，所以对参数子模型进行统计分析时，应该考虑与之对应的m个原始样本的平均值。令xi表示与参数序列中第i个参数估计向量θ^i相对应的解释变量的数值。那么X(i)=1m∑mj=1X(i+j-1)(7-46)式中，X(i)=［x1(i)，x2(i)，…，xp(i)］Ｔ为系统解释变量所组成的向量，小括号内为时间序号。3.建立参数子模型并进行变量筛选

设参数估计值θ^i和解释变量Xi之间存在线性关系，有θ^i=B０＋BXi+η(7-47)式中，B０=［β１，β２，…，βp］Ｔ，B为第二代参数矩阵，且B=β１２β12…β1pβ21β22…β2p︙︙︙βp1βp2…βppη为第二代参数模型中的噪声向量，η=［η１，η２，…，ηp］Ｔ。从式(7-47)中取出第k行，并略去时间编号i，有

αk=βk+∑pj=1βkjxj+ηk，k=1,2,…，p(7-48)根据二次回归的原理，这里对式(7-48)感兴趣的不是求解参数βk和βkj，而是判断参数αk和解释变量xj，j=1，2，…，p的相关程度。利用小节讨论过的t检验方法判别式(7-48)中各解释变量与参数αk的相关性，筛去相关性不强的解释变量。4.建立二次回归方程

把经过变量筛选的式(7-48)代回式(7-41)，构成二次回归方程。用总长度为n的样本集合估计二次方程中的参数，得到研究对象的模型。如果有必要，可以在第三级、第四级参数子模型中筛选解释变量，得到相当于三阶、四阶泰勒级数的二次回归方程。应该注意，并不是参数子模型的级数越高模型越好，当参数子模型的级数过高时，计算误差可能会导致建模过程完全失败。7.4.3对二次回归分析的说明

二次回归分析的思想也可以扩展到其他的动态建模方法之中。原那么上讲，只要动态建模过程是分级或分层进行的，都有可能使用二次回归分析的方法。例如在预报—校正法中，如果预报方程和校正方程都是线性方程，那么可以在利用残差分析求出校正方程的结构后，把它代回预报方程，构成一个统一的模型，并重新估计它的参数。在很多情况下，这种数据处理方法可以进一步提高预报—校正法的精度。同样，在差分法中，也可以把参数的修正模型和系统模型合并在一起，统一估计模型参数。不过最后需要说明的是，并不主张任意地推广二次回归分析的概念，在很多情况下，分级分层建立模型仍然是非常必要的，如果硬要把它们合在一起，反而会弄巧成拙。即使在可以使用二次回归分析的案例中，也要权衡一下不同方法的利弊，切不可把原理到处生搬硬套。7.5*变参数ARMA模型7.5.1根本定义

7.5.2几个定理

7.5.3时变ARMA模型的描述方法7.5.1根本定义设｛x(t)，t∈T｝是随机过程，那么对于每一个t∈T，x(t)均为概率空间(Ω，F，P)上的一个随机变量。考虑在时间t上x(t)的一阶矩，有如下定义。定义7-2如果随机过程x(t)的均值函数是独立于时间t的常数，即E［x(t)］=C那么称之为一阶矩平稳随机过程。所有的一阶矩平稳随机过程组成的集合记为S１。如果随机过程y(t)不属于S1，那么称之为一阶非平稳随机过程，它们组成集合S１的补集S１。定义7-3假设随机过程的均值为有限值，即E［x(t)］＜∞那么称之为一阶矩过程。所有一阶矩过程组成集合L１L1=｛x(t)｜E［x(t)］＜∞｝(7-49)显然可以看出，L1⊃S1。类似地，还可以考虑随机过程的二阶矩。定义7-4设S(s,t)为随机过程x(t)协方差函数，假设B(s,t)=D(s-t)(7-50)那么称之为二阶矩平稳随机过程，或者简称为二阶平稳过程。7.5.2几个定理定理7-4设随机过程x(t)∈H，那么x(t)为零均值正交随机过程的充要条件是x(t)=a(t)ξ(t)(7-54)式中，a(t)为时间的函数；ξ(t)为白噪声，其均值为零，方差为σ２１。7.5.3时变ARMA模型的描述方法根据定理7-7可以建立能描述一类非平稳过程x(n)的时变ARMA模型，即∑pk=0ak(n)x(n-k)＝∑qk=0bk(n)ξ(n-k)(7-68)式中，｛ξ(n)｝为零均值白噪声序列；｛x(n)｝为零均值二阶矩过程，其协方差函数为B(n,m)。当｜n-m｜＞q时，B(n,m)的值为0。一般来讲，式(7-68)描述了一个二阶非平稳序列，所以称之为非平稳自回归滑动平均模型，记为NARMA。和ARMA模型类似，式(7-68)表达的NARMA也可以简化为NAR模型或者NMA模型。NAR模型的一般形式为∑ak(n)x(n-k)=η1(n)(7-69)式中，｛η１(n)｝为零均值正交过程。7.6*非齐次马尔科夫链和时变图表7.6.1非齐次马尔科夫链

7.6.2转移矩阵估计

7.6.3时变图表7.6.1非齐次马尔科夫链本书在第6章比较详细地讨论了与马尔科夫序列有关的问题，介绍了状态转移概率矩阵P的概念。并且指出，为了强调P是系统从t时刻状态过渡到t+1时刻的状态转移概率矩阵，把它记为tP式中，n表示状态空间的维数；pij(t)表示在t时刻到t+1时刻之间对象从状态i转向状态j的概率。当tP不随时间变化，即tP≡P时，就构成齐次马尔科夫链。齐次马尔科夫链的建模问题在第6章已经进行了比较详细的讨论。本小节主要介绍非齐次马尔科夫链的建模问题。7.6.2转移矩阵估计假设矩阵的每个元素都是时间的函数，那么称之为函数矩阵。因此，非齐次马尔科夫链的状态转移概率矩阵实际上是一个函数矩阵。在下面的表达中为了简单起见，有时用函数矩阵一词代替非齐次马尔科夫链的状态转移概率矩阵。但是应该注意本节所讨论的函数矩阵，其矩阵元素pij(t)(i,j=1,2,…，n)的解析式一般不能准确地给定，只能通过建模手段把它们模拟出来。如果系统的历史状态空间，那么就有可能估计函数矩阵的运动形态。7.6.3时变图表

图和表都是非解析的建模方法。本书的3.5节曾介绍过静态图表的概念。如果相同构造的表格在不同的时间被不同的数据填充，那么就构成了时变图表。时变图表的建模工作比非齐次马尔科夫链的建模工作还要困难，因为后者有一些关于马尔科夫过程的研究成果可以借鉴，而前者几乎没有成熟的数学理论可供参考。从原那么上讲，处理时变图表中的模型问题，首先要从图表中寻找可以量化的因素，建立由观测量构成的序列，然后再考虑对这个序列建立解析模型。在实际工作中所建立的图表，有一些已经把时间作为一个根本变量考虑在其中了，这种情况一般可以用前面介绍的建模方法直接建模。还有一些图表只表示某一个特定时刻或时段的情况，不能表现系统的动态特征，这种图表处理起来就比