線性控制系統(tǒng)的計(jì)算機(jī)輔助分析

第4幸

鐵供控制/髭的計(jì)算機(jī)輔助臺(tái)折

2015/8/301

系統(tǒng)的分析方法

■充分利用計(jì)算機(jī)對線性系統(tǒng)進(jìn)行分析

■更新系統(tǒng)分析的觀念

■求解傳統(tǒng)方法難以求解的問題

■離散系統(tǒng)穩(wěn)定性如何分析？

■Nyquist圖、Nichols圖沒有頻率信息，如何彌補(bǔ)?

■高階系統(tǒng)的根軌跡如何繪制？

2015/8/302

本章主要內(nèi)容

■線性系統(tǒng)定性分析

■線性系統(tǒng)時(shí)域響應(yīng)解析解法

■線性系統(tǒng)的數(shù)字仿真分析

■根軌跡分析

■線性系統(tǒng)頻域分析

2015/8/303

4.1線性系統(tǒng)定性分析

■主要內(nèi)容

■線性系統(tǒng)穩(wěn)定性分析

■線性反饋系統(tǒng)內(nèi)部穩(wěn)定性分析

■線性系統(tǒng)的相似變換

■線性系統(tǒng)可控性分析

■線性系統(tǒng)可觀測性分析

■Kalman分解

■系統(tǒng)狀態(tài)方程的標(biāo)準(zhǔn)型

2015/8/304

4.1.1線性系統(tǒng)的穩(wěn)定性分析

■給定線性系統(tǒng)模型，如何分析穩(wěn)定性？

_10/+50?+100?+1005+40

,-十+21,2+]84/+87()54+2384s3+3664s2+2496s

■由控制理論可知，用Routh表

可以判定該系統(tǒng)穩(wěn)定性

■EdwardJohnRouth(1831-1907)

■歷史局限性

2015/8/305

狀態(tài)方程系統(tǒng)的穩(wěn)定性

■連續(xù)線性狀態(tài)方程

x(t)-Ax(t)+Bu(t)

[y(t)=Cx(t)+Du(t)

■解析階

以f)=c'"-，。)以“))+IC^(f~r)Bu(T)dT

Jo

■穩(wěn)定性：A矩陣的特征根均有負(fù)實(shí)部

2015/8/306

離散系統(tǒng)的穩(wěn)定性

■離散系統(tǒng)狀態(tài)方程

{x\(k+l)rj=Fx(kT)+Gu(kT)

[y(kT)=Cx(kT)+Du(kT)

■離散系統(tǒng)時(shí)域響應(yīng)解析階

k—1

x(kT)二尸"/(())+￡Fj'Gug

z=0

■穩(wěn)定性判定：所有特征根均在單位圓內(nèi)

2015/8/307

Routh判據(jù)的歷史局限性

■Routh判據(jù)提出時(shí)，沒有求多項(xiàng)式根的方法

■現(xiàn)在求解矩陣特征根、求解多項(xiàng)式方程的

根輕而易舉，無需間接方法

■Routh判據(jù)只能得出是否穩(wěn)定，進(jìn)一步信息

得不出來，如系統(tǒng)是否振蕩

■離散系統(tǒng)無法由Routh方法直接判定，得借

助于Jury判據(jù)，更復(fù)雜

■穩(wěn)定性分析方法不統(tǒng)一

2015/8/308

基于MATLAB的穩(wěn)定性判定方法

■直接判定

■狀態(tài)方程模型

■由eig(A)可以求出所有特征根

■離散系統(tǒng)：abs(eig(A))

■傳遞函數(shù)模型：完全同樣方法

■圖解判定法

■連續(xù)系統(tǒng)：pzmap(G)

■離散系統(tǒng)：pzmap(G)?同時(shí)畫出單位圓

2015/8/309

例4-1高階系統(tǒng)穩(wěn)定性判定

I0,S+50-+100./+]00.9+40

G(s)二F----7-----------;------;------------

s1+2156+184s5+870s4+2384s3+3664?+2496s

■直接分析方法

6為、?num=[10,50,100,100,40];

den=[1,21,184,870,2384,3664,2496,0];

G=tf(num,den);GG=feedback(G,1);

pzmap(GG)

eig(GG)'

10(5+2)(5+I)(52+2s+2)

■零極點(diǎn)模型G(s)二

(5+6.922)(3+2.635)(3+0.01577)

(52+4.127s+747)(/+7.3s+18.62)

2015/8/3010

例4-2高階離散單位負(fù)反饋系統(tǒng)模型

6工2-().6z-0.12

+0.2512+0.25—0.125

()，6

Gc(z)=0.3""

CKz+0.8

■MATLAB求解

今秦?den=[l-10.250.25-6).125];

