控制系统仿真 第四章Matlab仿真实现1学习资料

自动控制系统仿真Thesimulationofautomaticcontrolsystem第四章控制系统数字仿真的实现第四章控制系统数字仿真的实现计算机控制系统的数字仿真龙格库塔法求传递函数系统时域解控制系统的结构及其拓扑描述环节的离散化与非线性系统的数字仿真龙格库塔法求状态方程时域解面向系统结构的数字仿真第8次1.使用Matlab函数：ode451）编写xfun.m文件，定义xfun函数()2）在命令窗口调用ode45函数2.编写龙格库塔算法的**.m文件龙格库塔法求状态方程时域解[t,x]=ode45(‘xfun’,[t0,tf],x0)[t,x]=ode45(‘xfun’,[t0,tf],x0,option,u)functionxdot=xfun(t,x,flag,u)或)()()，(

率

斜

四

第)2()2()2，2(

第三斜率)2()2，2(

第二斜率)()(

率

斜

第一3342231121hthhhththhhthhhtttkkkkkkk+++=++=+++=++=++=++=+==KXAKXfKKXAKXfKKXAKXfKBrXAX,fKkkkkkkkk：：：：+(tkBr)2hBrBrXk+1=Xk+h/6(K1+2K2+2K3+K4)编写龙格库塔算法的rk4.m文件%round(x):对x取最接近的整数，h为步长%求各次斜率K1~K4%求Xk+1=Xk+h/6(K1+2K2+2K3+K4)%yk+1=CXk+1%A,B,C为系统状态方程系数阵%x0为x的初值，t0为初始时刻，tf为终止时刻第四章控制系统数字仿真的实现计算机控制系统的数字仿真龙格库塔法求传递函数系统时域解控制系统的结构及其拓扑描述环节的离散化与非线性系统的数字仿真龙格库塔法求状态方程时域解面向系统结构的数字仿真龙格库塔法求传递函数系统时域解传递函数：化为能控标准型(1)(2)(3)龙格库塔法求传递函数系统时域解0.化传递函数为状态方程形式1.使用Matlab函数：ode45(‘xfun’,[t0,tf],x0)

或2.编写龙格库塔算法的**.m文件%传递函数分子系数向量%传递函数分母系数向量%化为状态方程标准型实现%fliplr(X):矩阵X左右翻转%zeros(m,n):产生m×n的零阵%rot90(A):A阵逆时针转90°%length(x)：求向量x的维数%eye(m,n):产生m×n的单位阵第四章控制系统数字仿真的实现计算机控制系统的数字仿真龙格库塔法求传递函数系统时域解控制系统的结构及其拓扑描述环节的离散化与非线性系统的数字仿真龙格库塔法求状态方程时域解面向系统结构的数字仿真第8次控制系统的结构及其拓扑描述1.方框图化简2.想办法表达环节之间的输入输出关系y1y2y3y4y5y6y8y9y10u1u2u3u4u5u6u7y7u8u10u9麻烦连接矩阵写出每个环节输入ui

受哪些环节输出yi的制约和影响。控制系统的连接矩阵：W

称为连接矩阵W0

称为输入连接矩阵u=[u1

，u2

，…,un]’y=[y1，y2，…,yn

]’连接矩阵W

中元素值直接反映各环节之间连接情况：

wij=0,环节j不与环节i相连；

wij≠0,环节j与环节i有连接关系；

wij>0,环节j与环节i直接相连(wij=1)

或通过比例系数相连(wij为任意正实数)；

wij<0,环节j与环节i直接负反馈相连(wij=–1)

或通过比例系数负反馈相连(wij为任意负实数)；wii

≠0,环节i单位自反馈(wii=1或wii=–1)

或通过比例系数自反馈(wii

为任意实数)；控制系统的连接矩阵：习题4-5写出连接矩阵y1y2y3y4y5y6y8y9y10u1u2u3u4u5u6u7y7u8u10u9

U=WY+W0y0

第四章控制系统数字仿真的实现计算机控制系统的数字仿真龙格库塔法求传递函数系统时域解控制系统的结构及其拓扑描述环节的离散化与非线性系统的数字仿真龙格库塔法求状态方程时域解面向系统结构的数字仿真面向系统结构的数字仿真典型闭环系统的数字仿真：控制系统最常见的典型闭环系统结构

系统的开环传递函数G(s)，可按照能控标准型写出其开环状态方程：G(S)的控制量

又因为

其中：面向系统结构的数字仿真1.使用Matlab函数：ode45(‘xfun’,[t0,tf],x0)

或2.编写龙格库塔算法的**.m文件化为了标准形式：G(S)的状态方程：

面向系统结构的数字仿真例4-1．求图4-8所示系统的阶跃响应y(t)数值解：

解：①设开环放大系数k=1，反馈系数v=1(单位反馈系统)，阶跃输入幅值r=1；②利用conv()卷积函数功能，将系统开环传递函数G(s)化为传递函数分母、分子多项式系数向量的形式：a=conv([100],conv([0.251],[0.251]))

③设系统状态向量初值X0=[x10,x20,…,xn0]均为零；④系统运行参数,t0=0,tf=10,h=0.25，系统阶次为4阶；⑤求G(s)的能控标准型实现A、b、C，求反馈闭环系数矩阵Ab;

⑥调用龙格库塔法面向系统结构的数字仿真2面向系统结构的数字仿真例4-1．求图所示系统的阶跃响应y(t)数值解：

例4-1>>a=[0.06250.5100];>>b=[21];%k=1>>x0=[0000];>>v=1;n=4;t0=0;tf=10;h=0.25;r=1;>>rk4(a,b,x0,v,t0,tf,h,r,n)%Ab为闭环系数阵例4-1复杂连接的闭环系统数字仿真：思路：在计算机程序中应能方便表示出实际系统各典型环节；

将环节之间的连接关系输入计算机，由计算机程序自动形成闭环状态方程；3.运用龙格库塔数值积分方法求解时域响应。面向系统结构的数字仿真比例环节

惯性环节惯性比例环节积分环节积分比例环节二阶振荡环节高阶线性环节复杂连接的闭环系统数字仿真：常见环节完全可用一个通用一阶环节统一表示

uiyi模型参数阵：

于是系统中所有环节输出、输入关系统一用矩阵表示如下：

(A+Bs)Y=(C+Ds)U

面向系统结构的数字仿真复杂连接的闭环系统数字仿真：1.写成统一各环节数学模型复杂连接的闭环系统数字仿真：面向系统结构的数字仿真整理为矩阵形式：

U=WY+W0y0

其中：2.求连接矩阵各环节输入ui与输出yi有以下关系：

(A+Bs)Y=(C+Ds)(WY+W0

y0)

整理，得：

(B–DW)sY

=(CW–A)Y+CW0

y0

+DW0

sy0

简洁表达为：

QsY

=RY+V1y0

+V2

sy0

其中：Q=B–DW

R=CW–A

V1

=CW0

V2

=DW0

若Q阵逆存在，则式(4-9)两边同时左乘Q

–1

，得：

sY=Q–1RY+Q–1

V1

y0

+Q–1

V2

sy0

两边反拉氏变换，求得系统闭环状态方程时域表达式：

Ab

=Q

–1

R;b1

=Q–1V

1

；b2=Q–1V2

为闭环系统的系数阵和输入阵面向系统结构的数字仿真(A+Bs)Y=(C+Ds)U将U=WY+W0y0代入3.求系统闭环状态方程建立该系统仿真模型中应注意两点：①保证Q阵有逆。

②

去掉项，使V2=0面向系统结构的数字仿真Ab

=Q

–1

R;b1

=Q–1V

1

；b2=Q–1V2Q=B–DWR=CW–AV1

=CW0V2

=DW0步长值的选择：为系统开环频率特性的剪切频率或

tr为系统阶跃响应的上升时间

ts为系统阶跃响应的调节时间（过渡过程时间）

h0的选取应小于系统中最小时间常数τ的两倍，即：h0＜2τ

面向系统结构的数字仿真仿真程序流程图(1)系统参数输入方法面向系统结构的数字仿真输入各环节参数阵PA=diag(P(:,1));%取第一列构成对角阵B=diag(P(:,2));C=diag(P(:,3));D=diag(P(:,4));P(:,1)=[a1a2…an]’(2)联接矩阵输入方法面向系统结构的数字仿真有大量0元素，因此只输入非零元素WIJ=[ijwij]m×3i为被作用环节号；j为作用环节号；wij为关系值；m为非0元素个数。输入对环节1的作用环节3对环节1反馈作用环节1对环节2直接作用环节2的自反馈环节2对环节3的直接作用面向系统结构的数字仿真(2)联接矩阵输入方法※例4-2%建立初始W、W0阵调用4阶龙格库塔法%调用4阶龙格库塔法Y习题4-5y1y2y3y4y5y6y8y9y10u1u2u3u4u5u6u7y7u8u10u9习题4-5写出连接矩阵y1y2y3y4y5y6y8y9y10u1u2u3u4u5u6u7y7u8u10u9

U=WY+W0y0

第四章控制系统数字仿真的实现计算机控制系统的数字仿真龙格库塔法求传递函数系统时域解控制系统的结构及其拓扑描述环节的离散化与非线性系统的数字仿真龙格库塔法求状态方程时域解面向系统结构的数字仿真第10次数值积分方法仿真的不足1、单步法计算各次斜率（K1~K4），步骤繁琐2、多步法要求存储前r次时刻数据3、隐式算法求解需若干次迭代，计算速度受影响4、对状态方程求解，不能对非线性环节进行仿真求解差分方程——离散（采样）控制系统数学模型连续系统的差分方程连续系统的离散化非线性系统的仿真另辟蹊径连续系统的离散化

其中为状态初始值.则其中：Φ(t)为状态转移矩阵当状态方程为线性定常时，或：

Φ(t)=

eAt

=L-1{(sI–A)–1}故当t=kT时，求解差分方程法

t=(k+1)T时，所以：X((k+1)T)=Φ(T)X(kT)+Φm

(T)U

(kT

)

其中：Φ(T)=eAT,为t=T时的状态转移矩阵数值求解递推公式：Xk+1=Φ(T)Xk+Φm

(T)U

k

连续系统的离散化若希望递推公式精度更高些，应该考虑到在两次采样时刻kT

、(k+1)T之间u(τ)一直在变化，用一阶保持器近似更为合理，如图将u(τ)表为随Δu

k

(τ)变化的函数：

u(τ)=u(kT

)+Δu

k(τ)而Δu

k(τ)又可用下式近似表达：于是代入式(4-11)中积分项重新推导其中：

连续系统的离散化

利用式(4-14)就可以编程进行对某连续系统的仿真运算了。事先应离线求取Φ(t)=eAt再令阵中t=T,立即得离散化矩阵Φ(T)、Φm

(T)、Φj

(T)连续系统的离散化Matlab控制工具箱——c2d1.c2d:变连续时间系统为离散时间系统完成将状态空间模型从连续时间到离散时间的转换，并假定系统输入采用零阶保持器按系统离散化，不能解决非线性环节Matlab控制工具箱——c2dm2.c2dm:变连续状态空间模型为离散空间模型可实现状态空间模型从连续时间域变换成离散时间域‘zoh’：零阶保持器‘foh’：一阶保持器‘tustin’采用双线性逼近导数‘prewarp’利用频率预变的双线性逼近‘matched’利用匹配零-极点方法将SISO系统变为离散时间系统缺省为‘zoh’典型环节状态方程的离散化其中：

典型环节状态方程的离散化①积分环节离散化：

A=0，B=1，C=K，D=0由于A，B，C均与①相同，故Φ(T)、Φm

(T)、Φj

(T)和Φc

与①完全相同，相应状态方程也完全相同。但因D≠0，只有Φd

不同，所以应注意输出方程成为：

yk+1

=xk+1+Kbuk+1

②积分比例环节离散化：

A=0，B=1，C=K，D=bK

③惯性环节离散化：

A=a，B=1，C=K，D=0④比例惯性环节离散化：

A=a，B=1，C=K，D=离散化环节参数表K/sK(1+bs)/sK/(s+a)K(s+b)/(s+a)ФФmФj1KTKT2/2e-aT(K/a)(1-e-aT)(K/a2)[aT-(1-e-aT)]Фc11b-aФd0Kb0KKabC/B00C/B0D/CC/BA/B0D/BA/BC/D典型环节状态方程的离散化连续系统的离散化模型：按环节离散化数字仿真程序与实现典型环节数据输入后，首先判断A是否为0，即可分出①、②和③、④两组，这两组对应的状态方程离散化系数Φ(T)、Φm

(T)和Φj

(T)求取方法各自相同，可以直接套用相同求解公式求取后存入相应单元。但由于对应输出方程各有不同，故又需判断D是否为0，从而对输出方程离散化系数Φc

、Φd加以修正后，也存入相应单元。各环节离散化系数求得后，结果存入相应数组单元FI(I)、FIM(I)、FIJ(I)、FIC(I)以及FID(I)，其中：I表示环节序号。仿真运行时从各环节相应单元取出，分别求取各环节状态与输出即可。环节的离散化与非线性系统的数字仿真连续系统的离散化模型：求

t=(k+1)T时刻的各环节状态

Xk+1

的递推计算式中要用到

u

k

、而求

t=(k+1)T时刻的各环节输出Yk+1

的递推计算式中还要用到

u

k+1

①u

k

可通过联接矩阵直接求得，即：Uk

=WYk+W0y0②利用近似表达式求取③

u

k+1利用上面已求得的Uk

、在一个步长h内按一阶保持近似关系求取。环节的离散化与非线性系统的数字仿真离散法仿真程序实现离散法仿真程序实现%积分和积分比例环节情况%积分比例环节情况%惯性和比例惯性环节情况%比例惯性环节情况离散法仿真程序实现%记录前一时刻的Uk值，Ub=Uk-1%求Uk导数%Uk为各环节的输入%求Uf=Uk+1饱和非线性非线性系统的数字仿真2.死区非线性非线性系统的数字仿真3.滞环非线性非线性系统的数字仿真%Ucb表示前一次输出值%Urb表示前一次输出值4.继电非线性非线性系统的数字仿真非线性特征的判断利用按环节离散化的仿真程序，在输入数据时，设立非线性标志向量：

Z=[z1,z2

,…,zn

]

Z(i)=0线性环节

Z(i)=1线性环节前有饱和非线性，应修正U(i)

Z(i)=2线性环节前有死区非线性，应修正U(i)

Z(i)=3线性环节前有滞环非线性，应修正U(i)

Z(i)=4线性环节后有饱和非线性，应修正Y(i)

Z(i)=5线性环节后有死区非线性，应修正Y(i)

Z(i)=6线性环节后有滞环非线性，应修正Y(i)非线性系统的数字仿真程序功能扩展

1.求得uk

后，由Z(i)判断各环节入口有否非线性。若有，则根据标志值确定类型，转相应处理程序，修正uk

值；2.求得yk+1后，再由Z(i)判断各环节出口有否非线性。若有，则根据标志值确定类型，转相应处理程序，修正yk+1值；3.各种非线性特性按前节给出的程序，自定义为函数形式，以函数文件格式存储起来，由主程序在运行时调用。非线性系统的数字仿真例4–2．控制系统如图4–29所示设输入阶跃函数幅值Y0=10,滞环非线性参数s4=1(滞环宽度)⑴不考虑非线性环节影响时，求解y(t)的阶跃响应；⑵考虑非线性环节影响，其余参数不变，求解y(t)并与线性情况所得结果比较。解⑴①先将环节编号标入图4-29中；②MATLAB命令窗口下(以下语句前符号“》”即表示MATLAB命令窗口环境)，按编号依次将环节参数输入P阵非线性系统的数字仿真③按各环节相对位置和联接关系，有联接矩阵如下：④输入各环节初始值图4–30例4–2控制系统输出y4响应曲线⑤由于不考虑非线性影响，则非线性标志向量和参数向量均应赋零值⑥输入运行参数；开环截止频率ωc

约为1，故计算步距h取经验公式值⑦运行sp4_4.m求解环节Ⅳ输出y4数值解数据和响应曲线非线性系统的数字仿真非线性仿真程序实现Z=[0000],S=[0000]%线性环节后有滞环非线性非线性仿真程序实现第四章控制系统数字仿真的实现计算机控制系统的数字仿真龙格库塔法求传递函数系统时域解控制系统的结构及其拓扑描述环节的离散化与非线性系统的数字仿真龙格库塔法求状态方程时域解面向系统结构的数字仿真第12次典型计算机控制系统结构图采样控制系统结构图计算机控制系统的数字仿真采样开关采样开关+采样保持器连续信号为x(t)，拉氏变换为：X(s)=L{x(t)}对其拉氏变换，得：令：z=eTs，则定义x(t)

的Z变换：

记作：X(z)=Z{x(t)}采样控制系统的数学描述对采样信号的描述——Z变换：x(t)

经采样后，由下式描述：设控制对象传递函数为G(s)，其脉冲响应函数为g(t)，则扩展Z变换：G(z,Δ)=Z{g(

kT+ΔT)}=Z{G(s)

eΔTs

}

0<Δ<1由上式继续推导可得：

G(z,m)=Z{g(k+m)T–T)}=Z{g(kT–(1–m)T)}

