倒立擺系統(tǒng)實驗設計報告

1、學生實驗報告課程名稱：倒立擺系統(tǒng)課程設計組號：7姓名：學號：郵箱：2010年11月11 1日目錄倒立擺系統(tǒng)的構(gòu)成 3單級倒立擺數(shù)學模型的建立 3傳遞函數(shù) 6狀態(tài)空間方程 6系統(tǒng) MATLAB 仿真和開環(huán)響應 7穩(wěn)定性與可控性分析 11控制器設計 12基于狀態(tài)反饋的控制算法設計與仿真 LQR 12極點配置法 16PID控制算法19實驗結(jié)果及與仿真結(jié)果的對比分析 29感想和建議 30倒立擺系統(tǒng)的構(gòu)成圖1倒立擺系統(tǒng)的組成框圖如圖1所示為倒立擺的結(jié)構(gòu)圖。系統(tǒng)包括計算機、 運動控制卡、伺服機構(gòu)、 倒立擺本體和光電碼盤幾大部分， 組成了一個閉環(huán)系統(tǒng)。 光電碼盤1將小車的位移、速度信號反饋給伺 服驅(qū)動器和運

2、動控制卡，擺桿的位置、速度信號由光電碼盤 2反饋回控制卡。計算機從運動 控制卡中讀取實時數(shù)據(jù)，確定控制決策（小車向哪個方向移動、移動速度、加速度等），并由運動控制卡來實現(xiàn)該控制決策，產(chǎn)生相應的控制量，使電機轉(zhuǎn)動，帶動小車運動，保持擺 桿平衡。單級倒立擺數(shù)學模型的建立在忽略了空氣流動，各種摩擦之后，可將倒立擺系統(tǒng)抽象成小車和勻質(zhì)桿組成的系統(tǒng), 如下圖2所示M小車質(zhì)量（本實驗系統(tǒng)0.5 Kg)m擺桿質(zhì)量（本實驗系統(tǒng)0.2 Kg)b小車摩擦系數(shù)（本實驗系統(tǒng)0.1 N/m/sec)圖2單級倒立擺模型示意圖那我們在本實驗中定義如下變量：I擺桿轉(zhuǎn)動軸心到桿質(zhì)心的長度（0.3 m）I擺桿慣量（0.006 k

3、g*m*m ）F 加在小車上的力x小車位置0擺桿與垂直向上方向的夾角0擺桿與垂直向下方向的夾角（考慮到擺桿初始位置為豎直向下）下面我們對這個系統(tǒng)作一下受力分析。下圖3是系統(tǒng)中小車和擺桿的受力分析圖。其中,N和P為小車與擺桿相互作用力的水平和垂直方向的分量。因而矢量方向定注意：在實際倒立擺系統(tǒng)中檢測和執(zhí)行裝置的正負方向已經(jīng)完全確定, 義如圖，圖示方向為矢量正方向。-blRicticr圖3倒立擺模型受力分析分析小車水平方向所受的合力，可以得到等式：應用Newton方法來建立系統(tǒng)的動力學方程過程如下：分析小車水平方向所受的合力，可以得到以下方程:Mx = F - bx - N由擺桿水平方向的受力進行

