控制系统CAD案例

控制系统CAD 案例2
【例1】 已知单位反馈系统的被控对象的开环单位阶跃响应数据（存在随机扰动）在文件SourseData.mat 中，tout 是采样时间点，yout 是对应于采样时间点处的单位阶跃响应数据。

试设计串联补偿器，使得：① 在单位斜坡信号t t r )(的作用下，系统的速度误差系数K v ≥30sec -1；② 系统校正后的截止频率ωc ≥2.3sec -1；③ 系统校正后的相位稳定裕度P m ≥40°。

设计步骤：
1 根据原始数据得到被控对象的单位阶跃响应 被控对象的开环单位阶跃响应数据在文件SourseData.mat 中。

在MATLAB 指令窗中键入：
>> load SourseData
>> plot(tout,yout)
>> xlabel('tout')
>> ylabel('yout')
>> title('原始响应曲线')
由所给数据可看出，仿真步长为0.05s，仿真时间为10s。

以时间tout为横坐标，输出yout为纵坐标，得到单位阶跃响应随时间的变化曲线如图1所示。

图1 被控对象的单位阶跃响应
2系统模型的辨识
由于被控对象模型结构未知，题中仅给出了被控对象的单位阶跃响应。

因此需要根据原始的数据分析得出原系统的基本结构。

至于结构以及参数的最终确定，将使用MATLAB 的系统辨识工具箱来实现。

（1） 预处理
由图1可以看出被控对象模型中含积分环节，因此，对响应输出求导后进行辨识，即去积分。

再对辨识出的模型积分，即得被控对象模型。

对给定输出进行求导处理：
t
t y t t y t y ∆-∆+≈)()(d d 在MATLAB 指令窗中键入：
>> dy=(yout(2:end,:)-yout(1:end-1,:))/0.05;
% 差分代替导数
>> t=tout(1:end-1,:);
>> plot(t,dy)
>> xlabel('t')
>> ylabel('dy/dt')
>> title('原始响应的导数曲线')
以时间t 为横坐标，输出dy 为纵坐标，得到单位阶跃响应的导数随时间的变化曲线如图2所示。

图2被控对象的单位阶跃响应的导数
（2）系统模型结构的估计
系统辨识工具箱提供的模型结构选择函数有struc、arxstruc和selstruc。

函数struc 生成arx结构参数，调用格式为：
nn＝struc (Na, Nb, Nk)
其中，Na、Nb 分别为arx模型多项式A(q)、B(q) 的阶次范围；Nk为arx模型纯滞后的大小范围；nn为模型结构参数集构成的矩阵。

函数arxstruc用来计算arx模型结构的损失函数，即归一化的输出预测误差平方和，调用格式为：
v＝arxstruc (ze, zv, nn)
其中，ze＝[y u]为模型辨识的I/O数据向量或矩阵。

zv＝[yr ur]为模型验证的I/O数据向量或矩阵。

nn为多个模型结构参数构成的矩阵，nn的每行都具有格式：nn＝[na nb nk]。

v的第一行为各个模型结构损失函数值，后面的各行为模型结构参数。

函数selstruc 用来在损失函数的基础上进行模型结构选择，调用格式为：
[nn, vmod]＝selstruc (v, c)
其中v 由函数arxstruc获得的输出矩阵，为各个模型结构的损失函数。

c为可选参数，用于指定模型结构选择的方式。

根据图2所示曲线的形状初步估测被控对象的模型应该为二阶系统或者更高阶系统。

并且可以看出，纯滞后为0，故Nk恒为零。

在MATLAB指令窗中键入：
>> u=ones(size(dy));
>> Z=[dy,u];
>> v2=arxstruc(Z,Z,struc(2,0:2,0));
>> nn2=selstruc(v2,0);
>> v3=arxstruc(Z,Z,struc(3,0:3,0));
>> nn3=selstruc(v3,0);
>> v4=arxstruc(Z,Z,struc(4,0:4,0));
>> nn4=selstruc(v4,0);
得到如下结果：
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，则可以通过下面的指令辨识出系统的数学模型：
T=arx([y,u],[ n,m,d])
该函数将直接显示辨识的结果，且所得的T为一个结构体，其中T.A和T.B分别表示辨识得到的分子和分母多项式。

由给定的观测数据建立系统数学模型后，还需要进行检验，看模
型是否适用，如果不适用，则要修改模型结构，重新进行参数估计等。

MATLAB的系统辨识工具箱中用于模型验证和仿真的函数主要有compare、resid、pe、predict 和idsim。

此次实验主要用的是函数compare对模型进行验证。

函数compare可将模型的预测输出与对象实际输出进行比较。

验证过程与结果如下所示。

①对二阶系统的验证
在MATLAB指令窗中键入：
>>Z=iddata(dy,u,0.05);
>> M=arx(Z,[2,2,0]);
>> compare(M,Z)
得到图3所示结果。

图3二阶模型的匹配结果②对三阶系统的验证
在MATLAB指令窗中键入：
>>Z=iddata(dy,u,0.05);
>> M=arx(Z,[3,3,0]);
>> compare(M,Z)
得到图4所示结果。

图4三阶模型的匹配结果③对四阶系统的验证
在MATLAB指令窗中键入：
>>Z=iddata(dy,u,0.05);
>> M=arx(Z,[4,2,0]);
>> compare(M,Z)
得到图5所示结果。

图5四阶模型的匹配结果
④模型阶次的选取
从图3~图5可以看出，三阶和四阶模型的拟合率已高达90%，并且相差不大。

因此，为简便起见，本次试验选用三阶系统模型。

在MATLAB指令窗中键入：
>> M=oe(Z,[3,2,0]);
>> H=tf(M)
得到如下结果：
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": 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":
0.02835 s^2 + 1.037 s + 48.5
----------------------------
s^2 + 14.68 s + 48.51
Input groups:
Name Channels
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('原始数据', '辨识模型数据')
得到图7所示结果。

图7原始数据与辨识模型数据的比较从图7可以看出，由oe模型得到的传递函数的阶跃响应曲线并不能与原始数据很好的拟合。

故考虑在该传递函数的基础上，对其系数进行修改。

①由于传递函数分子的高次项系数远小于其一次项系数，因此将其忽略。

取分子为一次多项式，经验证对响应曲线影响并不大，但可以大大简化系统结构。

在MATLAB指令窗中键入：
>> num=[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('原始数据', '辨识模型数据')
得到图8所示结果。

图8原始数据与去除分子高次项的辨识模型数据的比较
②调整系统的零点和极点。

由于系统零点的作用为减小峰值时间，使系统响应速度加快，并且零点越接近虚轴，这种作用越显著。

然而只对分子一次项的系数进行修改时，效果并不理想。

因此，综合考虑，调节系统各个参数。

经过多
次尝试，在MATLAB指令窗中键入：>> num=[1.1 50.55];
>> den=[0.98 15.06 50.51];
>> G=tf(num,den)
>> sim('analyse',9.99)
>> plot(t,dy,'r'),hold on
>> plot(ymdl.time,ymdl.signals.values) >> legend('原始数据', '辨识模型数据') 得到如下结果和图9。

Transfer function:
1.1 s + 50.55
--------------------------
0.98 s^2 + 15.06 s + 50.51
图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))
得到如下结果：
fit =
92.8932
可见，模型数据与原始数据的匹配率达92.8932%。

