三級倒立擺建模與控制器設(shè)計仿真

1、三級倒立擺摘要：從剛體運動學(xué)角度分析了三級倒立擺系統(tǒng)的運動，運用狀態(tài)空間法對三級倒立擺系統(tǒng)進(jìn)行數(shù)學(xué)模型的建立;應(yīng)用二次型最優(yōu)控理論實現(xiàn)三級倒立擺控制；最后根據(jù)計算結(jié)果，通過matlab程序仿真，得出上中下擺的角度和小車的位移曲線圖。關(guān)鍵詞：倒立擺；最優(yōu)控；仿真；matlab三級倒立擺系統(tǒng)模型 一個典型的三級直線倒立擺系統(tǒng)主要由機(jī)電裝置和控制裝置兩部分組成,機(jī)電裝置的結(jié)構(gòu)見圖!小車圖1系統(tǒng)分析示意圖 其主要由擺桿、電機(jī)、小車、導(dǎo)軌、傳動齒輪、齒條及各種傳感器等構(gòu)成假設(shè)系統(tǒng)中的每一根擺桿都是勻質(zhì)剛體，驅(qū)動力與放大器的輸入成正比且無延遲地直接作用于小車上，并且可以在忽略實驗中的庫侖摩擦和動摩擦的前提

2、下，圖1所示的倒立擺系統(tǒng)的物理分析：yf22三級倒立擺受力分析圖：f21nx1f11f12f f11f12 小車受力分析圖yf31y下擺受力分析f33f21x23xf33f22f31 中擺受力分析上擺受力分析下擺牛頓定律及動力矩方程：f12-f22=m1d2dt2(r+l1sin1)m1g-f11+f21=m1d2dt2(l1cos1)j1=f11l1sin1+f21l1-l1sin1-f12l1cos1-f22(l1-l1)cos1中擺受力分析： f22-f31=m2d2dt2(r+l1sin1+l2sin2)m2g-f21+f33=m2d2dt2(l1cos1+l2cos2)j2=f21l

3、2sin2+f33l2-l2sin2-f22l2cos2-f31(l2-l2)cos2丄擺受力分析：f31=m3d2dt2(r+l1sin1+l2sin2+l3sin3)m3g-f33=m3d2dt2(l1cos1+l2cos2+l3cos3)j3=f33l3sin3-f31l3cos3 注：l1 l2 l3下，中，上擺質(zhì)心到鉸鏈處距離； l1 l2 下，中擺的長度； 1 2 3-上，中，下擺偏離豎直方向的角度； m1 m2 m3-下，中，上擺質(zhì)量;其動力學(xué)方程如下：m(1 ,2 ,3)r123=f(1,2 ,3,1,2 ,3)r123+n1 ,2 ,3r123+gu m(1 ,2 ,3)=m

4、0+m1+m2+m3 m1d1cos1+m2d4cos1+m3d4cos1m2d2cos2+m3d5cos2m3d3cos3 m1d1cos1+m2d4cos1+m3d4cos1j1+m1d12+m2d42+m3d42m2d2d4cos(1-1)+m3d4d5cos(2-1)m3d3d4cos(3-1) m2d2cos2+m3d5cos2m2d2d4cos(1-1)+m3d4d5cos(2-1)j2+m2d22+m3d52m3d3d5cos(3-2)m3d3cos3m3d3d4cos(3-1)m3d3d5cos(3-2)j3+m3d32f(1,2 ,3,1,2 ,3)=-f0000 m1d1s

5、in1+m2d4sin1+m3d4sin11-f1+f2f2-m2d2d4sin(2-1)+m3d4d5sin(2-1)-m3d3d4sin(3-1) m2d2sin2+m3d5sin22f2+m2d2d4+m3d4d52sin2-1-f3+f2f3-m3d3d52sin3-2m3d3sin33m3d3d4sin(3-1)f3+m3d3d53sin3-2-f3n(1,2 ,3)=0(m1d1+m2d4+m3d4)sin1(m2d2+m3d5)gsin2m3d3gsin3g=g0000t將此非線性模型在系統(tǒng)的平衡點r=1=2=3=0,r= 1=2=3=0處線性化得到系統(tǒng)的狀態(tài)方程: x=ax+b

6、uy=cx+du a=04*4i4*4a21a22 , b=04*1b2其中： 將參數(shù)帶入方程計算：a=00000000000000-3.655145.74-49.782-5.159000.12698-80.216169.9-114.49 010001000.44281-8.3647-4.78172.26900-16.80957.098-62.1436.4401000.016897-0.562720.98145-0.53453 000010-0.00564470.32956-0.659640.5167301-0.0032595-0.0721790.32957-0.76437 b=00008.7

7、2-29.62132.238-3.341tc=10000100001000010000000000000000 d=0000t控制器設(shè)計根據(jù)由圖所示的倒立擺系統(tǒng)模型的推導(dǎo)可知，該系統(tǒng)在一定的前提下，可以轉(zhuǎn)變?yōu)榫€性系統(tǒng)#對于線性系統(tǒng)的控制，需要經(jīng)過系統(tǒng)的穩(wěn)定性、能控性和能觀性判斷后才能選擇其最優(yōu)的控制策略#圖!顯示了式（$）所示的線性系統(tǒng)外控制器的設(shè)計過程狀態(tài)方程的重構(gòu)過程，實質(zhì)上就是本文從式式轉(zhuǎn)化為式式的過程，在matlab中，可以用如下代碼表示: 系統(tǒng)狀態(tài)方程的重構(gòu) 狀態(tài)變量 (x1,x2,x3,x4,x5,x6,x7,x8)=(r,zeta1,zeta2,zeta3,r,zeta1,zet

