直線(xiàn)一級(jí)倒立擺系統(tǒng)的建模及仿真

計(jì)算機(jī)控制技術(shù)課程設(shè)計(jì)實(shí)驗(yàn)：直線(xiàn)一級(jí)倒立擺系統(tǒng)的建模及仿真、已知條件:擺桿角度為輸出，小車(chē)的位移為輸入。導(dǎo)軌中點(diǎn)為坐標(biāo)軸的中心即零點(diǎn)，右向?yàn)樽鴺?biāo)值增加的方向，桿偏移其瞬時(shí)平衡位置右側(cè)的角度為正值。二、任務(wù)要求：總體任務(wù)通過(guò)調(diào)節(jié)PID參數(shù)，設(shè)計(jì)PID控制器實(shí)現(xiàn)擺桿在受到干擾的情況下，依然能恢復(fù)平衡。具體包括以下幾部分：理論推導(dǎo)包括倒立擺系統(tǒng)的動(dòng)力學(xué)模型，傳遞函數(shù)，離散傳遞函數(shù)，狀態(tài)空間或差分方程，穩(wěn)定性分析，PID控制器設(shè)計(jì)程序?qū)崿F(xiàn)實(shí)現(xiàn)內(nèi)容：倒立擺系統(tǒng)模型，控制器以及仿真結(jié)果的顯示。開(kāi)發(fā)語(yǔ)言和工具：Matlabm文件或C/C++（工具：VC++或其它）PID控制參數(shù)設(shè)定及仿真結(jié)果。分別列出不同桿長(zhǎng)的仿真結(jié)果（例如：L=0.25和L=0.5）。將理論推導(dǎo)、程序?qū)崿F(xiàn)、仿真結(jié)果寫(xiě)成實(shí)驗(yàn)報(bào)告。具體求解過(guò)程如下：一，倒立擺系統(tǒng)動(dòng)力學(xué)模型的建立1）圖1擺桿的受力分析圖以擺桿為研究對(duì)象，對(duì)其進(jìn)行受力分析，如圖11）FxCP+mgxCP-M=0IC I式中，F(xiàn)為桿的慣性力，表達(dá)式為F=ma=mCa+at+an)=F+Ft+FnIC IC C PCPCPIPicpicp，m為桿的質(zhì)量，g為重力加速度，M為桿的慣性力偶。慣性力及慣性力偶的大小分別為d2x d20 1 d202）F=ma=m,Ft=ma=m,M=Ja=_mL22）IPPdt2ICP dt2Ic3 dt2式中，a為桿的角加速度，a為小車(chē)的加速度，2L為桿的長(zhǎng)度，0為桿偏離中P心位置的角度，X偏離軌道中心的位移。對(duì)（2）式代入（1）式，并整理可得TOC\o"1-5"\h\zL^0-gsin0=-cos0 （3）3dt2 dt2當(dāng)擺動(dòng)較小時(shí)，可以進(jìn)行近似處理sin0q0,cos0沁1。故（3）式可化為4d20 d2xL -g0=- （4）3dt2 dt2對(duì)（4）式進(jìn)行拉普拉斯變換得3Ls20（s）-g0（s）=-s2X（s） （5）

則系統(tǒng)的開(kāi)環(huán)傳遞函數(shù)為()0(s) -s2 -as2GW_ _ _—X(s)4Ls2-g s2-ag3式中，a_-4L二，離散傳遞函數(shù)，差分方程、狀態(tài)空間及穩(wěn)定性分析1，離散傳遞函數(shù)將連續(xù)系統(tǒng)離散化，根據(jù)連續(xù)傳遞函數(shù)G(s)可求得相應(yīng)的脈沖傳遞函數(shù)為1-bz1-bz-i+z-27)1-e-TsG(s)_(1-z-1)Z-asss2-agG(z)=Z式中，b=eagT+e—agT。將參數(shù)L二1,g二10,T二0.01代入(7)式得該參數(shù)下的脈沖傳遞函數(shù)為G(zG(z)_—0.75+1.50z-i-0.7503z-1-2.001z-i+z-28)2，根據(jù)離散傳遞函數(shù)求系統(tǒng)差分方程由陽(yáng))?(z)-a+(a+ab/2)z-1-abz-2/2Glz丿_ _1-bz-1+z-1-bz-1+z-2可得z-1x(z)-吳-2x(z)弘)-bz啟)+z-1x(z)-吳-2x(z)進(jìn)行反變換即可得到對(duì)應(yīng)的差分方程6(k)6(k)-b6(k-1)+e(k-2)=一ax(k)'ab'

