第4章時域離散相似法仿真Matlab編程4.1連續(xù)模型的離散化線性連續(xù)狀態(tài)空間模型為(4.1-1)其采用零階保持器的離散相似模型為(4.1-2)其中(4.1-3)1第4章時域離散相似法仿真Matlab編程4.1連續(xù)模型的用于離散相似的Matlab函數(shù)連續(xù)模型轉(zhuǎn)換為離散模型：[G,H]=c2d(A,B,T)T為采樣周期可選參數(shù)離散相似：[G,H,C,D]=c2dm(A,B,C,D,T,’選項’)[numd,dend]=c2dm(num,den,T,’選項’)2用于離散相似的Matlab函數(shù)連續(xù)模型轉(zhuǎn)換為離散模型：[G,【例4.1.1】對下面的連續(xù)系統(tǒng)，在采樣周期T=0.1時進(jìn)行離散化>>A=[0,1;-2,-3];>>B=[0;1];>>[G,H]=c2d(A,B,0.1);G=0.99090.0861-0.17220.7326H=0.00450.0861經(jīng)計算，可得當(dāng)T=0.1時，計算結(jié)果與c2d函數(shù)得到的結(jié)果一致。3【例4.1.1】對下面的連續(xù)系統(tǒng)，在采樣周期T=0.1時進(jìn)行4.2針對離散狀態(tài)空間模型的仿真程序離散狀態(tài)空間模型為對此模型編寫仿真程序，只需要按模型方程迭代即可。function[t,y]=w_DiscreteSimu(tstart,tstop,h,x0,u0,G,H,C,D)%函數(shù)功能：對線性離散系統(tǒng)x(k+1)=G*X(k)+H*u(k),y(k)=C*x(k)+D*u(k)進(jìn)行仿真%輸入?yún)?shù)：tstart,tstop,h分別是起始時間、結(jié)束時間和仿真步%長，是標(biāo)量%x0,u0是狀態(tài)、輸入的初值，都是列向量%G,H,C,D是線性離散狀態(tài)空間模型的系數(shù)矩陣%輸出參數(shù)：t是仿真結(jié)果的時間序列%y是仿真結(jié)果系統(tǒng)的輸出序列44.2針對離散狀態(tài)空間模型的仿真程序離散狀態(tài)空間模型為對此function[t,y]=w_DiscreteSimu(tstart,tstop,h,x0,u0,G,H,C,D)t=[tstart:h:tstop];%t數(shù)一個行序列cntt=size(t,2);%返回列數(shù)cnty=size(C,1);%返回y的維數(shù)y=zeros(cnty,cntt);%構(gòu)造一個空矩陣，用來存儲結(jié)果y0=C*x0+D*u0;y(:,1)=y0;%將y0作為輸出的第1列curx=x0;%當(dāng)前一步的xcuru=u0;%當(dāng)前一步的ucury=y0;%當(dāng)前一步的yfori=1:1:cntt-1curx=G*curx+H*curu;%計算下一步的狀態(tài)cury=C*curx+D*curu;%計算輸出y(:,i+1)=cury;%將輸出加入到輸出序列里end離散狀態(tài)空間模型仿真程序5function[t,y]=w_DiscreteSimu(【例4.2.1】對下面的連續(xù)系統(tǒng)，采用數(shù)值積分法和離散相似法仿真。解：該系統(tǒng)連續(xù)狀態(tài)空間模型矩陣為其離散狀態(tài)空間模型矩陣為也可以采用c2d函數(shù)直接計算離散模型矩陣6【例4.2.1】對下面的連續(xù)系統(tǒng)，采用數(shù)值積分法和離散相似法仿真程序:用c2d進(jìn)行模型變換，將數(shù)值積分與離散相似2種方法的結(jié)果對比。數(shù)值積分仿真的結(jié)果用plot函數(shù)繪圖，表現(xiàn)連續(xù)的效果；離散模型仿真結(jié)果用stairs繪圖，各數(shù)據(jù)點之間是零階保持。functionsimu4_2tstart=0; tstop=6; 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);%得到離散模型矩陣%數(shù)值積分仿真[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*');holdon;stairs(t,y2,':b+');holdoff;axis([0,6,0,1.2]);xlabel('x');ylabel('y');legend('數(shù)值積分','離散相似');title('數(shù)值積分與離散相似仿真對比')7仿真程序:用c2d進(jìn)行模型變換，將數(shù)值積分與離散相似2種方法在步長h=0.5時，兩種方法的仿真結(jié)果及其接近。離散模型的每個計算點都在連續(xù)模型上當(dāng)步長h=1時，兩種方法的結(jié)果出現(xiàn)偏差。離散模型的點不一定在連續(xù)模型上。8在步長h=0.5時，兩種方法的仿真結(jié)果及其接近。離散模型的每4.3面向環(huán)節(jié)離散化仿真(無非線性環(huán)節(jié))【例4.3.1】將下面的系統(tǒng)按環(huán)節(jié)離散化仿真解：將系統(tǒng)分為2個典型環(huán)節(jié)，如圖所示，開環(huán)部分傳函為1】環(huán)節(jié)之間的連接方程94.3面向環(huán)節(jié)離散化仿真(無非線性環(huán)節(jié))【例4.3.1】將2】環(huán)節(jié)內(nèi)部離散模型（1）積分環(huán)節(jié)（2）慣性環(huán)節(jié)環(huán)節(jié)內(nèi)部動態(tài)方程環(huán)節(jié)輸出方程102】環(huán)節(jié)內(nèi)部離散模型（1）積分環(huán)節(jié)（2）慣性環(huán)節(jié)環(huán)節(jié)內(nèi)部動態(tài)3】系統(tǒng)輸出方程面向環(huán)節(jié)離散化仿真時應(yīng)該注意的問題應(yīng)該依次計算各環(huán)節(jié)的輸出，即計算出第1個環(huán)節(jié)，再利用第1個環(huán)節(jié)的輸出計算第2個環(huán)節(jié)，依次下去。而不要一次計算出各環(huán)節(jié)的輸入，再算出各環(huán)節(jié)的輸出，這樣相當(dāng)于在每個環(huán)節(jié)之間存在一個步長的延遲環(huán)節(jié)，模型的準(zhǔn)確性會大大降低。下面，我們將對這兩種方式做仿真對比，進(jìn)而說明此問題。113】系統(tǒng)輸出方程面向環(huán)節(jié)離散化仿真時應(yīng)該注意的問題11function[t,p]=w_DisNodesSimu2(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q)%函數(shù)功能：面向環(huán)節(jié)離散化仿真，%環(huán)節(jié)輸入方程：u(k)=W*x(k)+W0*r(k)%環(huán)節(jié)迭代方程：x(k+1)=G*x(k)+H*u(k)%環(huán)節(jié)輸出方程y(k)=C*x(k)+D*u(k);%整個系統(tǒng)的輸出方程p(k)=Q*y(k);%輸入?yún)?shù)：tstart,tstop,h分別是起始時間、結(jié)束時間和仿真步長%x0,y0是每個環(huán)節(jié)狀態(tài)和輸出的初值，%r0為系統(tǒng)外部輸入%輸出參數(shù)：t是仿真結(jié)果的時間序列%p是仿真結(jié)果系統(tǒng)的輸出序列1、整體刷新計算的按環(huán)節(jié)仿真函數(shù)整體刷新即根據(jù)環(huán)節(jié)連接方程計算每個環(huán)節(jié)的輸入，再用計算出來的輸入向量根據(jù)各環(huán)節(jié)的離散模型計算各環(huán)節(jié)的輸出。各環(huán)節(jié)之間相當(dāng)于并行運(yùn)算。12function[t,p]=w_DisNodesSimu2function[t,p]=w_DisNodesSimu2(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q)t=[tstart:h:tstop];%t數(shù)一個行序列cntt=size(t,2);%返回列數(shù)cnty=size(Q,1);%返回y的維數(shù)p=zeros(cnty,cntt);%構(gòu)造一個空矩陣，用來存儲結(jié)果p(:,1)=Q*y0;%fori=1:1:cntt-1u=W*y0+W0*r0;x0=G*x0+H*u;y0=C*x0+D*u;p(:,i+1)=Q*y0;%將y0作為輸出的第1列end13function[t,p]=w_DisNodesSimu22、逐個環(huán)節(jié)刷新的仿真程序該函數(shù)所有參數(shù)與w_DisNodesSimu2的一樣。逐個環(huán)節(jié)刷新是逐個根據(jù)環(huán)節(jié)的編號順序，計算環(huán)節(jié)的輸入，再計算輸出，再計算下一個環(huán)節(jié)的輸入和輸出，依次進(jìn)行，相當(dāng)于環(huán)節(jié)之間串連運(yùn)算。function[t,p]=w_DisNodesSimu(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q)t=[tstart:h:tstop];%t數(shù)一個行序列cntnodes=size(W,1);%得到環(huán)節(jié)的個數(shù)cntt=size(t,2);%返回列數(shù)cnty=size(Q,1);%返回y的維數(shù)p=zeros(cnty,cntt);%構(gòu)造一個空矩陣，用來存儲結(jié)果p(:,1)=Q*y0;%fori=1:1:cntt-1forj=1:1:cntnodesu(j)=W(j,:)*y0+W0(j,:)*r0;x0(j)=G(j,j)*x0(j)+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列end142、逐個環(huán)節(jié)刷新的仿真程序function[t,p]=w_functionsimu4_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);%離散相似法仿真%面向環(huán)節(jié)離散化仿真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’);holdon; stairs(t,y2,':k+'); stairs(t,y3,'-b*'); holdoff;axis([0,10,0,1.2]); xlabel('x');ylabel('y');legend('數(shù)值積分','整體離散','按環(huán)節(jié)離散');title('各環(huán)節(jié)整體刷新（h=0.5）');仿真調(diào)用函數(shù)數(shù)值積分的結(jié)果用plot繪圖，離散模型仿真結(jié)果用stairs繪圖。15functionsimu4_3仿真調(diào)用函數(shù)數(shù)值積分的結(jié)果用整體刷新時相當(dāng)于前后環(huán)節(jié)之間存在一個步長的延遲，導(dǎo)致模型不精確。環(huán)節(jié)逐個刷新符合系統(tǒng)實際運(yùn)算意義，模型比較精確16整體刷新時相當(dāng)于前后環(huán)節(jié)之間存在一個步長的延遲，導(dǎo)致模型不精4.4面向環(huán)節(jié)離散化仿真（有非線性環(huán)節(jié)）環(huán)節(jié)劃分原則1】線性環(huán)節(jié)都處理為1階；2】非線性環(huán)節(jié)附加在線性環(huán)節(jié)的前面或后面，環(huán)節(jié)個數(shù)等于線性環(huán)節(jié)的個數(shù)。一個環(huán)節(jié)的計算順序1】根據(jù)環(huán)節(jié)之間的連接方程計算該環(huán)節(jié)的輸入；2】根據(jù)前置非線性環(huán)節(jié)計算非線性環(huán)節(jié)的輸出；3】前置非線性環(huán)節(jié)的輸出作為線性環(huán)節(jié)的輸入，計算線性環(huán)節(jié)的輸出；4】線性環(huán)節(jié)的輸出作為后置非線性環(huán)節(jié)的輸入，計算整個環(huán)節(jié)的輸出。174.4面向環(huán)節(jié)離散化仿真（有非線性環(huán)節(jié)）環(huán)節(jié)劃分原則一個環(huán)模型描述1】環(huán)節(jié)之間連接方程，寫出各環(huán)節(jié)的輸入與其他環(huán)節(jié)的輸出和系統(tǒng)輸入之間的連續(xù)方程。2】線性環(huán)節(jié)離散模型3】其中18模型描述1】環(huán)節(jié)之間連接方程，寫出各環(huán)節(jié)的輸入與其他環(huán)節(jié)的輸4】非線性環(huán)節(jié)模型描述根據(jù)非線性環(huán)節(jié)類型不同，分別用不同的參數(shù)描述，包括非線性環(huán)節(jié)類型和必要的參數(shù)。非線性環(huán)節(jié)作為線性環(huán)節(jié)的前置或后置環(huán)節(jié)。5】整個系統(tǒng)的輸出方程仿真計算順序根據(jù)環(huán)節(jié)編號順序，計算一個環(huán)節(jié)，再計算下一個環(huán)節(jié)。不能整體計算每個環(huán)節(jié)的輸入，再計算每個環(huán)節(jié)的輸出，這樣每個環(huán)節(jié)之間將產(chǎn)生滯后。194】非線性環(huán)節(jié)模型描述5】整個系統(tǒng)的輸出方程仿真計算順序1【例4.4.1】面向環(huán)節(jié)的非線性系統(tǒng)仿真右圖所示系統(tǒng)有一個飽和非線性環(huán)節(jié)和一個線性環(huán)節(jié)。飽和非線性不能用微分方程表示，因而不能對整個系統(tǒng)采用數(shù)值積分法仿真。采用按環(huán)節(jié)離散化仿真，首先將線性環(huán)節(jié)分解為兩個基本環(huán)節(jié)，將非線性環(huán)節(jié)作為第1個線性環(huán)節(jié)的前置環(huán)節(jié)。仿真中，取c=20，K＝1，r=2520【例4.4.1】面向環(huán)節(jié)的非線性系統(tǒng)仿真右圖所示系統(tǒng)有一個1】各線性環(huán)節(jié)內(nèi)部離散方程（1）積分環(huán)節(jié)（3）慣性環(huán)節(jié)用矩陣表示為211】各線性環(huán)節(jié)內(nèi)部離散方程（1）積分環(huán)節(jié)（3）慣性環(huán)節(jié)用矩陣2】各環(huán)節(jié)之間連接方程3】系統(tǒng)輸出為4】飽和非線性環(huán)節(jié)作為第1個環(huán)節(jié)的前置環(huán)節(jié)222】各環(huán)節(jié)之間連接方程3】系統(tǒng)輸出為4】飽和非線性環(huán)節(jié)作為第function[t,p]=w_DisNodesSimuNonL(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q,NLBefore,NLBeforeParams,NLAfter,NLAfterParams)%函數(shù)功能：帶有非線性環(huán)節(jié)的面向環(huán)節(jié)離散化仿真，%環(huán)節(jié)輸入方程：u(k)=W*x(k)+W0*r(k)%環(huán)節(jié)迭代方程：x(k+1)=G*x(k)+H*u(k)%環(huán)節(jié)輸出方程y(k)=C*x(k)+D*u(k);%整個系統(tǒng)的輸出方程p(k)=Q*y(k);%輸入?yún)?shù)：tstart,tstop,h分別是起始時間、結(jié)束時間和仿真步長，是標(biāo)量%x0,y0是每個環(huán)節(jié)狀態(tài)和輸出的初值，%r0為系統(tǒng)外部輸入%NLBefore是前置非線性環(huán)節(jié)的名稱，如果沒有就寫'None',NLBeforeParams是前置非線性環(huán)節(jié)的參數(shù)，每個環(huán)節(jié)%的參數(shù)占一行。%NLAfter是后置非線性環(huán)節(jié)的名稱，如果沒有就寫'None',NLAfterParams是后置非線性環(huán)節(jié)的參數(shù)，每個環(huán)節(jié)%的參數(shù)占一行。%輸出參數(shù)：t是仿真結(jié)果的時間序列%p是仿真結(jié)果系統(tǒng)的輸出序列1、帶非線性環(huán)節(jié)的按環(huán)節(jié)離散化仿真程序原型23function[t,p]=w_DisNodesSimuNfunction[t,p]=w_DisNodesSimuNonL(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q,NLBefore,NLBeforeParams,NLAfter,NLAfterParams)t=[tstart:h:tstop];%t數(shù)一個行序列cntnodes=size(W,1);%得到環(huán)節(jié)的個數(shù)cntt=size(t,2);%返回列數(shù)cnty=size(Q,1);%返回y的維數(shù)p=zeros(cnty,cntt);%構(gòu)造一個空矩陣，用來存儲結(jié)果p(:,1)=Q*y0;fori=1:1:cntt-1forj=1:1:cntnodesu(j)=W(j,:)*y0+W0(j,:)*r0;%根據(jù)環(huán)節(jié)之間的連接計算，u(j)=NonlinearNode(NLBefore(j),u(j),NLBeforeParams(j,:));%計算第j個環(huán)節(jié)的前置非線性環(huán)節(jié)的輸出x0(j)=G(j,j)*x0(j)+H(j,j)*u(j);y0(j)=C(j,j)*x0(j)+D(j,j)*u(j);%根據(jù)線性部分計算線性部分的輸出y0(j)=NonlinearNode(NLAfter(j),y0(j),NLAfterParams(j,:));%計算第j個環(huán)節(jié)的后置非線性環(huán)節(jié)的輸出end;p(:,i+1)=Q*y0;%將y0作為輸出的第1列end24function[t,p]=w_DisNodesSimuN2、典型非線性環(huán)節(jié)函數(shù)functionu_out=NonlinearNode(NodeType,u_in,params)%函數(shù)功能：典型非線性環(huán)節(jié)的計算%輸入?yún)?shù)：NodeType節(jié)點類型，可以取'Saturation','Deadzone','Relay','RelayDeadzone',%'None'表示沒有非線性環(huán)節(jié)%u_in,環(huán)節(jié)的輸入；%params,環(huán)節(jié)的參數(shù)行向量%輸出結(jié)果：u_out，環(huán)節(jié)的輸出252、典型非線性環(huán)節(jié)函數(shù)functionu_out=NonlswitchNodeTypecase'Saturation'c=params(1);k=params(2);if(k*u_in)<-cu_out=-c;elseif(k*u_in)>cu_out=c;elseu_out=k*u_in;end;case'Deadzone'c=params(1);k=params(2);if(u_in<-c)|(u_in>c)u_out=k*u_in;elseu_out=0;end;case'Relay'c=params(1);ifu_in<0u_out=-c;elseifu_in>0u_out=c;elseu_out=0;end;case'RelayDeadzone'c=params(1);h=params(2);ifu_in<-hu_out=-c;elseifu_in>cu_out=c;elseu_out=0;end;otherwiseu_out=u_in;end;函數(shù)NonlinearNode的主體26switchNodeTypecase'Relay'函functionsimu_NonNodeststart=0;tstop=20;h=0.5;K=1;c=10;W=[0,-1;1,0];W0=[1,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*');xlabel('Time');ylabel('y');title('含非線性環(huán)節(jié)的面向環(huán)節(jié)仿真');3、仿真調(diào)用函數(shù)27functionsimu_NonNodes3、仿真調(diào)用函仿真結(jié)果：同樣的結(jié)果，用plot與stairs繪圖的對比。仿真中，取c=20，K＝1，r=2528仿真結(jié)果：同樣的結(jié)果，用plot與stairs繪圖的對比。仿4.5利用Simulink進(jìn)行面向環(huán)節(jié)仿真【例4.5.1】下圖所示線性與非線性混合系統(tǒng)，當(dāng)回環(huán)非線性特性c=1時，系統(tǒng)有自持振蕩，振幅約為10.8~9.4，振蕩周期約為6s.對該系統(tǒng)采用Simulink建模，進(jìn)行仿真294.5利用Simulink進(jìn)行面向環(huán)節(jié)仿真【例4.5.1】Simulink將此系統(tǒng)看做有6個連續(xù)環(huán)節(jié)。將仿真輸出用Scope繪圖，同時將結(jié)果數(shù)據(jù)輸出到Workspace里，以便用plot繪圖。從曲線，可以很明顯地看到輸出震蕩。30Simulink將此系統(tǒng)看做有6個連續(xù)環(huán)節(jié)。將仿真輸出用Sc【例4.5.2】用Simulik對下圖所示的系統(tǒng)進(jìn)行仿真仿真中，取c=10，K＝1，r=25，并且在t=1時才開始階躍輸入。31【例4.5.2】用Simulik對下圖所示的系統(tǒng)進(jìn)行仿真仿真仿真結(jié)果曲線：很明顯，飽和環(huán)節(jié)的輸出受限制。32仿真結(jié)果曲線：很明顯，飽和環(huán)節(jié)的輸出受限制。32【例4.5.3】對【例4.3.1】所示的系統(tǒng)采用Simulink仿真（1）積分環(huán)節(jié)（2）慣性環(huán)節(jié)當(dāng)T=0.5時33【例4.5.3】對【例4.3.1】所示的系統(tǒng)采用Simul本例采用離散化環(huán)節(jié)進(jìn)行建模34本例采用離散化環(huán)節(jié)進(jìn)行建模34第4章時域離散相似法仿真Matlab編程4.1連續(xù)模型的離散化線性連續(xù)狀態(tài)空間模型為(4.1-1)其采用零階保持器的離散相似模型為(4.1-2)其中(4.1-3)35第4章時域離散相似法仿真Matlab編程4.1連續(xù)模型的用于離散相似的Matlab函數(shù)連續(xù)模型轉(zhuǎn)換為離散模型：[G,H]=c2d(A,B,T)T為采樣周期可選參數(shù)離散相似：[G,H,C,D]=c2dm(A,B,C,D,T,’選項’)[numd,dend]=c2dm(num,den,T,’選項’)36用于離散相似的Matlab函數(shù)連續(xù)模型轉(zhuǎn)換為離散模型：[G,【例4.1.1】對下面的連續(xù)系統(tǒng)，在采樣周期T=0.1時進(jìn)行離散化>>A=[0,1;-2,-3];>>B=[0;1];>>[G,H]=c2d(A,B,0.1);G=0.99090.0861-0.17220.7326H=0.00450.0861經(jīng)計算，可得當(dāng)T=0.1時，計算結(jié)果與c2d函數(shù)得到的結(jié)果一致。37【例4.1.1】對下面的連續(xù)系統(tǒng)，在采樣周期T=0.1時進(jìn)行4.2針對離散狀態(tài)空間模型的仿真程序離散狀態(tài)空間模型為對此模型編寫仿真程序，只需要按模型方程迭代即可。function[t,y]=w_DiscreteSimu(tstart,tstop,h,x0,u0,G,H,C,D)%函數(shù)功能：對線性離散系統(tǒng)x(k+1)=G*X(k)+H*u(k),y(k)=C*x(k)+D*u(k)進(jìn)行仿真%輸入?yún)?shù)：tstart,tstop,h分別是起始時間、結(jié)束時間和仿真步%長，是標(biāo)量%x0,u0是狀態(tài)、輸入的初值，都是列向量%G,H,C,D是線性離散狀態(tài)空間模型的系數(shù)矩陣%輸出參數(shù)：t是仿真結(jié)果的時間序列%y是仿真結(jié)果系統(tǒng)的輸出序列384.2針對離散狀態(tài)空間模型的仿真程序離散狀態(tài)空間模型為對此function[t,y]=w_DiscreteSimu(tstart,tstop,h,x0,u0,G,H,C,D)t=[tstart:h:tstop];%t數(shù)一個行序列cntt=size(t,2);%返回列數(shù)cnty=size(C,1);%返回y的維數(shù)y=zeros(cnty,cntt);%構(gòu)造一個空矩陣，用來存儲結(jié)果y0=C*x0+D*u0;y(:,1)=y0;%將y0作為輸出的第1列curx=x0;%當(dāng)前一步的xcuru=u0;%當(dāng)前一步的ucury=y0;%當(dāng)前一步的yfori=1:1:cntt-1curx=G*curx+H*curu;%計算下一步的狀態(tài)cury=C*curx+D*curu;%計算輸出y(:,i+1)=cury;%將輸出加入到輸出序列里end離散狀態(tài)空間模型仿真程序39function[t,y]=w_DiscreteSimu(【例4.2.1】對下面的連續(xù)系統(tǒng)，采用數(shù)值積分法和離散相似法仿真。解：該系統(tǒng)連續(xù)狀態(tài)空間模型矩陣為其離散狀態(tài)空間模型矩陣為也可以采用c2d函數(shù)直接計算離散模型矩陣40【例4.2.1】對下面的連續(xù)系統(tǒng)，采用數(shù)值積分法和離散相似法仿真程序:用c2d進(jìn)行模型變換，將數(shù)值積分與離散相似2種方法的結(jié)果對比。數(shù)值積分仿真的結(jié)果用plot函數(shù)繪圖，表現(xiàn)連續(xù)的效果；離散模型仿真結(jié)果用stairs繪圖，各數(shù)據(jù)點之間是零階保持。functionsimu4_2tstart=0; tstop=6; 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);%得到離散模型矩陣%數(shù)值積分仿真[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*');holdon;stairs(t,y2,':b+');holdoff;axis([0,6,0,1.2]);xlabel('x');ylabel('y');legend('數(shù)值積分','離散相似');title('數(shù)值積分與離散相似仿真對比')41仿真程序:用c2d進(jìn)行模型變換，將數(shù)值積分與離散相似2種方法在步長h=0.5時，兩種方法的仿真結(jié)果及其接近。離散模型的每個計算點都在連續(xù)模型上當(dāng)步長h=1時，兩種方法的結(jié)果出現(xiàn)偏差。離散模型的點不一定在連續(xù)模型上。42在步長h=0.5時，兩種方法的仿真結(jié)果及其接近。離散模型的每4.3面向環(huán)節(jié)離散化仿真(無非線性環(huán)節(jié))【例4.3.1】將下面的系統(tǒng)按環(huán)節(jié)離散化仿真解：將系統(tǒng)分為2個典型環(huán)節(jié)，如圖所示，開環(huán)部分傳函為1】環(huán)節(jié)之間的連接方程434.3面向環(huán)節(jié)離散化仿真(無非線性環(huán)節(jié))【例4.3.1】將2】環(huán)節(jié)內(nèi)部離散模型（1）積分環(huán)節(jié)（2）慣性環(huán)節(jié)環(huán)節(jié)內(nèi)部動態(tài)方程環(huán)節(jié)輸出方程442】環(huán)節(jié)內(nèi)部離散模型（1）積分環(huán)節(jié)（2）慣性環(huán)節(jié)環(huán)節(jié)內(nèi)部動態(tài)3】系統(tǒng)輸出方程面向環(huán)節(jié)離散化仿真時應(yīng)該注意的問題應(yīng)該依次計算各環(huán)節(jié)的輸出，即計算出第1個環(huán)節(jié)，再利用第1個環(huán)節(jié)的輸出計算第2個環(huán)節(jié)，依次下去。而不要一次計算出各環(huán)節(jié)的輸入，再算出各環(huán)節(jié)的輸出，這樣相當(dāng)于在每個環(huán)節(jié)之間存在一個步長的延遲環(huán)節(jié)，模型的準(zhǔn)確性會大大降低。下面，我們將對這兩種方式做仿真對比，進(jìn)而說明此問題。453】系統(tǒng)輸出方程面向環(huán)節(jié)離散化仿真時應(yīng)該注意的問題11function[t,p]=w_DisNodesSimu2(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q)%函數(shù)功能：面向環(huán)節(jié)離散化仿真，%環(huán)節(jié)輸入方程：u(k)=W*x(k)+W0*r(k)%環(huán)節(jié)迭代方程：x(k+1)=G*x(k)+H*u(k)%環(huán)節(jié)輸出方程y(k)=C*x(k)+D*u(k);%整個系統(tǒng)的輸出方程p(k)=Q*y(k);%輸入?yún)?shù)：tstart,tstop,h分別是起始時間、結(jié)束時間和仿真步長%x0,y0是每個環(huán)節(jié)狀態(tài)和輸出的初值，%r0為系統(tǒng)外部輸入%輸出參數(shù)：t是仿真結(jié)果的時間序列%p是仿真結(jié)果系統(tǒng)的輸出序列1、整體刷新計算的按環(huán)節(jié)仿真函數(shù)整體刷新即根據(jù)環(huán)節(jié)連接方程計算每個環(huán)節(jié)的輸入，再用計算出來的輸入向量根據(jù)各環(huán)節(jié)的離散模型計算各環(huán)節(jié)的輸出。各環(huán)節(jié)之間相當(dāng)于并行運(yùn)算。46function[t,p]=w_DisNodesSimu2function[t,p]=w_DisNodesSimu2(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q)t=[tstart:h:tstop];%t數(shù)一個行序列cntt=size(t,2);%返回列數(shù)cnty=size(Q,1);%返回y的維數(shù)p=zeros(cnty,cntt);%構(gòu)造一個空矩陣，用來存儲結(jié)果p(:,1)=Q*y0;%fori=1:1:cntt-1u=W*y0+W0*r0;x0=G*x0+H*u;y0=C*x0+D*u;p(:,i+1)=Q*y0;%將y0作為輸出的第1列end47function[t,p]=w_DisNodesSimu22、逐個環(huán)節(jié)刷新的仿真程序該函數(shù)所有參數(shù)與w_DisNodesSimu2的一樣。逐個環(huán)節(jié)刷新是逐個根據(jù)環(huán)節(jié)的編號順序，計算環(huán)節(jié)的輸入，再計算輸出，再計算下一個環(huán)節(jié)的輸入和輸出，依次進(jìn)行，相當(dāng)于環(huán)節(jié)之間串連運(yùn)算。function[t,p]=w_DisNodesSimu(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q)t=[tstart:h:tstop];%t數(shù)一個行序列cntnodes=size(W,1);%得到環(huán)節(jié)的個數(shù)cntt=size(t,2);%返回列數(shù)cnty=size(Q,1);%返回y的維數(shù)p=zeros(cnty,cntt);%構(gòu)造一個空矩陣，用來存儲結(jié)果p(:,1)=Q*y0;%fori=1:1:cntt-1forj=1:1:cntnodesu(j)=W(j,:)*y0+W0(j,:)*r0;x0(j)=G(j,j)*x0(j)+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列end482、逐個環(huán)節(jié)刷新的仿真程序function[t,p]=w_functionsimu4_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);%離散相似法仿真%面向環(huán)節(jié)離散化仿真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’);holdon; stairs(t,y2,':k+'); stairs(t,y3,'-b*'); holdoff;axis([0,10,0,1.2]); xlabel('x');ylabel('y');legend('數(shù)值積分','整體離散','按環(huán)節(jié)離散');title('各環(huán)節(jié)整體刷新（h=0.5）');仿真調(diào)用函數(shù)數(shù)值積分的結(jié)果用plot繪圖，離散模型仿真結(jié)果用stairs繪圖。49functionsimu4_3仿真調(diào)用函數(shù)數(shù)值積分的結(jié)果用整體刷新時相當(dāng)于前后環(huán)節(jié)之間存在一個步長的延遲，導(dǎo)致模型不精確。環(huán)節(jié)逐個刷新符合系統(tǒng)實際運(yùn)算意義，模型比較精確50整體刷新時相當(dāng)于前后環(huán)節(jié)之間存在一個步長的延遲，導(dǎo)致模型不精4.4面向環(huán)節(jié)離散化仿真（有非線性環(huán)節(jié)）環(huán)節(jié)劃分原則1】線性環(huán)節(jié)都處理為1階；2】非線性環(huán)節(jié)附加在線性環(huán)節(jié)的前面或后面，環(huán)節(jié)個數(shù)等于線性環(huán)節(jié)的個數(shù)。一個環(huán)節(jié)的計算順序1】根據(jù)環(huán)節(jié)之間的連接方程計算該環(huán)節(jié)的輸入；2】根據(jù)前置非線性環(huán)節(jié)計算非線性環(huán)節(jié)的輸出；3】前置非線性環(huán)節(jié)的輸出作為線性環(huán)節(jié)的輸入，計算線性環(huán)節(jié)的輸出；4】線性環(huán)節(jié)的輸出作為后置非線性環(huán)節(jié)的輸入，計算整個環(huán)節(jié)的輸出。514.4面向環(huán)節(jié)離散化仿真（有非線性環(huán)節(jié)）環(huán)節(jié)劃分原則一個環(huán)模型描述1】環(huán)節(jié)之間連接方程，寫出各環(huán)節(jié)的輸入與其他環(huán)節(jié)的輸出和系統(tǒng)輸入之間的連續(xù)方程。2】線性環(huán)節(jié)離散模型3】其中52模型描述1】環(huán)節(jié)之間連接方程，寫出各環(huán)節(jié)的輸入與其他環(huán)節(jié)的輸4】非線性環(huán)節(jié)模型描述根據(jù)非線性環(huán)節(jié)類型不同，分別用不同的參數(shù)描述，包括非線性環(huán)節(jié)類型和必要的參數(shù)。非線性環(huán)節(jié)作為線性環(huán)節(jié)的前置或后置環(huán)節(jié)。5】整個系統(tǒng)的輸出方程仿真計算順序根據(jù)環(huán)節(jié)編號順序，計算一個環(huán)節(jié)，再計算下一個環(huán)節(jié)。不能整體計算每個環(huán)節(jié)的輸入，再計算每個環(huán)節(jié)的輸出，這樣每個環(huán)節(jié)之間將產(chǎn)生滯后。534】非線性環(huán)節(jié)模型描述5】整個系統(tǒng)的輸出方程仿真計算順序1【例4.4.1】面向環(huán)節(jié)的非線性系統(tǒng)仿真右圖所示系統(tǒng)有一個飽和非線性環(huán)節(jié)和一個線性環(huán)節(jié)。飽和非線性不能用微分方程表示，因而不能對整個系統(tǒng)采用數(shù)值積分法仿真。采用按環(huán)節(jié)離散化仿真，首先將線性環(huán)節(jié)分解為兩個基本環(huán)節(jié)，將非線性環(huán)節(jié)作為第1個線性環(huán)節(jié)的前置環(huán)節(jié)。仿真中，取c=20，K＝1，r=2554【例4.4.1】面向環(huán)節(jié)的非線性系統(tǒng)仿真右圖所示系統(tǒng)有一個1】各線性環(huán)節(jié)內(nèi)部離散方程（1）積分環(huán)節(jié)（3）慣性環(huán)節(jié)用矩陣表示為551】各線性環(huán)節(jié)內(nèi)部離散方程（1）積分環(huán)節(jié)（3）慣性環(huán)節(jié)用矩陣2】各環(huán)節(jié)之間連接方程3】系統(tǒng)輸出為4】飽和非線性環(huán)節(jié)作為第1個環(huán)節(jié)的前置環(huán)節(jié)562】各環(huán)節(jié)之間連接方程3】系統(tǒng)輸出為4】飽和非線性環(huán)節(jié)作為第function[t,p]=w_DisNodesSimuNonL(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q,NLBefore,NLBeforeParams,NLAfter,NLAfterParams)%函數(shù)功能：帶有非線性環(huán)節(jié)的面向環(huán)節(jié)離散化仿真，%環(huán)節(jié)輸入方程：u(k)=W*x(k)+W0*r(k)%環(huán)節(jié)迭代方程：x(k+1)=G*x(k)+H*u(k)%環(huán)節(jié)輸出方程y(k)=C*x(k)+D*u(k);%整個系統(tǒng)的輸出方程p(k)=Q*y(k);%輸入?yún)?shù)：tstart,tstop,h分別是起始時間、結(jié)束時間和仿真步長，是標(biāo)量%x0,y0是每個環(huán)節(jié)狀態(tài)和輸出的初值，%r0為系統(tǒng)外部輸入%NLBefore是前置非線性環(huán)節(jié)的名稱，如果沒有就寫'None',NLBeforeParams是前置非線性環(huán)節(jié)的參數(shù)，每個環(huán)節(jié)%的參數(shù)占一行。%NLAfter是后置非線性環(huán)節(jié)的名稱，如果沒有就寫'None',NLAfterParams是后置非線性環(huán)節(jié)的參數(shù)，每個環(huán)節(jié)%的參數(shù)占一行。%輸出參數(shù)：t是仿真結(jié)果的時間序列%p是仿真結(jié)果系統(tǒng)的輸出序列1、帶非線性環(huán)節(jié)的按環(huán)節(jié)離散化仿真程序原型57function[t,p]=w_DisNodesSimuNfunction[t,p]=w_DisNodesSimuNonL(tstart,tstop,h,x0,y0,r0,W,W0,G,H,C,D,Q,NLBefore,NLBeforeParams,NLAfter,NLAfterParams)t=[tstart:h:tstop];%t數(shù)一個行序列cntnodes=size(W,1);%得到環(huán)節(jié)的個數(shù)cntt=size(t,2);%返回列數(shù)cnty=size(Q,1);%返回y的維數(shù)p=zeros(cnty,cntt);%構(gòu)造一個空矩陣，用來存儲結(jié)果p(:,1)=Q*y0;fori=1:1:cntt-1forj=1:1:cntnodesu(j)=W(j,:)*y0+W0(j,:)*r0;%根據(jù)環(huán)節(jié)之間的連接計算，u(j)=NonlinearNode(NLBefore(j),u(j),NLBeforeParams(j,:));%計算第j個環(huán)節(jié)的前置非線性環(huán)節(jié)的輸出x0(j)=G(j,j)*x0(j)+H(j,j)*u(j);y0(j)=C(j,j)*x0(j)+D(j,j)*u(j);%根據(jù)線性部分計算線性部分的輸出y0(j)=NonlinearNode(NLAfter(j),y0(j),NLAfterParams(j,:));%計算第j個環(huán)節(jié)的后置非線性環(huán)節(jié)的輸出end;p(:,i+1)=Q*y0;%將y0作為輸出的第1列end58function[t,p]=w_DisNodesSimuN2、典型非線性環(huán)節(jié)函數(shù)functionu_out=NonlinearNode(NodeType,u_in,params)%函數(shù)功能：典型非線性環(huán)節(jié)的計算%輸入?yún)?shù)：NodeType節(jié)點類型，可以取'Saturation','Deadzone','Relay','RelayDeadzone',%'None'表示沒有非線性環(huán)節(jié)%u_in,環(huán)節(jié)的輸入；%params,環(huán)節(jié)的參數(shù)行向量%輸出結(jié)果：u_out，環(huán)節(jié)的輸出592、典型非線性環(huán)節(jié)函數(shù)functionu_out=NonlswitchNodeTypecase'Saturation'c=params(1);k=params(2);if(k*u_in)<-cu_out=-c;elseif(k*u_in)>c