系統(tǒng)辨識(shí)經(jīng)典辨識(shí)方法

1、系統(tǒng)辨識(shí) 經(jīng)典辨識(shí)方法經(jīng)典辨識(shí)方法報(bào)告1. 面積法1.1 辨識(shí)原理1.1.1 分子多項(xiàng)式為1的系統(tǒng) （1.1）由于系統(tǒng)的傳遞函數(shù)與微分方程存在著一一對(duì)應(yīng)的關(guān)系，因此，可以通過求取微分方程的系數(shù)來辨識(shí)系統(tǒng)的傳遞函數(shù)。在求得系統(tǒng)的放大倍數(shù)K后，要先得到無因次階躍響應(yīng)y(t)（設(shè)=0）。大多數(shù)自衡的工業(yè)過程對(duì)象的y(t)可以用下式描述來近似 （1.2）面積法原則上可以求出n為任意階的各系數(shù)。以n=3為例，注意到 （1.3）將式（2.1.2）的y(t)項(xiàng)移至右邊，在0，t上積分，得 （1.4）定義 （1.5）則由式（2.1.3）給出的條件可知，在t （1.6）將式a1y(t)移到等式右邊，定義 （1.

2、7）利用初始條件（2.1.3）當(dāng)t時(shí) （1.8）同理有a3=F3()以此類推，若n2，有an=Fn()1.1.2 分子、分母分別為m階和n階多項(xiàng)式的系統(tǒng)當(dāng)傳遞函數(shù)的形式如下所示時(shí) (1.9)定義 (1.10)由于 (1.11)則的Laplace變換為：(1.12)定義一階面積為：(1.13)令 (1.14)定義二階面積為：(1.15)同理，令 (1.16)定義階面積為。由此可得：(1.17)上式可寫成如下形式：(1.18)(1.19)通過該系數(shù)矩陣，即可求出傳遞函數(shù)分子分母系數(shù)的值。1.2 程序設(shè)計(jì)1.2.1 傳遞函數(shù)形式如式1.1的系統(tǒng)取系統(tǒng)傳遞函數(shù)如下：Gs=13s3+2s2+5.2s+1

3、MATLAB程序如下：clc %清空工作區(qū)cleardt=0.01; %設(shè)置采樣時(shí)間t=0:dt:50; %設(shè)置時(shí)間長(zhǎng)度num=1; %此系統(tǒng)分子為1den=3 2 5.2 1; %分母多項(xiàng)式系數(shù)%繪制原傳遞函數(shù)階躍響應(yīng)曲線fprintf('原系統(tǒng)傳遞函數(shù)為:')G=tf(num,den)y=step(num,den,t);Length=length(y); %數(shù)據(jù)長(zhǎng)度plot(t,y);grid;xlabel('t/s');ylabel('y(t)');%進(jìn)行辨識(shí)設(shè)計(jì)fprintf('辨識(shí)參數(shù)結(jié)果：');%求a1sum1=0;f

