范文一:跃响应曲线的系统传递函数辨识方法
z)1(
’Chin纛a。it!nce。and。Tec烹hno焉)gvReview-蛆I
sc
科学论坛
10口■■
1
基于阶跃响应曲线的系统传递函数辨识方法
张东宝t张庚云2
(1.长安大学工程机械学院陕西西安710064:
2.煤炭科学研究总院太原研究院
山西太原030006)
[摘要]介绍了一种利用拉式变换极限定理,由试验阶跃响应曲线辨识系统传递函数的方法,并以二阶系统为仿真实例,在MATLAB环境下编写M文件,运行得到系统模型结构和参数结果。结果表明,该方法误差小,对于低阶系统模型是一种有效的辨识方法.
[关键词]阶跃响应中图分类号:TUl92
拉式变换
参数辨识MATLAB文献标识码:A
文章编号:1009—914X(2009)36—0207—02
1引曹
系统是由相互联系、相互制约、相互作用的各部分组成,是具有一定整体功能和综合行为的统一体。系统的数学模型是系统本质特征的数学抽象,即研究系统内部因素(状态参数)之间以及系统与外部因素(外作用)之间的比较准确的定量关系,以便更深刻地了解系统的内在规律和特性,这些定量关系就是系统的数学模型。
建立系统模型的方法有两种:理论法和实验法。理论方法是从已知的定律、定理和原理出发,通过机理分析找出系统内在的运动规律,推导出系统中
各状态参数和外作用之间的解析关系式——数学模型,也称之为理论建模。
因为这类问题的基本规律是已知的,又称为“白箱”问题。实验方法是直接从系统运行和试验数据(外作用(输入)和系统响应(输出)数据),应用辨识方法,建立系统的数学模型,也称“辨识”建模。这种方法适用于系统的客观规律不清楚的情况,又称为“黑箱”问题。
本文利用Laplace极限定理法由单位阶跃响应曲线辨识常用的二阶定常数线性系统的传递函数,并在MATLAB提供的M文件编辑器中编写源程序,确定传递函数的参数。通过MATLAB仿真可以看出,对于低阶系统模型,该方法误差小,程序简单,是一种有效的辨识方法。
2
Laplace变换极限定理法原理
假设被辨识的线性定常系统传递函数的结构形式为
∞净孑at瓦≥?r-I}丽
当这个系统的输入为单位阶跃函数l(f)时,其输出为单位阶跃响应函数y(O,如图1(a)所示。
(a)1p)—一Go)卜—一y(t)
@lp)—一鱼盟卜一咒(f)
(c)lp)—1百万卜苁(f)
(r)lp)—1互p
i
K(f)
图1应用拉式变换极限定理的系统示意图
根据Lapl8ce变换的终值定理可知:
lit.J《r)=lim毋’(D
--)∞
J--tO
…
根据传递函数的定义,在此情况下系统的传递函数为:
㈣=锱=半
所以,】,(J):!G(s)
(2)
将(2)式代入(1)式中,得到
燃川)2炽G(,)
(3)
此即为Laplace变换的极限定理。为了表示上方便(与下面比较),则有
:。。32时,参数K的计算将比较复杂,还由于积分会造成较大的误差,但对于低阶模型的辨识,仍是一种有效的辨识方法。工程中类似的问题,都可以借鉴本文的方法解决。参考文献
[1]吴广玉.系统辨识与自适应控制[M].哈尔滨:哈尔滨工业大学出版社,
1987.
[2]楼顺天等.TvIATLAB5.X程序设计语言[M].西安:西安电子科技大学出版社,2000.
基于阶跃响应曲线的系统传递函数辨识方法
作者:作者单位:张东宝, 张庚云
张东宝(长安大学工程机械学院,陕西,西安,710064), 张庚云(煤炭科学研究总院太原研究院,山西,太原,030006)
刊名:中国科技博览
英文刊名:CHINA SCIENCE AND TECHNOLOGY REVIEW年,卷(期):2009,(36)被引用次数:
0次
参考文献(2条)
1. 楼顺天 MATLAB5.X程序设计语言 20002. 吴广玉 系统辨识与自适应控制 1987
本文链接:http://d.wanfangdata.com.cn/Periodical_zgbzkjbl200936193.aspx授权使用:山东大学(sddx),授权号:87401759-2df0-4eb6-8a05-9ef700a5c8e0
下载时间:2011年6月3日
范文二:线性系统传递函数的阶数辨识
1997 年 4 月 Journa l of Nan j in g Un iver s ity of Sc ien ce an d Techn o logy A p r. 1997
Ξ α 线性系统传递函数的阶数辨识
ΞΞ 张重雄任继舜
()南京理工大学电子工程与光电技术学院, 南京 210094
摘要 在实际系统的参数辨识中, 系统的结构往往并不知道。 因此, 要想正确辨识
系统参数, 就必须同时要辨识系统的阶数。 该文根据线性系统的频域分析原理, 讨
论了包含有积分环节的线性连续时间系统的阶数辨识方法, 并详细描述了线性系
统阶数辨识的实现过程以及用 检验法来辨识模型结构的步骤与程序设计方法。 F
实验结果表明, 这种方法具有一定的辨识精度和工程应用价值。
关键词 线性系统, 传递函数, 参数识别, 频率响应
分类号 217 T P
在进行系统参数辨识时, 都是假定模型的阶是已知的, 即对于下式: 2 m a 0 + a 1 s + a 2 s + + am s () H s= k 2 n ()s1 + bs + bs++ bs1 2 n
在知道系统的分子、分母阶数 m 、n 以及积分环节数 k 的前提下, 再根据实测的系统频率特性
() () am , b1 , b2 , , , bn 的值。然而, 在实际的系统辨识 的实部 P Ξi 和虚部 Q Ξi 来估计 a 0 , a 1 ,
问题中, 往往并不知道 k、m 、n 的值, 因而要想正确进行系统参数辨识, 就必须同时要辨识系
()() () 、m 、n 值的方法, 并统的阶数。本文讨论了一种根据实测数据 Ξi、P Ξi 、Q Ξi 辨识 H s的 k 给出了实验结果。
1 阶数辨识原理
最小二乘算法是在给定阶的条件下极小化二次型目标函数推导出来的。 参数的最终确 定是使得评价函数 J 为最小 的那一组参数。利用目标函数的极小化确定参数的概念可推广 k
到从一模型集中确定模型的阶数。 R 1 2 () () () ( ) H 如果假设 J k n = [ P Ξi + jQ Ξi -jΞi 作为参数辨识的评价函数, R 为 ||? R i= 1
( ) () () () 采样次数, n 表示 H jΞ的选定阶次 分子分母同为 n 阶, 设 J k n 1 和 J k n 2 分别表示模型 ( ) () ()H jΞ阶为 n 1 和 n 2 的评价函数, 设 n 1 > n 2 > n , 残差假定是正态分布的, 则 J k n 2 和 J k n 1
2 J (n ) 为具有 x 分布的独立随机变量, 自由度分别为 R - 2n , 2 (n - n )。设检测指标为 k 2 2 2 1 -
α 本文于 1996 年 6 月 1 日收到
Ξ 国防科技预研行业基金项目
ΞΞ 张重雄 男 41 岁 副教授
总第 92 期张重雄 任继舜 线性系统传递函数的阶数辨识 107
() () J n - J n R - 2n k 1 k 2 2 ) (r - n , R t = - 2n ]。F 检验又, 则 t 服从 F 分布, t, F 2 n 2 1 2 )(J (n ) 2 n 2 - n 1 k 2
() 称 显著性检验, 也就是说, 阶的改变能否明显地改变评价函数 J k n 。如增加阶, 评价函数将 减小, 但评价函数的减小也可能不是阶的增加引起的, 是由于某种偶然因素造成的。F 检验 就是在一定的置信度的前提下, 确定阶的变化会不会显著地影响评价函数。
用 F 检验法来辨识模型结构的步骤, 就是应用参数辨识方法估计出模型 (a + a s) ?(b 0 1 0
2 2 ) () () + bs, a + a s + a s?b+ bs + bs,, 下 相 应 模 型 参 数, 进 而 计 算 出 相 应 的1 0 1 2 0 1 2 () () J k 1, J k 2, , 最后进行 F 检验。设传递函数模型为 2 m a 0 + a 1 s + a 2 s + + am s () ()H s= 1 2 nb0 + bs + bs++ bs1 2 n 0, b= 0, 依此类1 如果系统中包含 1 个积分环节, 则表示 b0 = 0, 有 2 个积分环节, 表示 b0 = 推。考虑一般情况, 设 2 n ( ) ( ) ( ) a 0 + a 1 jΞ+ a 2 jΞ+ + a n jΞ ( ) ()H jΞ= 2 2 n- 1 nb+ b( jΞ) + b( jΞ) + + b( jΞ) + ( jΞ) 0 1 2 n - 1
考虑到拟合传递函数可能存在近似极点, 并使评价函数的物理意义明显, 从能量的观点定义 评价函数 R ( ) D jΞi L () () ()J = | [ P Ξ+ jQ Ξ-( ) 3 i i H jΞ] ? Ξ| r i i? ( ) D jΞi L - 1i= 1 n n- 1 ( ) jΞ, 下 ( ) ( ) ( ) ( ) b( jΞ) + 式中, D jΞi L ?D jΞi L - 1 为加权系数, D jΞ= b0 + b1 jΞ+ + n- 1
标 L , L - 1 分 别 代 表 迭 代 计 算 的 本 次 与 前 次 的 计 算 值。这 里 的 迭 代 计 算 是 指 在 给 出 ( ) () ( ) ( ) D jΞi L - 1 的任意初始条件下, 即当 L 时 D jΞi L - 1 = D jΞi 0 = 1。3式中, = 1
(Ξ+ Ξ) ?2 1 2 i = 1
) (? Ξi = Ξi+ 1 - Ξi- 1 2 i = , R - 1 2, 3, ? Ξ- ΞR R - 1i = R
J 5J 5J 55J5J 5J 0; 0, 0, 为了使得 J 为最小, 则有 == 0, , = = = 0, , = 0 。 将 5a 0 5a 1 5am 5b0 5b1 5bn- 1 () 3式代入上述各等式可得如下方程值:
5 r ( = - 7()4
() () 其中, 5 为 2n + 1× 2n + 1阶的系数矩阵 ,
V V V T S S 0 - T 3 - 0 0 - 2 4 1 2 4 0 S V - V - S - 4 - 0 2 4 0 2 T T 3 5
fi fi 5 = - U T 1 - S 2 - T 3 S 4 T 5 U 2 0 0 4 - S 0 - 4 S - T S 0 T U U 2 5 6 3 4 6
fi fi T() a 0 , a 1 , am , b0 , b1 , , b, ( = n- 1
() 7 是 5 矩阵中去掉的第 2n + 2 列向量 不包括该列中第 2n + 2 项元素, 5 的元素数值为 R ik i = 0, 2, , 2n V = Ξ- ? ΞΞr i kL k ? k = 1 R i() S = P ΞΞ? ΞΞr r r i = 0, 2, , 2n i k kL k k ? k = 1
南 京 理 工 大 学 学报第 21 卷第 2 期108
R i() T = Q ΞΞ? ΞΞi = 1, 3, , 2n - 1 r r r i k kL k k ? k = 1 R 2 2 i() () U = [ P Ξ+ Q Ξ] r Ξr ? Ξr Ξi = 2, 4, , 2n k i i i kL k ? k = 1 R 1 2() () () ( jΞ) ] |进行一定的迭代次数以后, 可得到 ( , 代入 J k n = | [ P Ξi + jQ Ξi -H i ? R i= 1
()() () () () 中, 就可求得 J k n 。在 2式中分别令 n = 1, 2, , 可得到 J k 1, J k 2, , 然后用 J k 1, () J k 2, , 进行 F 检验, 确定系统的阶数。
2 程序编制方法
F 检验程序的框图如图 1 所示。关于程序编写的几点说明如下:
() 1由于考虑到实际情况中一般的系统阶数不会超过 4 阶, 故在 F 检验中取 1, 6 阶作 为计算范围。
图 1 检验程序框图F
. 1 F igP ro g ram d iag ram o f F te st m e tho d R ( ) D jΞi L 1 () ( ) ( ) | 来估计D jΞi L 与D jΞi L - 1 的接2在 参数辨识中, 取 J p = | 1 -? ( ) D jΞi L - 1 R i= 1 近 程度基本是可行的, 因为输入的系统阶数是真实值。所以当确定一个拟合精度 e 后, 经过若 干次迭代后, 可以出现 J < e,="" 但做阶数辨识时,="" 所选择的阶不一定反映系统的真实模型,="" 因="" p="">
此如果仍用参数辨识的方法, 有时就会发生 J < e="" 始终满足不了的现象,="" 所以应预先设置一="" p="">
() () ( 个较大的迭代次数 本程序中设为 15 次, 在所有迭代中, 取 J n 为最小 对每一个固定的k
)() n 值。取 J n 为最小的一组参数作为参加 检验的数值, 这种方法经实验结果表明是可行 k F
的。 由于篇幅所限, 阶数辨识程序略。
总第 92 期张重雄 任继舜 线性系统传递函数的阶数辨识 109
3 仿真实例
2 100 + 100s + 10s () 对于 H , 取 Α=s= 0101, 频率范围取 0101, 100 H z, 取 80 个点, 2 ()s 100 + 10s + s
() () () () 1321279 1, J k 2= 211634 62, J k 3= J k n 随 n 的变化情况与辨识结果如下: J k 1=
- 6 - 5 () () ()31126 882 × 10, J k 4 = 61128 589 × 10; 分 子 系 数: a 0= 1001044 5, a 1=
- 5 () () () 82, a 3= 11541 112 × 10; 分母系数: b 0= - 21746 508 × 1001050 2, a 2= 91998
- 9 () () () 10, b 1= 1001044 5, b 2= 101002 07, b 3= 1 ; 系统阶数为 3, 误差 J = 31126 882× - 610。 从实例可以看到:
( 1) 当系统中存在积分环节时, 相应的 b、b等值与其它参数值相比将表现为特别小, 01
这样就可以估计出系统的积分环节数。
() 2系统传递函数分子、分母不一定是同阶的, 但 检验是在同阶情况下进行的。 如果 F
系统分子、分母多项式不同阶, 则多余参数项值将特别小, 这样可以去掉这些项。
在上述的判断中, 带有一定的人为主观性, 为了证实判断的正确性, 可以这样来进行检 验: 将最后人为确定的传递函数结构连同作 检验时用的数据一同输入到参数辨识程序 F
() () 中, 辨识参数, 计算残差平方和 J , 并考察 J 与 J n 的 2 个值的情况, 如果 J ν J n 或 J k k k k k k
() () < j="" k="" n="" 或接近于="" j="" k="" n="" ,="" 则认为所确定的传递函数结构是合理的。="">
参考文献
1 李白男 1 伪随机信号及相关辨识 1 北京: 科学出版社, 19871171, 173
2 张成乾, 张国强 1 系统辨识与参数估计 1 北京: 机械工业出版社, 1986. 158, 163
3 张重雄, 任继舜 1 线性系统传递函数的曲线拟合法测量 1 电子测量与仪器学报, 1992 ( 3) : 29
, 32
The O rder D iscernm en t of Tran sf er Fun c t ion f or
L in ear Sy stem
Zh an g C ho n gx io n g R en J ish u n
()Schoo l o f E lec t ro n ic E n g in ee r in g an d O p to e lec t ro n ic T ech no lo gy, N U ST , N an jin g 210094 A BSTRACT W h en a p rac t ica l sy stem is d isce rn ed, it s co n st ru c t io n is o f ten no t k now n.
, , T h e refo reif th e sy stem p a ram e te r s a re to d isce rn th e o rde r o f th e sy stem m u st b e d is2 ce rn ed a t th e sam e t im e. In th e p ap e r, acco rd in g to an a ly sis th eo ry o f lin ea r sy stem in th e
, 2, f requ en cy dom a in th e d isce rnm en t m e tho d o f lin ea r co n t in uo u st im e sy stem is d iscu ssed
, , in c lu d in g in teg ra l lin k an d th e im p lem en t p ro cedu re o f o rde r d isce rnm en t o f lin ea r sy stem th e step s o f d isce rn in g m o de l co n st ru c t io n an d th e w ay o f de sign in g p ro g ram w ith F te st
. m e tho dIt h a s b een m ade c lea r f rom th e exp e r im en t re su lt th a t th e m e tho d h a s a ce r ta in
.d isce rnm en t p rec isio n an d en g in ee r in g app lica t io n va lu e
, , , KEY W O RD S lin ea r sy stem t ran sfe r fu n c t io n p a ram e t r ic reco gn it io n f requ en cy re spo n se
范文三:基于频率响应的传递函数辨识
IDMEC/IST,
Technical University of Lisbon,
Avenida Rovisco Pais, 1049-001Lisboa, Portugal e-mail:dvalerio@dem.ist.utl.pt
Duarte Valério
UNINOVA and DEE of Faculdade de Ciênciase
Tecnologia of UNL, Campus da FCT da UNL,
Quinta da Torre, 2825-114Monte da Caparica,
Portugal
e-mail:mdortigueira@uninova.pt
Manuel Duarte Ortigueira
Identifying a Transfer Function From a Frequency Response
In this paper, the classic Levy identi?cationmethod is reviewed and reformulated using a complex representation. This new formulation addresses the well known bias of the clas-sic method at low frequencies. The formulation is generic, coping with both integer order and fractional order transfer functions. A new algorithm based on a stacked matrix and its pseudoinverse is proposed to accommodate the data over a wide range of frequencies. Several simulation results are presented, together with a real system identi?cation.This system is the Archimedes Wave Swing, a prototype of a device to convert the energy of sea waves into electricity. ?DOI:10.1115/1.2833906?
IDMEC/IST,
Technical University of Lisbon,
Avenida Rovisco Pais, 1049-001Lisboa, Portugal e-mail:sadacosta@dem.ist.utl.pt
JoséSáda Costa
Introduction
The identi?cationof linear systems is an interesting subject that deserved a lot of attention in the past ?1?. Identi?cationin the frequency domain is a particular case with great interest in appli-cations. Algorithms are traditionally based on Levy’swork ?2?, a least-squares based algorithm formulated in a real framework by separating the real and imaginary parts of the transfer function. This led to a formulation with lengthy expressions and to results not equally good at all frequencies ?3?. This frequency depen-dence has been faced ?3?by using an iterative method; another alternative, without iterations, was also proposed ?4?. Both adap-tations modify the basic algorithm introducing weights.
Here, we reformulate the original approach in a completely complex framework, leading to a set of normal equations from which we could remove the frequency dependency. Expressions obtained are shorter. Thus, we obtain two sets of linear equations for each frequency. The formulation is given for the general case of a fractional, commensurate transfer function.
The use of measurements corresponding to several frequencies has been addressed in the past by an average of the coef?cientscomputed from each frequency. This procedure is neither robust nor sensible to order changes. To obtain a more reliable algorithm, we propose a new algorithm with two steps:?rstly,collect the matrices corresponding to the different frequencies in a stacked matrix; secondly, compute the coef?cientsof the model by using the pseudoinverse. To illustrate the behavior of this new formula-tion, we present some simulation results and also a real practical example.
m kq
b 0+b 1s q +b 2s 2q +ˉ+b m s mq ?k =0b k s ?G ?s ?==n kq a 0+a 1s q +a 2s 2q +ˉ+a n s nq ?k =0a k s
?1?
where m and n are the preassigned orders of the numerator and
denominator, and q is the fractional derivative order. Without loosing generality, we set a 0=1.The frequency response of ?1?is given by
??j ??=G
m kq
?k =0b k ?j ??
1+
n kq ?k =1a k ?j ??
=
N ?j ??????+j ????
=
D ?j ??????+j ????
?2?
where N and D are complex valued and ?, ?, ?, and ??the real and imaginary parts thereof ?are real valued. From ?2?we see that
????=
????=
?b
k =0
m
k
Re ??j ??kq ?
?a
k =0
n
k
Re ??j ??kq ?=1+
?a
k =1
n
k
Re ??j ??kq ?
?3?
????=
????=
?b
k =0
m
k
Im ??j ??kq ?
?a
k =0
n
k
Im ??j ???=
kq
?a
k =1
n
k
Im ??j ??kq ?
Adapting Levy’sIdenti?cationMethod for Fractional Orders
Original Formulation. Let us suppose we have a plant de-scribed by a linear system with a transfer function G and a corre-sponding frequency response G ?j ??and that we want to model it using a transfer function
Contributed by the Design Engineering Division of ASME for publication in the J OURNAL OF C OMPUTATIONAL AND N ONLINEAR D YNAMICS . Manuscript received June 6, 2007; ?nalmanuscript received November 30, 2007; published online February 4, 2008. Review conducted by J. A. Tenreiro Machado. Paper presented at the ASME 2007Design Engineering Technical Conferences and Computers and Information in Engineering Conference ?DETC2007?, Las Vegas, NV , September 4–7,2007.
The error between model and plant, for a given frequency ?, will be
??j ??=G ?j ???
N ?j ??D ?j ??
?4?
Minimizing the error power would be an obvious but dif?cultway of adjusting the parameters of ?1?. Instead of this, Levy’smethod minimizes the square of the norm of
def
E ?j ??=??j ??D ?j ??=G ?j ??D ?j ???N ?j ???5?
This leads to a set of normal equations, easy to solve. Omitting the frequency argument ?to simplify the notation, we have
APRIL 2008, Vol. 3/021207-1
Journal of Computational and Nonlinear Dynamics
Copyright ?2008by ASME
E =GD ?N =?Re ?G ?+j Im ?G ????+j ?????+j ??
=?Re ?G ???Im ?G ?????+j ?Re ?G ??+Im ?G ?????
?6?
If we differentiate
?E ?2=?Re ?G ???Im ?G ?????2+?Re ?G ??+Im ?G ?????2
?7?
with respect to one of the coef?cientsb k ?k =0,1,..., m ?or a k ?k =1,2,..., n ?and equal the derivative to zero, we shall have
Alternative Formulation. Since ?E ?2=EE *, and since E *?j ??=G *?j ??D *?j ???N *?j ??, ?8?and ?9?might be given by
?E *?E *
E +E =0
?a k ?a k
?E *?E *
E +E =0
?b k ?b k
From ?2?, it can be seen that
k =1, ... , n
k =0, ... , m ?18?
??E ?2
=0??Re ?G ???Im ?G ?????Re ??j ??kq ??b k
+?Re ?G ??+Im ?G ?????Im ??j ??kq ?=0
?8?
?E
=G ????j ??kq ?a k
?E
=??j ??kq ?b k
?E *
=G *?????j ??kq ?a k
?E *
=???j ??kq ?b k
k =1, ... , n
k =0, ... , m ?19?
??E ?2
=0????Im ?G ??2+?Re ?G ??2?Re ??j ??kq ??a k
+???Im ?G ??+?Re ?G ???Im ??j ???+??Im ?G ?Im ??j ????Re ?G ?Re ??j ????+???Im ?G ?Re ??j ????Re ?G ?Im ??j ????=0
kq
kq
kq
kq
2
2
kq
Let us also de?ne,so as to simplify notation,
def
A ???=?G ????
def
?
?9?
????=arg ?G ????
?
G ???=A ???e j ????G *???=A ???e ?j ????
?
?20?
The m +1equations given by ?8?and the n equations given by ?9?form a linear system that may be solved to ?ndthe coef?cientsof ?1?:
Omitting the frequency argument ?to simplify the notation, we see that ?18?leads to
??????
A B
b a
C D
=e g
cq
?10?
A l , c =?Re ??j ??lq ?Re ??j ??cq ??Im ??j ??lq ?Im ??j ??cq ?
l =0, ... , m ∧c =0, ... , m
lq
lq
?11?
cq
B l , c =Re ??j ???Re ??j ???Re ?G ?j ???+Im ??j ???Re ??j ???
?Im ?G ?j ????Re ??j ??lq ?Im ??j ??cq ?Im ?G ?j ???+Im ??j ??lq ?Im ??j ??cq ?Re ?G ?j ???l =0, ... , m ∧c =1, ... , n
?12?
??
?G *D *?N *?G ?j ??kq +?GD ?N ?G *??j ??kq =0k =1, ... , n
??G *D *?N *??j ??kq ??GD ?N ?G *??j ??kq =0k =0, ... , m
G *G ?D *j kq +D ??j ?kq ??N *Gj kq ?NG *??j ?kq =0k =1, ... , n k =0, ... , m
G *D *j kq +GD ??j ?kq ??N *j kq ?N ??j ?kq =0
?
C l , c =?Re ??j ??lq ?Re ??j ??cq ?Re ?G ?j ???+Im ??j ??lq ?Re ??j ??cq ?
?Im ?G ?j ????Re ??j ??lq ?Im ??j ??cq ?
?Im ?G ?j ????Im ??j ??lq ?Im ??j ??cq ?Re ?G ?j ???l =1, ... , n ∧c =0, ... , m
2
2
lq
?13?
cq
D l , c =??Re ?G ?j ????+?Im ?G ?j ???????Re ??j ???Re ??j ???
+Im ??j ???Im ??j ????
lq
cq
l =1, ... , n ∧c =1, ... , n ?14?
e l ,1=?Re ??j ??lq ?Re ?G ?j ????Im ??j ??lq ?Im ?G ?j ???l =0, ... , m
?15?
g l ,1=?Re ??j ??lq ???Re ?G ?j ????2+?Im ?G ?j ????2?l =1, ... , n
?16?
The results of Levy’smethod are not equally good at all frequen-cies ?4?, but this last step helps address the problem. Since ?j =e ?j ??/2?, and since from Euler’sformula e jx =cosx +j sin x we can easily prove that e jx +e ?jx =2cos x , we will have
??j ?iq j kq +j iq ??j ?kq
=e ?j ??/2?iq e j ??/2?kq +e j ??/2?iq e ?j ??/2?kq =e j ??/2?q ?k ?i ?+e ?j ??/2?q ?k ?i ?=2cos
?
A
2
??a ???j ?
i
i =0
m
n
iq kq
j +j iq ??j ?kq ??iq ?
iq kq j ?
?A
??b ???j ?
i
i =0i ?j ?
j e +j iq ??j ?kq e ?j ???iq ?=0
k =1, ... , n A
??a ?e
i =0
m
n
??j ?iq j kq +e j ?j iq ??j ?kq ??iq ?
iq kq
?
??b ???j ?
i
i =0
j +j iq ??j ?kq ??iq ?=0
k =0, ... , m
?21?
b =?b 0ˉb m ?T a =?a 1ˉa n ?T
?17?
?Superscript T denotes transposition and ∧denotes the Boolean operator of conjunction. ?Further details may be found in the bib-liography ?5,6?.
021207-2/Vol. 3, APRIL 2008
?
?
q ?k ?i ?2
?
?22?
Transactions of the ASME
??j ?iq j kq e j ?+j iq ??j ?kq e ?j ?
=e ?j ??/2?iq e j ??/2?kq e j ?+e j ??/2?iq e ?j ??/2?kq e ?j ?=e j ?q ??/2??k ?i ?+??+e ?j ?q ??/2??k ?i ?+??
?
=2cos q ?k ?i ?+?
2e ?j ???j ?iq j kq +e j ?j iq ??j ?kq
??
??
?23?
=e ?j ?e ?j ??/2?iq e j ??/2?kq +e j ?e j ??/2?iq e ?j ??/2?kq =e j ?q ??/2??i ?k ?+??+e ?j ?q ??/2??i ?k ?+??
System ?25?is dealt with likewise.
????????
A 1B 1
b a
=e 1?b g 1e 2]
a
=A 1B 1
+
C 1D 1A 2B 2]
]
????
e 1
C 1D 1A 2B 2]
]
·
g 1e 2]
C 2D 2A f B f
g 2e f
C 2D 2A f B f
g 2e f
?28?
C f D f g f C f D f g f
?
=2cos q ?i ?k ?+?
2
?24?
Enhancing This Identi?cationMethod With Weights
The alternative formulation was developed so as not to need enhancement by weighting each contribution with a frequency-
Inserting ?22?–?24?into ?21?, we get A
?
i =0
n
????
a i cos q ?k ?i ?k =1, ... , n
?iq
??2
?????
m i =0
m i =0
b i cos ?+q ?k ?i ?
?iq
?2
????
?25?
=0A
dependent weight w ??p ?=w p . Nevertheless, it is possible to in-clude weights as well, if so desired.
When such weights are used, the original formulation in ?26?becomes
?=A
def
?
i =0
n
?
a i cos ?+q ?i ?k ??iq ?
2k =0, ... , m
?????
?
b i cos q ?k ?i ??iq
2
?w A
p p =1f
p p =1
f
p
?=B
?w B
p p =1f
f
p
?=C
?w C
p p =1f
f
p
=0
?=D
Coping With Many Frequencies. In principle, data from only
one frequency suf?ceto ?nda model, but in practice ?due to noise and unavoidable inaccuracies ?, the frequency behavior of the plant must be known in more than one frequency, lest the identi?edmodel should be poor. There are two ways of making use of data from f frequencies.
The ?rstway is summing the systems obtained for each fre-quency. If this approach is followed, matrices A , B , C , and D and vectors e and g in ?10?will be replaced by
?=A
?w D
p
?e =
?w e
p =1
p p
?g =
?w g
p =1
p p
?29?
Similar changes appear if matrices are stacked rather than
summed.
When the alternative formulation is used, system ?25?for a frequency ?p becomes
n
A p
?
i =0
?A
p =1f p =1
f
p
?=B ?e =
?B
p =1f
f
p
?=C ?g =
?C
p =1f
f
p
?=D
?D
p
?e
p =1
p
?g
p =1
p
?26?
A p
where A p , B p , C p , D p , e p , and g p are given by ?11?–?16?for a
particular frequency ?p . Similar summations shall appear in ?25?:
????????????
f
n
A p a i cos q ?k ?i ?
p =1i =0
m
?
i =0
?
b i cos ?p +q ?k ?i ??iq
2p
??
?iq
?2p
????????
?iq
?2p
????
?????
?????
?????
a i cos q ?k ?i ?
m
?iq
?w p 2p
?
b i cos ?p +q ?k ?i ?
i =0
?iq
?w p =02p
k =1, ... , n
n
a i cos ?p +q ?i ?k ?
m
i =0
?iq
?w p 2p
?
b i cos q ?k ?i ?
i =0
?iq
?w p =02p
k =0, ... , m
?30?
The resulting systems are then summed or stacked.
=0
k =1, ... , n
Vinagre’sWeights. One of the reasons why we may want to use weights is to counterbalance a known bias of Levy’smethod that often leads to models well ?ttedto high frequency data but poorly ?ttedto low frequency data. Weights that decrease with frequency can be used to ?xthis. One way to ?ndreasonable values for weights is ?4?
f n
A p a i cos ?p +q ?i ?k ?
p =1i =0
m
?
b i cos q ?k ?i ?
i =0
??
??
=0
k =0, ... , m
?27?
?iq
?2p
w p =
Here, A p =A ??p ?and ?p =???p ?.
The second way is stacking the several systems. Since we will end up with a system that is overde?ned,the pseudoinverse is used to ?tthe best possible solution. System ?10?becomes Journal of Computational and Nonlinear Dynamics
def def
?
?2??1
if p =1
2?21
?p +1??p ?1
if 1?p ?f
2?2p ?f ??f ?1
if p =f
2?2f
?31?
Sanathanan and Koerner’sIterations With Weights. An-other reason to use weights is that, while Levy’smethod relies on
APRIL 2008, Vol. 3/021207-3
Table 1Identi?cationresults for function G ?s …, ?34…, when no weights are used
G ?s ?, without noise, q =0.25,n =3,m =2Identi?edt.f.
?10?, matr. summed ?10?, matr. stacked ?25?, matr. summed ?25?, matr. stacked
?s 0.25?2.75??s 0.25+0.3916??s 0.25?2.75??s 0.25+1??s 0.25+0.3916?
?s 0.5?3s +8??s 0.25+1??s 0.5?3s +8??s 0.25?1.435??s 0.25+0.3937??s 0.25?1.435??s 0.25+1??s 0.25+0.3937?
?s 0.5?3s +8??s 0.25+1??s 0.5?3s +8?
9.7966?10?324.5143?10?291.4855?10?31
J 1.2967?10?29
G ?s ?, with noise, q =0.25,n =2,m =1Identi?edt.f. 0.90334?s 0.25?1.25??s 0.25?1.234??s 0.25+0.9392?0.88796?s 0.25?1.145??s 0.25?1.132??s 0.25+0.9097?0.90445?s 0.25?1.124??s 0.25?1.111??s 0.25+0.9416?0.90912?s 0.25?1.245??s 0.25?1.229??s 0.25+0.9453?
2.1802?10?42.7529?10?42.8829?10?4
J 2.2925?10?4
minimizing E instead of ?, it may be possible ?3?to ?ndthe result of minimizing E performing several iterations, minimizing, in each iteration L ,
E L =
GD ?N D L ?1
?32?essary poles and zeros cancel each other or have very large mag-nitudes ?and do not in?uencethe frequency behavior within the desired frequency range ?. As an example, the results obtained with q =0.25,n =3,m =2are given in Tables 1–3wherein n it is the number of iterations performed when ?33?is used, and J is the mean square error ?MSE ?performance index given by
1??j ???2J =?G ?j ???G f i =1
In the ?rstiteration, we make D 1=1to recover Levy’smethod. From then on, D L ?1will be the denominator found in the previous iteration. This corresponds to the use of a weight given by
w p =
1
?D L ?1??p ??2
?33?
?
f
?35?
If iterations converge, in the end D L ?1?D L , and that is why E L →?. This process is also expected to counterbalance the little in?uenceof low frequency data in the ?nalresult.
In what concerns the iterations proposed by Sanathanan and Ko-erner, they were implemented to stop a ?after 100iterations or b ?if the norm of the difference of the vectors of coef?cientsof two consecutive iterations is smaller than 10?3?log?n +m ?.
Test Data With Noise. Corrupting the data with noise presents a much harder test. Random values within ??1,+1?dB and ??1,+1?deg were added to each value of gain and phase of G ?s ?. Results for q =0.25,n =2,m =1given in Tables 1–3show that, even though this time the structure is much closer to the correct one, values of J are clearly higher, as was expected. The transfer functions obtained are close to G ?s ?, though now with some un-avoidable errors. It can be seen that stacking the matrices and solving the problem with the pseudoinverse presents an improve-ment over the more traditional approach of summation and linear system solving.
Archimedes Wave Swing. An application of the identi?cationmethods above to an engineering problem is provided by the Archimedes Wave Swing ?AWS ?, a device for producing electric-ity from sea waves ?see Fig. 1?.
The sea is an important source of renewable energy ?7–9?. Sev-
Identi?cationResults
These identi?cationmethods have been used with test data from a known function and with data from a wave energy con-verter.
Test Data Without Noise. The exact frequency response of
G ?s ?=
11+s 0.25
?34?
was reckoned at 0.1rad /s, 1rad /s, and 10rad /s, and the identi-?cationmethods above were used to reconstruct the function. When no noise is present, the methods usually reconstruct the original function, even if n and m have values larger than needed. Values of q of which 0.25is an integer multiple ?e.g., 0.125or 0.0625?also allow ?ndingG ?s ?back, if n is large enough. Unnec-
Table 2Identi?cationresults for function G ?s …, ?34…, when Vinagre’sweights are used
G ?s ?, without noise, q =0.25,n =3,m =2
Identi?edt.f.
?10?, matr. summed ?10?, matr. stacked ?25?, matr. summed ?25?, matr. stacked
?s 0.5?1.03s +1.346??s 0.25+1??s 0.5?1.03s +1.346?
?s 0.5?3s +8??s 0.25+1??s 0.5?3s +8??s 0.25+0.6875??s 0.25?0.5064??s 0.25?0.5064??s 0.25+0.6875??s 0.25+1?
?s 0.5?3s +8??s 0.25+1??s 0.5?3s +8?
2.2921?10?303.8929?10?281.2826?10?29
J 2.7188?10?28
G ?s ?, with noise, q =0.25,n =2,m =1
Identi?edt.f. 1.1648?s 0.25?1.439??s 0.25?1.463??s 0.25+1.302?1.1461?s 0.25?1.178??s 0.25+1.288??s 0.25?1.189?1.1364?s 0.25?1.079??s 0.25+1.278??s 0.25?1.087?1.1517?s 0.25?1.134??s 0.25+1.295??s 0.25?1.146?
9.6226?10?48.3611?10?49.2889?10?4
J 0.0011
021207-4/Vol. 3, APRIL 2008Transactions of the ASME
Table 3Identi?cationresults for function G ?s …, ?34…, when iterations with ?33…are used
G ?s ?, without noise, q =0.25,n =3,m =2Identi?edt.f.
?10?, matr. summed ?10?, matr. stacked ?25?, matr. summed ?25?, matr. stacked
?s 0.5?0.4453s +2.944??s 0.25+1??s 0.5?0.4453s +2.944?
?s 0.5?3s +8??s 0.25+1??s 0.5?3s +8??s 0.25?0.944??s 0.25+0.7302??s 0.25?0.944??s 0.25+1??s 0.25+0.7302?
?s 0.5?3s +8??s 0.25+1??s 0.5?3s +8?
1.5813?10?27
2
3.3615?10?30
100
1.0182?10?30
2
J 1.3238?10?31
n it 100
G ?s
?, with noise, q =0.25,n =2,m =1Identi?edt.f. 0.94251?s 0.25?1.423??s 0.25?1.404??s 0.25+1.003?0.88189?s 0.25?1.123??s 0.25?1.11??s 0.25+0.8972?0.91667?s 0.25?1.109??s 0.25?1.095??s 0.25+0.9662?0.91667?s 0.25?1.109??s 0.25?1.095??s 0.25+0.9662?
2.7827?10?4
100
2.7827?10?4
100
3.0943?10?4
100
J 1.8922?10?4
n it 57
eral different sources of energy related to the sea can be distin-guished, such as tidal energy, offshore wind energy, or wave en-ergy. In this last case, the energy present in sea undulation is converted into electricity. Sea waves are, in fact, a concentrated form of wind energy, which in turn is a concentrated form of solar energy, since it is the sun that causes the variations of atmospheri-cal temperature that originate the wind. Harnessing this energy, however, requires solving several engineering problems. The power of sea waves changes with time, and this means three things. First, the energy extracted will have seasonal variations. Second, since short-term variations of wave energy are important, there must be some means of ensuring that electricity produced will be synchronous and phase locked. Third, peak values of wave energy may be so high as to pose a danger to the device. Because of this, and because the sea is by its very nature an aggressive environment, devices for wave energy conversion ?WEC ?must have a very resistant design.
There are many possible con?gurationsfor WECs. Some are installed at the shoreline, some are deployed near the shore, others are deployed offshore; they may ?oator be submerged; and there are several possible working principles. To this day, several prom-ising prototypes have been built. Control engineering plays an important role in ongoing research to pass break-even point, to maximize the energy extracted and hence the rentability of the device.
The AWS is one of the several WECs currently under develop-ment. It is an offshore, underwater ?43m ?device. It is a point absorber, that is to say, its size is negligible compared to a typical wavelength. It consists mainly in a cylindrical, bottom-?xedchamber ?the silo ?, covered by a movable, cylindrical lid ?the ?oater?. Air is trapped between the silo and the ?oater.As the wave crests and troughs pass, the height of the water column above varies, and so does the pressure acting on the ?oater.Thus, and since the air trapped within is compressible, the ?oaterheaves ?see Fig. 1?. A linear generator converts this reciprocating move-ment into electricity.
From this description, it may be supposed that a mass-spring-damper system will provide a simple but reasonable model for the AWS. In reality, things are not so simple, since there are many nonlinearities. A very complete nonlinear simulator has been de-veloped for MATLAB ?10–12?. But, it is convenient, for control design purposes, at least in an initial stage, to have a linear model, even if less accurate. Plotting the simulated responses of the AWS to several sinusoidal waves of the same period, it is seen that the nonlinearities, though not negligible, are not so important as to make such a linear model useless. This linear model may be ob-tained from data provided by the nonlinear simulator. It was cho-sen to do so because, due to operational problems, very few ex-perimental data are available from the 2MW AWS prototype, which was tested at the Portuguese Northern Coast during 2004?and then decommissioned ?. Additionally, due to industrial se-crecy reasons, several parameters of the AWS nonlinear simulator have been modi?edto provide the data below.
To obtain the data, the nonlinear simulator was fed with several sinusoidal waves with an amplitude of 1.0m. Since the AWS is not linear, changing the amplitude of the waves will change ?even if only slightly ?the result obtained. This particular amplitude was chosen since it is the most frequent one in the location where the AWS prototype was deployed 1. The periods of the waves range from 4s to 14s. These are the wave periods expected to occur at
The AWS was tested 5km offshore Leix?es.The data on wave climate were obtained from the ONDATLAS software ?13?. The location for which ONDATLAS has wave climate data is called Leix?esbuoy ?41°12.2?N, 9°5.3?W ?.
1
Fig. 1The AWS before submersion and its working principle
Journal of Computational and Nonlinear Dynamics APRIL 2008, Vol. 3/021207-5
Table 4Data used in the identi?cation
Period ?s ?F exc ampl. ?kN ??amplitude ?m ??gain ?dB ??phase ?deg ?
431.880.1195?108.53?160.20
5108.571.1359?99.61?85.44
6202.641.1905?104.62?26.80
7290.461.1939?107.72?15.77
8365.291.2972?109.20?14.10
9423.351.3573?109.88?27.20
10467.481.5558?109.56?30.72
11501.301.4766?110.62?16.08
12527.621.4132?111.44?9.00
13548.441.3803?111.98?8.31
14565.101.4023?112.11?11.31
the AWS site. The input of the desired model is the wave excita-
tion force F exc , that is to say, the force that the waves would exert in the AWS if the device were stopped. The output is the position of the ?oater?.
Table 4shows the data to which the identi?cationmethods were applied. Several model structures ?that is, several values of q , n , and m ?were used. All fractional values of q led to unstable poles or zeros. Integer models with only one pole were clearly unsatis-factory, and most of the others had unstable poles or zeros. The exceptions were those with q =1,m =0,and n =2or n =3.The former case was chosen because it is simpler and because one further pole does not signi?cantlyimprove the result. This con-?rmsthe supposition above about the AWS being approximately a mass-spring-damper system ?and thus of integer order ?. The im-portance of fractional identi?cationin this case was to rule out the possible existence of fractional dynamics.
Figure 2shows how the model with a lower value of J in Table 5??rstline in the table, obtained with Sanathanan and Koerner’siterations ?is ?ttedto the data from which it was obtained. It can be seen that it is fairly accurate, save for the behavior around 6rad /s, that was neglected. This behavior seems to be a result of nonlinearities, since it could not be modeled even with a greater number of poles and zeros. Another indication that identi?cationresults are satisfactory is that the position of poles does not change signi?cantlywith the method used. It is possible to further validate this model by comparing its outputs with those of the nonlinear simulator; as an example, Fig. 2shows how close the two outputs are for a 100s period taken from an irregular wave similar to those expected to occur in January. See the bibliography ?14?for more details.
Conclusions
In this paper, we proposed a new formulation for the frequency domain identi?cationof linear systems. Two novel features are a complex form of Levy’smethodology and the use of a stacked matrix, leading to a system of linear equations that may be solved using the pseudoinverse. Some examples illustrate the behavior of the algorithms. These show that matrix stacking regularly leads to results better than those obtained summing the contributions of the several frequencies.
Fig. 2Top:Bode diagram of the chosen AWS model; the dots mark the points of Table 4; bottom:output of the chosen AWS model compared with the output of the nonlinear model of the AWS, the input being part of a wave typical for the month of January
Table 5Identi?cationresults for the AWS
No weights
Identi?edt.f.
?10?, matr. summed ?10?, matr. stacked ?10?, matr. summed ?10?, matr. stacked
2.259?10?6
0.6324s 2+0.1733s +12.232?10?6
0.6416s 2+0.1647s +12.369?10?6
0.6282s 2+0.2231s +12.249?10?6
0.6349s 2+0.1699s +1
J 4.6628?10?136.0520?10?139.8387?10?134.8735?10?13
Vinagre’sweights Identi?edt.f. 2.269?10?6
0.6316s 2+0.1766s +12.184?10?6
0.6662s 2+0.1667s +12.415?10?6
0.6185s 2+0.232s +12.253?10?6
0.6362s 2+0.1724s +1
J 4.6812?10?131.4088?10?121.0893?10?125.0636?10?13
Iterations with ?33?
Identi?edt.f. 2.306?10?6
0.6254s 2+0.1715s +12.169?10?6
0.6891s 2+0.1709s +12.394?10?6
0.6247s 2+0.1848s +12.394?10?6
0.6247s 2+0.1848s +1
J 4.2202?10?132.5257?10?124.4763?10?134.4763?10?13
n it 3744
021207-6/Vol. 3, APRIL 2008Transactions of the ASME
Acknowledgment
The authors acknowledge the contribution of Pedro Beir?o?In-stituto Superior de Engenharia de Coimbra, Portugal ?for obtain-ing the data concerning the AWS. The picture in Fig. 1is also his. Duarte Valériowas partially supported by Grant No. SFRH/BPD/20636/2004of FCT, funded by POCI 2010, POS C, FSE, and MCTES. Research for this paper was partially supported by Grant No. PTDC/EME-CRO/70341/2006of FCT, funded by POCI 2010, POS C, FSE, and MCTES.
References
?1?Ljung, L., 1999, System Identi?cation:Theory for the User , Prentice-Hall,
Upper Saddle River.
?2?Levy, E., 1959, “ComplexCurve Fitting,”IRE Transactions on Automatic
Control, 4, pp. 37–44.
?3?Sanathanan, C. K., and Koerner, J., 1963, “TransferFunction Synthesis as a
Ratio of Two Complex Polynomials,”IEEE Transactions on Automatic Con-trol, 8, pp. 56–58.
?4?Vinagre, B., 2001, “Modeladoy Control de Sistemas DinámicosCaracteriza-dos por Ecuaciones íntegro-Diferencialesde Orden Fraccional,”Ph.D. thesis, Universidad Nacional de Educacióna Distancia, Madrid.
?5?Valério,D., and Sáda Costa, J., 2005, “Levy’sIdenti?cationMethod Extended
to Commensurate Fractional Order Transfer Functions,”Fifth EUROMECH
Nonlinear Dynamics Conference , EUROMECH, pp. 1357–1366.
?6?Valério,D., and Sáda Costa, J., 2006, “Identi?cationof Fractional Models
From Frequency Data,”in Advances in Fractional Calculus:Theoretical De-velopments and Applications in Physics and Engineering , J. Tenreiro Machado, J. Sabatier, and O. Agrawal, eds., Springer-Verlag, New York. ?7?Falc?o,A., 2006, “TheHistory of and Progress in Wave Energy Conversion
Devices,”Ninth World Renewable Energy Conference .
?8?Boud, R., 2002, “Statusand Research and Development Priorities 2003:Wave
and Marine Current Energy,”UK DTI Technical Report No. FES-R-132; UK AEAT Report No. AEAT/ENV/1054.
?9?Pontes, M., and Falc?o,A., 2001, “OceanEnergies:Resources and Utilisa-tion,”18th World Energy Conference .
?10?Pinto, P., 2004, “TimeDomain Simulation of the AWS,”MS thesis, Technical
University of Lisbon, Portugal.
?11?Sáda Costa, J., Pinto, P., Sarmento, A., and Gardner, F., 2003, “Modellingand
Simulation of AWS:A Wave Energy Extractor,”Proceedings of the Fourth IMACS Symposium on Mathematical Modelling , Agersin-Verlag, pp. 161–170.?12?Sáda Costa, J., Sarmento, A., Gardner, F., Beir?o,P., and Brito-Melo, A.,
2005, “TimeDomain Model of the Archimedes Wave Swing Wave Energy Converter,”Proceedings of the Sixth European Wave and Tidal Energy Con-ference , pp. 91–97.
?13?Pontes, M. T., Aguiar, R., and Oliveira Pires, H., 2005, “ANearshore Wave
Energy Atlas for Portugal,”ASME J. Offshore Mech. Arct. Eng., 127, pp. 249–255.
?14?Beir?o,P., Valério,D., and Sáda Costa, J., 2007, “LinearModel Identi?cation
of the Archimedes Wave Swing,”IEEE International Conference on Power Engineering, Energy and Electrical Drives .
Journal of Computational and Nonlinear Dynamics APRIL 2008, Vol. 3/021207-7
范文四:传递函数频域辨识及M序列生成指南
传递函数频域辨识及M序列生成指南
12自动化
摘要:系统辨识、状态估计和控制理论是现代控制论中相互渗透的三个领域。控制理论的应
用离不开系统辨识技术,实际中,许多控制系统的模型在工作中是变化的,为了实现自适应
控制,需要系统辨识技术不断更新模型参数。通过学习使用MATLAB软件,初步体验系统
辨识方法。
关键字:系统辨识,控制理论,MATLAB。
Abstract:System identification, State estimation and The Principle of Automatic Control are three different disciplines of the modern control theory, which are interpenetrated with one another. In practice, the model of system is changing all the time. To control adaptively, the system model should be update its parameters, by the method of System identification. By learning the using of MATLAB, we are supposed to practice the method of system identification.
Key Words: System identification, System identification, MATLAB
目录
一、 引言 ......................................................................................... 3
1.1介绍 ......................................................................................... 3
1.2实验目的 ................................................................................. 3 二、 实验内容和方法 ...................................................................... 3
2.1实验内容 ................................................................................. 3
2.1.1实验一: ............................................................................ 3
2.1.2实验二: ............................................................................ 3
2.2实验步骤 ................................................................................. 4
2.2.1实验一: ......................................................................... 4
非周期激励信号 ...................................................................... 6
2.2.2实验二:M序列 ............................................................. 7
.......................................................................................................... 12 三、 实验分析分析 ........................................................................ 12
一、 引言
1.1介绍
在自然科学和社会科学的许多领域中,人们越来越重视对系统进行定量的系统分析、系统综合、仿真、控制和预测。将被控对象模型化,是开展工作的前提和基础。
基于传递函数的经典控制理论分析前需要确定模型的各个参数,忽略工作中的变化,适用于精度要求不高的情况下。而实际中许多控制系统内部参数是未知的,而且工作中有可能会变化,这就需要系统辨识来确定模型参数。通常为机理分析法和测试法相结合来辨识,对于系统机理已知的部分采用机理分析法,机理未知的部分采用测试法。 1.2实验目的
本次实验通过使用MATLAB 仿真系统辨识的过程,采用频域特性拟合的Levy方法,按要求完成传递函数的辨识。描述实验验证的数据准备、基本过程和实验结果。
结合课堂学习的M序列的特性和生成方法,使用matlab语言产生一个M序列并分析器自相关特性。
二、 实验内容和方法
2.1实验内容
2.1.1实验一:
1( 自己设定一个稳定系统,采用周期测试信号,测定系统的频率响应。
2( 对题1中的系统,采用非周期测试信号确定系统的频率响应,并与题1的结果对比。 3( 基于题1或题2 产生的频率响应数据,采用课堂讲授的频域特性拟合方法,辨识传递函数的参数。将
辨识结果与Matlab工具库中的等价的功能函数invfreqs产生的结果做对比。 2.1.2实验二:
根据最大长度现行反馈寄存器M序列生成机制,编写M序列生成的生成程序。
1( 自己设定移位寄存器的级数和初值,产生响应的M序列。
2( 绘制题1产生的M序列的自相关函数和功率谱密度图形。
2.2实验步骤
2.2.1实验一:
周期性信号:
周期性信号最典型的就是正弦周期函数,在不同频率下会产生不同的周期,进而有不同的
频率响应。将不同频率的正弦波组合起来作为输入信号。
U=Asin(,??+,)
Y=A′sin(,??+,) 所以: Am=A’/A
Phi= ,?,
设计系统函数:
num=[4 3];
den=[4 5 6];
sys=tf(num,den);
Transfer function:
4 s + 3
---------------
4 s^2 + 5 s + 6
根据这个原理,设计离散的正弦周期信号。
1、一个周期内取1000个采样点,采样是个周期。单位采样间隔二0.001*T
2、由于T=2*pi./w,所以每一次频率改变都黑改变周期,进而有不同的激励和响应图像。
3、在十个周期内,对第五和第六个采样周期进行分析,得到周期内的最大值和最大值的
坐标。
4、对比激励和响应的最大值的峰值比和相位差,进而得到不同频率的一系列频谱响应。
Matlab代码如下:
w=logspace(-1,2,1000); %w为0.1到100期间的1000的频率点 T=2*pi./w;
dt = 0.001*T;
tmax =round(10*T);
for i=1:1:1000
t = (0:dt(i):tmax(i))';
u1 = sin(w(i)*t);
u=lsim(sys,u1,t);%产生在激励下的相应
[Am(i) n]=max(u(5*1000:6*1000));%n为在5000~6000个点中峰值点的坐标。 Phi(i)=(n-T(i)/4/dt(i))*0.001*2*pi; %得到相频特性。
end
P=Am.*cos(-Phi);
Q=Am.*sin(-Phi);
H=P+1j*Q;
[n1 d1]=freq2levy(H,w,2,1); [n2 d2]=freq2levy(h,w,2,1); bode(tf(n2,d2),'r');hold on;bode(tf(n1,d1),'g');
结果如下:
Bode Diagram0
-10
-20
-30
-40Magnitude (dB)
-50
-6045
0
-45
Phase (deg)-90
-135-101231010101010
Frequency (rad/s) 绿色为辨识出来的系统的波特图,红色为用系统函数用matlab自带freqs函数画出来的图像。
可以看出这里二者的图像已经较为接近,如果继续增大w的密集程度,图像将更为接近。但计算结果较慢。
非周期激励信号
周期激励信号有很多种选择,如单位阶跃信号和脉冲信号,以及高斯信号和指数函数信号。在选择中我纠结了很长时间,开始选择用阶跃信号来做,结果不是很理想。后来想起了指数函数信号,决定拿来试一试。
非周期激励信号顾名思义就是不能用周期信号的周期性方法来解决这个问题。像正弦信号这种周期和频率密切相关的信号可以直接由激励响应得到关于不同频率的输出,而非周期信号如指数函数则不能这么做。要进行傅里叶变换将其变换到频域上去。
Matlab在矩阵运算方面有着得天独厚的优势,所以这里的傅里叶变换我采用了矩阵的计算方法。,代码非常简单容易理解。
y=y2'*exp(-1j*q'*w*dt);
u=u2*exp(-1j*q'*w*dt);
num=[32];
den=[234];
sys=tf(num,den);
dt=0.001;
t=0:0.001:30;
u2=exp(-t);
y2 = lsim(sys,u2,t);
w=[0:0.01:2]*2*pi;
n=length(w);
q=1:30001;
y=y2'*exp(-1j*q'*w*dt);
u=u2*exp(-1j*q'*w*dt);
l=y./u;
[n d]=freq2levy(l,w,2,1);
bode(tf(n,d),'r')
hold on
bode(tf(n2,d2),'g')
Bode Diagram10
0
-10
-20
Magnitude (dB)
-30
-4045
0
-45Phase (deg)
-90-2-10121010101010
Frequency (rad/s)
可以看出,辨识效果非常好,非常接近于系统自带函数freqs()的频响图。
2.2.2实验二:M序列
m序列码发生器是一种反馈移位型结构的电路,它由n位移位寄存器加异或反馈网络组成,其序列长度21nM ,,只有一个多余状态即全0状态,所以称为最大线性序列码发生器。由于其结构已定型,且反馈函数和连接形式都有一定的规律,因此利用查表的方式就设计出m序列码。列出部分m序列码的反馈函数F和移存器位数n的对应关系。如果给定一个序列信号长度M,则根据12, nM求出n,由n查表2-1便可以得到相应的反馈函数F。
根据m 序列的特征方程: 移位寄存器结构为
自相关特性:
1、定义:,: {0,1} , {1,-1}, 即,(0)=1, ,(1)=-1
2、M序列{a}与经τ迟延移位序列{a}分别表示: nn-τ
{a}={a,a,…,a}, {a}={a,a,…,a,a,…,a} n12pn-τ1 +τ2 +τp1τ
3、序列{,(a)}的自相关定义为 n
R(τ)={,(a),(a)+,(a),(a)+…+ ,(a),(a)+…+ ,(a),(a)}/N 11+τ22+τp-τppτp
={,(a,a)+,(a,a)+…+,(a,a)+…+,(a,a)}/N 11+τ22+τp-τppτp=(0的数目-1的数目)/N p
注意到,(a,a)=,(a),(a) nn+τnn+τ
4、根据移位相加特性,{a},{a}仍是M序列中的元素。M序列中一个周期中0 的数nn+τ
目与1 的数目之差。又由M序列的均衡性可知,0 比1 的个数少一个
R(0)= 1;R(τ)=-1/N p
生成M序列的代码如下:
function[new_sequence]=M_sequence(init_register,type_of_connec
tion)
N=length(init_register);%判断初始值的长度 length1=2^N-1;%可产生最大序列的长度 register=init_register;
temp_register=ones(1,N);
new_sequence=ones(1,length1); new_sequence(1)=init_register(N);%输出第一个M序列值
for i=2:1:length1
temp_register(1)=mod(sum(register.*type_of_connection),2);
for t=2:1:N
temp_register(t)=register(t-1);%移位
end
register=temp_register;
new_sequence(i)=register(N); end
end
如:
输入:M=M_sequence([0 0 0 0 1],[0 1 0 0 1]);
结果如下:
下面为自相关函数的代码:
function f=Autocorrelation(M) Np=length(M);
sum=0;
M1=ones(1,Np);
for k=1:1:150
sum=0;
k1=mod(k,Np);
for i=1:1:Np-k1
M1(i)=M(i+k1);
end
for j=1:1:k1
M1(Np-k1+j)=M(j);
end
for t=1:1:Np
sum=sum+juan(M(t),M1(t));
end
R(k)=sum/Np; end
h=1:1:150;
f=R;figure;
plot(h,R);
figure;
Hpsd = fft(R,Np); Pxx=abs(Hpsd);
plot_Pxx=10*log10(Pxx);
step(0:1:Np-1,Pxx); title('Gonglv') end
function f=juan(a,b) if a==b
f=1;
else
f=-1;
end
end
结果如下:
自相关函数:
1.21
0.80.60.40.20
-0.2050100150
功率谱:
功率谱密度
10
9
8
7
6
5
4
3
2
1
0051015202530 三、 实验分析分析
1、实验中发现频域法的频率取点非常重要,如果取的点不适合结果就差别很大。另外,我编写的辨识程序有一定的局限性,比如非周期的激励信号,在实验报告中的系统的辨识系统辨识效果非常好。如果我把系统函数中的各个参数都放大数倍,或者增加阶数,辨识效果可能就没有这么好了,在这个方面可能还需要改进,要达到系统自带函数的辨识度,还有一点差距。
2、有实验验证了M序列的确有优秀的自相关性质,加深了对理论知识的理解。
范文五:基于阶跃响应曲线的系统传递函数辨识方法
基于阶跃响应曲线的系统传递函数辨识方
法
基于阶跃响应曲线的系统传递函数辨识方法
张东宝1张庚云
(1.长安大学=】=_=程机械学院陕西西安710064:2.煤炭科学研究总院太原研究院山西太原030006)
科学论坛
啊I
[摘要]介绍了一种利用拉式变换极限定理,由试验阶跃响应曲线辨识系统传递函数的方法,并以二阶系统为仿真实例,在MATLAB环境下编写M文件,运行 得到系统模型结构和参数结果.结果表明,该方法误差小,对于低阶系统模型是一种有效的辨识方法.
[关键词]阶跃响应拉式变换参数辨识MATLAB
中图分类号:TU192文献标识码:A文章编号:1009—9】4x(2009)36—0207—02 1目l言
系统是由相互联系,相互制约,相互作用的各部分组成,是具有,定整 体功能和综合行为的统一体.系统的数学模型是系统本质特征的数学抽象, 即研究系统内部因素(状态参数)之间以及系统与外部因素(外作用)之间的比 较准确的定量关系,以便更深刻地了解系统的内在规律和特性,这些定量关系 就是系统的数学模型.
建立系统模型的方法有两种:理论法和实验法.理论方法是从已知的定
,推导出系统中 律,定理和原理出发.通过机理分析找出系统内在的运动规律
各状态参数和外作用之间的解析关系式——数学模型,也称之为理论建模. 因为这类问题的基本规律是已知的,又称为"自箱"问题.实验方法是直 接从系统运行和试验数据(外作用(输入)和系统响应(输出)数据),应用辨识方 法,建立系统的数学模型,也称"辨识"建模.这种方法适用于系统的客观 规律不清楚的情况,又称为"黑箱"问题.
本文利用Laplace极限定理法由单位阶跃响应曲线辨识常用的二阶定常
数线性系统的传递函数,并在MATLAB提供的M文件编辑器中编写源程序,确定
传递函数的参数.通过MATLAB仿真可以看出,对于低阶系统模型,该方法误
差小,程序简单,是一种有效的辨识方法
2Laplace变换极限定理法原理
假设被辨识的线性定常系统传递函数的结构形式为: G(j丽
当这个系统的输入为单位阶跃函数l(,)时,其输出为单位阶跃响应函数
(,),如图1(a)所示.
)l【f)—G)卜-_一r)
I(0鱼!卜—一(f)
(c)l)—_+[互(,)
(r)1(0—百万卜jIr(,)
网1应用拉式变换极限定理的系统示意图
根据Laplace变换的终值定理可知:
lim
~-i.at*
)'((1)
,—?O…
根据传递函数的定义,在此情况下系统的传递函数为: ?==
所以,】,):G((2)
J
将(2)式代入(i)式中,得到
州)=l
?
iraG((3)
此即为Lap]ace变换的极限定理.
为了表示上方便(与下面比较),则有
G(.):K(4)
现在,假设存在(设置)一个如图1(b)所示的系统,它的输入为单位阶跃函
数,传递函数为q),它的输出为(,),定义为: (,)【一jf)(5)
这时系统q(s)的放大倍数为: =
,
lira31(t)【一Jr)la(6) 根据前述的Lap1ace燹换的微限定理有: .,p)Gl(
由Laplace变换的积分性质和(5)式,可得 上)】=LLJ0tr)}
:三『一]
?::
7L-T--7-~
(8)
=
?【岛一G(】(9
由式(7)和式(9)可得:
q)=li—
m
L
,
[/~一G】
Kq(10)
进,步设置一个系统,如图1(c)所示,并定义和存在:
o);J:【一(r)(H)
=
f【墨一M(f)(12)
使用Laplace变换的极限定理得:
=q一(13)
下面研究一般情况,考虑系统q(),如图1(d),定义其输出为: o)=【一.(f)(14)
=
r【墨一.一Y,-j(r)r(15) 同样,采用前述的Laplace变换的极限定理,并用数学归纳法可得:
=
Kq一2+…+(一1)q
这样,可得到一线性方程组:
K=K
K=al
一q—
【:墨-r,一+…+(一1),-tq
由此方程组,可递推算出K,q,,…….
3仿真实例
MATLAB是一个功能强大的编程工具软件,同时它又是一种高效的编程语
言,文件是用~TLAB语言编写的程序代码文件.用Laplace极限定理法由
单位阶跃响应曲线辨识常用的二阶定常数线性系统的传递函数,假设给定的二
阶系统为:
…
7一
?)去1'+q十l3'+4+
(其中,q=4,=3,K=7.)
其单位阶跃响应函数为:
.
y(,)=3.5×(2+一一3ep)
科技博览l207
科学论坛
I?
在MATLAB提供的M文件编辑器中编写代码如下 clearall,clc
t=l:18:
y=7+3.5*exp(一t)一10.5*exp(一1/3*t): KO=7:KI=O:
gl=(KO+KO—Y(1))/2:
fori=2:18
KI=KI+(KO,Y(i)+KO—Y(卜1))/2: end
yl(1):(KO+KO—y(1))/2;
forj=2:18
yl(J)=yl(j-1)+(KO—Y(J)+KO—Y(j-1))/2: end
K2=(KI+K1一yl(1))/2:
fork=2:18
K2=K2+(K1一yl(k)+K1一yl(k-1))/2: end
K1
K2
al=K1/KO
a2=(Kl*al—K2)/KO.
Y=tf([5],[23l】):
YY:tf([KO],[a2,al,1]): step(Y)holdon:
step(YY):
grid
legend('给定系统阶跃响应曲线','辨识系统阶跃响应曲线') a=step(Y):
b=step(YY):
f0ri=1:100
e(i)=a(i):
d(i)=b(i):
end
R=corrcoef(c,d)
运行M文件,自动生成系统模型结构和参数结果:
啪【l?c】
图2给定系统与辨识系统的阶跃响应
数值积分参数:
=27.9254,=90.5004
参数估计:.
K=7,=3.9893,=2.9862
辨识结果:
7
G(s)丽丽
相关系数:R=0.9999
从仿真结果(图2)可以看出,用极限定理法由单位阶跃响应辨识二阶定常 线性系统的传递函数时,误差不大,最大相对误差(在t=2s)只有(1,2)%.造 成辨识误差的主要原因是计算(f:l,2,3?'?数值积分引起的误差.主要有 208l科技博览
两个因素:1)数值积分的方法:2)采样间隔的大小.本文采用梯形法求积分, 采样间隔时间At=1秒,采样总长度N=18.
结语
根据系统阶跃响应曲线及试验数据,利用Laplace变换的极限定理,借助 MATLAB强大的M文件编程功能,能够自动实现模型的传递函数和相关参数的求
取,以及所需曲线的自动生成.该方法克服了传统试探法的不足,但当i>2 时,参数的计算将比较复杂,还由于积分会造成较大的误差,但对于低阶模 型的辨识,仍是一种有效的辨识方法.工程中类似的问题,都可以借鉴本文的 方法解决.
参考文献
[1]吴广玉.系统辨识与自适应控制[M].哈尔滨:哈尔滨工业大学出版社, 1987.
[2]楼顺天等.MATLAB5.x程序设计语言[M].西安:西安电子科技大学出 版社,2000.
转载请注明出处范文大全网 » 跃响应曲线的系统传递函数辨识