=Z{g(kT–λT

)}=Z{g(t)e–λTs

}

0<λ<1扩展Z变换的求取方法与普通Z变换相同，典型函数的扩展Z变换通过查表方式求取更为方便。一般情况可用极点留数法：采样控制系统的数学描述对采样时刻之间信号变化的描述——扩展Z变换法：差分方程的一般形式表述如下：

y(k)+a1y(k-1)+…+any(k-n)=b0u(k)+b1u(k-1)+…+bmu(k-m)(n≥m)高阶差分方程也表为状态方程形式。单输入-单输出情况下，状态方程一般形式为：两边求零初始条件下的Z变换，可得：

Y(z)+a1

z-1

Y(z)+…+anz-nY(z)=b0

U(z)+b1

z-1

U(z)+…+bmz-mU(z)

(n≥

m)

采样控制系统的数学描述对离散(采样)信号相互作用的描述——差分方程：采样系统可以表示为：Z传递函数分子、分母系数向量形式：

(numd

,dend

)=([b0,b1,…,bm],[1,a1,…,an])零极点增益向量形式：

(Zd,Pd,Kd

)=([zd1,…,zdm],[pd1,…,pdn],Kd)部分分式向量形式：

(Rd,Pd,hd

)=([rd1,…,rdn],[pd1,…,pdn],hd)还有离散状态方程各系数矩阵形式：

