現(xiàn)代控制理論MATLAB上機(jī)實(shí)驗(yàn)指導(dǎo)1

PAGEPAGE10《現(xiàn)代控制理論》MATLAB實(shí)踐指導(dǎo)書(shū)劉紅軍 編寫(xiě)《現(xiàn)代控制理論》MATLAB實(shí)踐指導(dǎo)書(shū)概述LABoratory的縮寫(xiě)，早期主要用于現(xiàn)代控制中復(fù)雜的矩陣、向量的各種運(yùn)算。由于MTLAB工具包（toolbox，如控制系統(tǒng)工具包（controlsystemstoolbox；系統(tǒng)辨識(shí)工具包（systemidentificationtoolbo；信號(hào)處理工具包（signalprocessingtoolbo；魯棒控制工具包（robustcontroltoolbo；最優(yōu)化工具包(optimization等等。由于功能的不斷擴(kuò)展，所以現(xiàn)在的已不僅僅局限與現(xiàn)代控制系統(tǒng)分析和綜合應(yīng)用，它已是一種包羅眾多學(xué)科的功能強(qiáng)大的“技術(shù)計(jì)算語(yǔ)言（TheLanguageofTechnicalComputin。MathWorks公司于1992年推出了具有劃時(shí)代意義的4.0互式模型輸入與仿真系統(tǒng)，它使得控制系統(tǒng)的仿真與CAD應(yīng)用更加方便、快捷，用戶可以方便地在計(jì)算機(jī)上建模和仿真實(shí)驗(yàn)。1997MathWorks推出的5.0版允許了更多的數(shù)據(jù)結(jié)構(gòu)，1999的5.3版在很多方面又進(jìn)一步改進(jìn)了語(yǔ)言的功能。2000年底推出的6.0。MATLAB以矩陣作為基本編程單元，它提供了各種矩陣的運(yùn)算與操作，并有較強(qiáng)的繪圖功能。集科學(xué)計(jì)算、圖像處理、聲音處理于一身，是一個(gè)高度的集成系統(tǒng)，有良好的用戶界面，并有良好的幫助功能。工程、語(yǔ)音處理、圖像處理、信號(hào)分析、計(jì)算機(jī)技術(shù)等各行各業(yè)中都有極廣泛的應(yīng)用。如何獲得MATLAB幫助在MATLAB主窗口中鍵入help，即可獲得第一層幫助：help％加重型內(nèi)容為用戶鍵入的內(nèi)容，其它為執(zhí)行后顯示的內(nèi)容。HELPtopics:matlab\general Generalpurposecommands.matlab\ops Operatorsandspecialcharacters.matlab\lang Programminglanguageconstructs.matlab\elmat Elementarymatricesandmatrixmanipulation.matlab\elfun Elementarymathfunctions.matlab\specfun Specializedmathfunctions.matlab\matfun Matrixfunctions-numericallinearalgebra.simulink\simulink Simulinksimulink\blocks Simulinkblocklibrary.simulink\simdemos Simulink3demonstrationsandsamples.simulink\dee DifferentialEquationEditor(Notableofcontentsfile)toolbox\local Preferences.如果用戶對(duì)MATLAB的語(yǔ)言結(jié)構(gòu)lang感興趣，想進(jìn)一步了解，則鍵入：helplangProgramminglanguageconstructs.Controlflow.if Conditionallyexecutestatements.else IFstatementcondition.elseif IFstatementcondition.end TerminatescopeofFOR,WHILE,SWITCH,TRYandIFstatements.for Repeatstatementsaspecificnumberoftimes.while Repeatstatementsanindefinitenumberoftimes.如果想進(jìn)一步了解for語(yǔ)句，則鍵入：helpforFORRepeatstatementsaspecificnumberoftimes.ThegeneralformofaFORstatementis:FORvariable=expr,statement, ,statementENDThecolumnsoftheexpressionarestoredoneatatimeinthevariableandthenthefollowingstatements,uptotheEND,areexecuted.……Someexamples(assumeNhasalreadybeenassignedavalue).FORI=1:N,FORJ=1:N,A(I,J)=1/(I+J-1);ENDEND同樣，如果想了解MATLAB中有關(guān)矩陣的操作運(yùn)算函數(shù)，可以鍵入：helpmatfunMatrixfunctions-numericallinearalgebra.Matrixanalysis.norm Matrixorvectornorm.normest Estimatethematrix2-norm.rank Matrixrank.det Determinant.trace Sumofdiagonalelements.null Nullspace.orth Orthogonalization.rref Reducedrowechelonform.subspace--Anglebetweentwosubspaces.Eigenvaluesandsingularvalues.eig Eigenvaluesandeigenvectors.svd Singularvaluedecomposition.gsvd Generalizedingularvaluedecomposition.eigs Afeweigenvalues.svds Afewsingularvalues.poly Characteristicpolynomial.polyeig Polynomialeigenvalueproblem.condeig Conditionnumberwithrespecttoeigenvalues.hess Hessenbergform.qz QZfactorizationforgeneralizedeigenvalues.schur Schurdecomposition.Matrixfunctions.expm Matrixexponential.logm Matrixlogarithm.sqrtm Matrixsquareroot.funm Evaluategeneralmatrixfunction.上面所列的都是有關(guān)矩陣的操作函數(shù)。如eig(A)可求出A的特征根及其特征向量，具體執(zhí)行方法為：A=[01;-6-5]％輸入AA0 1-6 -5E=eig(A)%求出方陣A的特征根EE=-2-3AVAV=0.4472-0.8944D=-20

-0.31620.94870-3基本功能我們下面給出一些及其眾多TOOLBOX系統(tǒng)中用命令查閱其它功能。MATLAB的主要線性代數(shù)運(yùn)算helpmatfun常用線性代數(shù)函數(shù)矩陣轉(zhuǎn)置C=A+B矩陣相加C=A*B矩陣相乘C=A^k矩陣冪C=A.*B矩陣點(diǎn)乘，即兩維數(shù)相同的矩陣各對(duì)應(yīng)元素相乘expm(A)指數(shù)矩陣，也就是eAinv(A)矩陣的逆矩陣det(A)矩陣的行列式的值rank(A)計(jì)算矩陣的秩eig(A)矩陣的特征值[X,D]=eig(A)矩陣的特征向量X和以特征值為元素的對(duì)角陣Dp=poly(A)矩陣的特征多項(xiàng)式r=roots(p)特征多項(xiàng)式方程的根conv(p1,p2)兩多項(xiàng)式相乘常用的控制系統(tǒng)處理函數(shù)TF2SS將傳遞函數(shù)轉(zhuǎn)換到狀態(tài)空間表達(dá)式[A,B,C,D]=TF2SS(NUM,DEN)將系統(tǒng)：G(s)

NUM(s)

b smb sm1 bsbm m1 1 0DEN(s轉(zhuǎn)換成：AX YCX DU其中:

sna

n

sn1 a1

sa0NUM=[bm,bm-1,…,b1,b0]，DEN=[1,an-1,an-2,…,a1,a0]a

n1

a an2 1

a 0

1 1A 00

0 01 0 0 1

0 0 0

0B 00Cbn1

bn

b

Dbn0ZP2SS將零極點(diǎn)型傳遞函數(shù)轉(zhuǎn)換到狀態(tài)空間表達(dá)式0[A,B,C,D]=ZP2SS(Z,P,K)K(sz)(sz )(sz )除了Gs)

(sp1

1)(sp

2 m) (sp 2 n

以外，其它與TF2SS相同。SS2TF將狀態(tài)空間表達(dá)式轉(zhuǎn)換到傳遞函數(shù)[NUM,DEN]=SS2TF(A,B,C,D,iu)即求第iu個(gè)輸入信號(hào)對(duì)輸出y(t)的傳遞函數(shù)，即：G(s)

Y(s)

NUM(s) C(sIA)1B

NUM(s) U (s)iu

DEN(s)

sna

n1

sn1 a1

sa0SS2TF的調(diào)用返回值為G(s)的分子多項(xiàng)式的系數(shù)矩陣NUM和分母多項(xiàng)式的系數(shù)向量DEN。SS2ZP將狀態(tài)空間表達(dá)式到零極點(diǎn)形式的傳遞函數(shù)的轉(zhuǎn)換[Z,P,K]=SS2ZP(A,B,C,D,iu)（5）TF2ZP一般傳遞函數(shù)轉(zhuǎn)換到零極點(diǎn)型傳遞函數(shù)[Z,P,K]=TF2ZP（NUM,DEN）ZP2TF零極點(diǎn)型傳遞函數(shù)轉(zhuǎn)換到一般傳遞函數(shù)[NUM,DEN]=ZP2TF(Z,P,K)SS2SS狀態(tài)空間表達(dá)式的線性變換[A1,B1,C1,D1]=ss2ss(A,B,C,D,T)其中T為變換矩陣。注意變換方程為：X1=TX，而不是常見(jiàn)的X=TX1。所以要與用戶習(xí)慣的變換方程一致，我們必須用T的逆代入上式，即：[A1,B1,C1,D1]=ss2ss(A,B,C,D,inv(T))CANON求狀態(tài)空間表達(dá)式的對(duì)角標(biāo)準(zhǔn)型[As,Bs,Cs,Ds,Ts]=canon(A,B,C,D,'mod')其中Ts為變換矩陣，注意變換方程為：Xs=TsX。CTRBM=ctrb(A,B)OBSVN=obsv(A,C)IMPULSE[y,x]=impulse(A,B,C,D,in,t)[y,x]=impulse(num,den,t)STEP[y,x]=step(A,B,C,D,in,t)[y,x]=step(num,den,t)LSIM求系統(tǒng)對(duì)任意輸入函數(shù)u(t)[y,x]=lsim(A,B,C,D,u,t)[y,x]=lsim(num,den,u,t)C2D連續(xù)系統(tǒng)狀態(tài)方程轉(zhuǎn)換為離散狀態(tài)方程,T為采樣周期[G,H]=c2d(A,B,T)相關(guān)的函數(shù)還有D2C,D2D。LYAPAPPATQATPPAQ代入，即：P=Lya（’，。PLACE極點(diǎn)配置F=PLACE（A，B，P）PLOT繪圖函數(shù)plot(X,Y,’str’)可以用不同顏色、不同符號(hào)繪制曲線，其中‘str’可以是下列三組選項(xiàng)的任意組合。y yellowm magentac cyanr redg greenb bluew whitek black4

.pointo circlex x-mark+ plus* stars squared diamondv triangle(down)^ triangle(up)< triangle(left)> triangle(right)p pentagramh hexagram

-solid:dotted-.dashdot--dashed[例1]給定某控制系統(tǒng)的狀態(tài)空間描述為：0

1 0 x66y

11110

6x5

0u1試求對(duì)角規(guī)范型和變換矩陣W，并根據(jù)其對(duì)角規(guī)范型繪制系統(tǒng)在初值x10 15u0解：closeall;%狀態(tài)方程模型A=[01-1;-6-116;-6-115]; B=[0;0;1]; C=[100];%求取對(duì)角規(guī)范型[W,lamda]=eig(A); L=inv(W)*A*W; b=inv(W)*B; c=C*W;%顯示結(jié)果disp('TheDiagonalCanonicalFromofSystemis:');Lbcdisp('TransformationMatrix W%仿真數(shù)據(jù)初始化t=0:0.01:3; x0=[5;10;15];xx0=inv(W)*x0; x=zeros(3,n); xx=zeros(3,n);%求解狀態(tài)變量fori=1:nxx(:,i)=expm(L*t(i))*xx0;x(:,i)=W*xx(:,i);end%繪圖plot(t,x)title('SystemTimeResponsewithu=0');xlabel('Time-sec'); ylabel('Response-value');text(0.8,3.7,'x(1)');text(0.8,1.5,'x(3)');text(0.8,-0.4,'x(2)');>>TheDiagonalCanonicalFromofSystemis:L=-1.0000-0.0000-0.00000.0000-2.0000-0.0000-0.00000.0000-3.0000b=-2.8284-13.747710.8628c=0.7071-0.2182-0.0921TransformationMatrix is:W=0.7071 -0.2182 -0.09210.0000-0.4364-0.55230.7071-0.8729-0.8285的標(biāo)量的指數(shù)即可。在本例具體的求解的過(guò)程中，需要注意的是求得對(duì)角規(guī)范型的狀態(tài)xxWx，然后畫(huà)圖。從系統(tǒng)的響應(yīng)曲線來(lái)看，各狀態(tài)變量最終都趨于零。這一方面是因?yàn)檩斎肓縰=0，另一方面是因?yàn)槌踔挡粸榱?，但矩陣A的特征值均為零，系統(tǒng)的零輸入響應(yīng)呈衰減趨勢(shì)。但狀態(tài)變量x（3）下降幅度和速度最大，x（1）下降幅度和速度最小。這是因?yàn)闋顟B(tài)變量x（3）對(duì)應(yīng)的特征值是-3，而x（1）對(duì)應(yīng)的特征值是-1，其衰減速率不同。[例2]給定線性定常系統(tǒng)的狀態(tài)空間描述為：2 0 0 4 x

1 13 0x u0 0

4 1

1y

0 0 2 21 0 0試判斷系統(tǒng)的能控性。如果系統(tǒng)狀態(tài)完全能控，試求其能控標(biāo)準(zhǔn)型和變換矩陣。解：%系統(tǒng)狀態(tài)方程模型A=[2001;0413;0041;000B=[1;0;1;2];C=[1100];%系統(tǒng)階次n=length(A);%求解系統(tǒng)能控性矩陣Q=ctrb(A,B);%系統(tǒng)能控性矩陣的秩m=rank(Q);%判斷系統(tǒng)是否狀態(tài)完全能控，并求解能控標(biāo)準(zhǔn)型ifm==nAc1=inv(Q)*A*Q;elseend

Bc1=inv(Q)*B;Cc1=C*Q;disp('SystemisControllable.');disp('SystemFirstControllableCanonnicalFormAc1Bc1Cc1disp('TheTransformationMatrixis:');Qdisp('SystemStateVariablecannotbetotallyControlled.');disp('TherankofSystemControllableMatrixis:');m運(yùn)行結(jié)果：SystemisControllable.SystemFirstControllableCanonnicalFormis:Ac1=0 -0.0000 -0.0000 -64.00001.0000 -0.0000 -0.0000 96.00000 1.0000 -0.0000 -52.00000.00000.00001.000012.0000Bc1=1000Cc1=11158268TheTransformationMatrixis:Q= 1412320746236162812024816例3]求其能觀標(biāo)準(zhǔn)型和變換矩陣。解：%系統(tǒng)狀態(tài)方程模型A=[2001;0413;0041;000B=[1;0;1;2];C=[1100];%系統(tǒng)階次n=length(A);P=obsv(A,C);%系統(tǒng)能觀性矩陣的秩m=rank(P);ifm==nAo1=P*A*inv(P);Bo1=P*B;elseend

Co1=C*inv(P);disp('SystemisObservable.');disp('SystemFirst ObservableCanonnicalFormis:');Ao1Bo1Co1disp('TheTransformationMatrixis:');inv(P)disp('SystemStateVariablecannotbetotallyObserved.');disp('TherankofSystemObservableMatrixis:');m運(yùn)行結(jié)果：ystemisObservable.SystemFirst ObservableCanonnicalFormis:Ao1=010000100001-64 96 -52 12Bo1=11158268Co1=1 0 0 0TheTransformationMatrixis:ans=-14.0000 16.0000 -5.3750 0.562515.0000-16.00005.3750-0.562501.0000-0.75000.1250-8.00008.0000-2.50000.2500[例4]某控制系統(tǒng)的狀態(tài)方程描述如下。試判斷其穩(wěn)定性并繪制其時(shí)間響應(yīng)來(lái)驗(yàn)證上述判斷。3 8 2 4 11 0 x

0 0x u0 1 0 0

0 0 0 0y0 1 解：z=-1p=-1.4737+2.2638i-1.4737-2.2638i-0.0263+0.7399i-0.0263-0.7399ik=1Systemisstable-0.0263，非?？拷撦S，在實(shí)際過(guò)程中，例5]含積分環(huán)節(jié)的伺服系統(tǒng)設(shè)計(jì)，設(shè)對(duì)象為G(s) 1 s(s1)(s2)設(shè)計(jì)控制器uKx解：%Poleplacementcloseall;a=[010;001;0-2-3];b=[0;0;1];c=[100];d=[0];

r，使閉環(huán)系統(tǒng)具有極點(diǎn)2j2 3,10。1disp('Therankofcontrollabilitymatrix')rc=rank(ctrb(a,b))p=[-2+2*sqrt(3)*j-2-2*sqrt(3)*j-10];k=place(a,b,p)a1=a-b*k;b1=b*k(1);c1=c;d1=d;figure(1)step(a1,b1,c1,d1)title('StepResponeseofDesignedServoSystem')figure(2)[y,x,t]=step(a1,b1,c1,d1);plot(t,x)title('StepResponeseCurvesforx1,x2,x3')gridonTherankofcontrollabilitymatrixrc=3k= 160.0000 54.0000 11.0000 [例6]開(kāi)環(huán)系統(tǒng)的狀態(tài)方程為0 1 0 00 0 1x

0u 6y

110

6

1設(shè)計(jì)全維狀態(tài)觀測(cè)器，使觀測(cè)器的閉環(huán)極點(diǎn)為2j2 3,5。解：狀態(tài)觀測(cè)器的結(jié)構(gòu)見(jiàn)教材圖5-5，為求出狀態(tài)觀測(cè)器的反饋增益陣K，先為原系統(tǒng)構(gòu)造一對(duì)偶系統(tǒng)apcqbp然后采用極點(diǎn)配置方法對(duì)對(duì)偶系統(tǒng)進(jìn)行閉環(huán)極點(diǎn)培植，得到反饋增益陣kkek。程序如下：%Designofafull-orderstateobservera=[010;001;-6-11-6];b=[0;0;1];c=[100];%Checkobservabilitydisp('Therankofobservabilitymatrix')r0=rank(obsv(a,c))%Constructthedaulsystem(a1,b1,c1)a1=a';b1=c';c1=b';%sloveKusingpolesplacementp=[-2+2*sqrt(3)*j-2-2*sqrt(3)*j-5];k=acker(a1,b1,p);ke=k'運(yùn)行結(jié)果：Therankofobservabilityr0=3ke=3.00007.0000-1.00005上機(jī)習(xí)題1假設(shè)系統(tǒng)由兩個(gè)1

s)和G2s)串聯(lián)連接而成，已知G(s) s1 且

s3s5 )1 s23s4 2

s44s33s22s1若想求出總系統(tǒng)的狀態(tài)方程模型，請(qǐng)?jiān)贛ATLAB下比較下面兩種方法有何不同結(jié)果：將兩個(gè)傳遞函數(shù)模型進(jìn)行串聯(lián)連接，然后求出整個(gè)系統(tǒng)的狀態(tài)方程模型。求出兩個(gè)模型的狀態(tài)方程表示，然后求出整個(gè)系統(tǒng)的狀態(tài)方程模型。某位置控制系統(tǒng)方塊圖如下圖所示，求其狀態(tài)空間表達(dá)式。5-3試將5-2題所得動(dòng)態(tài)方程轉(zhuǎn)換成Jordan2=54語(yǔ)言求解。5-4k=2，k1=－0.05，k2=0.16，k3=.245-2所示系統(tǒng)的單位階躍響應(yīng)。5-5判斷