8、a2,zeta3,) 輸出變量y=x a=zeros(8,8); a(1:4,5:8)=eye(4,4); a(5:8,1:4)=t0m0 0 0-1n0 0 0 t0-1; a(5:8,5:8)= t0m0 0 0-1-f0 0 0 0 0 0 t0-1; 不用invese而用blackslash是由于blackslash有更好的錯誤檢查機(jī)制并且少用時間和內(nèi)存。 b=zeros(8,1); b(5:8)= t0m0 0 0-1g0 0 0 0t; c=zeros(4,8) c(1:4,1:4)=eye(4,4); d=zeros(4,1);系統(tǒng)的穩(wěn)定性、能控性及能慣性判斷 對于如式（2）式所

9、示的線性定常系統(tǒng)而言，其穩(wěn)定性可通過計算系統(tǒng)矩陣a 的特征值來確定(李雅普諾夫線性系統(tǒng)穩(wěn)定判據(jù))，如果特征值均處于s復(fù)平面左半平面，則系統(tǒng)穩(wěn)定。 matlab中，用函數(shù)eig（a）來計算系統(tǒng)矩陣的特征值，經(jīng)過計算，式式所示的線性定常系統(tǒng)的特征值為： r1 =0 13.2715 8.3919 3.1656 -18.2950 -2.9940 -8.5547 -13.7810因此，該系統(tǒng)為不穩(wěn)定系統(tǒng)，需加外部控制器形如式（2）所示的線性定常系統(tǒng)能控的充分必要條件是： rankb aban-1b=n;matlab中，用“ctrb(a,b)”來求系統(tǒng)能控矩陣pc=b aban-1b的秩，經(jīng)過計算rank

10、(pc)=8,系統(tǒng)能控。形如式（2）所示的線性定常系統(tǒng)能觀的充分必要條件為：rankccacan-1=nmatlab中，用“obsv(a,c)”來求系統(tǒng)能觀陣qo=ccacan-1的秩，計算的，rank(qo)=8,系統(tǒng)能觀。綜上所述，式（2）所示的線性定常系統(tǒng)為不穩(wěn)定的能控、能觀系統(tǒng)，可加外控制器實現(xiàn)系統(tǒng)的穩(wěn)定。線性最有調(diào)節(jié)器的設(shè)計 定義狀態(tài)方程的二次性能指標(biāo)函數(shù)為： j=120(xtqx+utru)dt 正定r陣的作用在于限制控制器的輸出不會太大而導(dǎo)致難于實現(xiàn)，所以不失一般性，可以將其設(shè)為1，lqr方法中，q陣可以選為非負(fù)的對角陣，q=diag（ q11，q22，q33，q44，q55，q

11、66，q77，q88）.各主對角系數(shù)的相對大小對系統(tǒng)的整體控制性能有很大的影響，目前普遍采用試探實驗的方法來確定加權(quán)矩陣q的系數(shù)大小。觀察知：q11大小的變化對小車的位移影響很大，因為小車位移大小受直線導(dǎo)軌長度的物理限制，所以，q11取值可以固定為一個適當(dāng)大小的數(shù)，使小車位移不至過大； q22對系統(tǒng)控制的影響最為明顯，q33次之，q44的影響較小。取值時可以去q22>q33>q44,取q=diag(500,500,300,100,0,0,0,0),r=1.求解riccati代數(shù)方程 pa+atp-pbr-1btp+q=0 可以得到狀態(tài)反饋陣 kt=r-1btp=-22.3607-8

12、0.103528.5577-408.4518-30.1003-53.7813-44.7459-53.7370此時各擺角的動態(tài)和穩(wěn)態(tài)指標(biāo)較為理想。初始條件為r0=0，10=0.05>0，20=-0.05<0，30=0.05>0。仿真結(jié)果程序：>>a=0,0,0,0,1,0,0,0;0,0,0,0,0,1,0,0;0,0,0,0,0,0,1,0;0,0,0,0,0,0,0,1;0,-3.6551,0.12698,0.44281,-16.809,0.016897,-0.0056447,-0.0032595;0,45.74,-80.216,-8.3647,57.098,-

13、0.56272,0.32956,-0.072179;0,-49.782,169.9,-4.781,-62.143,0.98145,-0.65964,0.32957;0,-5.159,-114.49,72.269,6.4401,-0.53453,0.51673,-0.76437;>>b=0;0;0;0;8.72;-29.621;32.238;-3.341;>> c=1,0,0,0,0,0,0,0;0,1,0,0,0,0,0,0;0,0,1,0,0,0,0,0;0,0,0,1,0,0,0,0;>> d=0;0;0;0;>> r1=eig(a)>

14、> q1=500;q2=500;q3=300;q4=100;q5=0;q6=0;q7=0;q8=0;>> q=q1 0 0 0 0 0 0 0;0 q2 0 0 0 0 0 0;0 0 q3 0 0 0 0 0;0 0 0 q4 0 0 0 0;0 0 0 0 q5 0 0 0;0 0 0 0 0 q6 0 0;0 0 0 0 0 0 q7 0;0 0 0 0 0 0 0 q8;r=1;>> k,p,r2=lqr(a,b,q,r)>> ac=(a-b*k);bc=b;cc=c;dc=d;>> t=0:0.02:20;u=zeros(size

15、(t);>> x0=0;0.05;-0.05;0.05;0;0;0;0;>> y,x=lsim(ac,bc,cc,dc,u,t,x0);>> figure(1);plot(t,y(:,1);>> xlabel('time/sec');ylabel('x(小車)/m'); title('小車位移');grid;>> figure(2);plot(t,y(:,2);>> xlabel('time/sec');ylabel('01/rad');title('下擺角度變化曲線');grid;>> fig