现代控制工程及测试技术(卓迅佳)

《现代控制工程及测试技术》作业

班级：硕911

姓名：卓迅佳

学号：3109009028

1. 用MATLAB 求解微分方程的不同命令求解如下微分方程。

5.02d d 3

d d 2

d d 2

2

33

=+++y t

y t

y t

y ，2.0)0(,4.0)0(,0)0(-===y y y

1) 至少选用两种求解微分方程的命令；

2) 在同一幅图上，用不同属性、颜色的曲线表示)(t y 和)(t y

； 解：编写m 文件程序如下：

%------第一种方法采用函数ode23或ode45解---- clear all; close all; t0=0; tf=15;

y0=[0,0.4,-0.2]';

[t,y]=ode23('vdpl',t0,tf,y0); figure(1)

plot(t,y(:,1),'g-',t,y(:,2),'r--')

title('用ode23函数实现微分方程的数值解') xlabel('time/sec') ylabel('value') legend('y','y''') grid

%-------第二种方法采用dsolve 函数求解----- t1=0:0.05:15

y=dsolve('D3y+2*D2y+3*Dy+2*y=0.5','y(0)=0,Dy(0)=0.4,D2y(0)=-0.2') s=subs(y ,t1); dy=diff(y); s1=subs(dy ,t1); figure(2)

plot(t1,s,'g-',t1,s1,'r--')

title('用dsolve 函数实现微分方程的符号解') xlabel('time/sec') ylabel('value') legend('y','y''')

grid

在方法一中将高阶微分方程等效表达成一阶微分方程组的程序如下：function xdot=vdpl(t,x)

xdot=zeros(3,1);

xdot(1)=x(2);

xdot(2)=x(3);

xdot(3)=-2.*x(3)-3.*x(2)-2.*x(1)+0.5;

程序运行的结果及输出图形如图1.1，图1.2所示：

y=1/4-3/20*exp(-t)+2/35*7^(1/2)*exp(-1/2*t)*sin(1/2*7^(1/2)*t)-1/10*exp(-1/2*t)*c os(1/2*7^(1/2)*t)

dy=3/20*exp(-t)+3/140*7^(1/2)*exp(-1/2*t)*sin(1/2*7^(1/2)*t)+1/4*exp(-1/2*t)*cos (1/2*7^(1/2)*t)

图1.1 用ode23函数求微分方程的解及解的一阶导数

图1.2 用dsolve 函数及diff 函数求微分方程的解及解的一阶导数

2. 负反馈系统的前向通道和反馈通道传递函数分别为

20

5164)(2

3

++++=

s s

s s s G ；

)

5)(3()1(2)(+++=

s s s s H

1) 求闭环系统的标准传递函数模型，零极点增益模型，状态空间模型；并

将状态空间表达模型转换成可控标准型和可观测标准型。

2) 用传递函数模型求系统的单位阶跃响应；用零极点增益模型求单位斜坡

响应；用状态空间表达模型求单位脉冲响应。

解：编写m 文件程序如下：

%-------------------------第二题---------------------- clear all; close all;

sys1=tf([4,16],[1,1,5,20]); %前向通道传递函数 sys2=zpk([-1],[-3,-5],2); %反馈通道传递函数

%-------生成标准传递函数模型、零极点增益模型、状态空间模型 disp('闭环系统的零极点增益模型为:') zpksys=feedback(sys1,sys2) [num,den]=tfdata(zpksys,'v');

disp('闭环系统的标准传递函数模型为:') tfsys=tf(num,den)

disp('闭环系统的状态空间模型为:') [A,B,C,D]=ssdata(zpksys) abcdsys=ss(zpksys)

%-------下面将状态控制模型转换成可控标准型和可观测标准型 %------- 判断系统是否可控

M=ctrb(A,B);

r1=rank(M);

l1=length(A);

if r1

disp('系统是状态不完全可控的！');

else

disp('系统是状态完全可控的！');

disp('将状态空间模型转换为可控标准型：')

JA=poly(A);

a4=JA(2);a3=JA(3);a2=JA(4);a1=JA(5);a0=JA(6);

W=[a1 a2 a3 a4 1;

a2 a3 a4 1 0;

a3 a4 1 0 0;

a4 1 0 0 0;

1 0 0 0 0];

%计算变换矩阵T

T=M*W;

Ac=inv(T)*A*T

Bc=inv(T)*B

Cc=C*T

Dc=D

end

%---------判断系统是否可观

