东北大学自动化控制系统计算机辅助设计实验模板.docx

东北大学自动化控制系统计算机辅助设计实验

38

资料内容仅供参考，如有不当或者侵权，请联系本人改正或者删除。

word编辑

11
控制系统计算机辅助设计
——MATLAB语言与应用

控制系统计算机辅助设计
第一部分
1
function dx=lorenzeq(t,x)
dx=[-x(2)-x(3);
x(1)+0.2*x(2);
0.2+(x(1)-5.7)*x(3);]

>> x0=[0;0;0];
>> [t,y]=ode45('lorenzeq',[0,100],x0);
>> plot(t,y)

>> figure;plot3(y(:,1),y(:,2),y(:,3)),grid,

2
function [y,yeq]=f2a(x)
yeq=[];
y=4*x(1)^2+x(2)^2-4;

>> Aeq=[];Beq=[];A=[];B=[];
>> xm=[0;0]; xM=[];x0=[0;0];
>> f1=inline('x(1)^2-2*x(1)+x(2)');
>> [x,f]=fmincon(f1,x0,A,B,Aeq,Beq,xm,xM,'f2a');x',f
ans =
1.0000 0
f =
-1
3
( a)
> >s=tf('s');G=(s^3+4*s+2)/s^3/(s^2+2)/[(s^2+1)^3+2*s+5]
Transfer function:
s^3 + 4 s + 2
------------------------------------------------------
s^11 + 5 s^9 + 9 s^7 + 2 s^6 + 12 s^5 + 4 s^4 + 12 s^3
( b)
>> z=tf('z',0.1);
>> H=(z^2+0.568)/(z-1)/(z^2-0.2*z+0.99)
Transfer function:
z^2 + 0.568
-----------------------------
z^3 - 1.2 z^2 + 1.19 z - 0.99
Sampling time: 0.1

4
>> A=[0 1 0;0 0 1;-5 -4 -13];
>> B=[0;0;2];
>> C=[1 0 0;0 0 0;0 0 0];
>> D=[0];
>> G=ss(A,B,C,D);
>> G
a =
x1 x2 x3
x1 0 1 0
x2 0 0 1
x3 -5 -4 -13
b =
u1
x1 0
x2 0
x3 2
c =
x1 x2 x3
y1 1 0 0
y2 0 0 0
y3 0 0 0
d =
u1
y1 0
y2 0
y3 0
Continuous-time model.
>> G=tf(G)
Transfer function from input to output...
2
#1: ----------------------
s^3 + 13 s^2 + 4 s + 5
#2: 0
#3: 0
>> GG=zpk(G)
Zero/pole/gain from input to output...
2
#1: ----------------------------------
(s+12.72) (s^2 + 0.2836s + 0.3932)
#2: 0
#3: 0
根据微分方程也能够直接写出传递函数模型:
>> num=[2];
>> den=[1,13,4,5];
>> G=tf(num,den);
>> G
Transfer function:
2
----------------------
s^3 + 13 s^2 + 4 s + 5
>> GG=zpk(G)
Zero/pole/gain:
2
----------------------------------
(s+12.72) (s^2 + 0.2836s + 0.3932)

5
>> num=[1,2];
>> den=[1,1,0.16];
>> H=tf(num,den,'Ts',1);
>> H
Transfer function:
z + 2
--------------
z^2 + z + 0.16
Sampling time: 1

6
function H=feedback(G1,G2,key)
if nargin==2;key=-1;end,H=G1/(sym(1)-key*G1*G2);H=simple(H);

>> syms J Kp Ki s;
>> gc=(Kp*s+Ki)/s;
>> g=(s+1)/(J*s^2+2*s+5);
>> gg=feedback(g*gc,1)
>> gg=feedback(g*gc,1)
gg =
((Ki + Kp*s)*(s + 1))/(J*s^3 + (Kp + 2)*s^2 + (Ki + Kp + 5)*s + Ki)

7
( a)
>> s=tf('s');
>> G=(211.87*s+317.64)/(s+20)/(s+94.34)/(s+0.1684);
>> Gc=(169.6*s+400)/s/(s+4);
>> Hs=1/(0.01*s+1);
>> GG=feedback(G*Gc,Hs)
Transfer function:
359.3 s^3 + 3.732e004 s^2 + 1.399e005 s + 127056
-----------------------------------------------------------------
0.01 s^6 + 2.185 s^5 + 142.1 s^4 + 2444 s^3 + 4.389e004 s^2

