控制系統(tǒng)仿真

1、 控制系統(tǒng)仿真實(shí)驗(yàn)報(bào)告姓名：大蔥哥學(xué)號(hào)：班級(jí)：測(cè)控12022015.10.28實(shí)驗(yàn)一 經(jīng)典的連續(xù)系統(tǒng)仿真建模方法一 實(shí)驗(yàn)?zāi)康? 了解和掌握利用仿真技術(shù)對(duì)控制系統(tǒng)進(jìn)行分析的原理和步驟。2 掌握機(jī)理分析建模方法。3 深入理解一階常微分方程組數(shù)值積分解法的原理和程序結(jié)構(gòu)，學(xué)習(xí)用 Matlab 編寫數(shù)值積分法仿真程序。4 掌握和理解四階 Runge-Kutta 法，加深理解仿真步長(zhǎng)與算法穩(wěn)定性的關(guān)系。二 實(shí)驗(yàn)內(nèi)容1. 編寫四階 Runge_Kutta 公式的計(jì)算程序， 對(duì)非線性模型 （3） 式進(jìn)行仿真。仿真程序:主程序：clcclear allu=zeros(2,1);u(1)=0.5%穩(wěn)態(tài)%u(1)

2、=0.55;%閥位增大%10%u(1)=0.45;%閥位減小%10u(2)=0.15; h=zeros(1,2);h(1,1)=1.5;h(1,2)=1.4;hStep = 10;Hlevel = h;nCounter = 25;for t=0:hStep:(nCounter-1)*hStep h = Z06_SystemSimulation_Lab01_Nonlinear_RK4(hStep,t,h,u); Hlevel=Hlevel;h;endfigure(1)plot(0:hStep:nCounter*hStep,Hlevel)grid四階龍格庫塔算法：function RK4_Resu

3、lt = Z06_SystemSimulation_Lab01_Nonlinear_RK4(h,t0,x0,u0)K1=Z06_SystemSimulation_Lab01_Linear_dxCompute(t0,x0,u0);K2=Z06_SystemSimulation_Lab01_Linear_dxCompute(t0+h/2,x0+h*K1/2,u0);K3=Z06_SystemSimulation_Lab01_Linear_dxCompute(t0+h/2,x0+h*K2/2,u0);K4=Z06_SystemSimulation_Lab01_Linear_dxCompute(t0+

4、h,x0+h*K3,u0);RK4_Result = x0 + h*(K1+2*K2+2*K3+K4)/6;（1） 將閥位 u 增大 10和減小 10，觀察響應(yīng)曲線的形狀;穩(wěn)態(tài)仿真曲線：增大10%：減小10%：（2） 研究仿真步長(zhǎng)對(duì)穩(wěn)定性的影響，仿真步長(zhǎng)取多大時(shí) RK4 算法變得不穩(wěn)定？仿真步長(zhǎng)的選取范圍可以通過理論計(jì)算獲知，這里我們通過仿真得到當(dāng)步長(zhǎng)大于63時(shí)，曲線發(fā)散，系統(tǒng)不再穩(wěn)定，發(fā)散曲線如下圖：（3） 利用 MATLAB 中的 ode45()函數(shù)進(jìn)行求解， 比較與 （1） 中的仿真結(jié)果有何區(qū)別。2. 編寫四階 Runge_Kutta 公式的計(jì)算程序，對(duì)線性狀態(tài)方程（18）式進(jìn)行仿真仿真

5、程序clcclear allclose allu=zeros(2,1);u(1)=0.0;%穩(wěn)態(tài)% u(1)=0.05;%閥位增大%10% u(1)=-0.05;%閥位減小%10u(2)=0; h=zeros(1,2);h(1,1)=0;h(1,2)=0;hStep = 10;Hlevel = h;nCounter = 25;for t=0:hStep:(nCounter-1)*hStep h = Z06_SystemSimulation_Lab01_Linear_RK4(hStep,t,h,u); Hlevel=Hlevel;h;end Hlevel(:,1) = Hlevel(:,1)+1