(Ad,Bd,Cd,Dd

)几种形式之间均可利用MATLAB语言控制系统工具箱中的数学模型转换函数tf2ss()、tf2zp()、ss2tf()、ss2zp()、zp2tf()、zp2ss()以及residue()等作相互转换采样控制系统的数学描述采样系统数学描述及其相互转换：差分方程递推求解法：已知：

则：U(z)=D(z)E(z)当

D(z)形如：

则：(1+c1

z-1+…+clz-l

)U(z)=(d0+d1z-1+…+drz-r)E(z)Z反变换，并整理得到递推式：

uk

=-c1

uk-1

-…-cluk-l+d0ek

+d1

ek-1+…+drek-r

相当于一个多步法递推算式，只要uk

和ek

的前若干步值已知，就可以递推得到

uk

。采样控制系统的仿真方法已知：

容易求出：

则：yk

=-a

1

yk-1

-…-anyk-n+

b0uk

+b1

uk-1+…+bmuk-m

再考虑：ek=rk

-yk就能按信号传递过程，从参考输入rk

开始，逐步求得各部分解ek、uk和输出yk。采样控制系统的仿真方法连续部分按环节离散化方法：当系统连续部分较复杂时，不必去化简和求取G(z)，而按照连续系统环节离散化仿真方法，将连续部分中各环节离散化处理后，与采样部分一并考虑进行仿真。连续部分各环节之间虚设采样开关和保持器，按环节离散化方法建立模型，应取数量级小于采样周期T较多的仿真步距h，才能很好地反映出连续系统在离散信号每隔周期T作用下，各环节在T内的细微变化采样控制系统的仿真方法控制器设计为连续系统环节D(s)的仿真方法：采样系统的控制器有时设计为连续系统环节形式，其传递函数为D(s)。但其发生作用又是每相隔一个采样周期T才经采样开关传递给控制对象连续部分。在T时刻内，D(s)、G0