+ 1.399e005 s + 127056
>> zpk(GG)
Zero/pole/gain:
35933.152 (s+100) (s+2.358) (s+1.499)
--------------------------------------------------------------------------
(s^2 + 3.667s + 3.501) (s^2 + 11.73s + 339.1) (s^2 + 203.1s + 1.07e004)
( b)
>> z=tf('z');
>> G=(35786.7*z^-1+108444)/[(z^-1+4)*(z^-1+20)*(z^-1+74.04)];
>> Gc=1/(z^-1-1);
>> H=1/(0.5*z^-1-1);
>> GG=feedback(G*Gc,H)
Transfer function:
-108444 z^6 + 1.844e004 z^5 + 1.789e004 z^4
------------------------------------------------------------------
1.144e005 z^6 + 2.876e004 z^5 + 274.2 z^4 + 782.4 z^3 + 47.52 z^2

+ 0.5 z
Sampling time: unspecified
>> zpk(GG)
Zero/pole/gain:
-0.94821 z^4 (z-0.5) (z+0.33)
-------------------------------------------------------------
z (z+0.3035) (z+0.04438) (z+0.01355) (z^2 - 0.11z + 0.02396)
Sampling time: unspecified

8
>> s=tf('s');
c1=feedback(1/(s+1)*s/(s^2+2),(4*s+2)/(s+1)^2);
c2=feedback(1/s^2,50);
G=feedback(c1*c2,(s^2+2)/(s^3+14))
Transfer function:
s^6 + 2 s^5 + s^4 + 14 s^3 + 28 s^2 + 14 s
---------------------------------------------------------------------------------
s^10 + 3 s^9 + 55 s^8 + 175 s^7 + 300 s^6 + 1323 s^5 + 2656 s^4 + 3715 s^3
+ 7732 s^2 + 5602 s + 1400

9
>> s=tf('s');
>> G=(s+1)^2*(s^2+2*s+400)/(s+5)^2/(s^2+3*s+100)/(s^2+3*s+2500);
>> G1=c2d(G,0.01)
Transfer function:
4.716e-005 z^5 - 0.0001396 z^4 + 9.596e-005 z^3 + 8.18e
-005 z^2 - 0.0001289 z + 4.355e-005
-----------------------------------------------------------------
z^6 - 5.592 z^5 + 13.26 z^4 - 17.06 z^3 + 12.58 z^2 - 5.032 z
+ 0.8521
Sampling time: 0.01
>>step(G1)

>> G2=c2d(G,0.1)
Transfer function:
0.0003982 z^5 - 0.0003919 z^4 - 0.000336 z^3 + 0.0007842 z^2
- 0.000766 z + 0.0003214
-----------------------------------------------------------------
z^6 - 2.644 z^5 + 4.044 z^4 - 3.94 z^3 + 2.549 z^2 - 1.056 z
+ 0.
Sampling time: 0.1
>>step(G2)

>> G3=c2d(G,1)
Transfer function:
8.625e-005 z^5 - 4.48e-005 z^4 + 6.545e-006 z^3 + 1.211e-005 z^2
- 3.299e-006 z + 1.011e-007
-----------------------------------------------------------------
z^6 - 0.0419 z^5 - 0.07092 z^4 - 0.0004549 z^3 + 0.002495 z^2
- 3.347e-005 z + 1.125e-007

Sampling time: 1
>> step(G3)

10
( a)
>> G=1/(s^3+2*s^2+s+2);
>>pzmap(G)

系统极点均在虚轴左侧, 系统稳定

( b)
>> G=1/(6*s^4+3*s^3+2*s^2+s+1);
>> pzmap(G)

系统极点在虚轴右侧侧, 系统不稳定
( c)
>> G=1/(s^4+s^3-3*s^2-s+2);
>> pzmap(G)

系统极点在虚轴右侧侧, 系统不稳定
11
( a)
>> z=tf('z',0.1);
>> H=(-3*z+2)/(z^3-0.2*z^2-0.25*z+0.05);
>> pzmap(H)

系统所有极点均在单位圆内, 因此系统稳定
( b)
>> z=tf('z',0.1);
>> H=(3*z^2-0.39*z-0.09)/(z^4-1.7*z^3+1.04*z^2+0.268*z+0.024);
>> pzmap(H)

系统所有极点不全单位圆内, 因此系统不稳定