4、分析可以得到下面等式:N 訓寫(x lsin)dtN =mx - ml jcos v _ml / sin把這個等式代入上式中，就得到系統(tǒng)的第一個運動方程:(1)可以得到下 a anO(M m)x bx ml vcos J - mh sin 門-F為了推出系統(tǒng)的第二個運動方程，我們對擺桿垂直方向上的合力進行分析, 面方程:d 2P - mg = m 2(l cos即：P - mg - -ml T sin - ml T2 cos力矩平衡方程如下:-Pl sin v - Nl cos）- I 二注意：此方程中力矩的方向，由于 v -二 ： cos 二-cost, sin 二-sin二，故等式前面有負

5、號。合并這兩個方程，約去 P和N，由|二舟ml2得到第二個運動方程：4 ml J mgls in - -mlxcos（ 2）設（ 是擺桿與垂直向上方向之間的夾角），假設與1 （單位是弧度）相比很小，即1,則可以進行近似處理：cost - -1，si nV -，（）0。用u來代dt表被控對象的輸入力 F，線性化后兩個運動方程如下：(s)s2 - g：(s)二 X(s)s23(M m)X(s)s2 bX(s)s-ml：(s)s2 =U(s)注意：推導傳遞函數(shù)時假設初始條件為0。由于輸出為角度，求解方程組(4)的第一個方程，可以得到4 g _ X(s) = l - 篤G(s)3 s把上式代入方程組(

6、4)的第二個方程，得到(M m) 山:(s)s2 b 山耳 G(s)s_ml：(s)s2 =U(s)ml smls2整理后得到傳遞函數(shù)：ml 2 sU(s)二q41 bml2 3 (M m)mgl 2 bmglssssqqq其中q=(M m)(I ml2)(ml)2狀態(tài)空間方程系統(tǒng)狀態(tài)空間方程為X 二 AX Buy 二 CX Du方程組(3)對x解代數(shù)方程，得到解如下:-4b(4M m)3mg(4M m)4(4M m)也x3四3 u(4M m)l (4M m)l (4M m)l整理后得到系統(tǒng)狀態(tài)空間方程:x00100_j-0 1-4b3mg0(4M +m)(4M +m)+(4M + m)001

7、CU-3b3g(M +m)0(4M +m)l(4M +m)l1(4M +m)l _0 x0 *0u一0系統(tǒng)Matlab仿真和開環(huán)響應實際系統(tǒng)參數(shù)如下，求系統(tǒng)的傳遞函數(shù)、Matlab 仿真。M小車質(zhì)量m擺桿質(zhì)量b 小車摩擦系數(shù)l擺桿轉(zhuǎn)動軸心到桿質(zhì)心的長度I擺桿慣量F加在小車上的力x小車位置0擺桿與垂直方向的夾角T采樣頻率狀態(tài)空間方程，并進行脈沖響應和階躍響應的1.096 Kg0.109 Kg0 .1N/m/sec0.2 5m0.0034 kg*m*m0.005 秒注意：在進行實際系統(tǒng)的Matlab仿真時，請將采樣頻率改為實際系統(tǒng)的采樣頻率傳遞函數(shù):在 Matlab 中，拉普拉斯變換后得到的傳遞函

8、數(shù)可以通過計算并輸入分子和分母矩陣來實現(xiàn)。求系統(tǒng)傳遞函數(shù)的m-文件內(nèi)容如下：M = 1.096;m = 0.109;b = 0.1;I= 0.0034;g = 9.8;l = 0.25;q = (M+m)*(l+m*L2)-(m*l)A2;%simplifies in putnum = m*l/q 0den = 1 b*(I+m*lA2)/q -(M+m)*m*g*l/q -b*m*g*l/qt = 0 : 0.01 : 5;impulse ( num , den , t )axis ( 0 1.1 0 70 )執(zhí)行上面的文件，就可以求出系統(tǒng)傳遞函數(shù)的分子與分母多項式的Matlab 表示：nu

9、m =2.3566 0den =1.0000 0.0883 -27.8285 -2.3094可 以 得 到 系 統(tǒng) 開 環(huán) 脈 沖 響 應 的 曲 線 如 下 ：706050Impulse Resp onseedLIP m40302010000.20.40.6Time （sec）圖4系統(tǒng)開環(huán)脈沖響應曲線 狀態(tài)空間法：0.8 1狀態(tài)空間法可以進行單輸入多輸出系統(tǒng)設計，（從實驗二開始，我們將嘗試同時對擺桿角度和小車位置進行控制）。為了更具挑戰(zhàn)性，給小車一個階躍輸入信號，設計指標如下：小車位置X和擺桿角度位置X的上升時間小于擺桿角度的超調(diào)量小于0的穩(wěn)定時間小于5秒;0.5 秒；20度（0.35弧度）

10、。下面，我們用 Matlab求出系統(tǒng)的狀態(tài)空間方程各矩陣，并仿真系統(tǒng)的開環(huán)階躍響應。在這里同樣給出了一個 m-文件，執(zhí)行這個文件，Matlab將會給出系統(tǒng)狀態(tài)空間方程的A,B, C和D矩陣，并可以繪出在給定輸入為一個0.2 m的階躍信號時系統(tǒng)的響應曲線。M = 1.096; m = 0.109; b = 0.1;I= 0.0034; g = 9.8;l = 0.25;p = l*(M+m)+M*m*F2;%de nomin ator for the A and B matricies0;A = 010 -(l+m*L2)*b/p (mA2*g*lA2)/p0;0 0 0 1;0 -(m*l*b

11、)/pm*g*l*(M+m)/p 0B = 0;(l+m*lA2)/p;0;m*l/pC = 1 0 0 0;0 0 1 0D = 0;0T=0:0.005:10;U=0.2*ones(size(T);Y,X=lsim(A,B,C,D,U,T);plot(T,Y)axis(0 2.5 0 100)執(zhí)行該 m 文件，可以求出系統(tǒng)的狀態(tài)空間A、B、C、D 矩陣 ,得到開環(huán)系統(tǒng)階躍響應的曲線。A =01.0000000-0.08830.629300001.00000-0.235727.82850B =00.883202.3566C =10 0000 10D =10090807060504030201

12、0000.511.52.5圖5倒立擺狀態(tài)空間開環(huán)系統(tǒng)階躍響應曲線圖中，實線是擺桿角度響應曲線，虛線是小車位置響應曲線。穩(wěn)定性與可控性分析我們先看一看系統(tǒng)的穩(wěn)定性，將數(shù)據(jù)代入狀態(tài)方程中，利用matlab程序可以求出系統(tǒng)的零極點。源代碼如下：M = 1.096;m = 0.109;b = 0.1;I= 0.0034;l = 0.25;a b c d=wer_ss(M,m,b,l);%自己編寫的函數(shù)，建立模型之用，具體程序見下面sysc=ss(a,b,c,d);sysd=c2d(sysc,0.005);da db dc dd=ssdata(sysd);z p gain=ss2zp(da,db,dc,

13、dd,1)z =-0.9999-0.99991.0000 + 0.0000i1.0000 - 0.0000i1.02750.97331.00000.99961.02850.9723gain =1.0e-004 *0.11130.3338wer_ss源程序：fun cti ona b c d=wer_ss(M,m,b,l)0;00 0 1;0-3*b/(4*M+m)*l)a=01 0 0;0 -4*b/(4*M+m)3*m*9.8/(4*M+m)3*9.8*(M+m)/(4*M+m)*l) 0;b=0;4/(4*M+m);0;3/(4*M+m)*l);c=1 0 0 0;0 0 1 0;d=0;

