實(shí)驗(yàn)：控制系統(tǒng)數(shù)字仿真之?dāng)?shù)值積分法

實(shí)驗(yàn)：控制系統(tǒng)數(shù)字

仿真之?dāng)?shù)值積分法實(shí)驗(yàn)?zāi)康模簩W(xué)會(huì)并掌握數(shù)值積分法的基本原理和方法，了解歐拉法，梯型法，龍格一庫(kù)塔法的區(qū)別，并熟練地使用這些方法。觀察并分析整體離散法、分環(huán)節(jié)離散法、歐拉法、梯形法、龍格?庫(kù)塔法這幾種方法原理上的差別，分析他們各自的優(yōu)缺點(diǎn)。實(shí)驗(yàn)原理：歐拉法：歐拉法是最簡(jiǎn)單的單步法，它是一階的，精度較差。但由于公式簡(jiǎn)單，計(jì)算方便，也易于理解，所以在討論微分方程初值問(wèn)題的數(shù)值解時(shí)通常先討論歐拉法。梯形法：梯形法與歐拉法相比，梯形法的e要比歐拉法的e更接近實(shí)際值，它舍棄的部分更少，它在每一步中用了兩個(gè)點(diǎn)的輸入，使得計(jì)算更加精確。龍格?庫(kù)塔法：龍格一庫(kù)塔法是采用間接利用臺(tái)勞展開(kāi)式的思路，即用在n個(gè)點(diǎn)上的函數(shù)值的線性組合來(lái)代替的導(dǎo)數(shù)，然后按臺(tái)勞展開(kāi)式確定其中的系數(shù),以提高算法的階數(shù)。這樣既能避免計(jì)算函數(shù)的導(dǎo)數(shù)，同時(shí)乂保證了計(jì)算精度。由于龍格薦法具有許務(wù)優(yōu)點(diǎn)，故在許 IM：包中，它是?個(gè)最壟本的算法之一。實(shí)驗(yàn)過(guò)程：分環(huán)節(jié)離散法得出的響應(yīng)曲線：整體離散法得出的響應(yīng)曲線:用一階歐拉法得出的系統(tǒng)響應(yīng)曲線：歐拉法是求出當(dāng)前系統(tǒng)的斜率(變化規(guī)律)，假設(shè)這個(gè)變化規(guī)律在下一次變化前不改變。那么系統(tǒng)下一次值就能夠通過(guò)4.當(dāng)前值2.斜率3.步長(zhǎng)來(lái)確定。比如說(shuō)系統(tǒng)當(dāng)前值x(t),斜率x'(t),仿真步長(zhǎng)dt。那么x(t+dt)=x(t)+x' (t)*dt程序代碼：clc;closeall;clearall;sampleTime=0?l;simuTime=2000;t=sampleTime:sampleTime:simuTime;K=1?2;n=3;T=20;[kp,ki]=PID_Gain(l?20z3,0);x=zeros(lr4);fori=l:fix(simuTime/sampleTime)u(i)=l;endfori=l:fix(simuTime/sampleTime)e=ST_RK_l(X/u(i)fkprkirTzK,n);x=xfe*sampleTime;y(i)=x(4)；endplot(try);匸ext=Tvaiuel(y,sampleTime);legend(text);