12
>> A=[-0.2 0.5 0 0 0;0 -0.5 1.6 0 0;0 0 -14.3 85.8 0;0 0 0 -33.3 100;0 0 0 0 -10];
>> B=[0;0;0;0;30];
>> C=[0 0 0 0 0];
>> G=ss(A,B,C,0)

a =
x1 x2 x3 x4 x5
x1 -0.2 0.5 0 0 0
x2 0 -0.5 1.6 0 0
x3 0 0 -14.3 85.8 0
x4 0 0 0 -33.3 100
x5 0 0 0 0 -10
b =
u1
x1 0
x2 0
x3 0
x4 0
x5 30
c =
x1 x2 x3 x4 x5
y1 0 0 0 0 0
d =
u1
y1 0
Continuous-time model.
>> pzmap(G)

系统所有极点均在虚轴左侧, 因此系统稳定
或
>> eig(G)
ans =
-0.
-0.5000
-14.3000
-33.3000
-10.0000
极点均在左半轴, 因此系统稳定

13
数值解:
>> f=@(t,x)[-5*x(1)+2*x(2);-4*x(2);-3*x(1)+2*x(2)-4*x(3)-x(4);-3*x(1)+2*x(2)-4*x(4)];
>> t_final=10;
>> x0=[1 2 0 0];
>> [t,x]=ode45(f,[0,t_final],x0);
>> plot(t,x)

