第4-1章-时域离散相似法Matlab编程课件.ppt

1、1第第4章章 时域离散相似法仿真时域离散相似法仿真Matlab编程编程4.1 连续模型的离散化连续模型的离散化)()()()()()(tDutCxtytButAxtx)()()()()()1(kTDukTCxkTykTHukTGxTkx线性连续状态空间模型为(4.1-1)其采用零阶保持器的离散相似模型为(4.1-2)其中TTAATBdeHeG0)(4.1-3)2用于离散相似的用于离散相似的Matlab函数函数连续模型转换为离散模型：G,H=c2d(A,B,T)T为采样周期可选参数离散相似：G,H,C,D=c2dm(A,B,C,D,T,选项)numd,dend=c2dm(num,den,T,选项

2、)3【例4.1.1】对下面的连续系统，在采样周期T=0.1时进行离散化232)(2sssGA=0,1;-2,-3;B=0;1;G,H=c2d(A,B,0.1);G=0.9909 0.0861 -0.1722 0.7326H=0.0045 0.0861经计算，可得TTTTTTTTeeeeeeeeTG22222222)(TTTTeeeeTH222121)(当T=0.1时，计算结果与c2d函数得到的结果一致。44.2 针对离散状态空间模型的仿真程序针对离散状态空间模型的仿真程序)()()()()()1(kTDukTCxkTykTHukTGxTkx离散状态空间模型为对此模型编写仿真程序，只需要按模型方

3、程迭代即可。function t,y=w_DiscreteSimu(tstart,tstop,h,x0,u0,G,H,C,D)%函数功能：对线性离散系统x(k+1)=G*X(k)+H*u(k),y(k)=C*x(k)+D*u(k)进行仿真%输入参数：tstart,tstop,h 分别是起始时间、结束时间和仿真步%长，是标量%x0,u0是状态、输入的初值，都是列向量%G,H,C,D是线性离散状态空间模型的系数矩阵%输出参数：t是仿真结果的时间序列%y是仿真结果系统的输出序列5function t,y=w_DiscreteSimu(tstart,tstop,h,x0,u0,G,H,C,D)t=ts

4、tart:h:tstop;%t数一个行序列cntt=size(t,2);%返回列数cnty=size(C,1);%返回y的维数y=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果y0=C*x0+D*u0;y(:,1)=y0 ;%将y0作为输出的第1列curx=x0;%当前一步的xcuru=u0;%当前一步的ucury=y0;%当前一步的yfor i=1:1:cntt-1 curx=G*curx+H*curu;%计算下一步的状态 cury=C*curx+D*curu;%计算输出 y(:,i+1)=cury;%将输出加入到输出序列里end离散状态空间模型仿真程序6【例4.2.1】对

5、下面的连续系统，采用数值积分法和离散相似法仿真。232)(2sssG解：该系统连续状态空间模型矩阵为TTTTTTTTeeeeeeeeTG22222222)(TTTTeeeeTH222121)(02,10,3210CBA其离散状态空间模型矩阵为也可以采用c2d函数直接计算离散模型矩阵7仿真程序:用c2d进行模型变换，将数值积分与离散相似2种方法的结果对比。数值积分仿真的结果用plot函数绘图，表现连续的效果；离散模型仿真结果用stairs绘图，各数据点之间是零阶保持。function simu4_2tstart=0;tstop=6;A=0,1;-2,-3;B=0;1;C=2,0;D=0;x0=0

6、;0;u0=1;h=0.5;G,H=c2d(A,B,h);%得到离散模型矩阵%数值积分仿真t,y1=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,RK4);%离散相似法仿真t,y2=w_DiscreteSimu(tstart,tstop,h,x0,u0,G,H,C,D);plot(t,y1,-r*);hold on;stairs(t,y2,:b+);hold off;axis(0,6,0,1.2);xlabel(x);ylabel(y);legend(数值积分,离散相似);title(数值积分与离散相似仿真对比)8在步长h=0.5时，两种方法的仿真结果及其