自程序ST_RK_1代碼：functionE=ST_RK_1(xrufkpfkizTrKzn)E(l)=(u-x(4))*ki;E(2)=(x(l)+kp*E(l)/ki)*K/T-x(2)/T;E(3)=x(2)/T-x(3)/T;E(4)=x(3)/T-x(4)/T;end用梯形法得出系統(tǒng)響應(yīng)曲線:X用梯形法得出系統(tǒng)響應(yīng)曲線:X=e(r)e［(kH)T］e(kT)牙［ee［(kH)T］e(kT)牙［e(燈)+e［伙+1)門］X(kT)kT(k+l)T上若采用歐拉法，誤差為紅色曲線圍成的面積，而如果用梯形法，誤差減少為藍(lán)色曲線闈成的面積。同時(shí)，要求出藍(lán)色曲線闈成的面積，就要先出下一個(gè)點(diǎn)的值。因此增加了計(jì)算量。算法：先用歐拉法求出下一個(gè)點(diǎn)的值，用下一個(gè)點(diǎn)的值求這個(gè)點(diǎn)的斜率，接著就能求出梯形的面積。用新的面積(代表斜率)求出下一個(gè)點(diǎn)的值。實(shí)驗(yàn)程序代碼(與之前相同的部分沒(méi)有復(fù)制):u(ifl)=u(i);fori=l:fix(simuTime/sampleTime)xO=x;el=ST_RK_l(xru(i),kp,kirTzK);x=xO+el*sampleTime;e2=ST_RK_l(xzu(i+1)/kpfkirTrK);e=(el+e2)/2;x=xO+e*sampleTime;y(i)=x(4);end程序響應(yīng)曲線：用龍格一庫(kù)塔法得出的系統(tǒng)響應(yīng)曲線：當(dāng)采取龍格一庫(kù)塔法(四階)時(shí)，需要用到多個(gè)系統(tǒng)斜率。具體原理不在闡述。程序代碼：u(ifl)=u(i);fori=l:fix(simuTime/sampleTime)xO=x;el=ST_RK_l(xru(i)/kprkirTfK);x=xO+el*sampleTime/2;e2=ST_RK_l(xr(u(i+1)+u(i))/2/kpfkifT,K);x=xO+e2*sampleTime/2;e3=ST_RK_l(xr(u(i+1)+u(i))/2/kprkirTrK);x=xO+e3*sampleTime;e4=ST_RK_l(xzu(i+1)fkpfkifTrK);e=(el+e4)/6+(e2+e3)/3;x=xO+e*sampleTime;y(i)=x(4);end系統(tǒng)的響應(yīng)曲線：從之前的兒個(gè)曲線來(lái)看，除了分布離散法與step（視為標(biāo)準(zhǔn)響應(yīng)曲線）相差比較大之外，剩余的曲線與標(biāo)準(zhǔn)曲線的區(qū)別不大，誤差約為「分之一。接下來(lái)我們要考慮，隨著仿真步距增大，系統(tǒng)失真速度的快慢。系統(tǒng)程序代碼：clc;closeall;clearail;sampleTime=simuTime=500;匸=sampleTime:sampleTime:simuTime;K=1?2;n=3;T=20;[kp,ki]=PID_Gain(l?2Z20f3r0);x=zeros(lf4);fori=l:fix(simuTime/sampleTime)u(i)=l;endu(ifl)=u(i);fori=l:fix(simuTime/sampleTime)e=STRK1(X/u(i),kpfkirT.K,n);x=xfe*sampleTime;yl(i)=x(4)；endx=zeros(lr4);fori=l:fix(simuTime/sampleTime)x0=x;el=ST_RK_l(xru(i),kpfkizTfK);x=xO+el*sampleTime;e2=STRK1(xfu(i+1)fkpfkirTrK);e=(el+e2)/2;x=xO+e*sampleTime;y2(i)=x(4);endx=zeros(lr4);fori=l:fix(simuTime/sampleTime)x0=x;el=ST_RK_l(x,u(i)fkprkizT,K);x=xO+el*sampleTime/2;e2=ST_RK_l(xr(u(i+1)+u(i))/2,kp,ki,T,K);x=x0+e2*sampieTime/2;e3=ST_RK_l(xr(u(i+1)+u(i))/2,kpfkifT,K);x=xO+e3*sampleTime;e4=ST_RK_l(xru(i+1)rkp,ki,TrK);e=(el+e4)/6+(e2+e3)/3;x=xO+e*sampieTime;y3(i)=x(4);endholdon;plot(tryl);pJLot(try2);plot(try3);;text=strvcat(text,Tvalue1(ylFsampleTime));text=stfvuat(textrTvalue1(y2fsampleTime));text二strvcat(text,Tvalue1(y3/sampleTime));legend(text);tir丄e([■仿貞』寸間=Jsi:ir丄]);

