MATLAB控制系统仿真作业1
一、 控制系统的模型与转换
1. 请将下面的传递函数模型输入到matlab 环境。
]52)1)[(2(24)(322
33++++++=s s s s s s s G )
99.02.0)(1(568
.0)(22+--+=z z z z z H ,T=0.1s >> s=tf('s');
G=(s^3+4*s+2)/(s^3*(s^2+2)*((s^2+1)^3+2*s+5)); G
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
>> num=[1 0 0.56];
den=conv([1 -1],[1 -0.2 0.99]); H=tf(num,den,'Ts',0.1)
Transfer function: z^2 + 0.56 ----------------------------- z^3 - 1.2 z^2 + 1.19 z - 0.99
2. 请将下面的零极点模型输入到matlab 环境。请求出上述模型的零极点,并绘制其位置。
)1)(6)(5()1)(1(8)(22+++-+++=s s s s j s j s s G )
2.8()
6.2)(2.3()(1
511-++=----z z z z z H ,T=0.05s
>>z=[-1-j -1+j]; p=[0 0 -5 -6 -j j]; G=zpk(z,p,8)
Zero/pole/gain: 8 (s^2 + 2s + 2) -------------------------- s^2 (s+5) (s+6) (s^2 + 1)
>>pzmap(G)
Zero/pole/gain:
z^5 (z+0.3125) (z+0.3846) ------------------------- (z-0.122)
Sampling time: 0.05
>>pzmap (H )
二、 线性系统分析
1. 请分析下面传递函数模型的稳定性。
2
21)(23+++=
s s s s G 1
3)50600300(1
3)(22+++++=
s s s s s s G
>> num=[1];
den=[1 2 1 2]; G=tf(num,den); eig(G)' ans =
-2.0000 0.0000 - 1.0000i 0.0000 + 1.0000i
可见,系统有两个特征根在虚轴上,一个特征根在虚轴左侧,所以系统是临界稳定的。
G=tf(num,den); eig(G)' ans =
-1.9152 -0.1414
0.0283 - 0.1073i 0.0283 + 0.1073i
可见,有两个特征根在虚轴右侧,所以系统是不稳定的。
2. 请判定下面离散系统的稳定性。
)
05.025.02.0(2
3)(23+--+-=
z z z z z H
)
34039.804.10215.20368.791
.1576.1112.2)(1
2345
12-++--++=-------z z z z z z z z H
>> num=[-3 2];
den=[1 -0.2 -0.25 0.05]; H=tf(num,den,'Ts',0.1); [eig(H) abs(eig(H))] ans =
-0.5000 0.5000 0.5000 0.5000
0.2000 0.2000
可以看出,由于各个特征根的模均小于1,所以可以判定闭环系统是稳定的。
[eig(H) abs(eig(H))] ans =
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4.1724 4.1724 0.3755 + 0.1814i 0.4170 0.3755 - 0.1814i 0.4170 -0.5292 0.5292 -0.2716 0.2716
0.1193 0.1193
可以看出,由于4.1724这个特征根的模大于1,所以可以判定闭环系统是不稳定的。
3. 设描述系统的传递函数为
40320
10958411812467284224494536546364032018576022208812266436380598251418)(2345678
234567+++++++++++++++=s s s s s s s s s s s s s s s s G ,假定系统
具有零初始状态,请求出单位阶跃响应曲线和单位脉冲响应曲线。
>> num=[18 514 5982 36380 122664 22088 185760 40320]; den=[1 36 546 4536 22449 67284 118124 109584 40320]; G=tf(num,den)
Transfer function:
18 s^7 + 514 s^6 + 5982 s^5 + 36380 s^4 + 122664 s^3 + 22088 s^2 + 185760 s + 40320 -----------------------------------------------------------------------------------------
s^8 + 36 s^7 + 546 s^6 + 4536 s^5 + 22449 s^4 + 67284 s^3 + 118124 s^2 + 109584 s + 40320
>> step(G ,10) >> impulse(G ,10)
单位脉冲响应:
1. 请分析下面传递函数模型阶跃响应。
2
21)(2
3
+++=
s s s s G
利用Simulink 建模,建立系统仿真模型如下:
单击启动仿真按钮,双击示波器得到系统的阶跃响应如下:
2. 请分析下面离散系统的脉冲响应。
13)50600300(1
43)(2
223++++++=s s s s s s s G
利用Simulink 建模,建立系统仿真模型如下:
单击启动仿真按钮,双击示波器得到系统的脉冲响应如下:
3. 对离散采样系统进行分析,并求出其阶跃响应。
G(s)
零阶保持器u(t)
y(t)
-+
其中:2
373
)(2
3++++=
s s s s s G
利用Simulink 建模,建立系统仿真模型如下:
单击启动仿真按钮,双击示波器得到系统的阶跃响应如下:
4. 设计控制器,使得下列系统稳定。
)2.1)(2)(3()3.2)(1()(+-+++=
s s s s s s G
利用Simulink 建模,未连入控制器时,仿真模型和响应如下:
利用Simulink 建模,设计控制器:
从响应输出图形可以看出,连入控制器后系统稳定,性能明显提高。
四、 基于MATLAB 的PID 控制器设计
设计题目:
1. 应用Ziegler — Nichols 算法设计PID 控制器,实现系统的闭环稳定,并比较对各个系统的控制效果。
3
)1s (1G (s)+=
未连入PID 控制器时的系统仿真及其性能指标如下:
可见,未调节时的系统性能有待提高,需设计PID控制器连入。
输入:
>> num=1;
den=conv( [1,1],conv( [1,1],[1,1] ));
Step(num,den);
K=dcgain (num,den)
得出:
K =1
根据图形,得出:
L=1.86 T=4.4
利用自定义的Ziegler_std函数求出Kp、Ti、Td
输入:
>> K=1;
L=1.86;
T=4.4;
[num,den,Kp,Ti,Td]=Ziegler_std (3,[K,L,T])
得出:
num =
2.6400 2.8387 1.5262
den =
根据得出的Kp 、Ti 、Td 值,设计PID 控制器,并利用利用Simulink 仿真建模。 仿真模型及其响应如下:
可见,加入PID 控制器调节后,系统性能明显改善。
5
)1s (1G (s)+=
未连入PID 控制器时的系统仿真及其性能指标如下:
输入:
>> num=1;
den=conv([1,1],conv([1,1],…;
conv( [1,1],conv( [1,1],[1,1] ))));
Step(num,den);
K=dcgain (num,den)
得出:
K = 1
根据图形,得出:
L=3.4 T=6.8
利用自定义的Ziegler_std函数求出Kp、Ti、Td
输入:
>> K=1;
L=3.4;
T=6.8;
[num,den,Kp,Ti,Td]=Ziegler_std (3,[K,L,T])
得出:
num =
4.0800 2.4000 0.7059
den =
1 0
Kp =
2.4000
Ti =
6.8000
Td =
1.7000
根据得出的Kp、Ti、Td值,设计PID控制器,并利用利用Simulink仿真建模。仿真模型及其响应如下:
可见,加入PID控制器调节后,系统性能明显改善。
3)1 s(1
-1.5s G(s)
++
=
利用Simulink建模,未连入控制器时,仿真模型和响应如下:
输入:
>>num=[-1.5 1];
den=conv( [1,1],conv( [1,1],[1,1] ));
Step(num,den);
K=dcgain (num,den)
得出:
K =1
根据图形,得出:
L=1.8 T=5.7
利用自定义的Ziegler_std函数求出Kp、
Ti、Td
输入:
>> K=1;
L=1.8
T=5.7;
[num,den,Kp,Ti,Td]=Ziegler_std (3,[K,L,T])
得出:
num =
3.4200 3.8000 2.1111
den =
10
Kp =
3.8000
Ti =
3.6000
Td =
0.9000
根据得出的Kp、Ti、Td值,设计PID控制器,并利用利用Simulink仿真建模。仿真模型及其响应如下:
可见,加入PID控制器调节后,系统性能明显改善。
五、模糊控制器设计
设计任务:试设计一个模糊控制器,实现对室内温度的控制的模拟。
参考输入:
(1)温度18-40℃范围内分为七个论域,NB NM NS ZE PS PM PB;隶属度函数满足高斯分布;
(2)温度变化率-2 ~2℃范围内分为七个论域,NB NM NS ZE PS PM PB;隶属度函数满足高斯分布;
参考输出:
变频空调输出的控制信号。在一定范围内分为七个论域,NB NM NS ZE PS PM PB,隶属度函数为常数1。
模糊推理过程,output=输入隶属度函数值*输出论域的中心值。
注:本模糊程序采用PAM控制方式的压缩机,则其输出的转速范围为:0~10500转/分。
控制规则:
% % if input is NB and errorinput is PS, then output is NM; % % if input is NB and errorinput is PM, then output is NM; % % if input is NB and errorinput is PB, then output is NS; % % if input is NM and errorinput is NB, then output is NB; % % if input is NM and errorinput is NM, then output is NM; % % if input is NM and errorinput is NS, then output is NM; % % if input is NM and errorinput is ZE, then output is NM; % % if input is NM and errorinput is PS, then output is NM; % % if input is NM and errorinput is PM, then output is NS; % % if input is NM and errorinput is PB, then output is NS; % % if input is NS and errorinput is NB, then output is NM; % % if input is NS and errorinput is NM, then output is NS; % % if input is NS and errorinput is NS, then output is NS; % % if input is NS and errorinput is ZE, then output is NS; % % if input is NS and errorinput is PS, then output is NS; % % if input is NS and errorinput is PM, then output is ZE; % % if input is NS and errorinput is PB, then output is ZE; % % if input is ZE and errorinput is NB, then output is NS; % % if input is ZE and errorinput is NM, then output is ZE; % % if input is ZE and errorinput is NS, then output is ZE; % % if input is ZE and errorinput is ZE, then output is ZE; % % if input is ZE and errorinput is PS, then output is ZE;
% % if input is ZE and errorinput is PM, then output is PS; % % if input is ZE and errorinput is PB, then output is PS;
% % if input is PS and errorinput is NB, then output is ZE; % % if input is PS and errorinput is NM, then output is PS; % % if input is PS and errorinput is NS, then output is PS;
% % if input is PS and errorinput is ZE, then output is PS;
% % if input is PS and errorinput is PS, then output is PS;
% % if input is PS and errorinput is PM, then output is PM; % % if input is PS and errorinput is PB, then output is PM; % % if input is PM and errorinput is NB, then output is PS; % % if input is PM and errorinput is NM, then output is PS; % % if input is PM and errorinput is NS, then output is PM; % % if input is PM and errorinput is ZE, then output is PM; % % if input is PM and errorinput is PS, then output is PM; % % if input is PM and errorinput is PM, then output is PM; % % if input is PM and errorinput is PB, then output is PB; % % if input is PB and errorinput is NB, then output is PS; % % if input is PB and errorinput is NM, then output is PM; % % if input is PB and errorinput is NS, then output is PM; % % if input is PB and errorinput is ZE, then output is PM; % % if input is PB and errorinput is PS, then output is PB;
>> x1 = (18:0.1:40)';
y0 = gaussmf(x1, [1 18]);
y1 = gaussmf(x1, [1 21]);
y2 = gaussmf(x1, [1 25]);
y3 = gaussmf(x1, [1 29]);
y4 = gaussmf(x1, [1 33]);
y5 = gaussmf(x1, [1 37]);
y6 = gaussmf(x1, [1 40]);
plot(x1,[y0 y1 y2 y3 y4 y5 y6])
2.误差图:
程序为:
>> x1 = (-2:0.1:2)';
y0 = gaussmf(x1, [0.5 -2]);
y1 = gaussmf(x1, [0.5 -1.3]); y2 = gaussmf(x1, [0.5 -0.7]); y3 = gaussmf(x1, [0.5 0]);
y4 = gaussmf(x1, [0.5 0.7]);
y5 = gaussmf(x1, [0.5 1.3]);
y6 = gaussmf(x1, [0.5 2]);
plot(x1,[y0 y1 y2 y3 y4 y5 y6])
3.程序为;
x=35;
ex=-0.8;
% define input type in fuzzy zone
y0 = gaussmf(x, [1 18]);
y1 = gaussmf(x, [1 21]);
y2 = gaussmf(x, [1 25]);
y3 = gaussmf(x, [1 29]);
y4 = gaussmf(x, [1 33]);
y5 = gaussmf(x, [1 37]);
y6 = gaussmf(x, [1 40]);
a=[y0 y1 y2 y3 y4 y5 y6];
b=max(a);
% caculate input in fuzzy zone,get input_type and input_authorityvalue if x<=40 & x>=18
if b==a(1)
type='NB';
authorityvalue=y0;
elseif b==a(2)
type='NM';
authorityvalue=y1;
elseif b==a(3)
type='NS';
authorityvalue=y2;
elseif b==a(4)
type='ZE';
authorityvalue=y3;
elseif b==a(5)
type='PS';
authorityvalue=y4;
elseif b==a(6)
type='PM';
authorityvalue=y5;
elseif b==a(7)
type='PB';
authorityvalue=y6;
end
else
if x>40
type='PB';
authorityvalue=1;
elseif x<18
type='NB';
authorityvalue
%error calculate.
ey0 = gaussmf(x, [0.5 -2]);
ey1 = gaussmf(x, [0.5 -1.3]);
ey2 = gaussmf(x, [0.5 -0.7]);
ey3 = gaussmf(x, [0.5 0]);
ey4 = gaussmf(x, [0.5 0.7]);
ey5 = gaussmf(x, [0.5 1.3]);
ey6 = gaussmf(x, [0.5 2]);
a=[ey0 ey1 ey2 ey3 ey4 ey5 ey6];
b=max(a);
% caculate input in fuzzy zone,get input_type and input_authorityvalue if x<=2 & x>=-2
if b==a(1)
etype='NB';
eauthorityvalue=y0;
elseif b==a(2)
etype='NM';
eauthorityvalue=y1;
elseif b==a(3)
etype='NS';
eauthorityvalue=y2;
elseif b==a(4)
etype='ZE';
eauthorityvalue=y3;
elseif b==a(5)
etype='PS';
eauthorityvalue=y4;
elseif b==a(6)
etype='PM';
eauthorityvalue=y5;
elseif b==a(7)
etype='PB';
eauthorityvalue=y6;
end
else
if x>2
etype='PB';
eauthorityvalue=1;
elseif x<-2
etype='NB';
eauthorityvalue=1;
end
if type=='NB'& etype=='NB'
out=authorityvalue*0.001;
% if input is NB and errorinput is NB, then output is NB;
elseif type=='NB'& etype=='NM'
out=authorityvalue*0.001;
% if input is NB and errorinput is NM, then output is NB;
elseif type=='NB'& etype=='NS'
out=authorityvalue*0.001;
% if input is NB and errorinput is NS, then output is NB;
elseif type=='NB'& etype=='ZE'
out=authorityvalue*1750;
% if input is NB and errorinput is ZE, then output is NM;
elseif type=='NB'& etype=='PS'
out=authorityvalue*1750;
% if input is NB and errorinput is PS, then output is NM;
elseif type=='NB'& etype=='PM'
out=authorityvalue*1750;
% if input is NB and errorinput is PM, then output is NM;
elseif type=='NB'& etype=='PB'
out=authorityvalue*3500;
% if input is NB and errorinput is PB, then output is NS;
elseif type=='NM'& etype=='NB'
out=authorityvalue*0.001;
% if input is NM and errorinput is NB, then output is NB;
elseif type=='NM'& etype=='NM'
out=authorityvalue*1750;
% if input is NM and errorinput is NM, then output is NM;
elseif type=='NM'& etype=='NS'
out=authorityvalue*1750;
% if input is NM and errorinput is NS, then output is NM;
elseif type=='NM'& etype=='ZE'
out=authorityvalue*1750;
% if input is NM and errorinput is ZE, then output is NM;
elseif type=='NM'& etype=='PS'
out=authorityvalue*1750;
% if input is NM and errorinput is PS, then output is NM;
elseif type=='NM'& etype=='PM'
out=authorityvalue*3500;
% if input is NM and errorinput is PM, then output is NS;
elseif type=='NM'& etype=='PB'
out=authorityvalue*3500;
% if input is NM and errorinput is PB, then output is NS;
elseif type=='NS'& etype=='NB'