7、接近。离散模型的每个计算点都在连续模型上当步长h=1时，两种方法的结果出现偏差。离散模型的点不一定在连续模型上。012345600.20.40.60.81Timey数 值 积 分 与 离 散 相 似 仿 真 对 比(h=0.5)数 值 积 分离 散 相 似012345600.20.40.60.81Timey数 值 积 分 与 离 散 相 似 仿 真 对 比(h=1)数 值 积 分离 散 相 似94.3 面向环节离散化仿真面向环节离散化仿真(无非线性环节无非线性环节)【例4.3.1】将下面的系统按环节离散化仿真解：将系统分为2个典型环节，如图所示，开环部分传函为312)3(2)(sssssGK1

8、】环节之间的连接方程)()()()()(1221kykukykrku)(01)(0110)()()(0krkykrWkyWku102】环节内部离散模型（1）积分环节)()()(2)()1(11111kxkykTukxkx（2）惯性环节)()()(1)()1(222331232kxkykuekxekxTT)()1(002)(001)()()1(3313kueTkxekHukGxkxTT)(0000)(1001)()()(kukxkuDkxCky环节内部动态方程环节输出方程113】系统输出方程)()(2kykp)(10)()(kykyQkp面向环节离散化仿真时应该注意的问题面向环节离散化仿真时应该

9、注意的问题应该依次计算各环节的输出，即计算出第1个环节，再利用第1个环节的输出计算第2个环节，依次下去。而不要一次计算出各环节的输入，再算出各环节的输出，这样相当于在每个环节之间存在一个步长的延迟环节，模型的准确性会大大降低。下面，我们将对这两种方式做仿真对比，进而说明此问题。12function t,p=w_DisNodesSimu2(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q)%函数功能：面向环节离散化仿真，%环节输入方程：u(k)=W*x(k)+W0*r(k)%环节迭代方程：x(k+1)=G*x(k)+H*u(k)%环节输出方程 y(k)=C*x(k)+

10、D*u(k);%整个系统的输出方程 p(k)=Q*y(k);%输入参数：tstart,tstop,h 分别是起始时间、结束时间和仿真步长%x0,y0是每个环节状态和输出的初值，%r0为系统外部输入%输出参数：t是仿真结果的时间序列%p是仿真结果系统的输出序列1、整体刷新计算的按环节仿真函数、整体刷新计算的按环节仿真函数整体刷新即根据环节连接方程计算每个环节的输入，再用计算出来的输入向量根据各环节的离散模型计算各环节的输出。各环节之间相当于并行运算。13function t,p=w_DisNodesSimu2(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q)t=ts

11、tart:h:tstop;%t数一个行序列cntt=size(t,2);%返回列数cnty=size(Q,1);%返回y的维数p=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果p(:,1)=Q*y0;%for i=1:1:cntt-1 u=W*y0+W0*r0;x0=G*x0+H*u;y0=C*x0+D*u;p(:,i+1)=Q*y0 ;%将y0作为输出的第1列end142、逐个环节刷新的仿真程序、逐个环节刷新的仿真程序该函数所有参数与w_DisNodesSimu2的一样。逐个环节刷新是逐个根据环节的编号顺序，计算环节的输入，再计算输出，再计算下一个环节的输入和输出，依次进行

12、，相当于环节之间串连运算。function t,p=w_DisNodesSimu(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q)t=tstart:h:tstop;%t数一个行序列cntnodes=size(W,1);%得到环节的个数cntt=size(t,2);%返回列数cnty=size(Q,1);%返回y的维数p=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果p(:,1)=Q*y0;%for i=1:1:cntt-1 for j=1:1:cntnodes u(j)=W(j,:)*y0+W0(j,:)*r0;x0(j)=G(j,j)*x0(j

13、)+H(j,j)*u(j);y0(j)=C(j,j)*x0(j)+D(j,j)*u(j);end;p(:,i+1)=Q*y0 ;%将y0作为输出的第1列end15function simu4_3tstart=0;tstop=10;A=0,1;-2,-3;B=0;1;C=2,0;D=0;x0=0;0;u0=1;h=0.5;G,H=c2d(A,B,h);%得到离散模型矩阵t,y1=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D,RK4);t,y2=w_DiscreteSimu(tstart,tstop,h,x0,u0,G,H,C,D);%离散相似法仿真%面向环节

