设计步骤：
1 根据原始数据得到被控对象的单位阶跃响应
被控对象的开环单位阶跃响应数据在文件 SourseData.mat 中。
在 MATLAB 指令窗中键入： >> load SourseData >> plot(tout,yout) >> xlabel('tout') >> ylabel('yout') >> title('原始响应曲线')
10
Transfer function from input "u1" to output "y1":
0.02835 z^2 + 0.02835 z + 0.02835
---------------------------------
z^2 - 1.395 z + 0.4799
Transfer function from input "v@y1" to output "y1":
在 MATLAB 指令窗中键入：
>> num=[1.037 48.5]; >> den=[1 14.68 48.51];
13
>> sim('analyse',9.99) >> plot(t,dy,'r'),hold on >> plot(ymdl.time,ymdl.signals.values) >> legend('原始数据', '辨识模型数据')
T=arx([y,u],[ n,m,d]) 该函数将直接显示辨识的结果，且所得的T为一个结构体， 其中T.A和T.B分别表示辨识得到的分子和分母多项式。由给 定的观测数据建立系统数学模型后，还需要进行检验，看模
6
型是否适用，如果不适用，则要修改模型结构，重新进行参 数估计等。MATLAB的系统辨识工具箱中用于模型验证和仿 真的函数主要有compare、resid、pe、predict 和idsim。此 次实验主要用的是函数compare对模型进行验证。函数 compare可将模型的预测输出与对象实际输出进行比较。验 证过程与结果如下所示。
控制系统 CAD 案例 2
【例 1】 已知单位反馈系统的被控对象的开环单位阶跃 响应数据（存在随机扰动）在文件 SourseData.mat 中，tout 是采样时间点，yout 是对应于采样时间点处的单位阶跃响应 数据。
试设计串联补偿器，使得：① 在单位斜坡信号 r(t) t 的
作用下，系统的速度误差系数 Kv≥30sec-1；② 系统校正后 的截止频率ωc≥2.3sec-1；③ 系统校正后的相位稳定裕度 Pm ≥40°。
v＝arxstruc (ze, zv, nn)
其中，ze＝[y u]为模型辨识的 I/O 数据向量或矩阵。zv＝[yr ur]为模型验证的 I/O 数据向量或矩阵。nn 为多个模型结构 参数构成的矩阵，nn 的每行都具有格式：nn＝[na nb nk]。 v 的第一行为各个模型结构损失函数值，后面的各行为模型 结构参数。
得到如下结果： nn2 = 2 1 0 nn3 = 3 2 0 nn4 = 4 1 0
于是，去除积分环节后的模型阶数为：二阶系统[2 1 0]、三 阶模型[3 2 0]和四阶模型[4 1 0]。
（3） 系统模型结构的确定 为了确定模型结构以及参数，使用 MATLAB 的系统辨识 工具箱中已有的辨识函数 arx()。辨识函数 arx()的使用方法 是：如果一直输入信号的列向量 u，输出信号的列向量 y， 并选定了系统的分子多项式阶次 m-1，分母多项式阶次 n 及 系统的纯滞后 d，则可以通过下面的指令辨识出系统的数学 模型：
得到图 7 所示结果。
12
图 7 原始数据与辨识模型数据的比较
从图 7 可以看出，由 oe 模型得到的传递函数的阶跃响应 曲线并不能与原始数据很好的拟合。故考虑在该传递函数的 基础上，对其系数进行修改。
① 由于传递函数分子的高次项系数远小于其一次项系 数，因此将其忽略。取分子为一次多项式，经验证对响应曲 线影响并不大，但可以大大简化系统结构。
2
（1） 预处理
由图 1 可以看出被控对象模型中含积分环节，因此，对响
应输出求导后进行辨识，即去积分。再对辨识出的模型积分，
即得被控对象模型。
对给定输出进行求导处理：
dy y(t t) y(t)
dt
t
在 MATLAB 指令窗中键入：
>> dy=(yout(2:end,:)-yout(1:end-1,:))/0.05; % 差分代替导数
① 对二阶系统的验证 在 MATLAB 指令窗中键入： >> Z=iddata(dy,u,0.05); >> M=arx(Z,[2,2,0]); >> compare(M,Z) 得到图 3 所示结果。
7
图 3 二阶模型的匹配结果 ② 对三阶系统的验证
在 MATLAB 指令窗中键入： >> Z=iddata(dy,u,0.05); >> M=arx(Z,[3,3,0]); >> compare(M,Z)
0.02835 s^2 + 1.037 s + 48.5
----------------------------
s^2 + 14.68 s + 48.51
Input groups:
Name
Channels
ห้องสมุดไป่ตู้11
Measured
1
（4）传递函数模型系数的修正 得到模型的传递函数后，需要查看其开环单位阶跃响应， 并与原始数据进行对比，看看是否吻合。为此，建立仿真模 块analyse.mdl，如图6所示。
图 6 辨识模型的 Simulink 模块
并在 MATLAB 指令窗中键入：
>> num=[0.02835,1.037, 48.5]; >> den=[1 ,14.68 ,48.51]; >> sim('analyse',9.99) >> plot(t,dy,'r'),hold on >> plot(ymdl.time,ymdl.signals.values) >> legend('原始数据', '辨识模型数据')
函数 selstruc 用来在损失函数的基础上进行模型结构选 择，调用格式为：
[nn, vmod]＝selstruc (v, c)
其中v 由函数arxstruc获得的输出矩阵，为各个模型结构的 损失函数。c为可选参数，用于指定模型结构选择的方式。
根据图 2 所示曲线的形状初步估测被控对象的模型应该 为二阶系统或者更高阶系统。并且可以看出，纯滞后为 0， 故 Nk 恒为零。
1
由所给数据可看出，仿真步长为 0.05s，仿真时间为 10s。 以时间 tout 为横坐标，输出 yout 为纵坐标，得到单位阶跃 响应随时间的变化曲线如图 1 所示。
图 1 被控对象的单位阶跃响应
2 系统模型的辨识
由于被控对象模型结构未知，题中仅给出了被控对象的单 位阶跃响应。因此需要根据原始的数据分析得出原系统的基 本结构。至于结构以及参数的最终确定，将使用 MATLAB 的系统辨识工具箱来实现。
得到如下结果和图 9。 Transfer function:
1.1 s + 50.55 -------------------------0.98 s^2 + 15.06 s + 50.51
15
图 9 原始数据与调整后的辨识模型数据的比较 ③ 计算匹配率 根据公式： FIT = 100(1-norm(Y-YHAT)/norm(Y-mean(Y))) (in %) 在 MATLAB 指令窗中键入： >> yhat=ymdl.signals.values; >> ymean=zeros(size(dy)); >> ymean(:)=mean(dy); >> fit=100*(1-norm(dy-yhat,2)/norm(dy-ymean,2))
0.009445
Input groups:
Name
Channels
Measured
1
Noise
2
Sampling time: 0.05
上述模型为离散时间系统模型，为了得到连续时间系统模 型，在 MATLAB 指令窗中键入：
>> G=d2c(H(1))
得到如下结果：
Transfer function from input "u1" to output "y1":
nn＝struc (Na, Nb, Nk) 其中，Na、Nb 分别为 arx 模型多项式 A(q)、B(q) 的阶次 范围；Nk 为 arx 模型纯滞后的大小范围；nn 为模型结构参 数集构成的矩阵。
4
函数 arxstruc 用来计算 arx 模型结构的损失函数，即归一 化的输出预测误差平方和，调用格式为：
>> t=tout(1:end-1,:); >> plot(t,dy) >> xlabel('t') >> ylabel('dy/dt') >> title('原始响应的导数曲线')
以时间 t 为横坐标，输出 dy 为纵坐标，得到单位阶跃响 应的导数随时间的变化曲线如图 2 所示。
3