6、.5; Hlevel(:,2) = Hlevel(:,2)+1.4;figure(3)plot(0:hStep:nCounter*hStep,Hlevel)gridhold on（1） 將閥位增大 10和減小 10，觀察響應(yīng)曲線的形狀;穩(wěn)態(tài)仿真曲線：閥位增大%10：閥位減小%10：（2） 研究仿真步長(zhǎng)對(duì)穩(wěn)定性的影響，仿真步長(zhǎng)取多大時(shí) RK4 算法變得不穩(wěn)定？當(dāng)步長(zhǎng)大于62時(shí)，算法變得不再穩(wěn)定（3） 閥位增大 10和減小 10，利用 MATLAB 中的 ode45()函數(shù)進(jìn)行求解階躍響應(yīng)，比較與（1）中的仿真結(jié)果有何區(qū)別。三 思考題1 討論仿真步長(zhǎng)對(duì)穩(wěn)定性和仿真精度的影響。仿真步長(zhǎng)越長(zhǎng)，系統(tǒng)越傾

7、向于不穩(wěn)定，仿真步長(zhǎng)越短，系統(tǒng)越穩(wěn)定，但是需要的仿真時(shí)間越長(zhǎng)。2 你是怎樣實(shí)現(xiàn)閥位增大和減小10%的？對(duì)于非線性模型和線性模型方法一樣嗎？ 通過改變u(1)的初值可以實(shí)現(xiàn)對(duì)閥位的增大和減小，由于線性模型法是對(duì)變量的變化量進(jìn)行計(jì)算的，而非線性模型法是直接對(duì)變量進(jìn)行計(jì)算，所以兩種方法實(shí)現(xiàn)對(duì)閥位的增大和減小方面存在一些區(qū)別，不完全一致。四 實(shí)驗(yàn)總結(jié)通過本次實(shí)驗(yàn)我學(xué)會(huì)了Matlab的基本使用方法，對(duì)該軟件有了一定程度的認(rèn)識(shí)，同時(shí)對(duì)仿真步長(zhǎng)對(duì)系統(tǒng)的穩(wěn)定性和仿真精度的關(guān)系有了一定的理解，對(duì)一階常微分方程組數(shù)值積分解法的原理和程序結(jié)構(gòu)有了進(jìn)一步的掌握，深入了解了Runge-Kutta法，對(duì)MATLAB軟件的

8、仿真有了一定程度的認(rèn)識(shí)，為以后的學(xué)習(xí)奠定了基礎(chǔ)。實(shí)驗(yàn)二 面向結(jié)構(gòu)圖的仿真第一部分 線性系統(tǒng)仿真一 實(shí)驗(yàn)?zāi)康? 掌握理解控制系統(tǒng)閉環(huán)仿真技術(shù)。2 掌握理解面向結(jié)構(gòu)圖的離散相似法的原理和程序結(jié)構(gòu)。3掌握 MATLAB 中 C2D 函數(shù)的用法，掌握雙線性變換的原理。 二 實(shí)驗(yàn)內(nèi)容 1.編寫仿真程序，實(shí)現(xiàn)無擾動(dòng)時(shí)給定值階躍仿真實(shí)驗(yàn)仿真程序：%面向結(jié)構(gòu)圖的離散相似法,線性模型clcclear all%定義參數(shù)%A=2; %水箱橫截面積ku=0.1/0.5; %閥門流量系數(shù)H10=1.5; %水箱1的平衡高度H20=1.4; %水箱2的平衡高度 alpha12 = 0.25/sqrt(H10); %水箱1

9、流向水箱2流量系數(shù)alpha2 = 0.25/sqrt(H20); %水箱2流出水流量系數(shù)R12=2*sqrt(H10)/alpha12; %線性化數(shù)學(xué)模型中的參數(shù)R2=2*sqrt(H20)/alpha2; %線性化數(shù)學(xué)模型中的參數(shù)H1SpanLo=0; %水箱1量程下限H2SpanLo=0; %水箱2量程下限H1SpanHi=2.52; %水箱1量程上限H2SpanHi=2.52; %水箱2量程上限Kp=1.78; %PI控制器比例系數(shù)Ti=85; %PI控制器積分時(shí)間常數(shù)R12*AR12Kc=Kp/Ti; %式(4)中參數(shù)bc=Ti; %式(4)中參數(shù)Kd = 1/A; %式(5)中參數(shù)