14、0由得到的p (極點)可知，有的極點在單位圓外，所以可知原系統(tǒng)是不穩(wěn)定的。 同樣，我們可以利用 matlab來得到系統(tǒng)的能控性，源代碼如下：ud=ctrb(da,db);ran k(ud)ans =4由得到的rank( ud)的值可知，原系統(tǒng)的能控性矩陣為4，所以我們可知原系統(tǒng)是能控的?？刂破髟O計基于狀態(tài)反饋的控制算法設計與仿真LQR由理論分析知，可以設計基于最優(yōu)控制的狀態(tài)調(diào)節(jié)器，使系統(tǒng)閉環(huán)穩(wěn)定。 設狀態(tài)反饋調(diào)節(jié)律的形式為u(k)二-KrX(k)通過使性能指標函數(shù)1 1 N丄J Xt(N)PX(N)XT(k)Qr(k)X(k) uT (k)Rr (k)u(k)2 2心為最小，根據(jù)在附錄1中我們

15、所介紹的求得Kr -Rr HtP(K 1)HHTP(k 1)G其中p由下列黎卡提方程獲得P(k) =Qr GTP(k 1)G _GTP(k 1)HRr H Tp(k 1)HHTP(k 1)G其中Qr , Rr分別用來對狀態(tài)向量X(k)，控制向量U(k)引起的性能度量的相對重要性進行加權(quán)。在實際運算中我們運用Matlab控制系統(tǒng)工具箱中的dlqr”函數(shù)直接進行運算。利用dlqr函數(shù)，我們需要提供兩個權(quán)值矩陣：Q、R。通常我們?nèi)=1，而對于Q我們只能通過不斷的湊取來得到。源代碼如下： R=1; Q=10 0 0 0;0 0 0 0;0 0 1 0;0 0 0 0Q =10 0000 0000 0