V=[C' A'*C' A'*A'*C' (A')^3*C' (A')^4*C'];

r2=rank(V);

l2=size(A,1);

if r2

disp('系统是不完全可观的');

else

disp('系统是状态完全可观的');

disp('将状态空间模型转换为可观测标准型：')

%计算变换矩阵Q

Q=inv(W*V');

Ag=inv(Q)*A*Q

Bg=inv(Q)*B

Cg=C*Q

Dg=D

end

%-------求系统的单位阶跃响应，单位斜坡响应，单位脉冲响应

t1=0:0.2:5;

figure(1)

step(tfsys,t1) %传递函数模型求系统的单位阶跃响应title('传递函数模型求系统的单位阶跃响应');

grid

%-------零极点增益模型求单位斜坡响应

%-------转换为求zpksys与1/s乘积的单位阶跃响应

zpk2sys=zpk([],[0],1); %zpk2sys=1/s

G=series(zpksys,zpk2sys);

t2=0:0.2:5;

figure(2)

step(G,t2);

title('零极点增益模型求单位斜坡响应');

grid

%------用状态空间模型求单位脉冲响应

t3=0:0.2:5;

figure(3)

impulse(A,B,C,D,1,t3)

title('状态空间模型求单位脉冲响应');

grid

程序运行的结果如下：

1）闭环系统的零极点增益模型为:

Zero/pole/gain:

4 (s+3) (s+4) (s+5)

-----------------------------------------------------

(s+3.332) (s+5.153) (s+2.029) (s^2 - 1.514s + 9.532)

闭环系统的标准传递函数模型为:

Transfer function:

4 s^3 + 48 s^2 + 188 s + 240

-------------------------------------------

s^5 + 9 s^4 + 28 s^3 + 83 s^2 + 275 s + 332

闭环系统的状态空间模型为:

A =

0.7568 2.9931 -0.7916 1.4741 0

-2.9931 0.7568 -0.3522 0.6558 0

0 0 -5.1533 2.1479 0

0 0 0 -3.3316 2.0000

0 0 0 0 -2.0288

B =

4

C =

3.2269 0 -0.2685 0.5000 0

D =

系统是状态完全可控的！

将状态空间模型转换为可控标准型：

Ac =

-0.0000 1.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 1.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 1.0000 0.0000

0.0000 0.0000 0.0000 0.0000 1.0000 -332.0000 -275.0000 -83.0000 -28.0000 -9.0000 Bc =

0.0000

-0.0000

-0.0000

0.0000

1.0000

Cc =

240.0000 188.0000 48.0000 4.0000 0

Dc =

系统是状态完全可观的

将状态空间模型转换为可观测标准型：

Ag =

0.0000 -0.0000 0.0000 -0.0000 -332.0000

1.0000 -0.0000 0.0000 -0.0000 -275.0000

0.0000 1.0000 0.0000 0.0000 -83.0000

0.0000 -0.0000 1.0000 0.0000 -28.0000

0.0000 -0.0000 0.0000 1.0000 -9.0000 Bg =

240.0000

188.0000

48.0000

4.0000

-0.0000

Cg =

0.0000 -0.0000 0.0000 -0.0000 1.0000 Dg =

2）运行结果如图2.1、图2.2、图2.3所示：

图2.1 传递函数模型求系统的单位阶跃响应

图2.2 零极点增益模型求单位斜坡响应

图2.3 状体空间模型求单位脉冲响应

3. 设反馈系统的开环传递函数为

)

164)(1()

1()()(2

++-+=

s s

s s s K s H s G

1) 绘制系统的根轨迹。

2) 确定系统稳定时的K 值范围。 解：编写m 文件程序如下： %----------第三题------------ clear all close all num=[1 1];

den=conv(conv([1,0],[1,-1]),[1,4,16]); G=tf(num,den)

rlocus(G) %求系统的根轨迹 grid

xlabel('实轴'); ylabel('虚轴'); title('根轨迹图');

%-----下面求根轨迹与虚轴交点的增益及频率 [K,Wcg]=plzy(G)

程序运行结果如图3所示：

图3 系统的根轨迹图

K =

23.3163 35.6837 Wcg =

1.5616

2.5615

从图3中可以看出，使系统稳定的K 值的取值范围为: 23.316335.6873

K <<

此时对应的频率为1.5616 rad/s 和2.5615 rad/s 。

4. 单位反馈系统具有如下的开环传递函数。

)

18)(15)(11.0()

1(10)()(++++=

s s s s s H s G

1) 绘制系统的Nyquist 图，判别系统的稳定性。

2) 绘制系统的Bode 图。确定系统的幅值交界频率，相位交界频率，相位裕

量和幅值裕量。

解：编写m 文件程序如下：

%------------第四题-------------------- clear all close all

num=[10,10];

den=conv(conv([0.1,1],[5,1]),[8,1]); sys=tf(num,den);

