CN102667755B - 通过迭代线性子空间计算对时变、变参数和非线性系统进行实验建模的方法和系统 - Google Patents

通过迭代线性子空间计算对时变、变参数和非线性系统进行实验建模的方法和系统 Download PDF

Info

Publication number
CN102667755B
CN102667755B CN201080050251.0A CN201080050251A CN102667755B CN 102667755 B CN102667755 B CN 102667755B CN 201080050251 A CN201080050251 A CN 201080050251A CN 102667755 B CN102667755 B CN 102667755B
Authority
CN
China
Prior art keywords
equation
state
estimation
linear
lpv
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201080050251.0A
Other languages
English (en)
Other versions
CN102667755A (zh
Inventor
华莱士·E.·拉里莫尔
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Individual
Original Assignee
Individual
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Individual filed Critical Individual
Publication of CN102667755A publication Critical patent/CN102667755A/zh
Application granted granted Critical
Publication of CN102667755B publication Critical patent/CN102667755B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/18Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F9/00Arrangements for program control, e.g. control units
    • G06F9/06Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/46Multiprogramming arrangements
    • G06F9/48Program initiating; Program switching, e.g. by interrupt
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00Aspects of pattern recognition specially adapted for signal processing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Computing Systems (AREA)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Testing And Monitoring For Control Systems (AREA)
  • Feedback Control In General (AREA)

Abstract

一种用于估计微分或差分方程的方法和系统,该方程可以控制非线性、时变和变参数的动态过程或系统。所述用于估计所述方程方法和系统是通过对所述方程的观测到的输出和(当需要时的)输入数据进行估计来实现的。所述方法和系统可采用任何如下的系统或过程:该系统或过程能够被描述成非线性、时变和变参数的差分方程并且可用于在描述详细的系统或方法性能时自动提取该差分方程,该系统或方法性能用于系统控制、故障检测、状态估计以及系统或方法中有变化时它们的预测和修改。

Description

通过迭代线性子空间计算对时变、变参数和非线性系统进行实验建模的方法和系统
相互参考应用
本申请根据35U.S.C.§120主张2009年9月3日提交的申请号为61/239,745的美国临时专利申请的优先权,其内容通过引用整体并入本申请。
背景技术
根据测得的输出数据和可能的输入数据来对非线性和时变动态过程或系统进行建模是该技术的新兴领域。根据理论或应用的领域,这在统计学中称为时序分析,在工程学中称为系统识别,在心理学中称为纵向分析,且在财务分析中称为预测。
过去,已对子空间系统识别方法进行改革、对具有反馈的系统的最优方法进行大量改进和改良,以及对包括双线性系统和线性变参数(LPV)系统的非线性系统的方法进行研究。子空间方法可以避免不收敛的迭代非线性参数最优化,并且可以对高阶大规模的系统使用具有相当大价值(considerablevalue)的数值稳定性方法。
在时变和非线性系统领域,已经开展了一些工作,虽然没有达到期望的结果。这些工作在最好将线性子空间方法的直接展开(directextension)用于非线性系统建模方面代表了本领域的现状。这种方法将旧值和新值表示为旧值输入和输出的非线性函数的线性组合。该方法的一个结果是,在测得的输入、输出、状态和所使用的旧值滞后的数量方面,旧值和新值的维数成指数增长。当只使用每种这些变量中的几个时,旧值的维数计算为超过104或甚至多于106。对于典型的工业过程,旧值的维数很容易超过109或甚至1012。这些极大的数字最多导致利用效率低的结果。
其它方法使用估计模型中非线性项的迭代子空间方法,因而需要很少的计算。这种方法包括启发式算法,且在具有随机调度函数(即,具有白噪声特性)的LPV系统情况下已经用于高精度模型识别。然而,问题之一是在大部分LPV系统中,调度函数通常由特定应用确定且特征是通常是非随机的。已经实施了一些修改以试图在非随机调度函数的情况下提高精度,结果是,所尝试的修改在充分提高模型精度方面没有获得成功。
在更通常的内容中,识别非线性系统的通常问题被认为是通常的非线性规范变量分析(CVA)过程。该问题可用洛伦茨吸引子(Lorenzattractor),由简单的非线性差分方程描述的混沌非线性系统,来阐述。以此方式,旧值和新值的非线性函数被确定来描述过程的状态,该过程的状态接着被用于描述系统的非线性状态方程。这种方法的一个主要难题是找到可行的计算实施,因为众所周知的是,所需要的旧值和新值的非线性函数的数量是成指数增长的。在寻找一般非线性系统所存在的系统识别问题的解决方案时经常会遇到这个难题。
因此,在下面描述的一些示例性的实施例中,描述的方法和系统可以实现相当大的改进且还可以在观测的‘大量采样’是可用的情况下产生最优结果。并且,该方法不是‘特定的(adhoc)’,而是可以包括最优的统计方法。
发明内容
一个示例性的实施例描述了一种用于使用非线性、时变和变参数动态过程的方法。该方法可用于产生减小的具有时变项的系统模型。该方法可包括如下步骤:展开状态差分方程;将差分方程表示为关于输出和增量输入的线性、时不变系统;以及估计状态方程的系数。
另一示例性的实施例描述了一种用于估计一组控制非线性、时变和变参数过程的方程的系统。该系统具有第一输入、第二输入、反馈框和时延框。并且,在该系统中,第一输入和第二输入经过反馈框到达时延框,从而产生输出。
具体实施方式
在下面的说明和指向本发明特定实施例的相关附图中描述了本发明的一些方面。本领域技术人员将意识到可以设计替代的实施例,只要其不背离权利要求的精神或范围。并且,这里不详细描述本发明示例性的实施例的公知原理,或省略这些公知原理,以免使得本发明的相关详细内容变得模糊不清。
如这里所使用的,词组“示例性的”表示“作为例子、实例或例证”。这里描述的实施例不只限于示例性的实施例,而可以包括示例性的实施例以外的实施例。应该理解的是,所描述的实施例不必然被解释为对其他实施例是优选的或具有优势的。并且,术语“本发明的实施例”、“实施例”或“发明”不要求本发明的所有实施例都包括所描述的特征、优势或操作方式。
此外,许多实施例是根据将要由例如计算设备的元件执行的操作顺序来描述的。应该认识到的是,可以通过特定电路(例如,特定用途集成电路(ASICs)),通过由一个或多个处理器执行的程序指令,或由这两者的组合来执行这里描述的各种操作。并且,这里描述的这些操作顺序被认为整个包含在任何形式的计算机可读存储媒介内,该存储媒介具有存储在其内的一组对应的计算机指令,当执行这些指令时会使得相关的处理器执行这里描述的功能(functionality)。因此,本发明的各个方面可以通过多种不同的形式具体体现,所有这些形式已经预期落入权利要求的主题范围内。并且,对于这里描述的每个实施例,任何这样的实施例的对应形式都在这里描述为,例如“逻辑被配置来(logicconfiguredto)”执行所描述的操作。
通常参考示例性的图1-6,描述了用于对时变、变参数和非线性差分方程进行实验建模的方法和系统。所述方法和系统可以被实施和使用来产生各种结果,并且所述方法和系统可以有效地实施。
如示例性的图1所示,显示了根据一个示例性的实施例,对时变、变参数和非线性差分方程进行实验建模的一套方法的流程图。这里,在102处,使用一组时变、变参数和(如果需要的话)非线性状态空间差分方程。然后在104处,关于所选择的一组基函数来展开所述方程,例如非线性输入-输出方程可以以xt和ut的多项式展开。然后在106处,例如根据输出yt和增量(augmented)输入ut将差分方程表示为线性时不变形式,该线性形式可以包括输入ut和基函数,例如,输入ut、调度函数ρt和状态xt的多项式。
示例性的图2显示了示例性的流程图,其中使用差分方程的线性、变参数系统。在该实施例中在202处,使用一组线性、变参数状态空间方程,如下面显示的方程1和方程2。
xt+1=A0xt+B0ut+[A1ρt(1)+...+Asρt(s)]xt+[B1ρt(1)+...+Bsρt(s)]ut方程1
yt=C0xt+D0ut+[C1ρt(1)+...+Csρt(s)]xt+[D1ρt(1)+...+Dsρt(s)]ut方程2
然后在204处,可以用关于调度函数ρt、状态xt和输入ut的多项式展开状态空间差分方程,例如如下面的方程3和4所示。
x t + 1 = A 0 x t + B 0 u t + [ A 1 ... A s ] ( ρ t ⊗ x t ) + [ B 1 ... B s ] ( ρ t ⊗ u t ) 方程3
y t = C 0 x t + D 0 u t + [ C 1 ... C s ] ( ρ t ⊗ x t ) + [ D 1 ... D s ] ( ρ t ⊗ u t ) 方程4
接下来,在206处,可以以初始输出yt和增量输入的形式来表示差分方程,增量输入是关于ut、xt和ρt的函数。差分方程具有可被估计的线性非时变的未知系数(A0,[B0AρBρ],C0,[D0CρDρ]),差分方程如下面显示的方程5和6。
x t + 1 = A 0 x t + [ B 0 A ρ B ρ ] u t ρ t ⊗ x t ρ t ⊗ u t 方程5
y t = C 0 x t + [ D 0 C ρ D ρ ] u t ρ t ⊗ x t ρ t ⊗ u t 方程6
然后,如在208处所示,增量输入和(在一些示例性的实施例中特别是)可以包含未知的状态向量xt,所以可以采用或需要迭代。因此,在这些示例性的实施例中,正如下面更详细描述的,可以采用使用迭代算法的迭代。
示例性的图3显示了可被执行用于迭代子空间识别的迭代算法的流程图。在这个实施例中,如在302处看到的,非线性差分方程可以用增加的(additive)基函数来展开,并且可以用增量输入ut的线性时不变形式来表示该非线性差分方程。在一些例子中,这可以包括包含输出、状态和调度函数的非线性基函数。
然后,在304处,状态估计是未知的,其等于不在LPV模型中的(AρCρ)项,因此对应的项可以从一组增量输入中删除。然后,该迭代算法可以将增量输入用作输入来执行并且可以计算模型参数的估计然后,例如对估计的参数值使用卡尔曼滤波器,根据一步预估的新方法(innovations)可以计算出状态估计然后,可以估计似然函数。
接下来,在306,可以执行迭代k(k>2)。这里,状态估计对于所有t都是初始化的。该迭代算法可以将增量输入用作输入来执行并且可以计算参数的估计再次,例如,对估计的参数值使用卡尔曼滤波器,则可以计算出状态估计。同样可以执行一步预估的新方法并且同样可以估计似然函数
在308处可以验证收敛性。这里,可以比较log似然函数的值的变化以及迭代k-1和迭代k之间的状态阶(orders)。如果,在一些示例中,状态阶相同并且log似然函数的变化小于所选择的阈值ε。该阈值在一些示例中小于1,例如为0.01,然后迭代结束或停止。否则,似然函数变化值高于所选择的阈值ε,则回到上面的步骤306并且执行迭代k+1。接下来执行迭代k+1,然后再次验证收敛性。
在另一个实施例中,对于使用例如线性非时变系统所采用的子空间方法的系统中自相关的误差和反馈的情况,可以使用不同的方法直接和简单地获得未知参数的最优或期望的估计。这可以通过用不同的形式表述这个问题来实现,这些不同的形式导致需要在状态估计上进行迭代;然而,迭代的次数可能是非常低的,并且,为了进一步简化系统和其展开式(development),可以消除随机噪声。
例如,考虑到线性系统,其中系统矩阵是如以下方程所示形式的调度参数向量ρt=(ρt(1)ρt(2)...ρt(s))T的时变函数:
xt+1=A(ρt)xt+B(ρt)ut方程7
yt=C(ρt)xt+D(ρt)ut方程8
对于调度参数的仿射相关性(affinedependence),状态空间矩阵可具有以下方程9至12的形式。
A ( ρ t ) - A 0 + Σ j = 1 s ρ t ( i ) A i 方程9
B ( ρ t ) = B 0 + Σ i = 1 s ρ t ( i ) B i 方程10
C ( ρ t ) = C 0 + Σ i = 1 s ρ t ( i ) C i 方程11
D ( ρ t ) = D 0 + Σ i = 1 s ρ t ( i ) D i 方程12
在以上方程中,应该注意到,左边的矩阵表示为右边由调度参数ρt表示的线性组合(1ρt(1)ρt(2)...ρt(s))T,该线性组合被称为仿射形式。在进一步的论述中,该仿射形式的符号将被用于系统矩阵A=[A0A1...As]=[A0Aρ],且类似地可用于B、C和D。
在一些其它示例性的实施例中,LPV系统类的系统识别方法具有多种潜在的应用和经济价值。这样的系统可以包括,但不限于空气动力学的和流体动力学的交通工具(例如飞机和轮船)、机动引擎动力学、涡轮引擎动力学、化学处理(例如搅拌槽反应器和蒸馏塔)和其他系统。一个特点是,系统动力学在任何已知的操作点ρt都可被描述为线性系统。调度参数ρt是关于操作点变量(例如,但不限于速度、压力、温度、燃料流和类似的变量)的复杂的非线性函数,该操作点变量可以是用于在可能未知的常数矩阵A、B、C和D范围内表征该系统动态的已知的或精确测得的变量。还可假设ρt是可计算的或可从任何这样的操作点变量的信息中确定。例如,机动引擎的LPV模型可以包括LPV状态空间方程,该LPV状态空间方程明确地将向量ρt项表示为关于各种操作点变量的极其复杂的非线性函数。在这里描述的一些示例性的实施例中,当执行系统识别计算时,仅期望的是,调度参数ρt是可用的。这样放宽了实时使用或如实时控制或滤波这样的应用的要求。
为了简化该论述,可以通过使调度参数ρt与输入ut、状态xt相关联将LPV方程写成时不变形式,如
x t + 1 = A 0 x t + B 0 u t + [ A 1 ... A s ] ( ρ t ⊗ x t ) + [ B 1 ... B s ] ( ρ t ⊗ u t ) 方程13
y t = C 0 x t + D 0 u t + [ C 1 ... C s ] ( ρ t ⊗ x t ) + [ D 1 ... D s ] ( ρ t ⊗ u t ) 方程14
这里,可以表示克罗内克积对于由i、j块形成的作为隔开矩阵的任何矩阵M和N,被限定为其中M的i、j项表示为mij。并且,符号[M;N]=[MTNT]T还被使用来堆栈(stacking)向量或矩阵M和N。则上述方程13和14还可以写成如下所示的方程15和16的格式。
x t + 1 = A 0 x t + [ B 0 A ρ B ρ ] u t ρ t ⊗ x t ρ t ⊗ u t 方程15
y t = C 0 x t + [ D 0 C ρ D ρ ] u t ρ t ⊗ x t ρ t ⊗ u t 方程16
如下面更详细的论述,上述方程可以表示为具有非线性反馈的线性时不变(LTI)系统,其中状态xt和输入ut乘以时变调度参数ρt。现在反馈ft输入可被认为是该LTI系统的实际输入。如下面进一步详细的论述,LPV系统描述的矩阵[AρBρ;CρDρ]是描述LPV系统的LTI反馈表示的恰当量。
并且,上述方程现在可表示为如下面所示的方程17和18。
方程17
方程18
因此,对于测量输出和输入 u ~ t = u t T ( ρ t ⊗ x t ) T ( ρ t ⊗ u t ) T T , 时不变矩阵可以分别为并且,在中的xt不是已知或测得的量的情况下,先前的估计xt是可用的或可以被使用,并且可以使用迭代来获得更精确的或期望的估计xt
在更进一步的示例性的实施例中,LPV系统可表示为具有非线性内反馈的线性时不变系统,该反馈可以包括已知的变参数函数ρt(见《非线性系统识别》部分2.1:状态空间法,由VincentVerdult著,荷兰特文特大学2002年博士论文,其内容通过引用整体并入本申请中)。在该示例性的实施例中,对i(1≤i≤s)使用奇异值分解将具有秩ri的系统矩阵Pi分解为如方程19所示。
P i = A i B i C i D i = B f , i D f , i [ C z , i D z , i ] 方程19
则参量可被定义如下,如方程20-23所示。
Bf=[Bf,1Bf,2…Bf,s]方程20
Df=[Df,1Df,2…Df,s]方程21
C z T = [ C z , 1 T C z , 2 T ... C z , s T ] 方程22
D z T = [ D z , 1 T D z , 2 T ... D z , s T ] 方程23
接下来,LTI系统P0=[A0B0;C0D0]中的内反馈被假设具有输出zt=Czxt+Dzut和将zt输入其中的非线性反馈该非线性反馈通过矩阵Bf和Df输入到LTI系统P0状态方程。反馈系统的xt+1和yt的状态方程可以如下面的方程24至27所示。
xt=A0xt+B0ut+Bfft方程24
yt=C0xt+D0ut+Dfft方程25
zt=Czxt+Dzut方程26
f t = [ ρ t ⊗ z t ] 方程27
关于上述现在来参考示例性的图4,如图4中显示的反馈结构以及方程26中,如果假设对于任何时刻t反馈ft不会对输出zt产生任何效应,则线性分式变换(LFT)技术中不存在参数相关性(见K.Zhou,J.Doyle和K.Glover(1996)著,普伦蒂斯霍尔出版社出版的《鲁棒与最优控制》,部分10.2,其内容通过引用整体结合在本文中)。如图4中所示,存在具有两个方框(404和402)的系统400,其中方框404可以是较少存储容量的非线性反馈系统,而方框402可以是线性时不变系统。因此,在一些示例性的实施例中,参数相关性可变为仿射的。
示例性的图4可以是方程15和16的示意图。状态方程15涉及402中的上部的方框,而测量方程16涉及402中的下部的方框。ΔT422是采样持续的延时,其中方程15的右边在444处进入422,并且在离开422时方程15的左边等于xt+1。这是与方程14和15类似的递推公式,从而对于该图,时间指标在下一迭代开始之前从“t”变为“t+1”,继续迭代直到进入方框420、430和410。调度参数ρt406、输入ut408和输出yt446是变量。上部的四个方框表示从左到右分别乘以B0418、A0420、Bρ414和Aρ416。类似地,用C代替A和用D代替B,下部的四个方框表示分别乘以Dρ424、Cρ426、D0428和C0430。在反馈方框404中,ρt与xt的克罗内克积以及ρt与ut的克罗内克积分别在410和412中形成。在垂直方向对准的方框对从左到右分别乘以变量xt、ut另外,所示的压线的箭头表示加法。并且,ΔT422表示持续时间的延时块ΔT,对于该示例性的实施例,该延时类似日界线(adateline)那样作用。因此,示例性的图4中的箭头表示操作的时间流程(atimeflow)或实际顺序,该时间流程在406、408以及422的输出处开始,继续进行执行该示意图。当到达ΔT422时,已经执行了抽样t的所有操作,并且当通过ΔT422时,抽样时间t+1开始。因此,例如,当离开ΔT422时,保持同样的量,但是所有的时间标志都变为t+1,且贯穿示例性的图4所示的进程。
并且,如下面所示,这等价于方程15和16所示的LPV形式,其中状态方程关于调度参数向量p是线性的。
在进一步示例性的实施例中且现在使用定义为的反馈定义(如方程26和27所指出的),xt+1和yt的状态方程如下面的方程28和29所示。
x t + 1 = A 0 x t + B 0 u t + Σ i = 1 s B f , i ρ t ( i ) C z , i x t + Σ i = 1 s B f , i ρ t ( i ) D z , i u k 方程28
y t = C 0 x t + D 0 u t + Σ i = 1 s D f , i ρ t ( i ) C z , i x t + Σ i = 1 s D f , i ρ t ( i ) D z , i u k 方程29
然后,如果方程19被用于定义上述因子,则方程28和29与方程15和16相同。
接下来,在进一步示例性的实施例中,当Pi的秩不是已知的,则输出zt被设置为状态和输入接着通过静态的非线性被反馈,从而有然后,[Cz,iDz,i]=[Idimx0;0Idimu]以及方程30和31如下设置。
P = A i B i C i D i = B f , i D f , i [ C z , i D z , i ] 方程30
= B f , i D f , i I dim x 0 0 I dim x = B f , i D f , i 方程31
因此,从上述可知,LPV系数矩阵Pi=[AiBi;CiDi]是左边状态方程变量(xt+1;yt)对非线性项的向量[ρi,txt;ρi,tut]的回归矩阵。通常,LPV系数矩阵Pρ=[AρBρ;CρDρ]=[Bf;Df]是左边状态方程变量(xt+1;yt)对非线性反馈项的向量的回归矩阵。
LTI非线性反馈表述可以解决将现有的子空间识别算法应用到LPV系统的识别时产生的主要障碍,并且克服了先前的关于其它方法中非线性项的数量成指数增长的问题。例如,上述LTI非线性反馈表述可以清楚表明:非线性项可以理解为LTI非线性反馈系统的输入。因此,使用线性子空间方法来直接估计LTI系统状态空间方程的矩阵是可能的,这种估计对于具有输入和反馈的过程是精确的。该方法可直接包括LTI非线性反馈系统的输出yt和增量输入的使用。
在另一示例性的实施例中,可以在减少了非线性反馈系统的LTI子系统之后确定LTI系统矩阵和状态向量,该非线性反馈系统包括已知的调度函数和LTI子系统的状态。该实施例可以包括对LTI系统状态以及描述该LTI系统的LTI状态空间矩阵进行迭代测量。
一个例子可以是将多项式系统认为是关于x和u的线性系统,该系统还具有另外增加的更高阶积的项的输入项,那么,另外的输入是假设调度变量ρt实时是可用的,如操作点或测量的变量。如果状态xt的精确估计也是可用的,那么,问题就只是对系统识别直接应用迭代算法。由于变量xt在解决了系统识别之后才是可用的,因此,可以使用不同的方法。
因此,在示例性的第一步骤中,对状态向量进行初始估计。这里,对状态方程中的项,包括变量xt、ut但不包括变量进行系统识别。由此,可以获得系统线性时不变(LTI)部分的近似值,该系统给出A0、B0、C0、D0、Bρ和Dρ的估计以及状态向量的估计。
然后,在示例性的第二步骤中,对状态向量进行迭代估计。这里,状态向量可用作方程15和16中的项中xt的初始估计。然后,应用迭代算法来获得系统矩阵A、B、C和D的估计以及的精确估计(refinedestimate)。并且,接着重复该步骤直到达到期望的收敛。
因而,上述示例性的步骤1至3只进行几次迭代。因此,以示例性的方式,迭代算法被用于解决先前公知的关于LPV系统识别的问题。因此,接下来是使用迭代算法分别直接确定非线性差分方程函数f(xt,ut,vt)和h(xt,ut,vt)的增加的多项式展开的系数Fij和Hij的示例性论述。对于这样的非线性系统来说,这是非常简洁和精简的参数。因而,可以使用这里描述的线性时不变系统的迭代算法,其在计算要求方面只具有很少的增加。并且,示例性使用该迭代算法来直接处理状态方程中增加的非线性项,该状态方程包括作为另外的输入的状态向量,如由于状态xt不是初始已知的,因此这可有助于进行迭代来估计只以系统的线性时不变(LTI)部分A0、B0、C0和D0开始的状态序列。
如方程13和14之后的上述示例性论述,关于调度变量ρt是仿射的线性变参数系统可以用包括另外的输入变量的时不变形式表示。注意到的是,这包括关于xt和ut的非线性函数ρt。动态系统在变量的这些非线性函数中是线性时不变的。
在进一步示例性的实施例中,可以通过迭代算法来追踪另外的输入产生的效应(effect)。下面概述的示例性的步骤可进一步反应在示例性的图5所示的表格以及示例性的图6的流程图中。
使用示例性的图5中的单元502且在示例性的图6的步骤602中,拟合增加阶的ARX模型以观察关于输出yt和输入的迭代k,其中对于k=1,项不存在输入中且进一步可等于从LPV模型结构中移除关于(AρCρ)的项。并且,在图6的步骤602中,计算AIC来确定最优阶。ARX模型拟合可具有线性回归(linearregression)问题:除估计的最大ARX阶以外,该问题不对ARX阶做先前的假设。如果确定的ARX阶接近估计的最大值,则估计的最大ARX阶被加倍且再次计算该回归。
在示例性的步骤604中,可以计算校正的新值。可以使用ARX模型来确定新值输入对新值输出产生的效应,且从输出中减去该效应。该效应可以计算新值输出,如果在新值时刻t(inthefutureoftimet)不存在输入,则可从在时刻t的状态获得该新值输出。
在示例性的步骤606中,可以执行或计算规范变量分析(CVA)。这里,可以计算旧值和校正的新值之间的CVA。此外,还可计算旧值和校正新值之间的协方差矩阵。这类似于对联合的(joint)旧值-新值协方差矩阵进行SVD以获得ARX模型,旧值-新值协方差矩阵具有旧值协方差的阶。该步骤的结果是获得所谓‘存储值(memory)’的系统状态的估计。
在示例性的步骤608中,使用状态方程来执行回归。步骤606中的‘存储值’好像数据一样可用于状态方程中,并且可以获得状态空间矩阵和噪声过程的协方差矩阵产生的估计。这些估计是状态空间模型的渐近ML的估计,而未对ARX或SS模型的参数值进行先前的假设。
因此,使用上面所描述的方法,可基于迭代k(如示例性的图5的504所示)中假设的输出和输入,在一次计算中获得迭代算法中的ML解,在迭代k中不对假设的参数值进行迭代。该迭代是非线性项中状态估计的精确结果,该项是在迭代k(504)中假设的数据的一部分。
回来参考步骤602,ARX模型拟合,增量输出的维数从dimu增加到dimua=dimu+dimp(dimx+dimu)。为了拟合ARX模型,由于非线性输入项并且根据存在于输出和增量输入变量之中的统计学显著性动态,确定的ARX阶lagp实质上更高。计算包括对维数为dimua*lagp的数据协方差矩阵计算SVD,其中lagp是估计的(considered)最大ARX阶。用于SVD的计算具有阶60*(dimua*lagp)3,因此,该计算正比于(dimp(dimx+dimu)/dimu)3增加。
因此,系统输入增加非线性项的一个结果是使旧值增长了因子dimp(dimx+dimu)/dimu,且使计算增长了该因子的立方倍。根据ut、xt和ρt的特定维数,这是非常重要的,然而,在项的数量或计算时间上不成指数增长。本发明的LPV子空间算法仍然对应于线性时不变系统的子空间系统识别,并且此外,由于状态估计Xt包含的项的非线性,对系统状态的估计进行迭代直到收敛是可期望的或者(在一些选项中)所需求的。
另一需要考虑的因子是,示例性的步骤606中CVA的结果是状态序列xt的估计的计算值,表示为mt。符号‘m’用作单词‘存储值’来使其区别于实际状态向量和如下面讨论的EM算法中所使用的状态向量的各种估计。结合EM算法的极大值步骤的估计可被示出为线性时不变系统情况下的系统识别问题的最大似然解。在这种情况下,只使用该算法的一次迭代就可获得全局(global)ML解。这与EM和梯度搜索算法不同,这些算法基本上总是使用许多次迭代。CVA估计mt实际上是系统状态序列的估计,该系统具有等于LTI系统和较大采样大小(largesamplesize)情况下最大似然估计的参数。这与如下通常的概念不同:首先获得ML参数估计,然后使用卡尔曼滤波器来估计该状态。该算法在一次迭代中不仅获得了ML参数的最优状态序列,而且在该过程中还确定了最优状态阶。在通常的ML方法中,期望的是或(在一些选项中)需要的是获得状态阶的每种选择的ML估计,然后继续进行假设测试以确定最优状态阶。
在另一示例性的实施例中,描述了LPV情况的迭代算法的收敛性。在该实施例中,对非线性差分方程的其它形式采用基本类似的方法。
对于任何迭代算法,通常会产生两个问题:(1)该算法是否总是收敛;和(2)以什么速率收敛。在估计的计算实例中,对估计的状态进行的迭代k可以非常快,例如六个步骤。因此,这里显示的是,假设LPV系统是稳定的情况下,迭代算法与显示为总是收敛的EM算法类(theclassof)紧密相关。并且,收敛的速率被计算成是几何学的。后者的结果是唯一的,因为EM算法通常在前面进程非常快但在结束时非常慢。将更详细地描述LPV算法的快速末端收敛(terminalconvergence)的原因。下面将详细描述初始化、稳定性和收敛的问题。
另外,可在EM算法的背景下论述这里所描述的方法和系统,因为这两者之间存在一些平行性(parallelism)。为了显示LPV算法的收敛性,与适于LPV算法的各种变形一起来讨论Gibson和Ninness中的推导,下面表示为GN且通过引用并入本申请,(S.Gibson和B.Ninness,“多变量动态系统的鲁棒极大似然估计”,自动化,第41卷,第5号,第1667-1682页,2005)。Gibson和Ninness(2005)中的所有方程的数目将包括文件GN中的数量之后的GN,例如(23GN)。
在该示例性的实施例中,对GN公式进行如下替换以获得LPV算法:用LPV状态方程替换LTI状态方程,以及对于丢失的数据,用迭代算法中的‘存储值’向量mt替换卡尔曼平滑中的状态估计。这个结果是很重要的,因为对于线性系统,当在GN中,这里所描述的迭代算法可以在较大的采样中用一个步骤获得全局ML参数估计。另一方面,对于线性系统,EM算法需要进行多次迭代来获得ML解。
用LPV模型替换LTI模型就是用LPV方程(15)和(16)替换状态空间模型方程11GN。这会产生定理3.1的方程的多种变形,因为GN中的ut和(A,B,C,D)分别用替换。因此‘数据’包括向量其中状态向量xt不是可用的。
为了执行示例性的图5中504所示的迭代k,可以使用关于EM算法的最直接的方法(straightforwardapproach)。‘丢失的数据’包括作为丢失数据的状态向量xt。(见S.Gibson,A.Wills和B.Ninness(2005),“双线性系统的极大似然参数估计”,IEEETrans.自动控制,第50卷,第10号,第1581-2005页,其内容通过引用整体并入本申请)。然后,我们可以继续讨论双线性系统,该双线性系统包括与LPV算法的完全类似的项但是这种EM方法有时会导致在接近最大值时EM算法通常慢收敛的现象。
因此,相反,可以使用LPV算法,且该算法会产生快速收敛。因此,在这个示例性的实施例中,子空间方法不是将‘丢失的数据’指定为卡尔曼平滑的估计,而是将CVA状态估计mt或‘存储值’指定为丢失的数据。该存储值mt是通过校正新值和旧值之间的规范变量分析所获得的状态向量的估计,校正新值和旧值是在迭代算法的示例性的步骤606中使用图5中描述的输入和输出向量来获得的。这类似于对与迭代k处输出和输入数据相关的全局ML参数估计进行卡尔曼滤波器状态估计,而不是对最后估计的参数值进行卡尔曼平滑状态估计。这会产生一种区别,因为CVA方法根据基于旧值的校正新值来表示似然函数,从而存储值估计mt只取决于输出和输入数据,且它们的分配取决于与用于迭代k中所使用的输出和输入数据相关的ML估计,而不是该状态的平滑估计。(26GN)的实际条件似然函数p0t/zt)是包括在所有的EM计算中的条件似然函数,其中这是与GN的定理3.3中的似然函数相同的似然函数,该似然函数包括在估计状态空间方程的参数的CVA算法的示例性的步骤608中。区别在于在CVA算法的示例性的步骤608中,期望是关于与迭代k处输出和输入数据相关的真实的全局ML估计,而GN估计是关于在先前的迭代中获得的参数值的期望。
并且,在该算法的期望步骤中,使用LPV模型和将存储值mt选作丢失的数据会产生GN中的如下结果。不需要修改部分2,即GN的最大期望(EM)算法中的基础理论。在GN的部分3中,丢失的数据被认为是基于图5的504中迭代k处输入和输出量的CVA状态估计因此,方程(22GN)至(28GN)中的xt可以用替换。
因此,保留GN的定理3.1,但是还可在一个步骤中实现与示例性的图5的输入-输出向量相关的全局ML估计。GN的定理3.2可以用迭代算法替换以获得存储值估计GN的定理3.3是与在该迭代算法中获得的结果相同的结果。另一个步骤可用于根据估计Θ[k]和LPV状态方程给出的线性时变状态方程来计算该步骤是期望获得用于开始下一迭代k1的状态估计该步骤可进一步提供迭代算法以产生Q(θ,θ′)≥Q(θ′,θ′),当且仅当时取等号。
然后,原理上,存储值可用来替换状态估计实际上,在方程(43GN)和(44GN)中状态方程的递归结构上输入(projected)的可以产生渐近ML状态估计和对应的最优滤波状态估计因此,在计算时将(而不是)用作‘输入数据’的一部分会产生基本相同的结果,除了结果中的一些‘噪声’外。但是可计算的噪声可以通过的另外的小的计算来避免。
并且,在一些示例性的实施例中,不需要对这里描述的迭代算法进行‘鲁棒’改进,因为使用主要奇异值分解计算已被推导为是稳健的(robust)并且十多年来也被证实为是如此。如果病态(ill-conditioned)的LPV动态系统将被解决以实现高精度,可能的例外是滤波的状态估计的计算,该计算可以使用Bierman的平方根方法(G.J.Bierman,离散序列估计的因子分解方法,学术出版社(1977);Dover再版,纽约(2006),其内容通过引用整体并入本申请)来执行。
在另一示例性的实施例中,可证实的是,LPV系统识别算法可以在接近极大似然解处以几何速率收敛。这里,可以推导出线性系统的结果。LPV系统可以使用同样的方法,但是下面的描述是取决于时间的且具有更大的复杂性。
并且,为了简化这里的推导,方程17和18中的时不变反馈可以通过用替换来估计。因此,具有噪声vt的方程17和18用新形式写成如下面的方程32和33所示。
xt+1=Axt+But+Kvt方程32
yt=Cxt+Dut+vt方程33
因此,从上述内容可以看出,在测量输出yt和输入 u t = u ~ t T ( ρ t ⊗ x t ) T ( ρ t ⊗ u t ) T T 时,时不变矩阵分别为(A,B,C,D)=(A0,[B0,Aρ,Bρ],C0,[B0,Cρ,Dρ])。然后,求解方程33中的vt并将其代入方程32中,则产生如下方程34。
xt+1=Axt+But+K(yt-Cxt+Dut)=(A-KC)xt+(B-KD)ut+Kyt方程34
接下来,通过递归将xt代入方程34的右边,可以产生如下方程35。
x t = Σ i = 1 ∞ ( A - K C ) i - 1 [ ( B - K D ) u t - i + Ky t - i ] = Jp t ( y t , [ u ~ t ( ρ t ⊗ x t ) , ( ρ t ⊗ u ~ t ) ] ) 方程35
在方程35中,J包括ARX系数,而ρt是旧值输出yt并且输入这里,是初始输入,从而ut包括非线性克罗内克积的项。
在进一步示例性的实施例中,渐近地对于较大的采样,进行足够大的迭代k时,J任意地接近常数,从而
Δx t [ k ] = x t [ k + 1 ] - x t [ k ] = Jp t [ 0 , 0 , ( ρ r ⊗ ( x t [ k ] - x t [ k - 1 ] ) ) 0 ] = L ρ t - 1 ⊗ Δx t - 1 [ k - 1 ] . . . ρ t - log ⊗ Δx t - log [ k - 1 ] = L t Δx t - 1 [ k - 1 ] . . . Δx t - log [ k - 1 ] ,
其中,Lt是随时间变化的且具有与项L组合在一起的时变调度参数ρt。并且,如果对于所有t,ρt都是有界的,则Lt类似地也是有界的。
现在,当用块向量的形式将状态方程写成时,其中vec操作是将矩阵的列项(columns)从左边列开始堆栈(stack)在上方,上述结果表明其中→表示LN-k延伸到右边。使用M来表示上三角矩阵可以为连续的迭代k和k-1处的状态序列之间的差值提供基础表述,如下面的方程36。
Δ X lag + 2 : N i = MΔ X lag + 1 : N - 1 i - 1 方程36
根据方程36可以看出,LPV情况下迭代的末端收敛速率由M控制,特别是M的最大奇异值。因而M的各种块可通过 ( A - KC ) i - 1 ( B - KD ) ρ ⊗ x ρ N - 1 ( 1 ) I dim x . . . ρ N - i ( dim p ) I dim x 来计算,其中下标表示选择B-KD的子矩阵,具有与中相应行对应的列。
在一些其它示例性的实施例中,迭代线性子空间计算的收敛性受LPV差分方程稳定性的影响,且更具体地,是受这里描述的LPV线性子空间系统识别的稳定性的影响。因为当且仅当状态转移矩阵的所有的特征值是稳定的,例如,所有特征值都小于1时,一组时不变线性状态空间差分方程才可能是稳定的。LPV情况更复杂,但是为了该示例性的实施例的目的,对于状态向量xt的每个特征向量的分量来说,每采样时间的增长或收缩的速率可以根据方程9由LPV状态转换矩阵的相应的特征值给出,现在如下面的方程37所示。
A ( ρ t ) = A 0 + Σ i = 1 s ρ t ( i ) A i 方程37
在方程37中,假设矩阵A=(A0Aρ)=(A0Al…As)是未知的常数矩阵。因此,明显的是,转换矩阵A(ρt)是矩阵Ai(0≤i≤s)的线性组合[1;ρt]。因此,对于矩阵Ai(0≤i≤s)的任何选择,存在可能的值ρt,该值ρt可在特定的采样时间t产生不稳定的奇异值。如果或当值ρt只是零星地产生和/或只是在有限个连续观测时间产生,则不会出现问题或产生误差。
在一些其它示例性的实施例中,例如在子空间算法的第k次迭代计算中,只可以使用状态序列的估计值和矩阵因此,如果这些数量的估计中的较大误差是允许的或可接受的,则在很大程度上存在不稳定的现象。例如,存在更有意义的重要性的应用领域,如包括飞机机翼颤振动力学识别,其中颤振动力学是边缘稳定的或甚至不稳定的,可以通过从传感器到机翼控制表面的主动控制反馈(activecontrolfeedback)以使颤振稳定。然后,在一些其它的应用中,可能保证的是:并不是调度参数ρt、矩阵Ai的值和不确定的任何组合都会产生不稳定。进一步地关于这方面,转换矩阵的稳定性可以提供问题的潜在可能,因此需要进一步的估计。
在一些其它示例性的实施例中,不稳定性可以产生这样的现象:预测的系统响应可快速地大量增长。例如,当奇异值都是有界的,且都小于1,则预测的系统响应也是有界的,但是如果估计误差很大,例如一段时间的一个奇异值等于2,则预测的响应在该时间内大约加倍每个采样。因此,如果1030=(103)10=(210)10=2100,则在如100个这样少的采样中,在计算时会损失30位(digits)的精度,这会产生毫无无意义的结果。
因此,根据上述论述,存在这样的情况,在该情况下,与如下时间段(periods)有关的数值精度会产生相当大的损失:在该时间段中,LPV转换矩阵是不稳定的,例如具有延长的时间间隔。并且,困难在于算法初始化和采样大小,因为在具有较大采样的最优解处,该算法是稳定的且收敛。如果用无限的精度进行计算是可能的,则可以避免病态的问题;但是,具有小数点后15或30位的精度,例如,飞机机翼颤振的一些实际数据组,该计算就可以从进一步的估计中获益。
因此,一些示例性的实施例可以如下方式进行处理:校正或者减少由算法不稳定导致的任何不需要的效应。例如,当根据先前提到的末端收敛结果时,如果状态序列充分接近最优值,则该迭代算法是稳定的(如果假设LPV系统是稳定的)。因此,在一些例子中,估计中大的初始误差会导致不稳定的计算。
在其它的例子中,如在时间间隔期间,当调度参数值导致A(ρt)中存在相当大的不稳定时,ft|qt的一些项(components)会变得非常大,如被校正用于新值输入qt产生的效应的新值输出ft。因此,在ft|qt的协方差矩阵计算中,这样的大值会导致数值精度存在相当大的病态和损失。
在下面的例子中,迭代算法和一些其它子空间方法允许在产生不稳定时任意编辑时间,对不稳定的区域执行异常值(outlier)编辑。例如,在ARX协方差计算中,计算校正的新值,可确定具有相当大的不稳定的时间间隔,并且可以从计算中删除该时间间隔。
在处理实验设计的例子中,如果操作点变量的轨线是指定的,则调度参数被安排来以几种方式增强系统识别。选择能避免计算不稳定的初始区域来获得LPV参数(A,B,C,D)的足够精确的初始估计。然后,在其它区域可使用该初始估计对算法的状态估计进行初始化,该算法具有不会出现计算不稳定的足够精确度。
在上述例子中,在每次迭代时删除不稳定的异常值是最通常和稳健的步骤。当状态序列和LPV参数的估计的值随着更多迭代改进时,可以期望减少异常值的数量直到具有快速末端收敛。该期望的反例是,当调度参数ρt时程(timehistory)的开始具有较小的变化,从而该数据这部分的LPV模型对于该数据这部分是有利的,然而却是差的全局模型。然后,在时程的后面部分,ρt中可存在估计变化,从而产生不稳定的现象。
并且,在这里描述的方法和系统的许多潜在的示例性的应用中,期望提议的算法要比目前工业问题方面不可行的现有方法表现得更好。并且,在许多情况下,期望设计实验来以尽可能少的时间和资源获得操作空间的指定全局区域的期望保真度的结果。由于迭代算法的线性子空间方法是基于极大似然的方法,因此可以改进设计用于LTI系统识别。并且,当用估计的干扰模型确定随机模型时,该干扰模型包括如动态频响函数的置信水平,可以根据干扰过程上的较少的先前信息来推导需要的采样大小和系统输入激励。
在一些进一步示例性的实施例中,这里描述的LPV方法和系统可以延伸到非线性系统。例如,示出的是,许多复杂的非线性系统可以用近似LPV的形式表示,该LPV形式充分应用了这里描述的LPV子空间系统识别方法。
在一个例子中,通常的非线性、时变、变参数动态系统可以用非线性状态方程的系统描述,如方程38和39所示的这些方程。
xt+1=f(xt,ut,ρt,vt)方程38
yt=h(xt,ut,ρt,vt)方程39
在方程38和39中,xt是状态向量,ut是输入向量,yt是输出向量且vt是白噪声测量向量。在一些示例性的实施例中,为了处理‘变参数’系统,时变参数的‘调度’变量ρt可以描述系统的当前操作点。非常普遍的函数类f(·)和h(·)可以用增加的(additive)不需要连续的波莱尔函数(borelfunctions)来表示。
以简化的方式,采用如Rugh,部分6.3(W.J.Rugh著,非线性系统原理:Volterra/Wiener方法,由马里兰州巴尔的摩的约翰·霍普金斯大学出版社1981出版,其内容通过引用整体并入本申请)中泰勒展开的函数的情况可以提供如下方程,如方程40和41所示,其中不存在ρt和vt
x t + 1 = Σ i = 0 I Σ j = 0 J F i j x t ( i ) u t j 方程40
y t = Σ i = 0 I Σ j = 0 J H i j x t ( i ) u t j 方程41
在方程38和39中,符号可以递归地定义为且类似地定义
因此,方程40和41是非线性函数f(·)和h(·)的多项式展开。注意到的是,非线性方程可以包括相对简单的形式的非线性函数,如只包括容易计算的用于各种目的积之和的近似多项式方程。然而,对于存在自相关误差时系数的经验估计,对于低维数的y、u和x,问题会变得困难,即使使用了子空间方法。对于子空间方法,矩阵的维数会随着可以使用的‘旧值’的维数成指数增长。这发生在展开方程40中,重复将方程40左边的xt代入状态方程40右边的中。因此,用t-1替换t的方程40的整个右边可以提高到幂(power)lagp,即通常选作估计的ARX模型阶的旧值的阶。对于相对低阶的旧值以及xt和ut的低维数,增加的项的数量成指数增长。
然而,在进一步示例性的实施例和下面的Rugh中,方程39和40可通过卡莱曼双线性转换成双线性向量差分方程,其中状态变量如方程42所示:
方程42
且输入幂和乘积变量如方程43所示。
方程43
在进一步示例性的实施例中,表示为状态-仿射的形式的方程40可写成如下面的方程44和45所示。
x t ⊗ = A ( x t ⊗ ⊗ u t ⊗ ) + Bu t ⊗ 方程44
x t ⊗ = C ( x t ⊗ ⊗ u t ⊗ ) + Du t ⊗ 方程45
接下来,假设系数矩阵(A,B,C,D)是调度参数ρt的线性函数,表示为(A(ρt),B(ρt),C(ρt)D(ρt))且是如方程9至12的状态-仿射的形式。调度参数ρt是操作点或者其它已知或精确测量的变量的非线性函数。例如,由于输入是乘法的且可实时假设为已知的或精确测量的,因此它们可用于调度参数ρt中,因而可能会减少它们的维数。双线性方程变成下面的方程46和47所示。
x t ⊗ = A ( ρ t ) x t ⊗ + B ( ρ t ) u t ⊗ 方程46
x t ⊗ = C ( ρ t ) x t ⊗ + D ( ρ t ) u t ⊗ 方程47
方程46和47是明显的LPV形式。因此,对于充分平滑的函数f(xt,uu,ρt,vt)和h(xt,uu,ρt,vt),这些函数在变量附近(neighborhood)存在明确定义的幂级数展开,可存在过程的LPV近似。
因此,将输入代入调度参数ρt中,方程可以是线性时变的且在一些例子中可以具有更高的维数ρt。这可以显示输入和调度参数的双功能,以及它们如何以各种方式相互交换,例如,它们的差值是随时间变化的速率(rate)。因此调度的知识可以从基础物理、化学或其它基础知识中获得。因此,当灰匣子(graybox)模型在如下范围内:它们能够结合关于过程的行为的相当大的全局信息,该过程可以并入过程动态如何随操作点改变的模型中时,有时可以估计LPV模型。
前面的描述和附图阐述了本发明的原理、优选的实施例和操作模式。然而,本发明不应认为限于上面描述的特定实施例。本领域技术人员将会知道上述实施例的另外的变形。
因此,上述实施例应该被认为是例证性的而不是限制性的。因此,应该理解的是,本领域技术人员可以对这些实施例进行变形,只要其不背离下面的权利要求限定的本发明的范围。

