f
=
−
1 2
ρv 2 ACd
(2.1)
其中 v 为乒乓球运动速度空气密度 ρ =1.19kg/m3,乒乓球截面积 A = 0.001256m2 ，
空气阻力系数 Cd = 0.1。
2
(1) 乒乓球向上运动时， ma = F − mg − f 。
(2) 乒乓球向下运动时， ma = F − mg + f 。
PLC
电压u与高度 H之间的非线
性关系
电压u
1
Ts +1
设定高度h
图 2-9 整体模型方框图
第三章 总结
本文 Airball Demo 建模涉及机理建模与使用自相关辨识的暗箱建模。M 序 列自相关辨识法使用简单，实践也证明其有很好的复现性，仿真的结果接近实际 模型。但是，选择合适的 M 序列又因实际模型不同而相差交大，有时甚至需要 多次测量系统的固有特性参数。机理建模理论上能跟实际系统模型完全吻合，但 由于存在随机的，不确定的因素，往往有一定的偏差。无论如何，系统模型辨识 永远无法反映真实的系统，所以单单依据辨识得到的模型进行控制仿真得到的控 制方法需要再实际系统中运行、多次改进之后才有实用性。现在出现了很多更为 先进的控制方法，例如把系统辨识融合到控制算法中，不断改变控制参数以适应 系统模型的改变，对与模型的准确性要求不高。
因为乒乓球运动的速度 v 较小，最高不会超过 0.1m/s。则空气阻力
f max
=
1 2
× 0.12
× 0.001256× 0.1 =
7.4 ×10−7
N
(2.2)
由式（2.2）看出，空气阻力与 F 或者 G 相比相差太大，可以忽略。则乒乓球向
上和向下运动的过程可以统一为 ma = F − mg ，微分与积分方程分别如式(2.3)与
N −1
∑ RMz (k) ≈ g( j)RMM (k − j)Δt j=0
M 序列的响应函数可写成
⎧
∑ ⎪
⎪
RMz (k )
⎨
∑ ⎪⎪⎩RMM (k)
= =
1 Np 1 Np
N p −1
M (i
i=0 N p −1
M (i
i=0
− k)z(i) − k)M (i)
从式（2.11）和（2.12），可推导出 M 序列的自相关函数为
1
电压u
Ts +1
电压u与高度H 之间的非线性
对应
运算处理
图 2-8 红外高度传感器子系统方框图
2.4 系统整体模型
到目前为止，三个子系统的模型都已建立。系统的整体模型如图 2-9 所示。
7
电机输入电压U U与F的响应函数 风扇作用力F+
1
g（k）
-
m
mg
乒乓球自身重力G
1
1 乒乓球高度H
s
s
运算处理
2.1 乒乓球子系统机理建模
在整个模型中，机理最为明确的的是乒乓球的受力过程。乒乓球受到三个力 作用分别是自身重力 G，风力 F 和空气阻力 f。乒乓球上升和下降过程受力分别 如图 2-1、2-2 所示。
F
Ff
v
v
Gf 图2-1 乒乓球上升过程受力图
G 图2-2 乒乓球下降过程受力图
根据牛顿第二定理， F合 = ma 。空气阻力计算公式[1]，
RMM
⎧ (k )⎪⎩⎪⎨−
a2 , a2 , Np
k = 0, N p ,2N p，L , k ≠ 0, N p ,2N p，L
系统输入输出的互相关函数为
RMz
(k)
=
(N
p
+ 1)a 2Δt Np
g(k)
−
c
∑ 式中
c
=
a 2 Δt
N p −1
g (i)
。对于稳定的系统，
c
是有界常数，而且很小。
(2.4)所示。
F
−
mg
=
m
d 2v dt 2
(2.3)
v
=
1 m
∫∫
(F
−
mg
)dt
(2.4)
乒乓球子系统方框图如图 2-3 所示。
风扇作用力F+ -
1
1
1 乒乓球高度H
m
s
s
mg
乒乓球自身重力G
图 2-3 乒乓球子系统方框图
2.2 M 序列自相关法辨识
2.2.1M 序列产生 M 序列属于伪随即信号，其自相关函数和功率谱在一定条件下接近于白噪声
第一章 模型简介
Airball Demo 实验装置由贝加莱公司提供，该系统由控制器、风扇、玻璃管、 乒乓球组成，如图 1-1 所示。各部件描述如下。 1. 底座：系统的支撑部分，主要是固定管道部分和乒乓球。此外，
还在底部安装了吹动乒乓球的直流风扇，与检测乒乓球 高度的红外传感器。 2. 管道：与底座相连，是乒乓球竖直运行的通道。管道外直径为 50mm，壁厚为 4mm， 可保证乒乓球在管道内自由（无 摩擦）运行。管道有效长度：400mm。 3. 乒乓球：标准乒乓球，直径 40mm，重量 2.7g。 接通风扇电源，风扇转动形成的气流将对乒乓球作用一个向 上的力，通过控制风扇的转速，可以调节对乒乓球的作用力，使 之在管道内上下运行并稳定在某一设定高度。 总的来说，整个系统的输入为电机的电压 U，输出为乒乓球 的高度 H。其控制方框图如图 1-2 所示。
的确定性、周期性的信号。M 序列是由 N 级串联的移位寄存器，按照一定的反 馈联结而产生的周期性方波。四级移位寄存器产生 M 序列结构图如图 2-4[2]所示。
图 2-4 四级移位寄存器产生 M 序列结构图
3
图中，⊕ 表示异或运算，C1 、C2 、C3 、C4 表示四个移位寄存器，x1 、x2 、 x3 、 x4 分别表示各移位寄存器的输出， x4 输出为 M 序列。四级移位寄
H
0
2 34 5
图 2-7 传感器阶跃响应曲线
可以使用一阶惯性环节近似，如下
C(s) = 1 Ts + 1
(2.16)
其中，T 为传感器的时间常数，T 越小，传感器器响应越快。可以通过测得
传感器阶跃响应曲线，类似图 2-7，得到时间常数T 。红外高度传感器子系统方
框图如图 2-8 所示。
PLC
高度H
图 2-6 高度与输出电压曲线 传感器存在惯性，它的输出不能立即复现输入信号，而是从零开始，按指数 规律上升，最终达到稳态值。使用阶跃相应法测传感器惯性的传递函数。使 t0 时
6
刻传感器输出 0 高度（经过非线性转换）， t1 时刻传感器输入 h 的高度并同时记 录传感器的经过转换的输出高度 H 的时间响应曲线，其一般形式如图 2-7 所示。
存器的连接方式可用式(2.5)表示。
⎧ x2 (k +1) = x1(k)
⎪⎪ ⎨ ⎪
x3 (k +1) = x2 (k) x3(k +1) = x3(k)
⎪⎩x1(k +1) = x3 (k) ⊕ x4 (k)
(2.5)
设初始状态为 1010，在移位脉冲 CP 作用下，寄存器各级状态的变化如表 2.1[2]
2.2.2 相关辨识原理
相关辨识的理论较为复杂，涉及使用变分法求泛函最小值函数问题，此处作 简要说明。相关辨识的基础是 Wiener – Hopf 方程，如式(2.6)所示。
∫ Rxy (τ ) =
∞ 0
g (t ) R xx
(t
−
τ
)dt
ቤተ መጻሕፍቲ ባይዱ
(2.6)
使用 M 序列相关辨识过程中，式(2.6)中的 g(t) 为系统输出 y(t) 与输入 x(t) 的比
可很大程度上简化计算过程,此时 Wiener – Hopf 变换为
∫ RMz (τ ) ≈
N 0
p Δt
g
(t
)RMM
(t
−
τ
)dt
(2.9)
式中，N p 是 M 序列的循环长度；Δt 是 M 序列移位脉冲的周期。当 N p 充分大时，
RMz (τ ) ≈ a2 g(τ )
当数据采集时间和 Δt 相等时，式（2.9）的离散形式为
N p i=0
到目前为止，可以得到 M 序列自相关辨识法求响应函数为
(2.10) (2.11) (2.12) (2.13) (2.14)
5
g(k)
=
(Np
Np + 1)a 2 Δt
[RMz
(k)
+
c]
(2.15)
2.2.3 电机输入电压 U 与风扇输出作用力 F 关系辨识
由于不能输出负电压给电机，M 序列输出的两极电压可取 0 与 24V。使用 PLC 编程实现 M 序列输出，通过 PWM 模块向电机供电。由于风扇对乒乓球的 作用力 F 无法直接测得，可以通过测乒乓球的高度 H，通过式(2.3)转换为 F。此 时得到 M 序列与 F 的离散值（注意数据采集时间和 Δt 相等才能使得两个数组同 维），根据式(2.15)可以推导出 U 与 F 之间的关系，其方框图如图 2-5 所示。
当前值的依赖关系的统计测度。对于离散系统，
Rxx
(τ
)
=
Lim
N →∞
1 N
[x(k)x(k
−τ
)]
(2.8)
由式(2.6)可知，如果知道输入的自相关函数 Rxx 和输入与输出的自相关函数 Rxy ，可以求得输入与输出之间的函数关系 g(t) 。但存在的问题是，由于存在积 分关系， g(t) 的求解过程并不简单。使用 M 序列伪随机信号作为激励输入，则
电机输入电压U U与F的响应函数 风扇对乒乓球作用力F g（k）