控制系统模型与转换

1、控制系统模型与转换第1页，共65页，2022年，5月20日，15点10分，星期三1控制系统的数学模型在控制系统的研究中有着相当重要的地位，要对系统进行仿真处理，首先应当知道系统的数学模型，然后才可以对系统进行模拟。同样，如果知道了系统的模型，才可以在此基础上设计一个合适的控制器，使得系统响应达到预期的效果，从而符合工程实际的需要。在线性系统理论中，一般常用的数学模型形式有：传递函数模型（系统的外部模型）、状态方程模型（系统的内部模型）、零极点增益模型和部分分式模型等。这些模型之间都有着内在的联系，可以相互进行转换。第2页，共65页，2022年，5月20日，15点10分，星期三2系统的分类 按系

2、统性能分：线性系统和非线性系统；连续系统和离散系统；定常系统和时变系统；确定系统和不确定系统。1、线性连续系统：用线性微分方程式来描述，如果微分方程的系数为常数，则为定常系统；如果系数随时间而变化，则为时变系统。今后我们所讨论的系统主要以线性定常连续系统为主。2、线性定常离散系统：离散系统指系统的某处或多处的信号为脉冲序列或数码形式。这类系统用差分方程来描述。3、非线性系统：系统中有一个元部件的输入输出特性为非线性的系统。第3页，共65页，2022年，5月20日，15点10分，星期三33.1 连续线性系统的数学模型3.11 线性系统的传递函数一、连续系统的传递函数模型连续系统的传递函数如下：

3、对物理可实现系统来说，必须满足mn,这种情况下又称系统为正则（proper）的。若ms=tf(s); G=(s3+4*s2+3*s+2)/(s2*(s+1)*(s+4)2+4)第8页，共65页，2022年，5月20日，15点10分，星期三8使用set命令可以获得tf( )函数的详细信息第9页，共65页，2022年，5月20日，15点10分，星期三9若有传递函数G，也可用tfdata( )函数提取系统的分子和分母多项式。 如：上例中的G当系统为多变量系统时，可用下面语句提取第1路输入和第1路输出的传递函数num=G.num1,1;den=G.den1,1第10页，共65页，2022年，5月20日

4、，15点10分，星期三103.1.2 线性系统的状态方程模型 线性时不变系统可以用一组一阶微分方程来描述，其矩阵形式即为现代控制理论中常用的状态空间表示法。状态方程与输出方程的组合称为状态空间表达式，又称为动态方程，经典控制理论用传递函数将输入输出关系表达出来，而现代控制理论则用状态方程和输出方程来表达输入输出关系，揭示了系统内部状态对系统性能的影响。在MATLAB中，系统状态空间用（A,B,C,D)矩阵组表示: G=ss(A, B, C, D)x为状态向量；u为输入向量；y为输出向量A为系统矩阵；B为输入矩阵；C为输出矩阵；D为前馈矩阵；第11页，共65页，2022年，5月20日，15点10

5、分，星期三11例4：系统为一个两输入两输出系统A=1 6 9 10; 3 12 6 8; 4 7 9 11; 5 12 13 14;B=4 6; 2 4; 2 2; 1 0;C=0 0 2 1; 8 0 2 2; D=zeros(2,2);G=ss(A, B, C, D)第12页，共65页，2022年，5月20日，15点10分，星期三12获取状态方程参数可使用ssdata( )函数。带有时间延迟的状态方程可使用语句：G=ss(A, B, C, D, ioDelay, )第13页，共65页，2022年，5月20日，15点10分，星期三133.1.3 线性系统的零极点模型 零极点模型实际上是传递函

6、数模型的另一种表现形式，其原理是分别对原系统传递函数的分子、分母进行分解因式处理，以获得系统的零点和极点的表示形式。其中，K为系统增益，zi为零点，pj为极点。在MATLAB中零极点增益模型zpk( )函数表示。即： Z=z1;z2;zm P=p1;p2;.;pn G=zpk(Z, P, K)系统G的零极点可以直接调用pzmap(G)函数显示。第14页，共65页，2022年，5月20日，15点10分，星期三14例5：输入该零极点模型：第15页，共65页，2022年，5月20日，15点10分，星期三15零极点增益模型：num=1,11,30,0;den=1,9,45,87,50; z,p,k=t