+a+——xI2丿(k-])一#x(k-2)10)將參數(shù)L二1,g二10,T二0.01代入(10)式得該參數(shù)下的差分方程為9(k)—2.0010(k-l)+e(k-2)=—0.75x(k)+1.50x(k-1)—0.7503x(k-2)(11)3，根據(jù)離散傳遞函數(shù)求系統(tǒng)狀態(tài)空間表達(dá)式根據(jù)G(z)=?(zG(z)=?(z)X(z)(a+ab/2)z-1-abz-2/21-bz-1+z-2=—a+'ab) 7a———z—1—bz—2I2丿1—bz—1+z—2設(shè)W(z)= (12)1—bz—1+z—2W(z)=bz-1W(z)-z-2W(z)+X(z)?(z)=-aX(z)+a—?dú)q丿z-1W(z)+”—字丿丿—2?(z)=-aX(z)+a—?dú)q丿z-1W(z)+”—字丿丿—2W(z)13)選取狀態(tài)變量X(z)=z-1W(z),X(z)=z-2W(z)=z-1X(z)12114)將W(z)代入X(z),X(z)，取z反變換，可得狀態(tài)方程12x(k+1)x(k+1)2rx(汀Tx1(k)2+_0_-10x(k)15)(k)=-ax(k)x(k)1x(k)2將參數(shù)L=1,g=10,T=0.01代入(15)(k+1)「=(k+1)_=0(k)=-0.75x(k)-0.0003x(k)-0.0003x(k)124，穩(wěn)定性分析對(duì)于開(kāi)環(huán)系統(tǒng)，由傳遞函數(shù)可得系統(tǒng)的特征方程為A(z)=z2-2.001z+1式得該參數(shù)下的狀態(tài)方程為x1x22.00081-10x(k)

x1(k)2x(k)16)17)特征方程的根為z=1.032,z=0.969。12由于特征根中有一個(gè)大于1，位于單位圓外，故系統(tǒng)是不穩(wěn)定的。四，PID控制器設(shè)計(jì)由以上分析可知，系統(tǒng)是不穩(wěn)定的，為使倒立擺在受到干擾時(shí)能保持穩(wěn)定必須對(duì)系統(tǒng)進(jìn)行PID控制器的設(shè)計(jì)。

倒立擺計(jì)算機(jī)控制系統(tǒng)的框圖如下倒立擺控制系統(tǒng)結(jié)構(gòu)圖模擬PID控制器的基本算式為x(k)=Ke(k)+K丈e(j)+Ke(k)-e(k-1) (18)PIDj=0根據(jù)式(16)和(18)可求得在不同K,K,K下的擺角值。具體的仿真過(guò)