（5）积分得出最终辨识模型
考虑到原系统含有积分环节，对辨识出的模型进行积分。

由给定原始数据看出，积分初值为0，从而得到最终辨识模型及其Simulink 模块 (图10, result.mdl)：
)51.5006.1598.0(55.501.1)(20+++=s s s s s G
图10 最终辨识模型的Simulink 模块
比较辨识出的最终辨识模型的阶跃响应与原始数据的匹配率，在MATLAB 指令窗中键入：
>> num=[1.1,50.55];
>> den=[0.98,15.06,50.51];
>> sim('result',10)
>> plot(tout,yout,'r'),hold on
>> plot(ymdl.time,ymdl.signals.values)
>> legend('原始数据', '辨识模型数据')
>> yhat=ymdl.signals.values;
>> ymean=zeros(size(yout));
>> ymean(:)=mean(yout);
>> fit=100*(1-norm(yout-yhat,2)/norm(yout-ymean,2)) 得到图11和如下结果：
fit =
99.0656
图11最终辨识模型的阶跃响应与原始数据之比较
可见，二者之间的匹配率高达99.0656%，满足要求。

3 根据最终辨识模型设计控制器法1（经典方法）
（1） 确定开环增益K 0
根据控制理论，给定被控对象为I 型系统，单位斜坡响应的速度误差系数K v =K =K 0≥30sec -1，其中K 0是系统的开环增益。

取K 0=30sec -1，则被控对象的传递函数为
)51.5006.1598.0()55.501.1(30)(20+++=s s s s s G
在MATLAB 指令窗中键入：
>> clear
>> num=[30*1.1,30*50.55];
>> den=[0.98,15.06,50.51,0];
>> sys0=tf(num,den);
>> figure(1); margin(sys0),grid
>> figure(2); step(feedback(sys0,1)), grid
得到如图12所示的系统bode 图和如图13所示的闭环系统的单位阶跃响应。

图12未校正系统的bode图
根据计算可知未校正系统的频域性能指标为：对数幅值稳定裕度G m0=-2.27dB
-180°穿越频率ωg0=8.8 rad/sec
相位稳定裕度P m0=-4.91°
截止频率ωc0=9.92 rad/sec
图13 未校正系统的单位阶跃响应
显见，这样的系统根本无法工作。

这一点也可从如图13所示发散振荡的阶跃响应曲线看出，系统必须进行校正。

采用两种方法设计串联校正器。

（2） 采用bode 图法设计串联校正器
由于给定的开环截止频率ωc ≥2.3rad/sec ，远小于ωc0=9.92rad/sec ，可以通过压缩频带宽度来改善相位裕度，因此采用串联滞后校正是合理的。

令校正器的传递函数为
s T s T Ts Ts s G c 211111)(++=++=β
显然，应有β>1。

① 确定新的开环截止频率ωc
希望的相位稳定裕度P m ≥40°，所以根据滞后校正的设计方法，应有
c
m P ωωωϕπ=+=+)()5~2( 其中，（2°~5°）是附加相位补偿角。

取其等于5°，则有
)(540c ωϕπ+=+
于是，有
135)(-=c ωϕ
从图12的对数相频特性图上可以取得对应于-135°的ω为3.09rad/sec （见图14），即有
sec /rad 09.3=c ω
② 计算高频衰减率β
从图12的对数幅频特性图上可以取得对应于sec /rad 09.3=c ω的L 为17.9dB （见图14）。

于是，有
9.17lg 20)09.3(==βL
从而
85.7=β
图14 从未校正系统的bode 图上取相关值
③ 计算两个转折频率ω1和ω2
0394.01110
1309.01009.3)101~51(1212=======ωβ
βωωωT T c ）（取 所以，校正器的传递函数为
s s s G c 41.25124.31)(++=
校正后的开环传递函数为
)141.25)(51.5006.1598.0()124.3)(55.501.1(30)()()(20+++++=⋅=s s s s s s s G s G s G c
④ 校验校正后系统的频域性能是否满足要求
在MATLAB 指令窗中键入：
>> n=[3.24,1];d=[25.41,1];sys1=tf(n,d);
>> sys=sys0*sys1;
>> figure(1); margin(sys),grid
>> figure(2); step(feedback(sys,1)), grid
得到如图15所示的校正后系统的bode 图和如图16所示的闭环系统的单位阶跃响应。

图15 校正后系统的bode 图
图16 校正后系统的单位阶跃响应
根据计算可知校正后系统的频率性能指标为：
对数幅值稳定裕度G m=14.8dB
-180°穿越频率ωg=8.43 rad/sec
相位稳定裕度P m=40°
截止频率ωc=3.12 rad/sec
显然，系统校正后的相位稳定裕度和截止频率均满足给定要求。

4 根据最终辨识模型设计控制器法2（优化方法） 采用参考模型法进行设计。

首要工作是建立一个理想的参考模型，它是根据给定的系统性能指标要求仔细地设计出来的。

这个模型可能与所研究系统的模型完全不同，但它在给定的输入信号作用下能够产生期望的输出响应，即它能满足所研究系统的所有性能指标要求。

以该模型作为参考，将系统的输入信号同时加到实际系统和参考模型上，如图17所示。

然后，将参考模型的输出与实际系统的输出进行分析、比较，并根据它们之间的偏差对系统的参数进行调整，使得在某种意义下偏差尽可能小，即
⎰
→f t t t e 02
min d )( 式中 )()()(t y t y t e m -=
图17 采用参考模型法进行参数寻优
根据偏差)(t e 可以利用各种优化方法来调整系统参数，使得系统输出响应跟踪或逼近参考模型输出。

显然，系统响
应特性的好坏在很大程度上取决于参考模型的特性。

所以，选择一个理想的参考模型尤为重要。

那么，何为“理想”的参考模型呢？首先，它要满足系统的所有性能指标要求。

其次，它的结构应该简单，便于数学处理。

很明显，一个极为严密的，包括每一现象细节的模型的结构会过于复杂，而建立这样一个模型本身也是一件非常困难的事。

从工程应用的观点而言，需要的是一个足以反映所研究的系统本质并且尽可能简单的模型。

这里，选择如图18所示的典型二阶系统作为参考模型。

图18 典型二阶系统
并且考虑到其开环频率特性应当和给定的性能指标要求尽
可能一致，选择104.02==n ωζ，。

对应于图17的Simulink
模块如图19所示。

图19 参考模型法目标函数的Simulink模块（EXAM.mdl）利用自编的控制系统优化软件对控制器参数进行仿真寻优，相关模型参数初始值、优化参数值及寻优结果如20所示。

图20参数设置及仿真寻优结果
即校正器的传递函数为
s s s G c 8.144112.151)(++=
在MATLAB 指令窗中键入：
>> n=[15.12,1];d=[144.8,1];sys1=tf(n,d);
>> sys=sys0*sys1;
>> figure(1); margin(sys),grid
>> figure(2); step(feedback(sys,1)), grid
得到如图21所示的校正后系统的bode 图和如图22所示的闭环系统的单位阶跃响应。

图21 校正后系统的bode 图（参考模型法）
图22 校正后系统的单位阶跃响应（参考模型法）根据计算可知校正后系统的频率性能指标为：
对数幅值稳定裕度G m=17.2dB
-180°穿越频率ωg=8.72 rad/sec
相位稳定裕度P m=49.2°
截止频率ωc=2.68 rad/sec
显然，经过参数寻优后系统校正的相位稳定裕度和截止频率均满足给定要求。

对比图16和图22，可以看出经过参数寻优的控制器优于
常规方法设计出来的控制器。

5 控制器效果的最终仿真检验
控制器是根据辨识出的模型设计的，它能否真正应用到实际系统中，达到满意的控制效果，需要进行检验。

事实上，实际被控对象的开环传递函数为
)12.0)(11.0(1)(0++=s s s s G
原始给定数据是在存在随机扰动（均值为0，方差为0.1，信噪比为0.001）通过仿真的方法（采样间隔为0.05sec ，采样时间为10sec ）得到的。

系统采样的Simulink 模块如图23所示。

图23 系统采样的Simulink 模块
（1）不加随机扰动时的控制效果
不加随机扰动时的Simulink模块及仿真结果分别如图24和图25所示。

图24 不加随机扰动时的Simulink模块（GcTest1）
图25 不加随机扰动时的仿真结果（黄线是图24中上面结果）
显见，经过参数优化后的控制器效果更好一些。

（2）加入随机扰动后的控制效果
加入随机扰动后的Simulink模块及仿真结果分别如图26和图27所示。

图26 加入随机扰动后的Simulink模块（GcTest2）
图27 加入随机扰动后的仿真结果
可见，仍然是经过参数优化后的控制器效果更好一些，其抗干扰能力更强。

现代控制理论的基本设计方法是状态反馈控制器设计，包括状态反馈增益矩阵的设计和状态观测器增益矩阵的设计，如图28所示。

图28 状态反馈控制器
开环系统的状态空间模型为
⎩⎨⎧+=+=)()()()()()(t t t t t t Du Cx y Bu Ax x (1)
由图28，将)()()(t t t Kx v u -=代入上式中，则在状态反馈增益矩阵K 下，闭环系统的状态空间模型可写为
⎩
⎨⎧+-=+-=)()()()()()()()(t t t t t t Dv x DK C y Bv x BK A x (2) 如果开环系统完全可控，则选择合适的K 矩阵，可以将闭环系统矩阵（A -BK ）的特征值配置到任意地方（当然，还要满足共轭复数的约束）。

在极点配置过程中，假设了系统的所有状态都可被量测并用于反馈。

然而在实际情况中，并非所有状态都可被直接量
测，这就需要对不可直接量测的状态)(ˆt x
进行估计。

对不可量测的状态进行估计称为状态观测，观测状态变量的装置（或计算机程序）称为状态观测器。

根据现代控制理论，状态观测器满足的误差状态方程为
)(~)()(~t t x LC A x -= (3)
其中，)()(ˆ)(~t t t x x
x -=。

只要开环系统是完全可观的，可以将矩阵（A -LC ）的特征值配置到任意地方，从而可以使
观测出的状态)(ˆt x
可以逼近原系统的状态)(t x （只要将该矩阵的特征值全部配置在左半s 平面上）。

将（1）式代入(3)式，并将观测出的状态)(ˆt x
作为输出反馈量)(ˆt y
，则得观测器满足的状态空间模型为 ⎩
⎨⎧=++-=)(ˆ)(ˆ)()()(ˆ)()(ˆt t t t t t x y Ly Bu x LC A x （4） 令
[][]T t t )(ˆ)((t)ˆ，ˆˆˆy u u 0D I C L B B LC A A
====-=，，，
则有
⎪⎩
⎪⎨⎧+=+=)(ˆˆ)(ˆˆ)(ˆ)(ˆˆ)(ˆˆ)(ˆt t t t t t u D x C y u B x A
x （5）
MATLAB 的控制系统工具箱提供了一系列函数，可用于状态反馈控制器的辅助设计。

它们常用的调用格式为 ctrb(sys) obsv(sys)
K=place(A, B, p)
L=place(A',C', p)
说明
∙ 以上函数要求对象模型sys 为状态空间模型。

∙ 第1种格式求出系统的可控性矩阵。

∙ 第2种格式求出系统的可观性矩阵。

∙ 第3种格式求出状态反馈增益矩阵。

A 为系统的状态矩
阵；B 为系统的输入矩阵；p 为指定的闭环系统极点；输出宗量K 为状态反馈增益矩阵。

∙ 第4种格式求出状态观测器增益矩阵。

A'为系统的状态
矩阵A 的转置；C'为系统的观测矩阵C 的转置；p 为指定的闭环系统极点；输出宗量L 为状态观测器增益矩阵。

【例2】 已知某控制系统的开环传递函数为
)3)(2)(1(10)(+++=s s s s G
试判断系统的可控性并求状态反馈增益矩阵K ，使得系统的闭环极点为：101-=λ，22j 23,2±-=λ。

在系统的闭环极
点不变的情况下，计算状态观测器增益矩阵L。

以状态反馈和状态观测器构成闭环回路完成单位阶跃仿真。

（1）判别系统的可控性并求出状态反馈增益矩阵K
编写M脚本文件exam2_1.m如下：
% exam2_1
clear
k=10; z=[]; p=[-1;-2;-3]; s=zpk(z,p,k);
sys=ss(s);
A=sys.a; B=sys.b; C=sys.c; % 求出所需的系统矩阵[n,m]=size(A); Qc = ctrb(A,B); % 求出系统的可控性矩阵if rank(Qc)==n
disp('系统是可控的。

')
disp('系统的可控性矩阵：')
disp(Qc)
pp=[-10,-2+2*sqrt(2)*j, -2-2*sqrt(2)*j];
K=place(A,B,pp); % 求出状态反馈增益矩阵disp('状态反馈增益矩阵：')
disp(K)
else
disp('系统是不可控的。

')
end
运行exam2_1.m后，得到如下结果：
系统是可控的。

系统的可控性矩阵：
0 0 4
0 4 -20
4 -12 36
状态反馈增益矩阵：
20.2500 4.2500 2.0000
（2）判别系统的可观性并求出状态观测器增益矩阵L 编写M脚本文件exam2_2.m如下：
% exam9_17_2
Qb=obsv(A,C);
if rank(Qb)==n
disp('系统是可观的。

')
disp('系统的可观性矩阵：')
disp(Qb)
L=place(A',C',pp); % 求出状态观测器增益矩阵disp('状态观测器增益矩阵(转置)：')
disp(L)
% 计算状态观测器的系数矩阵
AA=A-L'*C
BB=[B,L']
CC=eye(3)
DD=zeros(3,2)
else
disp('系统是不可观的。

')
end
运行exam2_2.m后，得到如下结果：系统是可观的。

系统的可观性矩阵：
2.5000 0 0
-2.5000 2.5000 0
2.5000 -7.5000 2.5000
状态观测器增益矩阵(转置)：
3.2000 0.4000 25.2000
AA =
-9.0000 1.0000 0
-1.0000 -2.0000 1.0000
-63.0000 0 -3.0000
BB =
0 3.2000
0 0.4000
4.0000 2
5.2000
CC =
1 0 0
0 1 0
0 0 1
DD =
0 0
0 0
0 0
说明
AA ，BB ，CC 和DD 即分别为（5）式中D C B A ˆˆˆˆ和，，。

根据所得结果，构造Simulink 模型如图29所示。

图29具有状态观测器的状态反馈系统的Simulink 模型
启动仿真后，得到如图30所示仿真结果。

图30具有状态观测器的状态反馈系统的阶跃响应。