4、or(i=1:Length) sum1=sum1+(1-y(i)*dt; F(i)=sum1;end a1=sum1 %求a2 sum2=0; for(i=1:Length) sum2=sum2+(F(i)-a1*y(i)*dt; f(i)=sum2; enda2=sum2 %求a3 sum3=0; for(i=1:Length) sum3=sum3+(f(i)-a2*y(i)*dt; enda3=sum3%繪制辨識(shí)后的傳遞函數(shù)dt=0.01;t=0:dt:50;num2=1;den2=a3 a2 a1 1;fprintf('系統(tǒng)辨識(shí)后的傳遞函數(shù)為:')G=tf(num2,de

5、n2)h=step(num2,den2,t); %辨識(shí)所得傳遞函數(shù)階躍響應(yīng)plot(t,y,'black',t,h,'blue');legend('原傳遞函數(shù)','辨識(shí)所得傳遞函數(shù)');title('原傳遞函數(shù)與辨識(shí)所得傳遞函數(shù)的階躍響應(yīng)對(duì)比')grid;xlabel('t/s');ylabel('y(t)和h(t)');fprintf('相關(guān)系數(shù)：'); %求相關(guān)系數(shù)r=corrcoef(y,h)運(yùn)行以上程序得到結(jié)果如下：原系統(tǒng)傳遞函數(shù)為:G = 1 - 3 s3

6、+ 2 s2 + 5.2 s + 1 Continuous-time transfer function.辨識(shí)參數(shù)結(jié)果：a1 = 5.2048a2 = 2.0608a3 = 2.8388系統(tǒng)辨識(shí)后的傳遞函數(shù)為:G = 1 - 2.839 s3 + 2.061 s2 + 5.205 s + 1 Continuous-time transfer function.相關(guān)系數(shù)：r = 1.0000 0.99990.9999 1.0000此時(shí)原傳遞函數(shù)和辨識(shí)所得傳遞函數(shù)的階躍響應(yīng)對(duì)比如下圖：圖1.1 原傳遞函數(shù)和辨識(shí)所得傳遞函數(shù)的階躍響應(yīng)對(duì)比由上圖可以看出，辨識(shí)所得結(jié)果比較準(zhǔn)確。1.2.1 傳遞函數(shù)形式

7、如式(1.9)的系統(tǒng)（無噪聲）取系統(tǒng)傳遞函數(shù)如下：Gs=s2+3s+12s3+5.5s2+3s+1MATLAB程序如下：clc %清空工作區(qū)cleardt=0.01; %設(shè)置采樣時(shí)間t=0:dt:20; %設(shè)置時(shí)間長(zhǎng)度num=1 3 1; %分子多項(xiàng)式系數(shù)den=2 5.5 3 1; %分母多項(xiàng)式系數(shù)%繪制原傳遞函數(shù)階躍響應(yīng)曲線fprintf('原系統(tǒng)傳遞函數(shù)為:')G=tf(num,den)y=step(num,den,t);Length=length(y); plot(t,y);grid;xlabel('t/s');ylabel('y(t)'

8、);%辨識(shí)程序設(shè)計(jì)，此系統(tǒng)m+n=5,故應(yīng)計(jì)算A1-A5fprintf('A1-A5階面積分別為：')%求A1sum1=0;for(i=1:Length-1) sum1=sum1+(1-(y(i)+y(i+1)/2)*dt; A(i)=sum1;endA1=sum1%求A2sum2=0;for(i=1:Length-1) sum2=sum2+(A(i)-A1*(y(i)+y(i+1)/2)*dt; B(i)=sum2;endA2=sum2%求A3sum3=0;for(i=1:Length-1) sum3=sum3+(B(i)-A2*(y(i)+y(i+1)/2)*dt; C(i

9、)=sum3;endA3=sum3%求A4sum4=0;for(i=1:Length-1) sum4=sum4+(C(i)-A3*(y(i)+y(i+1)/2)*dt; D(i)=sum4;endA4=sum4%求A5sum5=0;for(i=1:Length-1) sum5=sum5+(D(i)-A4*(y(i)+y(i+1)/2)*dt;endA5=sum5%求分子系數(shù)b1,b2M=(-1)*(inv(A3,A2;A4,A3)*A4;A5;fprintf('分子多項(xiàng)式系數(shù)為：')b1=M(1,1)b2=M(2,1)%求分母系數(shù)a1,a2,a3N=1 0 0;A1 1 0;A

10、2 A1 1*b1;b2;0+A1;A2;A3;fprintf('分母多項(xiàng)式系數(shù)為：')a1=N(1,1)a2=N(2,1)a3=N(3,1)%求辨識(shí)所得傳遞函數(shù)num1=b2 b1 1;den1=a3 a2 a1 1;fprintf('辨識(shí)所得傳遞函數(shù)為:')G=tf(num1,den1)h=step(num1,den1,t);plot(t,y,'black',t,h,'blue');legend('原傳遞函數(shù)','辨識(shí)所得傳遞函數(shù)');title('原傳遞函數(shù)與辨識(shí)所得傳遞函數(shù)的階躍響應(yīng)

11、對(duì)比')grid;xlabel('t/s');ylabel('y(t)');fprintf('相關(guān)系數(shù)：'); %求相關(guān)系數(shù)r=corrcoef(y,h)運(yùn)行以上程序結(jié)果如下：原系統(tǒng)傳遞函數(shù)為:G = s2 + 3 s + 1 - 2 s3 + 5.5 s2 + 3 s + 1 Continuous-time transfer function.A1-A5階面積分別為：A1 = 0.0076A2 = 4.3364A3 = -9.6285A4 = 15.7297A5 = 7.1087分子多項(xiàng)式系數(shù)為：b1 = 7.4409b2 = 12.8

12、942分母多項(xiàng)式系數(shù)為：a1 = 7.4485a2 = 17.2869a3 = 22.7359辨識(shí)所得傳遞函數(shù)為:G = 12.89 s2 + 7.441 s + 1 - 22.74 s3 + 17.29 s2 + 7.448 s + 1 Continuous-time transfer function.相關(guān)系數(shù)：r = 1.0000 0.9996 0.9996 1.0000此時(shí)原傳遞函數(shù)和辨識(shí)所得傳遞函數(shù)的階躍響應(yīng)對(duì)比如下圖：圖1.2 原傳遞函數(shù)和辨識(shí)所得傳遞函數(shù)的階躍響應(yīng)對(duì)比由上圖可以看出，在未加入噪聲之前，采用面積法辨識(shí)結(jié)果很精確，并且，分子可以為階次低于分母的任意階次。1.2.2 加

13、入噪聲為便于對(duì)比，仍取系統(tǒng)傳遞函數(shù)如下：Gs=s2+3s+12s3+5.5s2+3s+1MATLAB程序如下：clc %清空工作區(qū)cleardt=0.01; %設(shè)置采樣時(shí)間t=0:dt:20; %設(shè)置時(shí)間長(zhǎng)度pj=0; fc=sqrt(0.01); %加入方差為0.01的白噪聲序列num=1 3 1; %分子多項(xiàng)式系數(shù)den=2 5.5 3 1; %分母多項(xiàng)式系數(shù)%繪制原傳遞函數(shù)階躍響應(yīng)曲線fprintf('原系統(tǒng)傳遞函數(shù)為:')G=tf(num,den)L=length(t);cs=randn(1,L); %產(chǎn)生噪聲cs=cs/std(cs); cs=cs-mean(cs);

14、 cs=pj+fc*cs;cs=cs'y=step(num,den,t)+cs; %加入噪聲Length=length(y); plot(t,y);grid;xlabel('t/s');ylabel('y(t)');%辨識(shí)程序設(shè)計(jì)，此系統(tǒng)m+n=5,故應(yīng)計(jì)算A1-A5fprintf('A1-A5階面積分別為：')%求A1sum1=0;for(i=1:Length-1) sum1=sum1+(1-(y(i)+y(i+1)/2)*dt; A(i)=sum1;endA1=sum1%求A2sum2=0;for(i=1:Length-1) sum2

15、=sum2+(A(i)-A1*(y(i)+y(i+1)/2)*dt; B(i)=sum2;endA2=sum2%求A3sum3=0;for(i=1:Length-1) sum3=sum3+(B(i)-A2*(y(i)+y(i+1)/2)*dt; C(i)=sum3;endA3=sum3%求A4sum4=0;for(i=1:Length-1) sum4=sum4+(C(i)-A3*(y(i)+y(i+1)/2)*dt; D(i)=sum4;endA4=sum4%求A5sum5=0;for(i=1:Length-1) sum5=sum5+(D(i)-A4*(y(i)+y(i+1)/2)*dt;en

16、dA5=sum5%求分子系數(shù)b1,b2M=(-1)*(inv(A3,A2;A4,A3)*A4;A5;fprintf('分子多項(xiàng)式系數(shù)為：')b1=M(1,1)b2=M(2,1)%求分母系數(shù)a1,a2,a3N=1 0 0;A1 1 0;A2 A1 1*b1;b2;0+A1;A2;A3;fprintf('分母多項(xiàng)式系數(shù)為：')a1=N(1,1)a2=N(2,1)a3=N(3,1)%求辨識(shí)所得傳遞函數(shù)num1=b2 b1 1;den1=a3 a2 a1 1;fprintf('辨識(shí)所得傳遞函數(shù)為:')G=tf(num1,den1)h=step(num1,

17、den1,t);plot(t,y,'black',t,h,'blue');legend('原傳遞函數(shù)','辨識(shí)所得傳遞函數(shù)');title('原傳遞函數(shù)與辨識(shí)所得傳遞函數(shù)的階躍響應(yīng)對(duì)比')grid;xlabel('t/s');ylabel('y(t)');fprintf('相關(guān)系數(shù)：'); %求相關(guān)系數(shù)r=corrcoef(y,h)運(yùn)行程序，結(jié)果如下：原系統(tǒng)傳遞函數(shù)為:G = s2 + 3 s + 1 - 2 s3 + 5.5 s2 + 3 s + 1 Continu

18、ous-time transfer function.A1-A5階面積分別為：A1 = 0.0081A2 = 4.7477A3 = -13.2541A4 = 39.7921A5 = -121.2501分子多項(xiàng)式系數(shù)為：b1 = 3.6415b2 = 1.7846分母多項(xiàng)式系數(shù)為：a1 = 3.6496a2 = 6.5619a3 = 4.0492辨識(shí)所得傳遞函數(shù)為:G = 1.785 s2 + 3.641 s + 1 - 4.049 s3 + 6.562 s2 + 3.65 s + 1 Continuous-time transfer function.相關(guān)系數(shù)：r = 1.0000 0.907

19、6 0.9076 1.0000>>此時(shí)原傳遞函數(shù)和辨識(shí)所得傳遞函數(shù)的階躍響應(yīng)對(duì)比如下圖：圖1.3 加入噪聲后原傳遞函數(shù)和辨識(shí)所得傳遞函數(shù)的階躍響應(yīng)對(duì)比由運(yùn)行結(jié)果和上圖可以看出，加入噪聲后，系統(tǒng)的辨識(shí)結(jié)果受到了影響，準(zhǔn)確度下降，但基本上可以接受。2 差分方程法2.1 辨識(shí)原理設(shè)過程傳遞函數(shù)為 (2.1)則當(dāng)特征方程有n個(gè)單根s1,s2,···sn時(shí)，傳遞函數(shù)可寫成 (2.2)當(dāng)特征方程有n-r個(gè)單根s1,s2,···sn-r，r階重根s0時(shí)，傳遞函數(shù)可寫成 (2.3)為了確定ci和si，從獲取的過程輸出脈沖響應(yīng)g(t)中，選

20、取前n+1個(gè)坐標(biāo)點(diǎn)，每個(gè)坐標(biāo)點(diǎn)間隔相同的采樣時(shí)間T0，組成一個(gè)自回歸模型 (2.4)其中1，2n為待定系數(shù)。如果特征方程 (2.5)有一個(gè)單根，則其必是AR模型的解，它們的線性組合也是AR模型的解。當(dāng)其特征方程有n個(gè)單根時(shí)，自回歸模型的解為 (2.6)當(dāng)其特征方程有n-r個(gè)單根r階重根時(shí)，自回歸模型的解為 (2.7)根據(jù)分析，無論特征方程是單根還是重根，都有 (2.8)顯然，一旦求出xi和i，便可的傳遞函數(shù)。2.1 程序設(shè)計(jì)選取系統(tǒng)傳遞函數(shù)如下：Gs=2s+1.05s3+5.5s2+3s+1MATLAB程序設(shè)計(jì)如下：clc %清空工作區(qū)cleardt=0.01; %設(shè)置采樣時(shí)間t=0:dt:5

21、0; %設(shè)置時(shí)間長(zhǎng)度num=2 1.05; %分子多項(xiàng)式系數(shù)den=1 5.5 3 1; %分母多項(xiàng)式系數(shù)%繪制原傳遞函數(shù)階躍響應(yīng)曲線fprintf('原系統(tǒng)傳遞函數(shù)為:')G=tf(num,den)y=impulse(num,den,t); %脈沖響應(yīng)曲線Length=length(y); plot(t,y);grid;xlabel('t/s');ylabel('y(t)');%時(shí)刻為0s,1s 2s 3s 4s 5s時(shí)的脈沖響應(yīng),取T0=1fprintf('時(shí)刻為0s,1s 2s 3s 4s 5s時(shí)的脈沖響應(yīng),取T0=1');

22、Y=y(1,:) y(101,:) y(201,:) y(301,:) y(401,:) y(501,:) A=Y(:,2:4);Y(:,3:5);Y(:,4:6); %構(gòu)造系數(shù)矩陣B=Y(:,1:3)' %構(gòu)造方程右邊的數(shù)值fprintf('系數(shù)：');a=(inv(A)*(-1)*B) %求得待定系數(shù),在這里用a表示p=a(3,1) a(2,1) a(1,1) 1; %關(guān)于x的方程系數(shù)x=roots(p); %求得x1,x2,x3fprintf('系統(tǒng)極點(diǎn)：');s=log(roots(p) %求得傳遞函數(shù)極點(diǎn)%構(gòu)造系數(shù)矩陣，并求得，這里用E表示C=

23、1 1 1;x(1,1) x(2,1) x(3,1);(x(1,1)2) (x(2,1)2 (x(3,1)2);D=0 Y(:,2) Y(:,3)'fprintf('系數(shù)');E=(inv(C)*D % 繪制辨識(shí)所得傳遞函數(shù)的脈沖響應(yīng)num1,den1=residue(E,s,);fprintf('辨識(shí)所得傳遞函數(shù)');G1=tf(num1,den1) %傳遞函數(shù)h=impulse(num1,den1,t); %脈沖響應(yīng)曲線plot(t,y,'black',t,h,'blue');legend('原傳遞函數(shù)'

24、;,'辨識(shí)所得傳遞函數(shù)');title('原系統(tǒng)與辨識(shí)所得系統(tǒng)脈沖響應(yīng)對(duì)比');grid;xlabel('t/s');ylabel('y(t)和h(t)');運(yùn)行以上程序，得到結(jié)果如下：原系統(tǒng)傳遞函數(shù)為:G = 2 s + 1.05 - s3 + 5.5 s2 + 3 s + 1 Continuous-time transfer function.時(shí)刻為0s,1s 2s 3s 4s 5s時(shí)的脈沖響應(yīng),取T0=1Y = 0 0.3685 0.2956 0.2076 0.1261 0.0607系數(shù)：a = -141.2815 348.0135 -244.6919系統(tǒng)極點(diǎn)：s = -0.2835 + 0.3498i -0.2835 - 0.3498i -4.9329 系數(shù)E = 0.2028 - 0.1637i 0.2028 + 0.1637i -0.4055