PID程如下：1,取K=-2.0,K=-0.3固定不變，K分別取0.5,1.0,1.5,2.0,倒立擺的輸出變I D P化曲線(xiàn)如下：KP=0.5,KI=-2.0,KD=-0.3(a1)KP=1.0,KI=-2.0,KD=-0.3t/s(b1)KP=1.5,KI=-2.0,KD=-0.3t/s(c1)由圖（al）、（bl）、（cl）、（dl）可知，當(dāng)K=0.5時(shí)，倒立擺在受到干擾時(shí)可以P達(dá)到穩(wěn)定；繼續(xù)增大K到1.0,系統(tǒng)的調(diào)節(jié)時(shí)間變短；當(dāng)K=1.5時(shí)，雖然倒立PP擺仍能穩(wěn)定，但調(diào)節(jié)時(shí)間變長(zhǎng)；當(dāng)K=2.0時(shí)，倒立擺已無(wú)法恢復(fù)穩(wěn)定。P2,取K=1.0,K=—0.3固定不變，K分別取-0.5,-1.5,-2.5,-3.5，倒立擺的輸出變PDI化曲線(xiàn)如下：KP=1.0,KI=-0.5,KD=-0.3(a2)KP=1.0,KI=-1.5,KD=-0.3t/s(b2)KP=1.0,KI=-2.5,KD=-0.3t/s(c2)（d2）由圖（a2）、（b2）、（c2）、（d2）可知，當(dāng)K=-0.5時(shí)，倒立擺在受到干擾時(shí)開(kāi)始I可以達(dá)到穩(wěn)定，但當(dāng)時(shí)間變長(zhǎng)后無(wú)法達(dá)到穩(wěn)定；繼續(xù)減小K到-1.5，系統(tǒng)可以達(dá)I到穩(wěn)定；當(dāng)K=-2.5時(shí)，系統(tǒng)的調(diào)節(jié)時(shí)間變短；當(dāng)K減小到-3.5時(shí)，倒立擺已II無(wú)法恢復(fù)穩(wěn)定。3,取K=1.0,K=-2.0固定不變，K分別取-0丄-0.3,-0.5,-0.7，倒立擺的輸出變PID化,曲線(xiàn)如下：KP=1.0,KI=-2.0,KD=-0.1(a3)KP=1.0,KI=-2.0,KD=-0.3t/s(b3)KP=1.0,KI=-2.0,KD=-0.5t/s(c3)（d3）由圖（a3）、（b3）、（c3）、（d3）可知，當(dāng)K從-0.1減小到-0.7時(shí)，倒立擺在受到I干擾時(shí)達(dá)到穩(wěn)定的調(diào)節(jié)時(shí)間先減小后增大，到-0.7時(shí)系統(tǒng)已無(wú)法恢復(fù)穩(wěn)定。五，改變桿長(zhǎng)后的仿真結(jié)果將桿長(zhǎng)由L=1改為L(zhǎng)=0.5，系統(tǒng)的仿真結(jié)果如下：1,取K=—1.2,K=—0.1固定不變，K分別取020.4,0.8,1.0，倒立擺的輸出變IDP化曲線(xiàn)如下：KP=0.2,KI=-1.2,KD=-0.1(a11)KP=0.4,KI=-1.2,KD=-0.1(b11)KP=0.8,KI=-1.2,KD=-0.1(c11)012345678910t/s（d11）由圖（all）、（bll）、（ell）、（dll）可知，當(dāng)K=0.2時(shí)，倒立擺在受到干擾時(shí)可P以達(dá)到穩(wěn)定；繼續(xù)增大K到0.4,系統(tǒng)的調(diào)節(jié)時(shí)間變短；當(dāng)K=0.8時(shí)，雖然倒PP立擺仍能穩(wěn)定，但調(diào)節(jié)時(shí)間變長(zhǎng)；當(dāng)K=l.0時(shí)，倒立擺已無(wú)法恢復(fù)穩(wěn)定。P2,取K=0.4,K=-0.l固定不變，K分別取-l.0,-l.3,-l.6,-l.9,倒立擺的輸出變PDI化曲線(xiàn)如下：KP=0.4,KI=-1.0,KD=-0.1(a12)KP=0.4,KI=-1.3,KD=-0.1(b12)KP=0.4,KI=-1.6,KD=-0.1(c12)（d12）由圖（al2）、（bl2）、（cl2）、（dl2）可知，當(dāng)K=—1.0時(shí)，倒立擺在受到干擾時(shí)I開(kāi)始可以達(dá)到穩(wěn)定，但當(dāng)時(shí)間變長(zhǎng)時(shí)會(huì)失去穩(wěn)定；繼續(xù)減小K到-1.3，系統(tǒng)可以I達(dá)到穩(wěn)定；當(dāng)K=—1.6時(shí)系統(tǒng)可以穩(wěn)定但調(diào)節(jié)時(shí)間變長(zhǎng)；當(dāng)K減小到-1.9時(shí)，II倒立擺已無(wú)法恢復(fù)穩(wěn)定。3，取K=0.4,K=—1.2固定不變，K分別取-0.01,-0.1,-0.2,-0.24，倒立擺的輸出PID變化,曲線(xiàn)如下：KP=0.4,KI=-1.2,KD=-0.01(a13)KP=0.4,KI=-1.2,KD=-0.1(b13)KP=0.4,KI=-1.2,KD=-0.2(c13)(d13)由圖(a13)、(bl3)、(cl3)、(dl3)可知，當(dāng)k【從-0.01減小到-0.24時(shí)，倒立擺在受到干擾時(shí)達(dá)到穩(wěn)定的調(diào)節(jié)時(shí)間先減小后增大，到-0.24時(shí)系統(tǒng)已無(wú)法恢復(fù)穩(wěn)定。附錄：1，PID設(shè)計(jì)過(guò)程的Matlab程序如下l=1;g0=l0;%重力加速度tt=0.0l;%采樣周期a=3/(4*1);%代數(shù)anumGs=[-a00];%連續(xù)傳遞函數(shù)的分子denGs=[l0-a*g0];%連續(xù)傳遞函數(shù)的分母sysGs=tf(numGs,denGs);%系統(tǒng)連續(xù)傳遞函數(shù)[F,G,C,D]=ssdata(c2d(sysGs,tt,'zoh'));%求得ZOH下以tt為采樣周期的狀態(tài)空間(離散)t(l)=0;x=zeros(l,l000);%一行1000列0矩陣theta(l)=0;e=zeros(l,l000);%—行1000列0矩陣s=0;KP=0.5;KI=-2.0;KD=-0.3;%設(shè)置PID參數(shù)y=[0;0];fork=2:l000;t(k)=t(k-1)+tt;ifk==20delta=0.1;elsedelta=0;endtheta(k-1)=theta(k-1)+delta;e(k)=-theta(k-l);%PID公式的第一項(xiàng)s=s+e(k);%PID公式的第二項(xiàng)x(k)=KP*e(k)+KI*s+KD*(e(k)-e(k-1));%PID公式theta(k)=C*y+D*x(k);%初始輸出y為零y=F*y+G*x(k);endfigure(1)plot(t,theta),gridtitle('KP=0.5,KI=-2.0,KD=-0.3')xlabel('t/s'),ylabel('\theta')2,為便于改變參數(shù)，可利用Matlab的GUI界面來(lái)方便的進(jìn)行PID參數(shù)的設(shè)計(jì)，所設(shè)計(jì)的界面如下圖所示：此時(shí)需將程序命令添加到Plot按鈕下，即functionplot_push_Callback(hObject,eventdata,handles)l=str2double(get(handles.L_edit,'String'));tt=str2double(get(handles.T_edit,'String'));g0=str2double(get(handles.G_edit,'String'));KP=str2double(get(handles.KP_edit,'String'));KI=str2double(get(handles.KI_edit,'String'));KD=str2double(get(handles.KD_edit,'String'));a=3/(4*l);numGs=[-a00];denGs=[10-a*g0];sysGs=tf(numGs,denGs);[F,G,C,D]=ssdata(c2d(sysGs,tt,'zoh'));t(1)=0;x=zeros(1,1000);theta(1)=0;e=zeros(1,1000);s=0;y=[0;0];fork=2:1000;t(k)=t(k-1)+tt;ifk==20delta=0.1