7、f2zp(num,den)z= 0 -6 -5p= -3.0000+4.0000i -2.0000 -1.0000k= 1结果表达式：第16页，共65页，2022年，5月20日，15点10分，星期三163.1.4 多变量系统的传递函数矩阵模型 多变量系统的状态方程可由ss( )函数输入，也可直接给出传递函数矩阵第17页，共65页，2022年，5月20日，15点10分，星期三17例6：其输入命令如下或 g11=tf(0.1134,1.78,4.48,1); g12=tf(0.924,2.07,1); g21=tf(0.3378,0.361,1.09,1); g22=tf(-0.318,2.93,

8、1); G=g11,g12;g21,g22；G.ioDelay=0.72,0;0.3,1.29第18页，共65页，2022年，5月20日，15点10分，星期三183.2 离散系统模型3.2.1 离散系统传递函数模型可以用tf( )函数输入。 num=b0,b1,bn-1,bn; den=1,a1,a2,an; H=tf(num, den, Ts,T);其中，T为时既采样周期。也可定义算子z=tf(z,T),然后用数学表达式形式输入传递函数。第19页，共65页，2022年，5月20日，15点10分，星期三19例7：离散系统传递函数为其采样周期T=0.1s,则输入命令为第20页，共65页，2022

9、年，5月20日，15点10分，星期三20Notice:该模型还可用z=tf(z,0.1);H=(0.31*z3-0.57*z2+038*z-0.088)/ (z4-3.23*z3+3.98*z2-2.22*z+0.47)若将离散系统的传递函数分子和分母同时除以zn,则系统可变换为用q取代z-1,则有第21页，共65页，2022年，5月20日，15点10分，星期三21 对于这种多用于表示滤波器的形式，其输入命令有 num= bn, bn-1, , b1, b0; den=an, an-1, , a1, 1; H=tf(num, den, Ts,T, Variable, q);例8：输入采用如下语

10、句z=0.67+0.35i; 0.67-0.35i; 0.48;p=0.85+0.18i; 0.85-0.18i; 0.75+0.19i;0.75-0.19i;H=zpk(z, p, 1, Ts, 0.1)第22页，共65页，2022年，5月20日，15点10分，星期三223.2.2 离散状态方程模型离散系统的状态方程可用ss( )函数输入，即 H=ss (F, G, C, D, Ts, T)若系统带有时间延迟则其输入语句为 H=ss (F, G, C, D, Ts, T, ioDelay, m)第23页，共65页，2022年，5月20日，15点10分，星期三233.3 框图描述系统的化简3.

11、3.1 控制系统的典型连接机构串联系统模型由G=G2*G1求出。G1(s)G2(s)G2(s)G1(s)uuyy第24页，共65页，2022年，5月20日，15点10分，星期三24b) 并联系统模型由G=G2+G1求出。G1(s)G2(s)yG1(s)+G2(s)yuu第25页，共65页，2022年，5月20日，15点10分，星期三25c) 反馈 但这样得出的模型阶次可能高于实际的阶次，需要用minreal( )函数求取得出模型的最小实现形式。或：-G(s)K(s)uy+G(s)K(s)uy第26页，共65页，2022年，5月20日，15点10分，星期三26 a,b,c,d=feedback(

12、a1,b1,c1,d1,a2,b2,c2,d2) 将两个系统按反馈方式连接，一般而言，系统1为对象，系统2为反馈控制器。 a,b,c,d=feedback(a1,b1,c1,d1,a2,b2,c2,d2,sign)系统1的所有输出连接到系统2的输入，系统2的所有输出连接到系统1的输入，sign用来指示系统2输出到系统1输入的连接符号，sign缺省时，默认为负，即sign= -1。总系统的输入/输出数等同于系统1。系统2系统1+y1u2y2u1第27页，共65页，2022年，5月20日，15点10分，星期三27num,den=feedback(num1,den1,num2,den2,sign)

13、可以得到类似的连接，只是子系统和闭环系统均以传递函数的形式表示。sign的含义与前述相同。G=feedback(G1,G2)第28页，共65页，2022年，5月20日，15点10分，星期三28例9：考虑下图的系统框图，试求系统的总模型H(S)Gc(S)+R(S)G(S)Y(S)第29页，共65页，2022年，5月20日，15点10分，星期三29例10：上例的反馈系统，假设受控对象模型为双输入、双输出的状态方程模型控制器为传递函数矩阵，为对角矩阵，其子传递函数为g11(s)=(6s+2)/s，g22 (s)=(3s+5)/s，反馈环节为单位矩阵，试求系统模型。第30页，共65页，2022年，5月