16、100 000 T=0.005; syse k=wer_lqr(da,db,dc,dd,Q,R,T);%wer_lqr 是自己定義的函數(shù)，具體見下面程序 x0=0.05;0;0.0175;0; t=0:0.005:10; y x1=in itial(syse,x0,t);plot(t,y(:,1),red,t,y(:,2),blue)wer_lqr源程序：fun cti onsysresult k=wer_lqr(da,db,dc,dd,Q,R,T);%sysresult k=wer_lqr(da,db,dc,dd,Q,R,T);k S e=dlqr(da,db,Q,R);G=da-db*k;

17、sysresult=ss(G,db,dc,dd,T);我們已開始的 Q為：Q1=10 0 0 0;0 0 0 0;0 0 1 0;0 0 0 0；結(jié)果得到的圖為圖6；Q2取為：Q=100 0 0 0;0 0 0 0;0 0 10 0;0 0 0 0;結(jié)果得至曲勺圖為圖 7 ；通過比較，我們發(fā)現(xiàn)當Q11、Q33比值一定時，取大的值時系統(tǒng)的響應速度加快，但是超調(diào)加大；反之則響應變慢但是超調(diào)減小。0.03 -小車位移-0.01 Lc1J012345678910圖6 Q1響應圖0 12345678910圖7 Q2響應圖 在左右權(quán)衡之間，我們最終選取了：Q=300 0 0 0;0 0 0 0;0 0 3

18、0 0;0 0 0 0此時的響應曲線如圖8，而k值為：-16.6147-12.422656.590910.24440.050.040.03小車位移0.020.01-0.01 匚08910圖8最優(yōu)的響應曲線此時的單位階躍響應曲線為圖9:0.020.010-0.01-0.02-0.03-0.04-0.05-0.06-0.0710擺角小車位移圖9單位階躍響應曲線 從仿真效果來看，零狀態(tài)響應和單位階躍響應都符合要求。極點配置法采用極點配置法設計多輸出的倒立擺系統(tǒng)的控制方案。 制擺桿和小車的位置。圖10是控制系統(tǒng)的示意圖?？梢杂猛耆珷顟B(tài)反饋來解決,圖10控制系統(tǒng)框圖假定所有的狀態(tài)變量都可以測量和反饋，

19、可以證明：若所研究的系統(tǒng)是狀態(tài)完全可控的， 那么，利用狀態(tài)反饋的方法， 經(jīng)過適當?shù)臓顟B(tài)反饋增益矩陣， 就可以把閉環(huán)系統(tǒng)的極點配置 到任何期望的位置。設開環(huán)控制系統(tǒng)的離散狀態(tài)方程為：x(k+1) = Gx(k) + Hu(k)其中，假設系統(tǒng)是狀態(tài)完全可控的x(k)為在第k次采樣時刻的狀態(tài)矢量(n維矢量)u(k)為在第k次采樣時刻的控制信號(標量)G = n n矩陣H = n 1矩陣設極點配置的控制律形式為u(k) KpX(k)式中Kp是狀態(tài)反饋增益矩陣(1 n矩陣)，于是該系統(tǒng)就成為一個閉環(huán)控制系統(tǒng)。其閉環(huán)狀態(tài)方程為gXdHKpXk)注意，G 一 HK p的特征值就是所要求的閉環(huán)極點我們希望利用

20、狀態(tài)反饋u(k) = -K pX(k)把閉環(huán)極點布置在z八i，z八2,Z n。即要求特征方程為：zl _G +HK p =(z _ 氣)(z _巴)(Z _)=zn a1zn4 a2znJ 梟 ;: an4z an = 0根據(jù)Cayley_hamiton定理，經(jīng)過推導(此略)可以得到Kp = 0 00 1H GHGn4H廣(G)其中上式給出所要求的狀態(tài)反饋增益矩陣 曼公式。狀態(tài)反饋增益矩陣按這樣的方法確定，降到零。注意，對于一個給定的系統(tǒng)，矩陣Kp。矩陣Kp的這種特殊表達式就是常說的阿克即可使誤差（由擾動所引起的） 以足夠快的速度K p并不是唯一的，而是取決于所期望的閉環(huán)極點位置（它決定響應速