根據(jù)上圖，我們可以看出當(dāng)仿真步距比較小的時(shí)候，數(shù)值積分法常用的三種方法得出的實(shí)驗(yàn)精度差不多。但是當(dāng)仿真步距逐漸增大，計(jì)算量較小，解題階次較低的歐拉法最先出現(xiàn)失真，接著是梯形法。貼上一張整體離散法，仿真步距為20時(shí)：當(dāng)仿真步距為1時(shí)，歐拉法的穩(wěn)定時(shí)間達(dá)到199So當(dāng)仿真步距為5時(shí)，分環(huán)節(jié)離散法穩(wěn)定時(shí)間超過(guò)200so仿真步距為10和20時(shí)，整體離散法和梯形法穩(wěn)定時(shí)間均改變，梯形法改變量略低于整體離散法。要說(shuō)一句的是：這里的整體離散法用的n為2。當(dāng)n更大時(shí)，引入的計(jì)算會(huì)更多，但是計(jì)算更精確。當(dāng)系統(tǒng)的階次提高時(shí)，系統(tǒng)的響應(yīng)分別為：當(dāng)輸入為fori=1:fix(simuTime/sampleTime)u(i)=0.8+0.0002*i;endfori=1:fix(simuTime/sampleTime)u(i)=1+0.0002*i-0.0000002*i.A2;end響應(yīng)曲線:我們假設(shè)龍格一庫(kù)塔為標(biāo)準(zhǔn)的響應(yīng)曲線，可以看出歐拉法和梯形法與龍格法相差并不大。自己設(shè)計(jì)的方法：我認(rèn)為在大多數(shù)情況下，輸入為階越信號(hào)或者階越信號(hào)加上一點(diǎn)白噪聲，或者說(shuō)信號(hào)變化穩(wěn)定在一個(gè)較小范圍內(nèi)。因此系統(tǒng)變化的斜率一般會(huì)逐漸變小。如:

因?yàn)橛锰菪畏ㄋ愠龅闹荡蠖鄶?shù)情況下會(huì)小于實(shí)際值(除非斜率為一個(gè)定值),我設(shè)想利用高中的不等式變換，將e稍微的擴(kuò)大，使誤差相抵，減少一部分。編程只要把梯形法中的一句話e=(el+e2)/2;改成：e=sqrt((el.V+W)/2);實(shí)際效果圖：得出系統(tǒng)發(fā)散的結(jié)論，說(shuō)明這個(gè)方法不行。實(shí)驗(yàn)結(jié)果分析：其實(shí)數(shù)值積分法的要領(lǐng)是如何使求岀的e逼近實(shí)際值。用歐拉法會(huì)造成相對(duì)來(lái)說(shuō)十分大的誤差，梯形法次Z,龍格一庫(kù)塔法最接近實(shí)際的e值。因?yàn)樵黾臃抡娌骄鄷?huì)使得系統(tǒng)的準(zhǔn)確性下降，當(dāng)使用相同的較大的仿真步距時(shí)，系統(tǒng)反應(yīng)的越準(zhǔn)確，系統(tǒng)在相同的較小的仿真步距下也應(yīng)該更準(zhǔn)確。那么精確性排爼為：歐拉法小于分環(huán)節(jié)離散法小丁梯形法小J：n取2時(shí)整體離散法小于四階龍格?庫(kù)塔法。在使用數(shù)值積分法時(shí)，如果求e的時(shí)候出現(xiàn)的誤差較小時(shí)，對(duì)系統(tǒng)的影響較小。（比如說(shuō)龍格一庫(kù)塔法和歐拉法之間的誤差）如果是較大的誤差，如戶吾和字的誤差，可導(dǎo)致系統(tǒng)直接發(fā)散。論達(dá)到相同的精確度時(shí)需要的計(jì)算最，龍格庫(kù)塔法小于梯形法小于歐拉法。換句話說(shuō)，在次數(shù)足夠多時(shí)，做了相同次計(jì)算量的龍格庫(kù)塔法要比梯形法和歐拉法精確。（假設(shè)龍格庫(kù)塔每次循環(huán)需要4次計(jì)算，梯形法需要2次計(jì)算，歐拉法需要1次計(jì)算）實(shí)驗(yàn)中出現(xiàn)的問(wèn)題:問(wèn)題有兩個(gè)：梯形法和龍格庫(kù)塔法是因?yàn)樘幚砀唠A輸入信號(hào)時(shí)能夠更準(zhǔn)確，然而在我做的兩次實(shí)驗(yàn)中，當(dāng)輸入高階信號(hào)時(shí)，歐拉法和梯形法能夠跟得上龍格一庫(kù)塔法的變化規(guī)律，而且穩(wěn)定時(shí)間也只差一個(gè)仿真步距，為什么？是因?yàn)槲易龅膶?shí)驗(yàn)次數(shù)不夠多還是