%-------绘制Nyquist 图 figure(1)

nyquist(sys,'-b');

title('系统Nyquist 图'); v=[-5 10 -8 8];axis(v);

xlabel('实轴');ylabel('虚轴'); %-------画系统的Bode 图 figure(2) bode(sys) grid;

title('系统Bode 图')

%-------Wcm 为幅值交界频率，Wcg 为相位交界频率，Pm 为相位裕量，Gm 为幅值裕量

[Gm,Pm,Wcg,Wcm]=margin(sys) 运行结果如图4.1，图4.2所示：

图4.1 系统的Nyquist图

图4.2 系统的Bode图

从图4.1可以看出，极坐标图不包围（-1，j0）点，而系统开环传递函数没有右半s平面的极点，因此闭环系统是稳定的。

运行程序可知，系统的幅值交界频率Wcm = 0.5018，相位交界频率Wcg = Inf，相位裕量Pm =59.4929和幅值裕量Gm =Inf。

5. 如题图所示的系统。

题5图

在rltool环境设计一个PID控制器)(

c s

G，使闭环系统一对共轭复数极点位于3

1j

s±

-

=。分析系统给定值阶跃扰动时系统的瞬态响应过程和性能指标；并说

明极点3

1j

s±

-

=能否作为闭环系统的主导极点，为什么？

解：编写m程序如下：

%---------第五题程序----------

clear all

close all

num=1;

den=[1 0 1];

sys=tf(num,den)

产生的受控对象模块sys为：

Transfer function:

1

-------

s^2 + 1

在MATLAB的Command Window窗口键入rltool打开根轨迹设计的GUI窗口，导入受控对象模块sys。首先观察C(s)=K时的根轨迹图, 可知无论如何调整增益

K，根轨迹都不可能经过3

1j

s±

-

=的希望极点，可见单纯的P控制器不能满足控制要求，因此还需要添加零点或极点。在s平面上添加一实数零点，用鼠标拖动零点观察此时根轨迹的变化趋势，当零点移至-1.49后，根轨迹就经过希望的

闭环极点3

1j

s±

-

=，如图5-1所示，此时采用理想微分作用的PD控制器，控

制器的传递函数为()3(0.671)

C s s

=+，闭环系统的零点为-1.49，故系统的闭环传递函数为：

()

G s=。

打开LTI Viewer for SISO Design Tool窗口,显示闭环阶跃响应，如图5.2所示。性能参数：超调量为41.2%，峰值时间为1.06S，调整时间为3.72S。

由于此时的零点为-1.49，极点为3

1j

s±

-

=，在极点附近有闭环零点，故3

1j

s±

-

=不能作为闭环主极点。

图5.1

()C s k

=时的根轨迹图

图5.2 采用PD 控制器的阶跃响应

6. 控制系统具有如下控制对象传递函数

)

3)(2)(1(10

)

()(+++=

s s s s U s Y

定义状态变量

y

x =1，12x

x =，23x x =

利用状态反馈Kx -=u ，把闭环极点配置到

3

221j s +-=，3

222

j s --=，10

3

-=s

试用MATLAB 求所需的状态反馈增益矩阵K 。 解：根据系统的传递函数可以得到系统的线型微分方程为：

u

y y y

y

106116.

)

2()

3(=+++

整理得：???

????+---===u y y y x x x x x 106116)

2(.

.33.

22.

1

写成矩阵向量形式即为：B u Ax X .+=。式中：

????

??????=321x x x X ，??????????---=611

6100

010A ，??

??

?

?????=1000B 。 编写m 程序如下：

%---------------------第六题----------------------- %----依题意可知系统的状态方程模型------- A=[0 1 0;0 0 1;-6 -11 -6]; B=[0;0;10];

P=[-2+j*2*1.732 -2-j*2*1.732 -10]; K=acker(A,B,P)

程序运行结果如下：

K =

15.3993 4.4999 0.8000

7. 题图所示为电站锅炉三冲量给水控制系统典型示意图。图中

D

为蒸汽流量信号；

d

m 为蒸汽流量测量变送装置放大系数，

????

?

?=h

t %

2.0d m 。 W

为给水流量信号； w m 为给水流量测量变送装置放大系数，

???

? ??=h

t %

2.0w m 。

H

为汽包水位信号；

h

m 为汽包水位测量变送装置放大系数，

??

?

??=mm %5.0h m 。

)

111(3.0)

()()(+=

=

s s s W s H s G w 为给水流量对汽包水位传递函数。

s

s s D s H s G d 035.01

87.0)

()()(-

+=

=

为蒸汽流量对汽包水位传递函数。

?

?

? ??=%h t 8.1zf

H

为执行器构传递函数。

)(),(c2c1s G s G 分别为主、副控制器。

4.0,6.0==w d n n 分别为蒸汽、给水流量分流系数。

题7图 锅炉三冲量给水控制系统

试用SIMULINK 及LTI Viewer 对控制器参数进行设计和整定。基本要求如下：

1) 系统的超调量不大于%20；

2) 主控制器选择PI 或PID 控制器；副控制器选择P,PD 或PI 控制器； 3) 选择不同的主、副控制器类型组合，进行仿真和参数阐述整定研究，画出阶跃响应曲线，综合评价比较性能指标。