21、度）的選擇。選擇期望的閉環(huán)極點或期望的特征方程是在誤差矢量響 應的快速性與對擾動和測量噪聲敏感型之間的一個折衷方案。也就是說，如果我們使誤差響應的速度提高，那么擾動和測量噪聲的有害影響往往也會增強。在確定給定系統(tǒng)的狀態(tài)反饋增益矩陣Kp時，通常是通過比較按不同的期望閉環(huán)極點或期望特征方程得到的矩陣Kp，K并從中選出使整個系統(tǒng)達到的特性最好的那個矩陣p。在實際設計時，我們運用Matlab控制系統(tǒng)工具箱中的“ place ”函數(shù)直接進行仿真和運算。先在連續(xù)域中進行計算，然后再轉(zhuǎn)到離散域。根據(jù)系統(tǒng)的性能要求，我們可取 =0.707，Wn =3，這時完全滿足題中的性能要求。那么我們就可以進行離散域設計了

22、: z仁exp(-3*2A0.5)/2+(-3*2A0.5)*j/2)*0.005) z1 =0.9894 - 0.0105i z2=exp(-3*2A0.5)/2-(-3*2A0.5)*j/2)*0.005) z2 =0.9894 + 0.0105i z3=exp(-10*0.005)z3 =0.9512 z4=exp(-12*0.005)z4 =0.9418 p=z1 z2 z3 z4;K=place(da,db,p)K =-38.6579-25.5096103.324717.9041 G=da-db*K; syse=ss(G,db,dc,dd,0.005); t=0:0.005:10;

23、x0=0.05; 0; 0.0175; 0; y1,x=in itial(syse,x0,t); plot(t,y1(:,1),red,t,y1(:,2),blue)0.060.050.040.03小車位移0.020.01-0.01 L0910圖11極點配置圖零輸入響應而它的單位階躍響應是:u=on es(1,le ngth(t); y,x=lsim(syse,u,t)0.0150.010.005擺角0-0.005-0.01-0.015-0.02小車位移-0.025-0.03圖12極點配置單位階躍響應 從仿真效果來看，也是基本上達到了系統(tǒng)的要求。PID控制算法目的：設計PID控制器，使得當在小

24、車上施加1N的脈沖信號時，閉環(huán)系統(tǒng)的響應指標為:1、穩(wěn)定時間小于5秒2、 穩(wěn)態(tài)時擺桿與垂直方向的夾角變化小于0.1弧度分析：系統(tǒng)輸出量為擺桿的位置，它的初始位置為垂直向上，我們給系統(tǒng)施加一個擾動，觀察擺桿的響應。系統(tǒng)框圖如下：ControllerKD(j)圖13 PID系統(tǒng)框圖圖中KD (s)是控制器傳遞函數(shù)，G(s)是被控對象傳遞函數(shù)。考慮到輸入r(s) =0，結(jié)構(gòu)圖可以很容易的變換成他=PlantCor toilerKD(f)圖14 PID系統(tǒng)反饋控制框圖該系統(tǒng)的輸出為numy(s)=G(s)1 KD(s)G(s)F(s= (nunXxnum)(de nPID)(de n)F(s)num(

25、den PID)(denPID)(den) (num PID)(num)其中,num 被控對象傳遞函數(shù)的分子項被控對象傳遞函數(shù)的分母項numPID pid控制器傳遞函數(shù)的分子項denPID pid控制器傳遞函數(shù)的分母項被控對象的傳遞函數(shù)是G(s) _ml s qnumU(s) 3 鼬ml2 2(M m)mglbmgl dens +ssqqq其中q 二(M m)(l ml2) (ml)22KdS KpS Kisnu mPID den PIDPID控制器的傳遞函數(shù)為KKD(s) =KdS Kp 1 s調(diào)節(jié)PID控制器的各個參數(shù)，以得到滿意的控制效果。前面討論的輸出量只考慮了擺桿角度，那么，在我們施

26、加擾動的過程中，小車位置如何變化？考慮小車位置，得到改進的系統(tǒng)框圖如下:Cor trailerPlant 2HG心)4Plant%)圖15改進的PID系統(tǒng)控制框圖其中，Gi(s)是擺桿傳遞函數(shù)， G2(s)是小車傳遞函數(shù)。由于輸入信號r(s) =0，所以可以把結(jié)構(gòu)圖轉(zhuǎn)換成:圖16等價PID系統(tǒng)控制框圖其中，反饋環(huán)代表我們前面設計的擺桿的控制器。并沒有對小車位置進行控制。注:從此框圖我們可以看出此處只對擺桿角度進行了控制,小車位置輸出為:X(s)num2den，F(xiàn) (s)1 + (n umPID)( num,) (de nPID)(den,)(nu m2)(de nPID)(der,)G2(s)