14、离散化仿真W=0,-1;1,0;W0=1,0;Q=0,1;C=diag(1,1);D=zeros(2);G=diag(1,exp(-3*h);H=diag(2*h,(1-exp(-3*h)/3);y0=x0;t,y3=w_DisNodesSimu2(tstart,tstop,h,x0,y0,u0,W,W0,G,H,C,D,Q);plot(t,y1,-r);hold on;stairs(t,y2,:k+);stairs(t,y3,-b*);hold off;axis(0,10,0,1.2);xlabel(x);ylabel(y);legend(数值积分,整体离散,按环节离散);title(各环节

15、整体刷新（h=0.5）);仿真调用函数仿真调用函数数值积分的结果用plot绘图，离散模型仿真结果用stairs绘图。16整体刷新时相当于前后环节之间存在一个步长的延迟，导致模型不精确。环节逐个刷新符合系统实际运算意义，模型比较精确01234567891000.20.40.60.81xy各 环 节 整 体 刷 新（h=0.5）数 值 积 分整 体 离 散按 环 节 离 散01234567891000.20.40.60.81xy各 环 节 整 体 刷 新（h=0.5）数 值 积 分整 体 离 散按 环 节 离 散174.4 面向环节离散化仿真（有非线性环节）面向环节离散化仿真（有非线性环节）环节划

16、分原则环节划分原则1】线性环节都处理为1阶；2】非线性环节附加在线性环节的前面或后面，环节个数等于线性环节的个数。一个环节的计算顺序1】根据环节之间的连接方程计算该环节的输入；2】根据前置非线性环节计算非线性环节的输出；3】前置非线性环节的输出作为线性环节的输入，计算线性环节的输出；4】线性环节的输出作为后置非线性环节的输入，计算整个环节的输出。18模型描述模型描述1】环节之间连接方程，写出各环节的输入与其他环节的输出和系统输入之间的连续方程。rWyWu02】线性环节离散模型)()()()()()1(kuDkxCkykuHkxGkxnnnnddDccChhHggG000000,00000000

17、0000,00000011113】其中19yQp4】非线性环节模型描述根据非线性环节类型不同，分别用不同的参数描述，包括非线性环节类型和必要的参数。非线性环节作为线性环节的前置或后置环节。5】整个系统的输出方程 仿真计算顺序仿真计算顺序根据环节编号顺序，计算一个环节，再计算下一个环节。不能整体计算每个环节的输入，再计算每个环节的输出，这样每个环节之间将产生滞后。20【例4.4.1】面向环节的非线性系统仿真右图所示系统有一个饱和非线性环节和一个线性环节。饱和非线性不能用微分方程表示，因而不能对整个系统采用数值积分法仿真。采用按环节离散化仿真，首先将线性环节分解为两个基本环节，将非线性环节作为第1

18、个线性环节的前置环节。仿真中，取c=20，K1，r=25211】各线性环节内部离散方程（1）积分环节（3）惯性环节)()()()()1(11111kxkykTuKkxkx)()()(1)()1(22222kxkykuekxekxTT0000,1001,)1(00,001DCeTKHeGTT)()()()()()1(kDukCxkykHukGxkx用矩阵表示为222】各环节之间连接方程)()()()()(1221kykukykrku3】系统输出为)()(2kykp)(01)(0110)()()(0krkykrWkyWku)(10)()(kykyQkp4】饱和非线性环节作为第1个环节的前置环节ck

19、ucckukuckuckuo)()()()()(111123function t,p=w_DisNodesSimuNonL(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q,NLBefore,NLBeforeParams,NLAfter,NLAfterParams)%函数功能：带有非线性环节的面向环节离散化仿真，%环节输入方程：u(k)=W*x(k)+W0*r(k)%环节迭代方程：x(k+1)=G*x(k)+H*u(k)%环节输出方程 y(k)=C*x(k)+D*u(k);%整个系统的输出方程 p(k)=Q*y(k);%输入参数：tstart,tstop,h 分别是

20、起始时间、结束时间和仿真步长，是标量%x0,y0是每个环节状态和输出的初值，%r0为系统外部输入%NLBefore是前置非线性环节的名称，如果没有就写None,NLBeforeParams是前置非线性环节的参数，每个环节%的参数占一行。%NLAfter是后置非线性环节的名称，如果没有就写None,NLAfterParams是后置非线性环节的参数，每个环节%的参数占一行。%输出参数：t是仿真结果的时间序列%p是仿真结果系统的输出序列1、带非线性环节的按环节离散化仿真程序原型、带非线性环节的按环节离散化仿真程序原型24function t,p=w_DisNodesSimuNonL(tstart,t

21、stop,h,x0,y0,r0,W,W0,G,H,C,D,Q,NLBefore,NLBeforeParams,NLAfter,NLAfterParams)t=tstart:h:tstop;%t数一个行序列cntnodes=size(W,1);%得到环节的个数cntt=size(t,2);%返回列数cnty=size(Q,1);%返回y的维数p=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果p(:,1)=Q*y0;for i=1:1:cntt-1 for j=1:1:cntnodes u(j)=W(j,:)*y0+W0(j,:)*r0;%根据环节之间的连接计算，u(j)=Non

22、linearNode(NLBefore(j),u(j),NLBeforeParams(j,:);%计算第j个环节的前置非线性环节的输出 x0(j)=G(j,j)*x0(j)+H(j,j)*u(j);y0(j)=C(j,j)*x0(j)+D(j,j)*u(j);%根据线性部分计算线性部分的输出 y0(j)=NonlinearNode(NLAfter(j),y0(j),NLAfterParams(j,:);%计算第j个环节的后置非线性环节的输出 end;p(:,i+1)=Q*y0 ;%将y0作为输出的第1列end252、典型非线性环节函数、典型非线性环节函数function u_out=Nonli

23、nearNode(NodeType,u_in,params)%函数功能：典型非线性环节的计算%输入参数：NodeType 节点类型，可以取Saturation,Deadzone,Relay,RelayDeadzone,%None表示没有非线性环节%u_in,环节的输入；%params,环节的参数行向量%输出结果：u_out，环节的输出26switch NodeType case Saturation c=params(1);k=params(2);if(k*u_in)c u_out=c;else u_out=k*u_in;end;case Deadzone c=params(1);k=para

24、ms(2);if(u_inc)u_out=k*u_in;else u_out=0;end;case Relay c=params(1);if u_in0 u_out=c;else u_out=0;end;case RelayDeadzone c=params(1);h=params(2);if u_inc u_out=c;else u_out=0;end;otherwise u_out=u_in;end;函数函数NonlinearNode的主体的主体27function simu_NonNodeststart=0;tstop=20;h=0.5;K=1;c=10;W=0,-1;1,0;W0=1,

25、0;Q=0,1;C=diag(1,1);D=zeros(2);G=diag(1,exp(-h);H=diag(K*h,1-exp(-h);x0=0;0;u0=25;y0=x0;NLBefore=Saturation,None;NLBeforeParams=c,1;0,0;NLAfter=None,None;NLAfterParams=0,0;0,0;t,y=w_DisNodesSimuNonL(tstart,tstop,h,x0,y0,u0,W,W0,G,H,C,D,Q,NLBefore,NLBeforeParams,NLAfter,NLAfterParams);stairs(t,y,-r*)

26、;xlabel(Time);ylabel(y);title(含非线性环节的面向环节仿真);3、仿真调用函数仿真调用函数2802468101214161820051015202530Timey含 非 线 性 环 节 的 面 向 环 节 仿 真阶 梯 状连 续 状仿真结果：同样的结果，用plot与stairs绘图的对比。仿真中，取c=20，K1，r=25294.5 利用利用Simulink进行面向环节仿真进行面向环节仿真【例4.5.1】下图所示线性与非线性混合系统，当回环非线性特性c=1时，系统有自持振荡，振幅约为10.89.4，振荡周期约为6s.对该系统采用Simulink建模，进行仿真30Si

27、mulink将此系统看做有6个连续环节。将仿真输出用Scope绘图，同时将结果数据输出到Workspace里，以便用plot绘图。051015202530354000.20.40.60.811.21.41.6Timeoutput从曲线，可以很明显地看到输出震荡。31【例4.5.2】用Simulik对下图所示的系统进行仿真仿真中，取c=10，K1，r=25，并且在t=1时才开始阶跃输入。3202468101214161820-5051015202530系 统 输 出饱 和 环 节 的 输 出仿真结果曲线：很明显，饱和环节的输出受限制。33【例4.5.3】对【例4.3.1】所示的系统采用Simulink仿真（1）积分环节)()()(2)()1(11111kxkykTukxkx（2）惯性环节)()()(1)()1(222331232kxkykuekxekxTT当T=0.5时)()()()()1(11111kxkykukxkx)()()(2590.0)(2231.0)1(22222kxkykukxkx34本例采用离散化环节进行建模