解：（1） 搭建SIMULINK 下控制系统的仿真模型，如下图7.1所示：

图7.1 锅炉三冲量给水控制的Simulink 模型

（2）按照两步法来整定串级控制系统的参数：将主控制器设置成1

1p K =的工作状

态，用衰减曲线法整定副回路，当副回路75%衰减率时的比例增益2 2.5

p s

K

=,震荡

周期230s

T

s

=。然后置副控制器比例增益2 2.5

p s

K

=,求得当主回路75%衰减率时主控

制器的比例增益11p s

K

=,瞬态响应周期130s T s =。即：

2 2.5p s K =；230s T s =；11p s K =；130s T s

=

① 当副控制器采用P 作用，主控制器采用PI 作用时

2 2.5p K =；10.8310.83p K =?=，110.515i s T T s ==，1

11

0.055

p i i K K T =

=

此时系统单位阶跃响应曲线图如图7.2所示：

图7.2 2 2.5p K =，10.83p K =，10.055i K =时的阶跃响应

此时不满足超调量在20%范围内的要求，对控制器参数进行调整如下：

2 2.5p K =；10.1p K =；10.00001i K =

此时系统的单位阶跃响应如图7.3所示：

图7.3 2 2.5p K =，10.1p K =，10.00001i K =时的阶跃响应

超调量为6%，响应时间为117s ，满足要求。 ②当副控制器采用P 作用，主控制器采用PID 作用时

2 2.5p K =；1 1.251 1.25p K =?=；1

11

0.139p i i K K T =

=；111 3.75p d d K K T =?=。

此时系统的阶跃响应如图7.4所示：

图7.4 2 2.5p K =，1 1.25p K =，10.139i K =，1 3.75d K =时的阶跃响应 由图可知，超调量为57%，不符合要求。经调整后的PID 参数如下：

2 2.5p K =，1 1.1p K =，10.039i K =，19d K =

此时系统的阶跃响应如图7.5所示：

图7.5 2 2.5p K =，1 1.1p K =，10.039i K =，19d K =时的阶跃响应

超调量为17.5%，调整时间为100s ，满足要求。

③当副控制器采用PI 作用，主控制器采用PI 作用时

2 2.50.8

3 2.075p K =?=；22/0.138p i i K K T ==； 10.83p K =；10.055i K =

此时系统的阶跃响应如图7.6所示：

图7.6 2 2.075p K =,20.138i K =,10.83p K =,10.055i K =时的响应

此时系统发散，不符合要求。对控制器参数进行调整如下：

24p K =，20.0001i K =，10.08p K =，10.00001i K =

系统的阶跃响应如图7.7所示：

图7.7 24p K =，20.0001i K =，10.08p K =，10.00001i K =时的响应

此时响应曲线的超调量为7.5%，响应时间为105s 。满足要求。 ④当副控制器采用PI 作用，主控制器采用PID 作用时

2 2.075p K =；20.138i K =；1 1.25p K =；10.139i K =；1 3.75d K =

此时系统的阶跃响应如图7.8所示：

图7.8 2 2.075p K =，20.138i K =，1 1.25p K =，10.139i K =，1 3.75d K =时的响应

由图7.8可以看出，系统超调量超过20%，不符合要求，对控制器参数进行调整如下：

24p K =；20.0001i K =；10.225p K =；10.00001i K =；12d K =

此时系统的阶跃响应如图7.9所示：

图7.9 24p K =，20.0001i K =，10.225p K =，10.00001i K =，12d K =时的响应

此时响应曲线无超调，响应时间为50s 。满足要求。

（3）综合评价性能指标

当主控制器采用PI控制器时，能使系统稳定，经整定参数后，系统响应曲线的超调量较小，满足要求，但调整时间相对较长。

当主控制器采用PID控制器时，能使系统稳定，但超调量较大，进一步调整参数后，可使超调量满足要求，同时，与主控制器为PI控制器相比，系统的调整时间较短。

副控制器采用P控制器时，调整主控制器参数，也能满足要求，此时，调整时间虽然较长一些，但控制系统相对简单。