14、20日，15点10分，星期三30思考题：1）exp3_2.m 系统1为： 系统2为：求按串联、并联、正反馈、负反馈连接时的系统状态方程及系统1按单位负反馈连接时的状态方程。第31页，共65页，2022年，5月20日，15点10分，星期三313.3.3 节点移动的等效变换G(S)G(S)X1X2X2X2X1X2G(S)G(S)X2X1X11、分支点等效移动规则前向移动节点后向移动节点G(S)1/G(S)X1X2X1第32页，共65页，2022年，5月20日，15点10分，星期三32G(S)1/G(S)X1X2X3-G(S)X1X2X3-x2x3x1G(s)G(s)G(s)x1x2x3 (1) 相

15、加点前移2、相加点等效移动规则相加点前移，在移动支路中串入所越过的传递函数的倒数方框 (2) 相加点后移 相加点后移,在移动支路中串入所越过的传递函数方框。 第33页，共65页，2022年，5月20日，15点10分，星期三33G1G2G3G4G5G7G6-BAG1G2G3G4G4G5G7G6-第34页，共65页，2022年，5月20日，15点10分，星期三34G1G2G3G4G5G7G6/G4-第35页，共65页，2022年，5月20日，15点10分，星期三35注：以上语句要调用自定义feedback( )函数 function H=feedback(A,B,key) if nargin=2;

16、key=-1;end H=A/(sym(1)-key*A*B);H=simple(H);将其放置在work目录下的sym子目录，则可以直接处理符号模型的化简问题。第36页，共65页，2022年，5月20日，15点10分，星期三363.3.3 纯时间延迟环节的处理 对于带有时间延迟的模块G1(s)和G2(s)，若两模块串联，则系统模型是两个模块的延时相加。但当两模块并联时，要求其延迟必须相同。若不同，可先对模块进行近似，将纯时间延迟近似为高阶有理式。 作者编写了paderm( )函数（见附录1） num,den=paderm(T, r, m)其中，T为延迟时间常数，r，m分别为近似传递函数分子和