10、 ad = 1/(A*R12); %式(5)中參數(shù)a1 = 1/(A*R12); %式(6)中參數(shù)K1 = ku/A; %式(6)中參數(shù)a2 = 1/(A*R2); %式(7)中參數(shù)K2 = 1/(A*R12); %式(7)中參數(shù)%定義初值%uc(1)=0; %Gc環(huán)節(jié)(8)的輸入初值xc(1)=0; %Gc環(huán)節(jié)(8)的狀態(tài)初值yc(1)=0; %Gc環(huán)節(jié)(8)的輸出初值ud(1)=0; %Gd環(huán)節(jié)(9)的輸入初值xd(1)=0; %Gd環(huán)節(jié)(9)的狀態(tài)初值yd(1)=0; %Gd環(huán)節(jié)(9)的輸出初值 u1(1)=0; %G1環(huán)節(jié)(10)的輸入初值x1(1)=0; %G1環(huán)節(jié)(10)的狀態(tài)初值

11、y1(1)=0; %G1環(huán)節(jié)(10)的輸出初值u2(1)=0; %G2環(huán)節(jié)(11)的輸入初值x2(1)=0; %G2環(huán)節(jié)(11)的狀態(tài)初值y2(1)=0; %G2環(huán)節(jié)(11)的輸出初值%nCounter = 70; %迭代次數(shù)T=10; %步長(zhǎng)k=1;deltaQd=0; %擾動(dòng)量H20_percent=(H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100; %穩(wěn)態(tài)工作點(diǎn)的液位百分比H2=80; %給定液位百分比tend = nCounter*T;for t=T:T:tend k=k+1; uc(k)= (H2 - (y2(k-1)+H20-H2SpanLo)/(H2S

12、panHi-H2SpanLo)*100)/100; %系統(tǒng)(8)的輸入百分比 ud(k)=deltaQd; %系統(tǒng)(9)的輸入 u1(k)=yc(k-1); %系統(tǒng)(10)的輸入 u2(k)=y1(k-1); %系統(tǒng)(11)的輸入 xc(k) = xc(k-1) + Kc*T*uc(k-1); yc(k)=xc(k)+bc*Kc*uc(k); xd(k) = exp(-ad*T)*xd(k-1) + Kd/ad*(1-exp(-ad*T)*ud(k);yd(k)=xd(k); x1(k) = exp(-a1*T)*x1(k-1) + K1/a1*(1-exp(-a1*T)*u1(k);y1(k

13、)=x1(k); x2(k) = exp(-a2*T)*x2(k-1) + K2/a2*(1-exp(-a2*T)*u2(k);y2(k)=x2(k); endHlevel(:,1)=(y1+H10-H1SpanLo)/(H1SpanHi-H1SpanLo)*100;Hlevel(:,2)=(y2+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100;yc=(yc+0.5)*100;y2sp=H2*ones(size(y1);yv=yc; textPositionH1=max(Hlevel(:,1); textPositionH2=max(Hlevel(:,2); H2S

14、teady=Hlevel(size(Hlevel(:,1),1),1)*ones(size(y1); xmax=max(0:T:tend); xmin=0; ymax=110; ymin=50; scrsz = get(0,ScreenSize); gca=figure(Position,5 10 scrsz(3)-10 scrsz(4)-90); %gca=figure(Position,5 10 scrsz(3)/2 scrsz(4)/1.5) set(gca,Color,w); plot(0:T:tend,Hlevel(:,1),r,LineWidth,2) hold on plot(0

15、:T:tend,Hlevel(:,2),b,LineWidth,2) hold on plot(0:T:tend,yv,k,LineWidth,2) hold on plot(0:T:tend,y2sp,g,LineWidth,2) hold on plot(0:T:tend,H2Steady,y,LineWidth,2) line(tend/2 tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/10 (ymax-ymin)/2+ymin-(ymax-ymin)/10,Color,r,LineWidth,6) text(tend/2+27,(ymax-ymin)

16、/2+ymin-(ymax-ymin)/10, 第1個(gè)水箱的液位H1,FontSize,16) line(tend/2 tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/6 (ymax-ymin)/2+ymin-(ymax-ymin)/6,Color,b,LineWidth,6) text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/6, 第2個(gè)水箱的液位H2,FontSize,16) line(tend/2 tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/4.2 (ymax-ymin)/2

17、+ymin-(ymax-ymin)/4.2,Color,g,LineWidth,6) text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/4.2, 第2個(gè)水箱的液位給定值,FontSize,16) line(tend/2 tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/3.2 (ymax-ymin)/2+ymin-(ymax-ymin)/3.2,Color,k,LineWidth,6) text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/3.2, 閥位變化情況,FontSize,1

18、6) axis(xmin xmax ymin ymax); text(tend/5,ymax+1.5, 實(shí)驗(yàn)二 不考慮閥位飽和特性時(shí)的控制效果,FontSize,22) grid %hold off仿真實(shí)驗(yàn)曲線：由控制的效果可以看出，閥位的變化出現(xiàn)了超過100%的情況，這是控制系統(tǒng)中不允許出現(xiàn)的，說明控制功能需要改進(jìn)。三 思考題在未考慮調(diào)節(jié)閥飽和特性時(shí)， 討論一下兩個(gè)水箱液位的變化情況， 工業(yè)上是否允許？討論閥位的變化情況，工業(yè)上是否能實(shí)現(xiàn)？答：在未考慮調(diào)節(jié)閥飽和特性時(shí)，兩個(gè)水箱的液位在閥位一開始大開時(shí)，H1，H2液位上升迅速，很快就達(dá)到預(yù)期值。但是其控制不能在工業(yè)上實(shí)現(xiàn)，沒有實(shí)際可行度。除此

19、之外，閥位有其本身的最大最小的限制，在仿真中出現(xiàn)的超過100%的情況在現(xiàn)實(shí)生活中不可能出現(xiàn)，因此這一部分對(duì)應(yīng)的控制效果在工業(yè)上不能實(shí)現(xiàn)。四 實(shí)驗(yàn)總結(jié) 通過本次實(shí)驗(yàn)，我對(duì)控制系統(tǒng)閉環(huán)仿真技術(shù)有了進(jìn)一步的認(rèn)識(shí)，掌握了面向結(jié)構(gòu)圖的離散相似法的原理和程序結(jié)構(gòu)，了解了雙線性變換的原理，收獲頗豐。第二部分 含有非線性環(huán)節(jié)的控制系統(tǒng)仿真一 實(shí)驗(yàn)?zāi)康? 掌握理解控制系統(tǒng)閉環(huán)仿真技術(shù)。5 掌握理解面向結(jié)構(gòu)圖的離散相似法的原理和程序結(jié)構(gòu)。6 掌握理含有非線性環(huán)節(jié)的控制系統(tǒng)的仿真方法。二 實(shí)驗(yàn)內(nèi)容根據(jù)上面的各式，編寫仿真程序，實(shí)現(xiàn)無擾動(dòng)時(shí)給定值階躍仿真實(shí)驗(yàn)仿真程序：%面向結(jié)構(gòu)圖的離散相似法,非線性模型clcclea

20、r allA=2;ku=0.1/0.5;H10=1.5;H20=1.4;alpha12 = 0.25/sqrt(H10);alpha2 = 0.25/sqrt(H20);R12=2*sqrt(H10)/alpha12;R2=2*sqrt(H20)/alpha2;H1SpanLo=0;H2SpanLo=0;H1SpanHi=2.52;H2SpanHi=2.52; Kp=3.91/2.2;Ti=0.85*100;%Kp=3.21;%Ti=99999999999999;ad = 1/(A*R12);a1 = 1/(A*R12);a2 = 1/(A*R2);Kc=Kp/Ti;bc=Ti;Kd = 1/

21、A;K1 = ku/A;K2 = 1/(A*R12);uc(1)=0;uv(1)=0;ud(1)=0;u1(1)=0;u2(1)=0;xc(1)=0;xv(1)=0;xd(1)=0;x1(1)=0;x2(1)=0;yc(1)=0;yv(1)=0;yd(1)=0;y1(1)=0;y2(1)=0;nCounter = 70;T=10;k=1;deltaQd=0;c=0.5;H20_percent=(H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100;H2set_percent=82;tend = nCounter*T;for t=T:T:tend k=k+1; uc(k)=

22、 (H2set_percent - (y2(k-1)+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100)/100; uv(k)=yc(k-1); ud(k)=deltaQd; if uv(k)c yv(k)=c; end if uv(k)-c yv(k)=0; end if uv(k)=-c yv(k)=uv(k); end u1(k)=yv(k); u2(k)=y1(k-1); xc(k) = xc(k-1) + Kc*T*uc(k-1); yc(k)=xc(k)+bc*Kc*uc(k); xd(k) = exp(-ad*T)*xd(k-1) + Kd/ad*(1-

23、exp(-ad*T)*ud(k);yd(k)=xd(k); x1(k) = exp(-a1*T)*x1(k-1) + K1/a1*(1-exp(-a1*T)*u1(k);y1(k)=x1(k); x2(k) = exp(-a2*T)*x2(k-1) + K2/a2*(1-exp(-a2*T)*u2(k);y2(k)=x2(k);endHlevel(:,1)=(y1+H10-H1SpanLo)/(H1SpanHi-H1SpanLo)*100;Hlevel(:,2)=(y2+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100;yv=(yv+0.5)*100;y2sp=H2s

24、et_percent*ones(size(y1); textPositionH1=max(Hlevel(:,1); textPositionH2=max(Hlevel(:,2); H2Steady=Hlevel(size(Hlevel(:,1),1),1)*ones(size(y1); xmax=max(0:T:tend); xmin=0; ymax=110; ymin=50; scrsz = get(0,ScreenSize); gca=figure(Position,5 10 scrsz(3)-10 scrsz(4)-90) %gca=figure(Position,5 10 scrsz(

25、3)/2 scrsz(4)/1.5) set(gca,Color,w); plot(0:T:tend,Hlevel(:,1),r,LineWidth,2) hold on plot(0:T:tend,Hlevel(:,2),b,LineWidth,2) hold on plot(0:T:tend,yv,k,LineWidth,2) hold on plot(0:T:tend,y2sp,g,LineWidth,2) hold on plot(0:T:tend,H2Steady,y,LineWidth,2) line(tend/2 tend/2+27,(ymax-ymin)/2+ymin-(yma

26、x-ymin)/10 (ymax-ymin)/2+ymin-(ymax-ymin)/10,Color,r,LineWidth,6) text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/10, 第1個(gè)水箱的液位H1,FontSize,16) line(tend/2 tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/6 (ymax-ymin)/2+ymin-(ymax-ymin)/6,Color,b,LineWidth,6) text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/6, 第2

27、個(gè)水箱的液位H2,FontSize,16) line(tend/2 tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/4.2 (ymax-ymin)/2+ymin-(ymax-ymin)/4.2,Color,g,LineWidth,6) text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/4.2, 第2個(gè)水箱的液位給定值,FontSize,16) line(tend/2 tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/3.2 (ymax-ymin)/2+ymin-(ymax-ymin)/3

28、.2,Color,k,LineWidth,6) text(tend/2+27,(ymax-ymin)/2+ymin-(ymax-ymin)/3.2, 閥位變化情況,FontSize,16) axis(xmin xmax ymin ymax); text(tend/5,ymax+1.5, 實(shí)驗(yàn)三 考慮閥位飽和特性時(shí)的控制效果,FontSize,22) grid %hold off仿真實(shí)驗(yàn)曲線：三 思考題與實(shí)驗(yàn)三相比，考慮調(diào)節(jié)閥飽和特性前后，響應(yīng)有何不同？答：與實(shí)驗(yàn)三相比，調(diào)節(jié)飽和特性閥前后調(diào)進(jìn)行對(duì)比，H1、H2的液位在考慮飽和特性之后，響應(yīng)曲線比不考慮的時(shí)候略微平緩一些，閥位的調(diào)節(jié)沒有出現(xiàn)超過1

29、00%的情況。四 實(shí)驗(yàn)總結(jié)通過本次實(shí)驗(yàn)，我認(rèn)識(shí)了面向結(jié)構(gòu)圖的離散相似法的原理和程序結(jié)構(gòu)，掌握了理含有非線性環(huán)節(jié)的控制系統(tǒng)的仿真方法。實(shí)驗(yàn)三 采樣系統(tǒng)的仿真一 實(shí)驗(yàn)?zāi)康? 掌握理解數(shù)字控制系統(tǒng)的仿真技術(shù)。2 掌握理解增量式 PID 數(shù)字控制器的實(shí)現(xiàn)方法。二 實(shí)驗(yàn)內(nèi)容 1根據(jù)上面的各式，編寫仿真程序。仿真程序：clcclear allclose all;A=2;ku=0.1/0.5;H10=1.5;H20=1.4;alpha12 = 0.25/sqrt(H10);alpha2 = 0.25/sqrt(H20);R12=2*sqrt(H10)/alpha12;R2=2*sqrt(H20)/alpha

30、2;H1SpanLo=0;H2SpanLo=0;H1SpanHi=2.52;H2SpanHi=2.52;Kp=1.89;Ti=30;Td=10;ad = 1/(A*R12);a1 = 1/(A*R12);a2 = 1/(A*R2);Kc=Kp/Ti;bc=Ti;Kd = 1/A;K1 = ku/A;K2 = 1/(A*R12);beta1=1/(a1*a2);beta3=beta1*a1/(a2-a1);beta2=-beta1-beta3;u(1)=0;u(2)=0;u(3)=0;y(1)=0;y(2)=0;y(3)=0;nCounter = 70;T=10;k=2;deltuU=0;e(1

31、)=0;e(2)=0;e(3)=0;% uc(1)=0;H20_percent=(H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100;H2set_percent=80;tend = nCounter*T;for t=2*T:T:tend k=k+1; e(k)=(H2set_percent - (y(k-1)+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100)/100; deltaU=Kp*(e(k)-e(k-1)+Kp*T/Ti*e(k)+Kp*Td/T*e(k)-2*e(k-1)+e(k-2); e(k-2)=e(k-1); e(k-1)=

32、e(k); u(k)=u(k-1)+deltaU; y(k)=(exp(-a1*T)+exp(-a2*T)*y(k-1) - exp(-(a1+a2)*T)*y(k-2)+K1*K2*(beta1+beta2+beta3)*u(k)- . K1*K2*(beta1*(exp(-a1*T)+exp(-a2*T)+beta2*(1+exp(-a2*T) . +beta3*(1+exp(-a1*T)*u(k-1)+K1*K2*(beta1*exp(-(a1+a2)*T)+beta2*exp(-a2*T)+ . beta3*exp(-a1*T)*u(k-2) ; y(k-2)=y(k-1); y(k-1)=y(k); u(k-2)=u(k-1); u(k-1)=u(k); endy=(y+H20-H2SpanLo)/(H2SpanHi-H2SpanLo)*100;y2sp=H2set_percent*ones(size(y);u=(u+0.5)*100;% textPositionH1=max(Hlevel(:,1);% textPositionH2=max(Hlevel(:,2);% H2Steady=Hlevel(size(Hlevel(:,1),1),1)*ones(size(y1); xmax=max(0:T:tend); xmin=0; ymax