解析解
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)>1e-5)),
if (abs(dd(4))>1e-5),
Aa=[Aa dd(2)*Ba,dd(3)*Ba;...
zeros(2,length(Aa)),[dd(1),-dd(4);dd(4),dd(1)]];
Ca=[Ca dd(2)*D dd(3)*D];Xa=[Xa;1;0];Ba=[Ba;0;0];
else,
Aa=[Aa dd(2)*B;zeros(1,length(Aa)) dd(1)];
Ca=[Ca dd(2)*D];Xa=[Xa;1];Ba=[B;0];
end
end
if (length(cc)>0&sum(abs(cc))>1e-5),M=length(cc);
Aa=[Aa Ba zeros(length(Aa),M-1);zeros(M-1,length(Aa)+1)...
eye(M-1);zeros(1,length(Aa)+M)];
Ca=[Ca D zeros(1,M-1)];Xa=[Xa;cc(1)];ii=1;
for i=2:M,ii=ii*i;Xa(length(Aa)+i)=cc(i)*ii;
end,end
Ga=ss(Aa,zeros(size(Ca')),Ca,D);

>> cc=[2];dd=[-3,0,2,2];x0=[1;2;0;1];
>> A=[-5,2,0,0;0,-4,0,0;-3,2,-4,-1;-3,2,0,-4];B=[0;0;0;1];C=[2 1 0 0];D=0;
>> G=ss(A,B,C,D);
>> [Ga,xx0]=ss_augment(G,cc,dd,x0);Ga.a,xx0'
ans =
-5 2 0 0 0 0 0
0 -4 0 0 0 0 0
-3 2 -4 -1 0 0 0
-3 2 0 -4 0 2 1
0 0 0 0 -3 -2 0
0 0 0 0 2 -3 0
0 0 0 0 0 0 0
ans =
1 2 0 1 1 0 2
>> syms t;y=Ga.c*expm(Ga.a*t)*xx0;
>> latex(y);
>> latex(y)
ans =
-6\,{e^{-5\,t}}+10\,{e^{-4\,t}}

14
( a)
>> s=tf('s');
>> g=(s+6)*(s-6)/s/(s+3)/(s+4+4i)/(s+4-4i);
>> rlocus(g)

( b)
>> s=tf('s');
>> G=(s^2+2*s+2)/(s^4+s^3+14*s^2+8*s);
>> rlocus(G)

K<5.62

15
>> s=tf('s');
>> G=(s-1)/(s+1)^5;
>> G.ioDelay=2
Transfer function:
s - 1
exp(-2*s) * --------------------------------------
s^5 + 5 s^4 + 10 s^3 + 10 s^2
+ 5 s + 1
>> rlocus(pade(G,2))

16
( a)
>> s=tf('s');
>> G=8*(s+1)/s^2/(s+15)/(s^2+6*s+10);
>> nyquist(G),
>> nyquist(G),grid

>> bode(G)

>> figure;nichols(G),grid

>> [gm,pm,wg,wp]=margin(G)
gm =
30.4686
pm =
4.2340
wg =
1.5811
wp =
0.2336
>> GG=feedback(G,1)
Transfer function:
8 s + 8
------------------------------------------
s^5 + 21 s^4 + 100 s^3 + 150 s^2 + 8 s + 8
>> pzmap(GG)

经过以上的分析, 能够看出该系统是稳定的。
采用阶跃响应进行验证如下图:
>> pzmap(GG)
>> step(GG)

( b)
>> Z=[-1.31;-0.054;0.957];
>> P=[0;1;0.368;0.99];
>> G=zpk(Z,P,0.45,'Ts',0.1)
Zero/pole/gain:
0.45 (z+1.31) (z-0.957) (z+0.054)
---------------------------------
z (z-1) (z-0.99) (z-0.368)
Sampling time: 0.1
>> nyquist(G);grid

>> bode(G)

>> nichols(G),grid

>> [gm,pm,wg,wp]=margin(G)
gm =
0.9578
pm =
-1.7660
wg =
10.4641
wp =
10.7340
>> GG=feedback(G,1);
>> pzmap(GG)

>> step(GG)

17
>> z=[-2.5];p=[0;-0.5;-50];
>> G=zpk(z,p,100/2.5*0.5*50);
>> z=[-1;-2.5];p=[-0.5;-50]
>>Gc=zpk(z,p,1000);
>> GG=feedback(G*Gc,1)

Zero/pole/gain:
1000000 (s+1) (s+2.5)^2
------------------------------------------------------
(s+1) (s^2 + 4.99s + 6.239) (s^2 + 95.01s + 1.002e006)
>> bode(GG)

经过bode图能够看出该系统是稳定的。
验证如下:
>> step(GG)

第二部分
2

>> Y=dsolve('D4y+5*D3y+6*D2y+4*Dy+2*y=exp(-3*t)+exp(-5*t)*sin(4*t+pi/3)',...)
'y(0)=1','Dy(0)=1/2','D2y(0)=1/2','D3y(0)=0.2');
>> latex(Y)

3
Simulink仿真图形:

系统输出曲线:

系统误差曲线:

4
x1(t)=sin(x2(t)*exp(-2.3*x4(t)))
x2(t)=x1(t)
x3=sin(x2*exp(-2.3*x4))
x4=x3

5
Simulink仿真图形:

阶跃响应曲线:

6
>> s=tf('s');
>>g=210*(s+1.5)/((s+1.75)*(s+16)*(s+1.5+3*j)*(s+1.5-3*j));
>> gc=52.5*(s+1.5)/(s+14.86);
>> step(feedback(g,1))

>> step(feedback(g*gc,1))

经过两个图比较能够看出, 原系统超调量过大, 震荡严重, 加入控制器后, 系统变得不稳定。

7
>> A=[0 1 0 0;0 0 1 0;-3 1 2 3;2 1 0 0];
>> B=[1 0;2 1;3 2;4 3];
>> Q=diag([1,2,3,4]);
>> R=eye(2);
>> [K,S]=lqr(A,B,Q,R)
K =
-0.0978 1.2118 1.8767 0.7871
-3.8819 -0.4668 2.6713 1.0320
S =
5.4400 0.6152 -2.3163 0.0452
0.6152 1.8354 -0.0138 -0.7582
-2.3163 -0.0138 1.9214 -0.3859
0.0452 -0.7582 -0.3859 0.8540
>> eig(A-B*K)
ans =
-12.2563
-1.6786 + 0.9981i
-1.6786 - 0.9981i
-1.4627
>> step(ss(A-B*K,B,eye(4),zeros(4,2)))

8
>> A=[-0.2 0.5 0 0 0;0 -0.5 1.6 0 0;0 0 -14.3 85.8 0;0 0 0 -33.3 100;0 0 0 0 -10];
>> B=[0;0;0;0;30];
>> C=[1 0 0 0 0];
>> D=0;
>> G=ss(A,B,C,D);
>> pole(G)
ans =
-0.
-0.5000
-14.3000
-33.3000
-10.0000
>> zero(G)
ans =
Empty matrix: 0-by-1
>> P=[-1 -2 -3 -4 -5];
>> K=place(A,B,P)
K =
Columns 1 through 4
0.0004 0.0004 -0.0035 0.3946
Column 5
-1.4433
>> eig(A-B*K)
ans =
-5.0000
-4.0000
-3.0000
-2.0000
-1.0000
9
>> open_system(simulink)

PID控制参数如下:
[0.40725 0.19096 0.21986]
7.9397 0.9011 78.2377