17、分母的阶次。第37页，共65页，2022年，5月20日，15点10分，星期三37完全近似：分母近似：注：分母近似比完全近似的结果更加精确。第38页，共65页，2022年，5月20日，15点10分，星期三38例12：假设系统的开环传函为且整个系统为单位负反馈结构构成，试选择分子0阶、分母2阶的Pade近似，其闭环系统模型为第39页，共65页，2022年，5月20日，15点10分，星期三393.4系统模型的相互转换3.4.1 连续模型和离散模型的相互转换如果连续系统由传函给出，可用s=2(z-1)/(T(Z+1)代入连续系统的传函模型，获得离散系统的传递函数模型。在MATLAB中，可采用c2d(

18、)函数将系统离散化。格式：Gd=c2d(G, Ts)第40页，共65页，2022年，5月20日，15点10分，星期三40例13：多变量系统模型如下，假设采样周期T=0.1s，求离散化的状态方程。第41页，共65页，2022年，5月20日，15点10分，星期三41第42页，共65页，2022年，5月20日，15点10分，星期三42例14：假设连续系统的数学模型为采样周期T=0.1s，求离散化的状态方程。 G=tf(1,1 3 3 1,ioDelay,0.5); G1=c2d(G,0.1,zoh) %零阶保持器变换Transfer function: 0.0001547 z2 + 0.000574

19、 z + 0.0001331z(-5) * - z3 - 2.715 z2 + 2.456 z - 0.7408 Sampling time: 0.1第43页，共65页，2022年，5月20日，15点10分，星期三43 G2=c2d(G,0.1,foh) %一阶保持器变换 Transfer function: 3.925e-005 z3 + 0.0004067 z2 + 0.000383 z + 3.278e-005z(-5) * - z3 - 2.715 z2 + 2.456 z - 0.7408Sampling time: 0.1 G3=c2d(G,0.1,tustin) %Tustin变

20、换 Transfer function: 0.000108 z3 + 0.0003239 z2 + 0.0003239 z + 0.000108z(-5) * - z3 - 2.714 z2 + 2.456 z - 0.7406Sampling time: 0.1第44页，共65页，2022年，5月20日，15点10分，星期三44由离散系统变换连续系统，可代入 z=(1+sT/2)/(1-sT/2)。 在MATLAB中可采用d2c( )函数。例15：采用d2c( )函数还原例13中的连续系统。 A=0 1 0;0 0 1;-6 -11 -6; B=1 0;2 -1;0 2;C=1 -1 0;2

21、 1 -1; D=zeros(2,2); G=ss(A,B,C,D); T=0.1;Gd=c2d(G,T); G1=d2c(Gd) a = x1 x2 x3 x1 4.025e-015 1 7.98e-016 x2 -1.804e-015 -1.527e-015 1 x3 -6 -11 -6 b = u1 u2 x1 1 -5.204e-017 x2 2 -1 x3 -5.551e-016 2 c = x1 x2 x3 y1 1 -1 0 y2 2 1 -1 d = u1 u2 y1 0 0 y2 0 0Continuous-time model.第45页，共65页，2022年，5月20日，1

22、5点10分，星期三453.4.2 系统传递函数的获取已知系统的状态方程，可通过tf( )函数获得系统的传递函数。例16：已知例13中的系统状态方程，求其传递函数。 A=0 1 0;0 0 1;-6 -11 -6;B=1 0;2 -1;0 2; C=1 -1 0;2 1 -1; D=zeros(2,2);G=ss(A,B,C,D); G1=tf(G)Transfer function from input 1 to output. -s2 - 4 s + 29 #1: - s3 + 6 s2 + 11 s + 6 4 s2 + 56 s + 52 #2: - s3 + 6 s2 + 11 s +

23、 6Transfer function from input 2 to output. s2 + 3 s - 4 #1: - s3 + 6 s2 + 11 s + 6 -3 s2 - 17 s - 14 #2: - s3 + 6 s2 + 11 s + 6第46页，共65页，2022年，5月20日，15点10分，星期三463.4.3 控制系统的状态方程实现 系统的状态方程可由ss( )函数实现。例17：系统传函为试将其转换为状态方程形式。 g11=tf(0.1134,1.78,4.48,1); g12=tf(0.924,2.07,1); g21=tf(0.3378,0.361,1.09,1);

24、 g22=tf(-0.318,2.93,1); G=g11,g12;g21,g22；G.ioDelay=0.72,0;0.3,1.29G1=ss(G)第47页，共65页，2022年，5月20日，15点10分，星期三473.4.4 状态方程的最小实现 当系统存在相同的零极点时，可以进行对消，完全对消相同零极点的形式称为系统的最小实现。如： 最后得到的一阶模型G(s)=5/(s+4)，即为系统的最小实现。在MATLAB中，可采用minreal( )函数获得。第48页，共65页，2022年，5月20日，15点10分，星期三48A=-6 -1.5 2 4 9.5;-6 -2.5 2 5 12.5;-5

25、 0.25 -0.5 3.5 9.75;-1 0.5 0 -1 1.5;-2 -1 1 2 3;B=6 4;5 5;3 4;0 2;3 1;C=2 0.75 -0.5 -1.5 -2.75;0 -1.25 1.5 1.5 2.25;D=0 0;0 0;G=ss(A,B,C,D);G1=minreal(G)第49页，共65页，2022年，5月20日，15点10分，星期三493.5 线性系统的模型辨识3.5.1 连续系统的模型辨识 用数值分析的手段，根据系统输入输出的实测数据，构造其数学模型的方法称为系统辨识。 对给定的离散频率采样点 若已知系统的频率响应为Pi，Qi， ，而连续系统传递函数为第5

26、0页，共65页，2022年，5月20日，15点10分，星期三50其中， （1）则可以用下面几个表达式定义中间数值 （2）第51页，共65页，2022年，5月20日，15点10分，星期三51引入拟合的性能指标 式中 e(jk)=G(jk)-G(jk)为频率拟合的误差， G(jk)为辨识出来的模型计算出的频率响应数据,D(jk)为加权系数。若对性能指标取导数 （3） 则获得性能指标J的最小值。由式（3）可推导出线性代数方程 （4）第52页，共65页，2022年，5月20日，15点10分，星期三52式中 = 可见，有了式（2）中定义的各个参数后，待辨识系统的传函可通过求解方程（4）获得。第53页，共

27、65页，2022年，5月20日，15点10分，星期三53MATLAB信号处理箱中，提供了辨识系统传函模型的invfreqs( )函数： num,den=invfreqs(H,w,nn,nd)其中，H为频域响应数据，w为离散频率点构成的向量，nn和nd分别为待辨识系统的分子和分母阶次。如果给出系统的幅频mag和相频响应数据phase，则可以由下面的方式来调用invfreqs( )函数 num,den=invfreqs(mag.*exp(sqrt(-1)*phase),w,nn,nd)例18：系统的频域响应已知(见教材p100表3-1)，且传函模型的分子、分母阶次分别为3和4，求其系统的传函。w=

28、0.01,0.0374,72.79,100;H=0.9999-0.01083i,0.9998-0.01488i,0.0002989-0.009981i;B,A=invfreqs(H,w,3,4);G1=tf(B,A)第54页，共65页，2022年，5月20日，15点10分，星期三543.5.2 离散系统的模型辨识离散系统的传递函数可写为其对应的差分方程为对测出的各组输入和响应数据，可写为下面的矩阵形式即 YM=XM(5)第55页，共65页，2022年，5月20日，15点10分，星期三55 这样得出的系数矩阵一般情况下有Mn，故该方程为超定方程，可用最小二乘法求解此外，MATLAB的系统辨识工具

29、箱还提供了各种系统辨识函数，ARX（Auto-Regressive exogenous, 自回归遍历）模型，可采用arx( )函数实现 T=arx(y, u,m, n, d)其中，y、u分别为输出、输入信号列向量；m、n分别为系统分子、分母多项式阶次；d为纯滞后；T为一个结构体，T.A和T.B分别表示辨识出的分母和分子多项式。(6)第56页，共65页，2022年，5月20日，15点10分，星期三56下面分别用两种方法辨识系统参数。例19：系统的实测输入/输出数据（教材p102表3-2），系统分子和分母的阶次分别为3、4。利用式(5)、(6)辨识系统参数y=0, 0.1374, 0.2978，.

30、, 3.3669, 3.4131;u=0.4398, 0.3400, 0.3142,., 0.5527, 0.4001, 0.1988;Phi=0;y(1:end-1) 0;0;y(1:end-2) ,. 0;0;0;y(1:end-3) 0;0;0;0;y(1:end-4),. 0;u(1:end-1) 0;0;u(1:end-2),. 0;0;0;u(1:end-3) 0;0;0;0;u(1:end-4);T=Phi y; T;Gd=tf(T(5:8),1, -T(1:4), Ts, 0.1)第57页，共65页，2022年，5月20日，15点10分，星期三57试采用arx( )函数辨识系统

31、参数先输入y,u列向量，然后直接调用arx( )函数y=0, 0.1374, 0.2978,3.3669, 3.4134;u=0.4398, 0.34,0.4001, 0.1988;t1=arx(y,u,4,4,1)Notice:要确定准确的采样信息，可采用iddata( )函数创建数据对象。然后调用arx( )函数。U=iddata(y, u, 0.1);T=arx(U,4, 4, 1);H=tf(T); G=H(1) %H为双输入传函矩阵，H(2)是由误差信号(k)到输入的传函。第58页，共65页，2022年，5月20日，15点10分，星期三58另：输入和输出数据的精确性，对于辨识系统的准

32、确性影响很大。y=0, 0.1374, 0.2978, 1.576, 1.743, 1.979, 2.33, 2.734, 2.965, 3.109, 3.402, 3.508, 3.673, 3.921, 4.099, 4.327, 4.296, 4.205, 4.244, 4.171, , 3.146, 3.288, 3.367, 3.413;u=0.4398, 0.34,0.4001,0.1988;Transfer function: 0.312 z3 - 0.4151 z2 + 0.1892 z - 0.01747 - z4 - 2.723 z3 + 2.64 z2 - 1.014 z

33、 + 0.1068输入、输入数据保留小数点后4位时，辨识系统的传函为Transfer function: 0.3126 z3 - 0.581 z2 + 0.3955 z - 0.09168 - z4 - 3.254 z3 + 4.04 z2 - 2.266 z + 0.4854第59页，共65页，2022年，5月20日，15点10分，星期三593.5.3 离散系统辨识信号的生成 idinput( )函数可生成指定的输入信号。如：u=idinput(N,rgs ) %生成N个点的Random, Gaussian Signal。 例20：生成31个点的PRBS信号u=idinput(31,PRBS

34、)； %生成31个PRBS(pseudo-random binary sequence)信号。t=0:.1:3; stairs(u), set (gca, Xlim, 0,31,Ylim,-1.1 1.1) %gca获得当前坐标轴句柄第60页，共65页，2022年，5月20日，15点10分，星期三60利用上述信号，可获得指定系统的输入和输入数据。 num=0.3124 -0.5743 0.3879 -0.0889;den=1 -3.233 3.9869 -2.2209 0.4723;G=tf(num,den,Ts,.1); t=0:.1:3; u=idinput(31,PRBS)； y=lsi

35、m(G,u,t); % LTI(linear time invariable 线性时不变)系统在给定输入下的时间响应。T1=arx(y,u,4 4 1)Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + e(t) A(q) = 1 - 3.233 q-1 + 3.987 q-2 - 2.221 q-3 + 0.4723 q-4B(q) = 0.3124 q-1 - 0.5743 q-2 + 0.3879 q-3 - 0.0889 q-4 Estimated using ARXLoss function 3.2187e-030 and FPE 5.45779e-030Sampling interval: 1第61页，共65页，2022年，5月20日，15点10分，星期三613.5.4 多变量离散系统的辨识 p路输入、q路输出的多变量系统的数学模型的差分方程为其中，A(z-1)和B(z