(s)同时都在按照自己规律连续变化对这类系统仿真，要顾及到两部分连续系统各环节D(s)、G0

(s)各自的实际变化过程，对两部分都应正确地通过仿真得到准确的结果，并在每采样周期T到来时，将变化结果及时传递到相应环节去采样控制系统的仿真方法高阶差分方程的仿真程序实现：若采样系统直接给出输入—输出闭环Z传递函数GB(z)形式，即：

相应差分方程：

y(k)+a1y(k

-1)+…+any(k

-n)=b0u(k)+b1u(k-1)+…+bmu(k-m)需求解一组高阶差分方程：可见，以上算法最终都归结为求解高阶差分方程问题。在计算机上实现高阶差分方程的求解程序，应当注意以下问题：1.建立向量存储单元,保存和记忆输入uk

、输出yk

前若干时刻的值uk-1,…,uk-m

和yk-1

,…,yk-n

。2.每运算一个时刻值后，要及时刷新和摒弃相应的存储单元内容。即：只保留uk

及其前m个时刻的数值，yk

及其前n个时刻的数值，因此，要安排相应的平移操作程序采样控制系统仿真程序实现采样控制系统仿真程序实现纯滞后环节的处理方法及仿真程序实现：采样控制系统中常见一些控制对象包含有滞后(延迟)环节。其数学模型：即：y(t)=u(t-τ)将y(t)离散化成为y(kT

)，并将滞后时间常数τ表为T的函数，则：

τ=(M1

+M

2

)T其中：M1为整数，M