num=[6-?.6-0.12];H=tf(num,den,'Ts',0.1);

z=,'Ts',?.1);Gc=0.3*(z-0.6)/(z+?.8);

GG=feedback(H''Gc,1);

pzmap(GG),absfeig(GG)')

2015/8/3011

4.1.2線性反饋系統(tǒng)的內(nèi)部穩(wěn)定性

■輸入、輸出穩(wěn)定是不夠的，因?yàn)槿魞?nèi)部信

號(hào)可能過大，對系統(tǒng)的硬件破壞

■應(yīng)該引入內(nèi)部穩(wěn)定性概念，保證內(nèi)部信號(hào)

也是穩(wěn)定的。

2015/8/3012

■由給定穩(wěn)定輸入/.，d,n到內(nèi)部信號(hào)工］,、2,.口

都穩(wěn)定的系統(tǒng)稱為內(nèi)部穩(wěn)定系統(tǒng)

■傳遞函數(shù)矩陣

I—G(s)H(s)—"(s)

_1

Gc(s)1—Gc(s)"(s)d

一用(s)

G(s)Gc(s)G(s)1n

其中例(s)=1+G(s)Gc(s)〃(s)

■逐一判定每個(gè)子傳遞函數(shù)的穩(wěn)定性很煩瑣

■內(nèi)部穩(wěn)定性定理

2015/8/3013

內(nèi)部穩(wěn)定性定理

■閉環(huán)系統(tǒng)內(nèi)部穩(wěn)定的充要條件為

■1+H(s)G(s)Gc(s)沒有不穩(wěn)定零點(diǎn)

■〃(.s)G(s)Gc(s)沒有不穩(wěn)定零極點(diǎn)對消

■第一個(gè)條件等效于輸入輸出穩(wěn)定性

■判定第2條件即可

■可以編寫MATLAB函數(shù)判定內(nèi)部穩(wěn)定性

key=intstable(G,,〃)

2015/8/3014

■判定的MATLAB函數(shù)

functionkey=intstable(G,Gc,H)

GG=minreal(feedback(G*Gc,H));

Go=H*G*Gc;Gol=minreal(Go);

p=eig(GG);z?=eig(Go);zl=eig(Gol);

zz=setdiff(zO,zl);

if(G.Ts>1),%禺散系統(tǒng)劃定

key=any(abs(p)>l);

ifkey==Q,key=2*any(abs(zz)>l);end

else,%連續(xù)系統(tǒng)判定

key=any(real(p)>?);

ifkey=0,key=2*any(real(zz)>?);end

end

2015/8/3015

4.1.3線性系統(tǒng)的線性相似變換

■系統(tǒng)的狀態(tài)方程表示稱為系統(tǒng)實(shí)現(xiàn)

■不同狀態(tài)選擇下，狀態(tài)方程不惟一

■相似變換

■非奇異矩陣T

■狀態(tài)變換N=T~]x

■新/狀態(tài)方程模型

2⑺=Az(t)十=1

1,且z(())=T—l力(0)

y⑴=GN⑺+D{u(t)

2015/8/3016

■狀態(tài)變換公式

l1

A{=T-AT,Bt=T-B

Ct=CT,Dt=D

■MATLAB求解方法

G!=SS2SS(G,T)I

2015/8/3017

例4-3已知系統(tǒng)和轉(zhuǎn)換矩陣

0I00

00I0

由⑺二

()00I

-24-50-35-IO

「V⑺=|24247

■MATLAB求解

?A=[0l00;0010;??01;-24-50-35-10];

Gl=ss(A,,[242471],0);

T=fliplr(eye(4));G2=ss2ss(Gl,T)

2015/8/3018

■變換結(jié)果

-10-35-50-241

I0000

2⑺=N⑺+Z)

0100()

00100

172424]z(t)

■可見，相似變換能改變系統(tǒng)的結(jié)構(gòu)