27、1 KD(s)G(s) F -(denPID)(den,)(den2)(numPID)(nugXderO F(s)其中，nug，deg， num?，den分別代表被控對象1和被控對象2傳遞函數(shù)的分子和分母。numPID 和 denPID 代表PID控制器傳遞函數(shù)的分子和分母。下面我們來求G2(s)，根據(jù)前面實驗二的推導，有可以推出小車位置的傳遞函數(shù)為4 23ml 2 s X(s)qG2 (s)4U (s)33 bml 2 (M + m)mgl bmglmglqs 十s s q其中q = (M m)(I ml2) -(ml)2可以看出，den =derb=den，小車的算式可以簡化成:(nu m

28、D(de nPID)X(s) 一 (denPID)(den ) k(numPID)(num,)F(s)Matlab 仿真 實際系統(tǒng)參數(shù)如下：M小車質(zhì)量1.096 Kgm擺桿質(zhì)量0.109 Kgb小車摩擦系數(shù)0 .1N/m/secl擺桿轉(zhuǎn)動軸心到桿質(zhì)心的長度0.2 5mI擺桿慣量0.0034 kg*m*mF加在小車上的力x小車位置0擺桿與垂直方向的夾角T采樣頻率0.005 秒注意：在進行實際系統(tǒng)的 Matlab 仿真時，請將采樣頻率改為實際系統(tǒng)的采樣頻率。1 參照上例的處理方法， 把實際系統(tǒng)參數(shù)代入， 進行 PID 控制仿真， 并找到合適的參數(shù)；2 將小車推到導軌正中間位置，并且使擺桿處于自由下

29、垂的靜止狀態(tài)；3 給計算機和電控箱通電；4. 打開計算機，在 DOS操作系統(tǒng)下，鍵入 Pend”啟動程序，并按 S”使系統(tǒng)處于準備狀 態(tài)；5. 按“V鍵起擺，等擺桿立起來并穩(wěn)定下來之后（此時為LQR控制），選擇”控制器”菜單中的“ PID”，輸入?yún)?shù)，觀察小車和擺桿的運動；（注意由于控制器只對擺桿進行了控制，所以在 PID 中輸入?yún)?shù)后小車可能向一個方向運動，此時需用手輕輕扶一 下擺桿）6. 按“T”停止擺桿，當擺桿處于自由下垂的靜止狀態(tài)，并在”控制器”菜單中選擇“PID”，并參照仿真結(jié)果，輸入PID控制器參數(shù),觀察小車和擺桿的運動；（通過調(diào)整參數(shù)可以控制擺桿擺起并能夠豎直向上，此時可能需用手

30、輕輕扶一下擺桿，以避免 小車“撞墻” ）。7. 如果控制效果不理想，調(diào)整控制器參數(shù)，直到獲得較好的控制效果；8. 運動曲線文件保存在當前路徑下，其中 pos_x.dat 為小車位置， pos_t.dat 為擺桿角 度， con_u.dat 為控制量；9. 到 MatLab 中運行如下指令， MatLab 會將保存的曲線重新繪制出來。S=load（ 路徑 文件名 ）plot(S) 1擺桿角度控制首先，需要把輸出為擺桿角度勺寸系統(tǒng)的傳遞函數(shù)用Matlab表示出來，建立一個m-文件，將下面幾行表示傳遞函數(shù)的語句拷貝進去，其中k代表比例系數(shù) Kp :M = 1.096;m = 0.109;b = 0.

