系统辨识方法

n p
新疆大学电气工程学院
3/13/2010
M序列的性质 n-----周期 N = 2 − 1 N −1 N +1 逻辑1状态 2 ，逻辑0状态 2 “游程”-----若干个状态连在一起：总数（周 期）、各不同游程所占比例。 所有M序列都具有“移位相加”的性质。 M序列的相关函数在原点处为最大
n p
Ruz (τ ) =
+∞
+∞
∫
0
g (t ) Ruu (τ − t )dt
∫
0
f (t )δ (t − τ )dt = f (τ )
Ruy (τ ) = Kg (t )
新疆大学电气工程学院
3/13/2010
白噪声的自相关是冲激函数，因此，理论上应该用白 噪声信号作为输入信号，但是实际应用时，理想的白噪 音信号不容易产生，而且由于其随机的特性，一般的执 行器很难按它的变化工作，这些都限制了它在实际中的 应用。实际辨识中，常用伪随机信号代替白噪声，最常 用随机信号是指二进制伪随机序列（PRBS），二进制伪 随机信号常见的有M序列、逆M序列等。
z(k)=-a1z(k-1)-...-a n z(k-n) +b1u (k − 1) + b2u (k − 2) + ... + bnu (k − n) + e(k ) = hT (k )θ + e(k )
h(k ) = [−z(k-1) ...-z(k-n) u (k − 1) ... u (k − n)]T
p
p
新疆大学电气工程学院
3/13/2010
C1
CP
C2
C3
C4
M
4级M序列
新疆大学电气工程学院
3/13/2010
最小二乘批处理方法
差分方程：
z(k)+a1z(k-1)+...+a n z(k-n) = b1u (k − 1) + b2u (k − 2) + ... + bn u (k − n) + e(k )
1.0
0 y(t ) = t −τ − (1 − e T )
0 ≤ t <τ t ≥τ
0.632
τ
τ +T
t
新疆大学电气工程学院
3/13/2010
两点法：
y (t )
1.0
0.63
y (t1 ) = 1 − e t −τ −2 y (t ) = 1 − e T 2
新疆大学电气工程学院
τ
t1
t2
t
3/13/2010
相关辨识法： 经典辨识的主要方法有阶跃响应法、脉冲响应法、频率响应法 等，一般先是通过获取相应的非参数模型然后再转化为参数模 型。这类方法的主要缺点是要求无噪声或噪声很小，一般是离 线进行。还有一种方法是相关辨识法，它具有很好的抗干扰能 力和在线辨识能力， 相关法辨识原理基于Wiener-Hopf方程：
2 w
新疆大学电气工程学院
3/13/2010
PRBS
Pseudo—random—binary series 伪随机二进制序列
最大长度PRBS-------M序列 （1）PRBS可由移位寄存器产生 （2）模2加法规则 （3）一个移位寄存器的输出序列的最大周期 n 级 ------- N = 2 − 1，M序列 n = 4 ------- N p = 15
3/13/2010
Ke G (s) = s (Ts + 1)
新疆大学电气工程学院
−τ s
放大倍数的确定：
y (∞) 自衡对象： K = u y′ ( ∞ ) 非自衡对象： = K u
新疆大学电气工程学院
3/13/2010
一阶惯性滞后环节的辨识
y (t )
Ke G (s) = Ts + 1
−τ s
0
-45 -90 -135 -180 -1 10 10
0
频频 (rad/sec)
10
1
10
2
阻尼比为1.067，肌肉关节的截止频率1.6HZ，自然角频率为3.26。
新疆大学电气工程学院
3/13/2010
10
8
50阶阶阶阶 实实实实
50阶阶阶阶 仿仿实实
6 角角/ 0 / 4 2
30阶阶阶阶 仿仿实实
Ruz (τ ) =
新疆大学电气工程学院
+∞
∫ g (t ) R
0
uu
(τ − t )dt
3/13/2010
Wiener-Hopf方程表明：当系统的输入是原输入信号的自相 关函数时，系统的输出是原输入信号和它对应的输出信号的互相 关函数。 Wiener-Hopf方程是一个积分方程，要想解出脉冲响应的是很困 难的，但是当Ruu等于冲激函数时， Wiener-Hopf化简为：
3/13/2010
这里利用MATLAB提供的赤池检验法函数aic()来判定系统阶 次。其主要MATLAB程序代码： TAIC=zeros(4,4); sys=iddata(y,u,0.005)；%u、y为输入输出序列，0.005采 样周期 for n=1:4 for m=1:4 G=arx(sys,[n,m,0]); TAIC(n,m)=aic(G)； end end
3/13/2010
例1、下列数据是退火温度T与黄铜延展性Y 的试验结果 T Y 300 40 400 500 600 700 50 55 60 67 800 70
1画出散点图 2求出Y对T的线性回归方程
新疆大学电气工程学院
3/13/2010
例2：已知：y=a0+a1*x1+a2*x2 x1 x2 y 0.2 0.1 0.17 0.5 0.3 0.26 0.6 0.4 0.28 0.8 1.0 0.9 1.1 0.23 0.27 1.1 1.4 0.24
t −τ −1 T
0.39
t2 − t1 T = ln(1 − y (t )) − ln(1 − y (t )) 1 2 τ = t2 ln(1 − y (t1 )) − t1 ln(1 − y (t2 )) ln(1 − y (t1 )) − ln(1 − y (t2 ))
30阶阶阶阶 实实实实
0 0
2
4
6
8 时时/s
10
12
14
16
新疆大学电气工程学院
3/13/2010
20
15
角角 角角/ 0
10
5
130的阶阶阶阶 仿仿实实
110阶阶阶阶 仿仿实实
130的阶阶阶阶 实Biblioteka Baidu实实
110阶阶阶阶 实实实实
0 0
2
4
6
8
10 时时/s
12
14
16
18
20
新疆大学电气工程学院
1.0
0.632
τ
τ +T
t
0 y(t ) = t −τ − Ku (1 − e T )
0 ≤ t <τ t ≥τ
新疆大学电气工程学院
3/13/2010
切线法：
0 y(t ) = t −τ − Ku (1 − e T )
归一化：
0 ≤ t <τ t ≥τ
y (t )
系统辨识
张宏立
新疆大学电气工程学院
3/13/2010
自衡对象：阶跃输入下，系统的输出能渐进于一 个新的稳定状态
自衡对象：
−τ s
Ke G (s) = Ts + 1
非自衡对象：
Ke −τ s G ( s) = (T1s + 1)(T2 s + 1)
Ke −τ s G ( s) = s (T1s + 1)(T2 s + 1)
θ = [a1 ...a n b1 ... b n ]T
新疆大学电气工程学院
3/13/2010
z (n + 1) − z (n) z (n + 2) − z (n + 1) = ... ... z (n + N ) − z (n + N − 1)
2
4
6 时时/s
8
10
12
14
新疆大学电气工程学院
3/13/2010
174.9 16.46 × 3.26 2 G(s) = 2 = 2 s + 6.957 s + 10.63 s + 2 × 1.067 × 3.26 s + 3.26 2
Bode Diagram 40 20 幅角 (dB) 相相 (deg) 0 -20 -40
... − z (1) ... − z (2) ... ...
u ( n) u (n + 1) ...
... ... ...
... − z ( N ) u (n + N − 1) ...
Z = Φθ + e T −1 T θ LS = (Φ Φ ) Φ Z
新疆大学电气工程学院
a1 u (1) ... e(n + 1) a e( n + 2) u (2) n + ... b1 ... u ( N ) ... e( n + N ) bn
新疆大学电气工程学院
3/13/2010
实例介绍
新疆大学电气工程学院
3/13/2010
新疆大学电气工程学院
3/13/2010
15 10 5 角角/ 0
10 8 6
角角 / 0
初初输输1.25/正正 初初输输1.25/反正 初初输输1.75/正正 初初输输1.75/反正
0 -5 -10
4 2 0 预预输输1V 预预输输1.25V 预预输输1.5V 预预输输1.75V
-15 -1
-0.8
-0.6
-0.4
-0.2
0 输输/v
0.2
0.4
0.6
0.8
1
-2 0
1
2
3
4
5 时时/s
6
7
8
9
10
控制电压与输出角度关系图
不同预紧压力的阶跃响应
新疆大学电气工程学院
3/13/2010
6 4 2 输输角角/ 0
0 -2 -4
-6 -0.4
-0.3
-0.2
-0.1
0 输输输输/v
新疆大学电气工程学院
3/13/2010
m n
m=1 875.17 143,84 64.53 25.18
m=2 253.89 11.52 28.61 49.22
m=3 174.26 40.88 30.66 19.57
m=4 53.72 70.55 643.32 1045.57
n=1 n=2 n=3 n=4
新疆大学电气工程学院
3/13/2010
白噪声 白噪声是一种均值为零、谱密度为非零常数 的平稳随机过程。 ∞ τ = 0 R (τ ) = σ δ (τ ) = 自相关函数 τ ≠0 0 谱密度 S w ( w) = σ 2 以白噪声为输入，最小二乘辨识是无偏的 有色噪声可利用白噪声通过一个成形滤波器获 得白噪声序列的产生方法
1000
广义二进制伪随机序列
广义二进制伪随机序列自相关特性
P [u (t ) = −u (t − 1) ] = pw P [u (t ) = u (t − 1) ] = 1 − pw
新疆大学电气工程学院
3/13/2010
4 3 2 角角 角角/ 0 1 0 -1 -2 -3 0
0.1
0.2
0.3
0.4
关节系统的李萨育图形
新疆大学电气工程学院
3/13/2010
0.4 0.3 0.2
25 20 15
0.1 幅幅/v 0 -0.1 -0.2
10 5 0
-0.3 -0.4 0 50 100 采采采 150 200 250
-5 0
100
200
300
400
500
600
700
800
900
−1 −d
A( z −1 ) = 1 + a1 z −1 + a2 z −2 + ... + ana z − na
ARX模型：
A( z ) y (k ) = z B( z )u (k ) + ξ (k )
−1 −d −1
A( z −1 ) = 1 + a1 z −1 + a2 z −2 + ... + ana z − na B( z −1 ) = b1 z −1 + b2 z −2 + ... + bnb z − nb
新疆大学电气工程学院
3/13/2010
例4：仿真对象模型：
z (k ) − 1.5 z (k − 1) + 0.7 z (k − 2) = u (k − 1) + 0.5u (k − 2) + ξ (k )
辨识模型：
z (k ) − a1 z (k − 1) + a2 z (k − 2) = b1u (k − 1) + b2u (k − 2) + ξ (k )
A( z −1 ) = 0.00175 + 0.002149 z −1 + 0.001075 z −2
B( z −1 ) = 1 − 1.966 z −1 + .0.9658 z −2
新疆大学电气工程学院
3/13/2010
MATLAB辨识工具箱
AR模型：
A( z ) y (k ) = z u (k ) + ξ (k )
新疆大学电气工程学院
3/13/2010
例3、下表是美国旧车价格调查资料，x表示车的使用 βx 年限，y表示相应价格(美元)，若已知 y = α e 试确定 α , β 的值。
x y 1 2 3 4 1087 5 765 6 7 8 9 10
2651 1943 1494
538 484 290 226 204