■引入相似變換矩陣，可以將已知系統(tǒng)轉(zhuǎn)換

成其他的形式

2015/8/3019

4.1.4線性系統(tǒng)的可控性分析

■可控性定義

■假設(shè)系統(tǒng)由狀態(tài)方程(A3,。,。)給出,對

任意的初始時(shí)刻如如果狀態(tài)空間中任一狀

態(tài)勺⑺可以從初始狀態(tài)犬偵)處，由有界的輸

入信號(hào)u(t)的驅(qū)動(dòng)下,在有限時(shí)間年內(nèi)能夠到

達(dá)任意預(yù)先指定的狀態(tài)??(")，則稱此狀態(tài)是可控

的。如果系統(tǒng)中所有的狀態(tài)都是可控的，則稱該

系統(tǒng)為兌全可控的系統(tǒng)。

■系統(tǒng)的可控性就是指系統(tǒng)內(nèi)部的狀態(tài)是不是可以由外

部輸出信號(hào)控制的性質(zhì)，

2015/8/3020

線性系統(tǒng)的可控性判定

■可控性判定矩陣

2

Tc=[B,AB,AB.

■若矩陣Tc是滿秩，則系統(tǒng)完全可控。

■基于MATLAB的判定方法

rank(T)

■構(gòu)造可控性判定矩陣7^=ctrb(A,B)

2015/8/3021

例4-4離散狀態(tài)方程的可控性

-2.2-0.71.5—169

0.2-6.36-1.546

ce[(A+l)T]=力(A7)+u(kT)

0.6-0.9-2-0.544

1.4-0.1—1-3.584

■MATLAB求解

___

〉〉A(chǔ)=[_2.2,0.7,Q.2,6.3,6,l.5;.

0.6,-0.9,-2,-0.5;1.4,-?.1,-1,-3.51;

B=[6,9;4,6;4,4；8,4]；Tc=ctrb(A,B)

rank(Tc)

2015/8/3022

■判定矩陣

69-18-225452-162-118

46-12-183658-108-202

44-12-103626-108-74

84-24-6722-21634

■判定矩陣構(gòu)造方法

Tcl=[B,A-B,A-2*B,AT*B];

■這樣的判定方法同樣適合于連續(xù)系統(tǒng)和離

散系統(tǒng)。也適用于多變量模型

2015/8/3023

由Gram矩陣判定可控性

■引入可控Gram矩陣

「8T

At7A

Lc=Ie-BBe-'dt

■該矩陣滿足Lyapunov方程

TT

ALC+LCA二-BB

■MATLAB求解|&二

■矩陣構(gòu)造|Gc二gram(G,'c')

2015/8/3024

例4-5求Gram矩陣

6?-0.6z-0.12

Hg=

J—N3+().253+O.25Z—0.125

■MATLAB命令

?num=[6-0.6-?.12];

den=[1-1?,25?.25-0.125];

H=tf(num,den,,Q.1)

Lc=gram(ss(H),，c，)

■Gram矩陣10.76515.7547.35180

15.75443.06131.5083.6759

IJc—

7.351831.50843.0617.8769

03.67597.87692.6913

2015/8/3025

可控性階梯分解

■對于不完全可控的系統(tǒng)階梯分解

■階梯標(biāo)準(zhǔn)型

——_

Ac二才，8c二E,Cc—\CQ,CCI

A21AcBe

■MATLAB函數(shù)調(diào)用

[Ac,Bc,Cc2]=ctrbf(48,C)|

■若原系統(tǒng)狀態(tài)方程完全可控，則不必分解

2015/8/3026

例4-6不完全可控系統(tǒng)

-2.2-0.71.569

0.2—6.36-1.54

x[(k+\)T\=x(kT)+.:u(kT)

0.6-0.9-2-0.54

1.4-0.1—1-3.584

今｛>>A-[-2.2,-0.7,Q.2,-6.3,6,-1.5;...

0?6,-。.9,-2,-Q.5;1.4,—0.1,—1,-3.5J]

B=[6,9;4,6;4,4;8,4];C=口2341;

[Ac,Bc,Cc,Tc]=ctrbf(A,B,C);

-400000

x[(h-\)T]=-4.638-3.823-0.5145-0.1277"2.754-2(575l，U/)

-3.6370.1827-3.492-0.1215

-4.114-1.8881.275-2.685-11.15-11.93

2015/8/3027

4.1.5線性系統(tǒng)的可觀測性分析

■可觀測性定義

■假設(shè)系統(tǒng)由狀態(tài)方理(AB,C.。)給出，對

任意的初始時(shí)刻外，如果狀態(tài)空間「I?任一狀

態(tài)芍⑺在任意有限時(shí)刻"的狀態(tài)勺(")可以由

輸出信號(hào)在這一時(shí)間區(qū)間內(nèi)te[z0Jf]的值精確地

確定出來，則稱此狀態(tài)是可觀測的。如果系統(tǒng)中

所rr的狀態(tài)都是可觀測的，則稱該系統(tǒng)為完全e

觀測的系統(tǒng)。

■系統(tǒng)的可觀測性就是指系統(tǒng)內(nèi)部的狀態(tài)是不是可以

由系統(tǒng)輸出信號(hào)重建起來的性質(zhì)

2015/8/3028

可觀測性判定

■判定矩陣[c-

CA

2

To=CA

.*

CAn-1

■等同于(AT,CT)系統(tǒng)可控性判定

■Gram矩陣上。=A％丁@-4山

Jo

■MATLAB求解|gramdjc77

2015/8/3029

■Gram矩陣滿足Lyapunov方程

TT

A￡0+L0A=-CC

■對偶問題Q4,8,C)

(ATCT,BT)

2015/8/3030

4.1.6Kalman規(guī)范分解

■Kalman規(guī)范分解

4,64,200

00；0

之⑺=-*****??N⑺+

<A.IA3.2人,643,4

0A4I，2J0/c,o

u?)=I0a,。0a,o]z⑺

2015/8/3031

■子空間

子空間（二叮）,0,0為既不可控，又不可觀測的

子空叫（上0,0,a⑻大可整但不可觀測的子空

1川，（Ac,6,Be、。，0）和（Ac,o,Bc,o,Cc,O）

■示意圖

2015/8/3032

4.1.6系統(tǒng)狀態(tài)方程標(biāo)準(zhǔn)型的

MATLAB求解

■常用標(biāo)準(zhǔn)型

■單變量系統(tǒng)的標(biāo)準(zhǔn)型

■MATLAB默認(rèn)的標(biāo)準(zhǔn)型

■可控標(biāo)準(zhǔn)型實(shí)現(xiàn)

■可觀測標(biāo)準(zhǔn)型實(shí)現(xiàn)

■和Jordan標(biāo)準(zhǔn)型實(shí)現(xiàn)

■多變量系統(tǒng)Leunberge標(biāo)準(zhǔn)型

■側(cè)重點(diǎn)：如何用MATLAB直接獲取標(biāo)準(zhǔn)型

2015/8/3033

單變量系統(tǒng)的標(biāo)準(zhǔn)型

■可控標(biāo)準(zhǔn)型0I??.00

00??-00

*??*

*??*

(x=Acx+BciiX-■*??X+:u

\y=C(：x+Dcii00?*?I0

??

一4I-6/2,一Q〃I

■可觀測標(biāo)準(zhǔn)型y=I??E.

???

000-a\h‘I

I0???0一。2b?

0I0

x-Aox+B0uX-一。3X+u

??*

v=Coi+。0〃??

00???I一〃〃Ei

???

V二[0,0\\x

2015/8/3034

■可控可觀測標(biāo)準(zhǔn)型轉(zhuǎn)換

functionGs=sscanform(G,type)

switchtype

case,Ctrl5

G=tf(G);Gs二口；

G.num{l}=G.num{l}/G.den{l}(l);%

G.den{1}=G.den{1}/G.den{1}(1);d=G.num{1}(1);

G1=G;Gl.ioDelay=0;Gl=Gl-d;

num=Gl.num{1};den=Gl.den{1};n=length(G.den{1})-1;

A=[zeros(n-1,1)eye(n-l);-den(end:-1:2)];

B=[zeros(n-1,1);1];C=num(end:-1:2);D=d;

Gs=ss(A,B,C,D,'Ts',G.Ts,'ioDelay',G.ioDelay);

2015/8/3035

case'obsv'

Gc=sscanform(G,'Ctrl5);

Gs=ss(Gc.a',Gc.c',Gc.b',Gc.d','Ts',G.Ts,?..

'ioDelay',G.ioDelay);

otherwise

error('Onlyoptions''Ctrl''and''obsv''.

areapplicable.J)

end

■可控標(biāo)準(zhǔn)型和可觀測標(biāo)準(zhǔn)型，對偶關(guān)系

人二與，8-〃Cc=BlD-從

2015/8/3036

Jordan標(biāo)準(zhǔn)型

■假設(shè)系統(tǒng)矩陣A的狂征根為4],心，…，*，

第i個(gè)特征根力對應(yīng)特征向量為3,則

Avj=AjVj,i=1,2,???,

■矩陣A對應(yīng)的模態(tài)矩陣A定義為

\J\1

A=TlAT=2*■

Jk_

■MATLAB變換[G1,T]二canon(G,'modal')

2015/8/3037

多變量系統(tǒng)的Leunberge標(biāo)準(zhǔn)型

■由可控性判定矩陣

S=僅1，4加,???，4丁均也,???,

4廠2-1偽，???，4々廠1環(huán)）

■構(gòu)造矩陣

一提取此行

一提取此行

2015/8/3038

■得出Leunberge變換矩陣

■編寫leunberge.m函數(shù)

T二leunberge(A,B)

2015/8/3039

■MATLAB函數(shù)清單

functionT=leunberge(A,B)

n=size(A,1);p=size(B,2);S二口；sigmas二口；k=l;

fori=l:p

forj=?:n-l

ifrank(S)=二k,k=k+1;

else,sigmas(i)=j-l;S=S(:,1:end-1);break;end,

end

ifk>n,break;end

end

k=k-l;%如果不是完全可控，則川隨機(jī)數(shù)補(bǔ)足

2015/8/3040

ifk<n

whilerank(S)~=n,S(:,k+l:n)=rand(n,n-k);end

end

L=inv(S);iT=[];

fori=l:p

forj=0:sigmas(i)

iT=[iT;L(i+sum(sigmas(l:i)),:M'Kj];

end,

end

ifk<n,iT(k+1:n,:)=L(k+1:end,:);end

T二inv(iT);%構(gòu)造變換矩陣

2015/8/3041

標(biāo)準(zhǔn)型的變換方法總結(jié)

J5

■可控標(biāo)準(zhǔn)型Gs=sscanform(G,Ctrl)

■可觀測標(biāo)準(zhǔn)型Gs=sscanformCG,JobsvJ)

■Jordan標(biāo)準(zhǔn)型Gs二canon(GJmodal，)

■Leunberge標(biāo)準(zhǔn)型

T=leunberge(A,B),G『ss2ss(G,T)

2015/8/3042

6s4+2s2+8s+10

例4一7G(s)=

54+6"+4s+8

■求解可觀測標(biāo)準(zhǔn)型

X?num=[602810];den=[20648];

G二tf(num,den);Gs=sscanform(G,'obsv')

000

■標(biāo)準(zhǔn)型

■/\100

N⑺=o1o

001

)，(/)二|oool]z(f)+3〃⑺

2015/8/3043

例4-8已知模型

156-12933

4148—422

密⑺=力⑺+

2410-2-2-2

96-121539

?A=[15,6,-12,9;4,14,8,-4;2,4,l?,-2；9,6,-12,15]；

B=[3,3;2,2;-2,-2;3,9];T=leunberge(A,B)

Al=inv(T)*A*T,Bl=inv(T)-B

1861.2010000

—48-79.29-1443()-57.69.610

T=,z(t)N⑺+u(r)

4843.2-2000100

18-46.8900-108240I

2015/8/3044

4.2線性系統(tǒng)時(shí)域響應(yīng)解析解法

■給線性系統(tǒng)一個(gè)激勵(lì)信號(hào)，輸出是什么?

■有兩大類方法

■解析解方法

■求解微分方程、差分方程解析解

■數(shù)值解方法

■主要內(nèi)容

■基于狀態(tài)方程的解析解方法

■基于傳遞函數(shù)部分方式展開的解析解方法

階系統(tǒng)的解析解方法

2015/8/3045

4.2.1基于狀態(tài)方程的解析解方法

■狀態(tài)方程模型(x(t)-Ax(t)+Bu{t}

\2/(0=Ci⑺+Du(t)

■解析解

x(t)=e""-，。)宓(“))+Ic'"—丁)8"(丁)”

Jo

■求解難點(diǎn)ClA(t-T)7D/xj

2015/8/3046

狀態(tài)增廣方法

■消除方矩陣，變成自治系統(tǒng)

■單位階躍信號(hào)〃⑺=1⑺，若假設(shè)有另外

一個(gè)狀態(tài)變量.%+1⑺=〃⑺，則其導(dǎo)數(shù)

為xn+1（0=0

■增廣狀態(tài)方程

宏⑺AB1（。

冊+1。）00入〃+1⑺

■自治系統(tǒng)演。二方（0）可以直接求解析解

2015/8/3047

一般輸入信號(hào)的系統(tǒng)增廣

■一般輸入信號(hào)模型

m

it(t)-〃|(7)+〃2(/)=c/+e"i["2COS(”4，)+“3sin("“)

/=()

■引入增廣狀態(tài)變量

xn+1=cos(d“)

x〃+2-sin"http://)

-，?+3=〃I⑺，.??，?，?+〃z+3=%⑺

2015/8/3048

■增廣狀態(tài)方程模型

而⑺=eAtx(O)

其中

AchBdBIf001

y以/)1(0)

cl\一“4

x〃+W)1

〃44工〃+2。)0

A=010,5?(z)=巧計(jì)3(，),x(0)=(,0

00X〃+4(J)

?

cnr

000?\〃+〃？+3(/)m

■解析解x(t)-c出方(0)

2015/8/3049

■MATLAB實(shí)現(xiàn)函數(shù)

function[Ga,Xa]=ss_augment(G,cc,dd,X)

G二ss(G);Aa=G.a;Ca=G.c;Xa=X;Ba=G.b;D=G.d;

if(length(dd)>0&sum(abs(dd))>le-5),

if(abs(dd(4))>le-5),

Aa=[Aadd(2)';-Ba,dd(3)"Ba;...

zeros(2,length(Aa)),[dd(l),-dd(4);dd(4),dd(l)]];

Ca=[Cadd(2)*Ddd(3)*D];Xa=[Xa;1;0];Ba=[Ba;0;0];

else,

Aa=[Aadd(2)*B;zeros(l,length(Aa))dd(l)];

Ca=[Cadd(2)*D];Xa=[Xa;1];Ba=[B;0];

end

end

2015/8/3050

if(length(cc)>?&sum(abs(cc))>le-5),M=length(cc);

Aa=[AaBazeros(length(Aa),M-1);zeros(M-l,length(Aa)+l)

eye(M-l);zeros(l,length(Aa)+M)];

Ca=[CaDzeros(l,M-l)];Xa=[Xa;cc(l)];ii=l;

fori-2:M,ii=ii*i;Xa(length(Aa)+i)=cc(i)*ii;

end,end

Ga=ss(Aa,zeros(size(CaJ)),Ca,D);

■調(diào)用格式

[G-,胡]/ss_augment(G,c,d,胡)|

■信號(hào)描述

c=[co,5，???且d=\d\,d^d^d^\

2015/8/3051

例4-10連續(xù)系統(tǒng)模型

-19—16—16-191

211617190

x(r)=1⑺+

201716201

-20—16—16-192

(y(r)=|2,1,0,0]x(t)

■初值/(())=2]

■輸入信號(hào)O=2+2e-3fsin(2z)

■求解析解

2015/8/3052

系統(tǒng)增廣

2、?cc=[2];dd=[-3,?,2,2];x?=[0;1;1;2];

A=[-19,-16,-16,-19;21,16,17,19;

20,17,16,20;-20,-16,-16,-19];

B=[1;0;1;2];C=[21??];D二位；G=ss(A,B,C,D);

[Ga,xx?]=ss_augment(G,cc,dd,x?);Ga.a,xx?'

增廣模型o2no

-19-16-16-19o0on

21161719

2n—

20171620

而

方zzo.

/42\2

L((!-

1")二—20—16—16—19\\z-

2/X

II1

0000-3X/

23zX0

()

0000x/-

z/)

fl2

0000k\

2015/8/3053

■解析解求解

3*?symst;y=Ga.c"expm(Ga.a"t)"xxO;

%求解系統(tǒng)的解析解

latex(y);

■解析解求解結(jié)果

127-119-

y(z)=-54+----re1+57e/+----e'+

?48

9_/13577

4re1-------e3tcos(21)+—esin(2

84

■穩(wěn)定性

2015/8/3054

4.2.2基于部分分式展開方法求解

■連續(xù)系統(tǒng)的解析解法

b1V"+b2MLi+...+歷〃$+片+T

G(s)=

n

S+4”，L1+42S，L2+???++an

■輸入信號(hào)的Laplace變換U(s)

■輸出信號(hào)的Laplace變換Y(s)=G(s)U(s)

■無重根時(shí)部分方式展開

y(s)=q+q+???+q

s—Pl5一〃2s_Pm

2015/8/3055

■由Laplace反變換求解析解

/),rit

Y(0=2-[Y(s)]=設(shè)〃"+r2匕〃”+???+rmc

■有重根時(shí)

_2_+叮+1+...十叮+〃1

s-Pj(s—p/)2(s-ppm

■相應(yīng)項(xiàng)的解析解為

rye"/+!/)+]re，/+???+

1(/7?-1)!

11

勺+[。+〃+＞一+("7_1F"cP/

2015/8/3056

■部分分式的MATLAB求解

[笑，p,K]=residue(num,den)

例4-10$3+7sl+3s+4

G(s)--7------------5----------

d+7s3+\7@+17s+6

輸入信號(hào)為階躍信號(hào)R(s)=1/5

■輸出信號(hào)計(jì)算

s'+7/+3s+4

y（s）=

4

+7ky+17$3+17s?+6s

2015/8/3057

■MATLAB求解

"?num=[1734];den=[1717176];

[R,P,K]=residue(num,[den,0]);

[R,P]

■解析解

y(n=2.5833e-3z-9e-2z+5.75e-z-3.5re-z+0.667

■解析解精確值(rat°)

31凸、_23-7-2

v(0=——e/-9e2」z+—e1——tei+-

12423

2015/8/3058

例4-11帶有復(fù)數(shù)極點(diǎn)的系統(tǒng)

s+3

G(S)——

54+2?+11?+18s+18

■階躍響應(yīng)解析解

氫?num=[1,3];den=[12111818];

[r,p,k]=residue(num,[den,0]);[r,p]

■解析解

)、⑺二(0.002+0.0255j)e3jz+(0.002-0.0255j)e-3jr

+(-0.0853+0.0088j)e(-1+j)r

+(-0.0853-0.0088j)e(-1-j)z+0.1667

2015/8/3059

解析解的進(jìn)一步化簡

■基于Euler公式的化簡

(4+〃j)c"廠由⑴"+(4-/%)《"一"""二Ac"sin(3f+。)

其中A=-2V”+，2,0=arctan(-〃/。)

■新MATLAB函數(shù)

［丁，p,K］=pfrac(num,den)

■若P(i)為實(shí)數(shù),則(R⑺,P(i))同residueO

■若P(i)為復(fù)數(shù),則［/?(/),/?(/+1)］對返回A和0

2015/8/3060

新MATLAB函數(shù)清單

function[R,P,K]=pfrac(num,den)

[R,P,K]=residue(num,den);

fori=l:length(R),

ifimag(P(i))>eps

a=real(R(i));b=imag(R(i));

R(i)=-2*sqrt(屋2+b八2);R(i+l)=-atan2(a,b);

elseifabs(imag(P(i)))<eps,R(i)=real(R(i));

end,end

2015/8/3061

例4-12仍考慮

s+3

G(s)-1--------------5--------------

54+2s3+1152+18s+18

■MATLAB求解

?num=[1,3];den=[12111818];

[r,p,k]=pfrac(num,[den,0]);[r,p]

■解析解

y(t)=-0.05llsin(3r-0.0768)-

().1715e-zsin。+1.4677)+().1667

2015/8/3062

基于Laplace變換的求解

■參附錄A

■步驟：

■定義符號(hào)變量

■描述原函數(shù)表達(dá)式

■調(diào)用laplace()函數(shù)或ilaplace()函數(shù)求解

■結(jié)果化簡，如simple()函數(shù)

■求解舉例

2015/8/3063

§3+7s?+3s+4

■例1G(s)

d+7s3+17s2+17s+6

■MATLAB求解

?symss;G=(s，'3+7*s，'2+3*s+4)/...

(sF+7*s八3+17/s3+17*s+6);

y=ilaplace(G/s)

latex(y)

■解析解

771/07\r

v(r)=^+—e—°，+I-7/21+—e-z-9e--z

?312\4

2015/8/3064

■例2°⑸\4+2S3+；1+A+I8

■MATLAB求解

宜0?symss;G=(s+3)/.?.

(s八4+2%人3+11*s-2+18*s+18);

ilaplace(G/s),latex(ans)

■解析解

13

=募cos(31)-sin(3f)+1/6-

255

293

e-zcos⑺-e-zsin(?)

170170

2015/8/3065

離散系統(tǒng)的解析解法

■Z變換

■無重根時(shí)2--^―=—甲

_z-l-PlP\P

■部分分式展開

1777

丫⑴二與+???H--------

+「—P21

-P\一Pm

■解析解

V（/7）=鴛一匕⑵]=上閆、20”———（―

P\\P\!Pl\P2/Pm\Pm

2015/8/3066

■考慮采樣周期

P\\P\!P2\P2)Pm\Pm)

zr(z—1/3)

例4-13('?)—(.]/2)(--1/4)G+1/5)

?D=conv([1-1/2],conv([1-1/4],...

conv([l1/5],[1-11)))；

N=[00conv([l-1/3],[1Q])];

N=N(end:-1:1);D=D(end:-1:1);

[R,P,K]=residue(N,D);[R,P,-R./P]

2015/8/3067

■輸出信號(hào)

7.05473.95063.8095-1.4815

丫⑴

—I—1

■解析解

n

)，(〃)=1.4109+1.4815

■Z變換求解步驟

■定義符號(hào)變量

■調(diào)用iztrans()函數(shù)求解

■化簡

2015/8/3068

■利用符號(hào)運(yùn)算工具箱求解

9大?symsz

G=(z-l/3)/((z-l/2)*(z-l/4)*(z-l/5))

iztrans(G"z/(z-l))

■求解結(jié)果

100204080

)'(/?)=4-(1/5)〃/7+--y(1/2)〃77——(1/4)

■方法更規(guī)范，結(jié)果更簡單

2015/8/3069

有重根問題的解析解

■部分分式表達(dá)式的Z反變換

q(-1)〃%

及I(〃+1)(〃+2)???(〃+/〃-1)

Q—1-p)〃7(/〃一1)!(一〃)"+'"

5^—7

例4J5G(z)=-------二^----

(z-1/2)3仁—1/3)

■部分分式展開

夕???*?D=conv([1-1/2],conv([l-1/2],conv([1,-1/2],...

conv([l,-1/3],[1,-1]))));

N=[0,0,0,5-2,0];

[R,P]=residue(N(end:-1:1),D(end:-1:1));[RP]

2015/8/3070

■部分分式展開

324-240-96192-36

-J—3+L—2+(廠—2)2—2)3+z~l-

■解析解

192/2/1V7

+(〃十1)(/7+2)+36

9

-108(—12/—60〃+72)+36

2015/8/3071

■符號(hào)運(yùn)算求解

11?symsz;G=(5*z-2)/(z-l/2”3/(z-l/3)

iztrans(G*z/(z-l))

■解析解

iV7/1/i\,?

-108-+72--60-77+36-12-n1

\2)\2/

■更直觀，不建議用前者求解，而直接采用Z

變換的符號(hào)運(yùn)算方法求解

2015/8/3072

時(shí)間延遲系統(tǒng)的解析解法

■連續(xù)系統(tǒng)模型G(s)e-八

■求解G(s)的解析解，用一L替代「即可

■離散系統(tǒng)傳遞函數(shù)H(”k

■求解”⑶的解析解，用n-k替代nBPnJ

5z-2

例4/6G(z)r5=

(11/2)3(11/3)

2015/8/3073

■無延遲解析解

)’(〃)二—108口+(；)(-12/72-60/?+72)+36

■有延遲解析解

/]5/15

y(/7)=-108(-+(5)?-12(/7-5)2-

60(/7-5)+72]+36x1(/?-5)

n-5—5

-1083+(-⑵2+60〃+72)

+36x1(77-5)

2015/8/3074

4.2.3二階系統(tǒng)的階躍響應(yīng)及

階躍響應(yīng)指標(biāo)

二咻2

■二階系統(tǒng)模型G°(s)

s(s+243n)

?

■閉環(huán)模型G(s)二

.尸0++必9

■記3d=.]一