31、1;I = 0.0034;g = 9.8;l = 0.25;q =(M+m)*(l+m*L2) -(m*l)A2;%simplifies in putnum = m*l/q 0 0den = 1 b*(l+m*L2)/q -(M+m)*m*g*l/q -b*m*g*l/q 0kd=1k=1ki=1n umPID= kd k ki ;den PID= 1 0 ;n umc= conv ( n um, den PID )denc= polyadd ( con v(de nPID, den ), conv( n umPID, num )t = 0 : 0.05 : 5;impulse ( n umc

32、 , dene , t )其中函數(shù)polyadd是求兩個多項式之和的函數(shù)，它不是Matlab工具，因此必須把它拷貝到polyadd.m文件中，并把該文件的目錄用addpath命令加到路徑中。Polyadd函數(shù)內(nèi)容如下：fun cti on poly=polyadd(poly1,poly2)if len gth(poly1)0poly=zeros(1,mz),short+l ong;elsepoly=lon g+short;end在這里我們假定比例、積分和微分控制都是必需的。m-文件中加入下面的語現(xiàn)在，就可以進行系統(tǒng)脈沖響應的 PID 控制仿真了。在前面的 句就可以得到系統(tǒng)的脈沖響應仿真結(jié)果：運

33、行結(jié)果及響應曲線如下：num =2.3566 0 0den =1.0000 0.0883 -27.8285 -2.3094 0kd =1k =1ki =1numc =2.3566 0 0 0denc =1.0000 2.4449 -25.4720 0.0471 0 0x 10Impulse Resp onse w ith P ID control:k=1,kd=1,ki=1dene =1.00002.4449207.82680.04719876543200.511.522.53Time (sec)3.544.55圖17初始PID參數(shù)擺角狀態(tài)圖系統(tǒng)響應是不穩(wěn)定的，不能滿足要求，需要調(diào)整參數(shù)Kp，

34、K D 和 K1，直到獲得滿意的控制結(jié)果。首先增加比例系數(shù)Kp，觀察它對響應的影響，取Kp=100，kd=1.系統(tǒng)響應如下：num=2.35660 0den=1.00000.0883-27.8285-2.30940kd =1k =100ki =1numc =2.35660 00105To0505o-0Impulse Responsew ith PID control k=1OO,ki=1,kd=1LLIiI II i-0.15-0.200.511.522.533.544.55Time (sec)圖18調(diào)節(jié)PID參數(shù)后擺角輸出圖系統(tǒng)穩(wěn)定時間約為4秒，滿足要求。由于此時穩(wěn)態(tài)誤差為 0,所以不需要改

35、變積分環(huán)節(jié)(你可以改變積分系數(shù)，觀察系統(tǒng)響應如何變化)；系統(tǒng)響應的超調(diào)量比較大，為了減小超調(diào)，增加微分系數(shù)Kd，取Kd=2o，響應結(jié)果和響應曲線如下：num =2.3566 00den =1.00000.0883 -27.8285 -2.30940kd =20k =100ki =1numc =2.3566dene =1.0000 47.2194 207.82680.047100 numc = conv(num2,denPID)0.10.090.080.070.060.050.040.030.020.010Impulse Responsew ith PID control k=100,ki=20

36、,kd=100.511.522.53Time (sec)3.544.55系統(tǒng)穩(wěn)定時間約為圖19微調(diào)PID參數(shù)后擺角狀態(tài)曲線1秒，超調(diào)約為0.04，響應滿足指標要求。2 小車位置變化仿真仿真小車位置變化的m-文件內(nèi)容如下:M = 1.096;m = 0.109;b = 0.1;I = 0.0034;g = 9.8;l = 0.25;q =(M+m)*(I+m*lA2) -(m*l)A2;%simplifies in putnum1 = m*l/q 0 0den1 = 1 b*(l+m*L2)/q -(M+m)*m*g*l/q -b*m*g*l/q 0num2 = -(I+m*lA2)/q 0 m

37、*g*l/qden2 = den1kd = 20k = 100ki = 1numPID = kd k ki;den PID = 1 0;denc = polyadd(conv(denPID,den2),conv(numPID,num1) t=0:0.05:5;impulse(numc,denc,t) 仿真結(jié)果如下：num1 =2.3566 0 0den1 =1.0000 0.0883 -27.8285 -2.3094 0num2 =-0.8832 0 23.0942den2 =1.0000 0.0883 -27.8285 -2.3094 0kd =20k =100ki =1numc =-0.8832 0 23.0942 0denc =1.0000 47.2194 207.8268 0.0471 0 00.60.5 _” _0.4 _e 0.3pm 0.2J0.1 L.J *0干-0.1 II :00.511.522.533.544.55Time (sec)圖20小車位置曲線實驗結(jié)果及與仿真結(jié)果的對比分析我們組于2010年10月前往倒立擺實驗室做實驗，我們共做了PID、LQR和極點配置法，由于事先已經(jīng)將參數(shù)設計好并經(jīng)過