Claims (2)

1.一种通过将子空间识别方法用于具有一般调度函数的线性变参数(LPV)和非线性变参数系统以使飞机机翼颤振稳定的方法,包括:
通过传感器检测飞机机翼颤振;
将检测到的机翼颤振数据存储于存储器中;
对所存储的机翼颤振数据执行增加阶的自回归(ARX)模型拟合,所述拟合包括:计算AIC以确定最优阶,和选择具有最优阶的ARX模型;
处理器执行校正的新值计算,所述校正的新值计算通过使用ARX模型来确定一个或者多个机翼颤振数据新值输入对新值输出的效应,和从新值输出中减去所述效应来确定一个或多个校正后的机翼颤振数据新值输出;
处理器通过执行校正后的新值输出和旧值之间的规范变量分析(CVA)来确定状态估计;
通过AIC来确定状态阶;
将状态估计代入一个或多个状态方程;
对所述一个或多个状态方程执行线性回归运算以确定状态方程的系数矩阵;以及
用具有上述系数矩阵的状态方程的形式提供飞机机翼颤振动态模型,以将子空间识别方法用于具有一般调度函数的LPV和非线性变参数系统。
2.根据权利要求1所述的方法,其中,ARX模型拟合具有假定估计的最大ARX阶的线性回归问题。
CN201080050251.0A 2009-09-03 2010-09-03 通过迭代线性子空间计算对时变、变参数和非线性系统进行实验建模的方法和系统 Active CN102667755B (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US23974509P 2009-09-03 2009-09-03
US61/239,745 2009-09-03
PCT/US2010/047832 WO2011029015A2 (en) 2009-09-03 2010-09-03 Method and system for empirical modeling of time-varying, parameter-varying, and nonlinear systems via iterative linear subspace computation

Publications (2)

Publication Number Publication Date
CN102667755A CN102667755A (zh) 2012-09-12
CN102667755B true CN102667755B (zh) 2016-08-03

Family

ID=43626135

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201080050251.0A Active CN102667755B (zh) 2009-09-03 2010-09-03 通过迭代线性子空间计算对时变、变参数和非线性系统进行实验建模的方法和系统

Country Status (7)

Country Link
US (2) US8898040B2 (zh)
EP (1) EP2465050A4 (zh)
JP (2) JP2013504133A (zh)
KR (1) KR101894288B1 (zh)
CN (1) CN102667755B (zh)
CA (1) CA2771583C (zh)
WO (1) WO2011029015A2 (zh)

Families Citing this family (37)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2465050A4 (en) * 2009-09-03 2014-01-29 Wallace E Larimore METHOD AND SYSTEM FOR THE EMPIRICAL MODELING OF NONLINEAR SYSTEMS WITH TIME AND PARAMETER VARIATIONS THROUGH REPEATED LINEAR SURFACE CALCULATION
US20120330717A1 (en) * 2011-06-24 2012-12-27 Oracle International Corporation Retail forecasting using parameter estimation
US20140278303A1 (en) * 2013-03-15 2014-09-18 Wallace LARIMORE Method and system of dynamic model identification for monitoring and control of dynamic machines with variable structure or variable operation conditions
AT512977B1 (de) * 2013-05-22 2014-12-15 Avl List Gmbh Methode zur Ermittlung eines Modells einer Ausgangsgröße eines technischen Systems
KR20160021209A (ko) * 2013-06-14 2016-02-24 월러스 이. 래리모어 동적 모델 식별을 모니터링 및 가변 구조 또는 가변 동작 조건을 가진 동적 기기의 제어 방법 및 시스템
JP6239294B2 (ja) * 2013-07-18 2017-11-29 株式会社日立ハイテクノロジーズ プラズマ処理装置及びプラズマ処理装置の運転方法
WO2016194025A1 (ja) * 2015-06-02 2016-12-08 日本電気株式会社 線形パラメータ変動モデル推定システム、方法およびプログラム
US10634580B2 (en) * 2015-06-04 2020-04-28 The Boeing Company Systems and methods for analyzing flutter test data using damped sine curve fitting with the closed form shape fit
US10552408B2 (en) * 2016-11-02 2020-02-04 Oracle International Corporation Automatic linearizability checking of operations on concurrent data structures
CN106845016B (zh) * 2017-02-24 2019-10-15 西北工业大学 一种基于事件驱动的量测调度方法
CN106934124B (zh) * 2017-02-24 2020-02-21 西北工业大学 一种基于量测变化检测的自适应变划窗方法
CN107122611A (zh) * 2017-04-28 2017-09-01 中国石油大学(华东) 青霉素发酵过程质量相关故障检测方法
CN107065557B (zh) * 2017-05-04 2020-04-28 南京理工大学 一种具有随机滤波增益变化的有限域滤波器设计方法
CN107168053B (zh) * 2017-05-04 2020-10-30 南京理工大学 一种具有随机滤波增益变化的有限域滤波器设计方法
CN107194110A (zh) * 2017-06-13 2017-09-22 哈尔滨工业大学 一种线性变参数双率系统的全局鲁棒参数辨识及输出估计方法
US10344615B2 (en) 2017-06-22 2019-07-09 General Electric Company Method and system for schedule predictive lead compensation
US20190271608A1 (en) * 2018-03-01 2019-09-05 GM Global Technology Operations LLC Method to estimate compressor inlet pressure for a turbocharger
CN108764680B (zh) * 2018-05-18 2022-04-19 秦皇岛港股份有限公司 一种用于港前待泊区的船舶调度方法
CN109033021B (zh) * 2018-07-20 2021-07-20 华南理工大学 一种基于变参收敛神经网络的线性方程求解器设计方法
WO2020118512A1 (zh) * 2018-12-11 2020-06-18 大连理工大学 一种基于lft的航空发动机传感器及执行机构故障诊断方法
CN109840069B (zh) * 2019-03-12 2021-04-09 烟台职业学院 一种改进的自适应快速迭代收敛解方法及系统
CN110298060B (zh) * 2019-04-30 2023-04-07 哈尔滨工程大学 一种基于改进自适应遗传算法的间冷燃气轮机状态空间模型辨识方法
CN110513198B (zh) * 2019-08-13 2021-07-06 大连理工大学 一种涡扇发动机控制系统主动容错控制方法
CN110597230B (zh) * 2019-09-24 2020-08-18 清华大学深圳国际研究生院 一种主动故障诊断方法、计算机可读存储介质、控制方法
JP7107339B2 (ja) * 2019-09-30 2022-07-27 株式会社村田製作所 非線形特性計算の方法、非線形特性計算プログラムとその使用方法、並びに記録媒体
KR102429079B1 (ko) 2019-12-23 2022-08-03 주식회사 히타치하이테크 플라스마 처리 방법 및 플라스마 처리에 이용하는 파장 선택 방법
CN111324852B (zh) * 2020-03-06 2020-11-24 常熟理工学院 基于状态滤波和参数估计的cstr反应器时延系统的方法
CN111538237B (zh) * 2020-03-20 2021-10-29 北京航空航天大学 一种倾转旋翼无人机非线性浅灰模型辨识与校正方法
CN112199856B (zh) * 2020-10-23 2023-03-17 中国核动力研究设计院 基于modelica核反应堆管路系统模型构建与强耦合方法及装置
CN112487618B (zh) * 2020-11-19 2024-06-21 华北电力大学 基于等值信息交换的分布式抗差状态估计方法
CN112610330B (zh) * 2020-12-08 2023-05-09 孚创动力控制技术(启东)有限公司 一种内燃机运行状态的监测及分析系统和方法
CN112613252B (zh) * 2020-12-29 2024-04-05 大唐环境产业集团股份有限公司 一种吸收塔搅拌器的节能运行方法
CN113626983B (zh) * 2021-07-06 2022-09-13 南京理工大学 基于状态方程递推预测高炮射弹脱靶量的方法
CN113486523B (zh) * 2021-07-09 2024-06-11 南京航空航天大学 一种线性变参数振动系统全局辨识方法
CN113641100B (zh) * 2021-07-14 2023-11-28 苏州国科医工科技发展(集团)有限公司 针对未知非线性系统的通用辩识方法
CN113704981B (zh) * 2021-08-11 2024-04-09 武汉理工大学 温升过程中高速轴承时变动力学行为分析方法
CN113742936B (zh) * 2021-09-14 2024-04-30 贵州大学 基于函数型状态空间模型的复杂制造过程建模与预测方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1577339A (zh) * 1999-04-07 2005-02-09 凯登丝设计系统公司 用于对随时间变化系统和非线性系统建模的方法和系统

Family Cites Families (29)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH0250701A (ja) * 1988-08-12 1990-02-20 Fuji Electric Co Ltd 適応制御方法
JPH0695707A (ja) * 1992-09-11 1994-04-08 Toshiba Corp モデル予測制御装置
US5465321A (en) * 1993-04-07 1995-11-07 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Hidden markov models for fault detection in dynamic systems
JPH0962652A (ja) * 1995-08-29 1997-03-07 Nippon Telegr & Teleph Corp <Ntt> 最尤推定方法
US5991525A (en) * 1997-08-22 1999-11-23 Voyan Technology Method for real-time nonlinear system state estimation and control
US6173240B1 (en) * 1998-11-02 2001-01-09 Ise Integrated Systems Engineering Ag Multidimensional uncertainty analysis
US7493240B1 (en) * 2000-05-15 2009-02-17 Cadence Design Systems, Inc. Method and apparatus for simulating quasi-periodic circuit operating conditions using a mixed frequency/time algorithm
US7035782B2 (en) * 2000-06-02 2006-04-25 Cadence Design Systems, Inc. Method and device for multi-interval collocation for efficient high accuracy circuit simulation
JP4310043B2 (ja) * 2000-12-05 2009-08-05 本田技研工業株式会社 フラッタ試験用模型
US20020129038A1 (en) * 2000-12-18 2002-09-12 Cunningham Scott Woodroofe Gaussian mixture models in a data mining system
JP4184028B2 (ja) * 2001-10-12 2008-11-19 シャープ株式会社 流体・構造連成数値モデルの作成方法
AU2003213369A1 (en) * 2002-03-28 2003-10-13 Yuichi Nagahara Random number generation method based on multivariate non-normal distribution, parameter estimation method thereof, and application to simulation of financial field and semiconductor ion implantation
US7487078B1 (en) * 2002-12-20 2009-02-03 Cadence Design Systems, Inc. Method and system for modeling distributed time invariant systems
US20050021319A1 (en) * 2003-06-03 2005-01-27 Peng Li Methods, systems, and computer program products for modeling nonlinear systems
JP2005242581A (ja) * 2004-02-25 2005-09-08 Osaka Prefecture パラメータ推定方法、状態監視方法、パラメータ推定装置、状態監視装置及びコンピュータプログラム
US7133048B2 (en) * 2004-06-30 2006-11-07 Mitsubishi Electric Research Laboratories, Inc. Variable multilinear models for facial synthesis
US7184938B1 (en) * 2004-09-01 2007-02-27 Alereon, Inc. Method and system for statistical filters and design of statistical filters
US20070028220A1 (en) * 2004-10-15 2007-02-01 Xerox Corporation Fault detection and root cause identification in complex systems
JP2007245768A (ja) * 2006-03-13 2007-09-27 Nissan Motor Co Ltd ステアリング装置、自動車、及びステアリング制御方法
US7603185B2 (en) * 2006-09-14 2009-10-13 Honeywell International Inc. System for gain scheduling control
US20080126028A1 (en) * 2006-09-26 2008-05-29 Chang Gung University Method of reducing a multiple-inputs multiple-outputs (MIMO) interconnect circuit system in a global lanczos algorithm
FR2911202B1 (fr) * 2007-01-05 2009-02-13 Airbus France Sas Procede d'optimisation de panneaux raidis sous contraintes '
US7881815B2 (en) * 2007-07-12 2011-02-01 Honeywell International Inc. Method and system for process control
WO2010056736A2 (en) * 2008-11-11 2010-05-20 Axis Network Technology Ltd. Resource efficient adaptive digital pre-distortion system
EP2465050A4 (en) * 2009-09-03 2014-01-29 Wallace E Larimore METHOD AND SYSTEM FOR THE EMPIRICAL MODELING OF NONLINEAR SYSTEMS WITH TIME AND PARAMETER VARIATIONS THROUGH REPEATED LINEAR SURFACE CALCULATION
US8346693B2 (en) * 2009-11-24 2013-01-01 King Fahd University Of Petroleum And Minerals Method for hammerstein modeling of steam generator plant
US8903685B2 (en) * 2010-10-27 2014-12-02 King Fahd University Of Petroleum And Minerals Variable step-size least mean square method for estimation in adaptive networks
US9170572B2 (en) * 2011-07-06 2015-10-27 Honeywell International Inc. Dynamic model generation for implementing hybrid linear/non-linear controller
US8649884B2 (en) * 2011-07-27 2014-02-11 Honeywell International Inc. Integrated linear/non-linear hybrid process controller

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1577339A (zh) * 1999-04-07 2005-02-09 凯登丝设计系统公司 用于对随时间变化系统和非线性系统建模的方法和系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
机翼颤振的非线性动力学和控制研究;丁千 等;《科技导报》;20090131;第53-61页 *

Also Published As

Publication number Publication date
EP2465050A4 (en) 2014-01-29
WO2011029015A2 (en) 2011-03-10
KR20120092588A (ko) 2012-08-21
JP2013504133A (ja) 2013-02-04
JP2015181041A (ja) 2015-10-15
WO2011029015A3 (en) 2011-07-21
EP2465050A2 (en) 2012-06-20
US20110054863A1 (en) 2011-03-03
US20150039280A1 (en) 2015-02-05
US8898040B2 (en) 2014-11-25
CA2771583C (en) 2018-10-30
KR101894288B1 (ko) 2018-09-03
CA2771583A1 (en) 2011-03-10
CN102667755A (zh) 2012-09-12

Similar Documents

Publication Publication Date Title
CN102667755B (zh) 通过迭代线性子空间计算对时变、变参数和非线性系统进行实验建模的方法和系统
Gan et al. A variable projection approach for efficient estimation of RBF-ARX model
Peherstorfer et al. Online adaptive model reduction for nonlinear systems via low-rank updates
Kühl et al. A real-time algorithm for moving horizon state and parameter estimation
CN105814499B (zh) 用于监测和控制具有可变结构或可变运行条件的动态机器的动态模型认证方法和系统
Jouffroy et al. Finite-time simultaneous parameter and state estimation using modulating functions
Surana et al. Koopman operator framework for constrained state estimation
Al-Matouq et al. Multiple window moving horizon estimation
Calderón et al. Koopman operator-based model predictive control with recursive online update
Chen et al. Nonlinear process identification in the presence of multiple correlated hidden scheduling variables with missing data
Pilipovic et al. Parameter estimation in nonlinear multivariate stochastic differential equations based on splitting schemes
Azam Linear discrete-time state space realization of a modified quadruple tank system with state estimation using Kalman filter
Yan et al. What Observables Are Needed for Precision Data-Enabled Learning of Inverse Operators?
Liu et al. Iterative state and parameter estimation algorithms for bilinear state-space systems by using the block matrix inversion and the hierarchical principle
Rezaian et al. Non-intrusive balancing transformation of highly stiff systems with lightly damped impulse response
Huang et al. Performance assessment of process with uncertainties based on minimum variance
Kubalinska Optimal interpolation-based model reduction
Liu et al. A data-driven method for skr identification and application to stability margin estimation
Fujimoto et al. Numerical solutions of Hamilton-Jacobi inequalities by constrained Gaussian process regression
Chu et al. A Deep Learning Approach For Epistemic Uncertainty Quantification Of Turbulent Flow Simulations
Liu et al. Two-stage gradient-based iterative algorithm for bilinear stochastic systems over the moving data window
CN113191082B (zh) 基于机器学习的模型参数获取方法、系统及可读介质
Li et al. Structured covariance estimation for state prediction
Wu Data filtering based multi‐innovation gradient identification for bilinear systems with colored noise
Huang et al. Subspace approach to MIMO feedback control performance assessment

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant