华中科技大学电气学院matlab选修课大作业__PI控制器.docx

2008 级

《MATAB语言与控制系统仿真》课程
大作业

姓 名 赖智鹏
学 号 u200811806
所在院系 电气与电子工程学院
班 级 电气0809
日 期 1月16日
同组人员

作业评分 评阅人
设计报告评分表

项 目
应包括的主要内容或考核要点
常见问题
扣分
得分
基本部分
63分
方案论证9分
性能指标分析；控制方法及实现方案
未分析指标
没有说明方案理由
未进行方案比较
有分析指标，未用程序分析
有方案说明，未用程序说明
-5
-1
-1
-1
-1

设计过程35分
控制器设计与参数计算
有设计过程，未运用程序
有简单程序应用，无算法
有程序，没有程序说明
未说明设计理由
-20
-5
-5
-2

结果分析10分
对设计结果的分析与核算，分析原因和改进
没有误差分析
没有对结果验证
有量化指标，但没有分析
有分析验证，没有运用程序
没有总结或总结太虚
-9
-3
-2
-2
-2

格式规范9分
重点考查完整性，图表，公式的规范性
图号图名等问题
没有对图表说明
中英文参杂，不一致
字体不一致
图形截屏
参考文献规范问题
参考文献的引用问题
-2
-2
-1
-1
-1
-1
-1

提高部分
27分
第1项9分
提出改进的性能指标，完成分析，设计并对结果予以验证
未提出其他方案
有，未用程序验证
有程序，未对程序说明
-9
-4
-2

第2项9分
考虑参数变化，干扰影响等其他因素，完成分析，设计并对结果予以验证
未对误差干扰进行分析
有，未用程序验证
有程序，未对程序说明
-9
-4
-2

第3项9分
提出其他更完善的性能指标，完成分析，设计并对结果予以验证
未提出其他指标要求
有，未用程序验证
有程序，未对程序说明
-9
-4
-2

报告得分
90分

合计

特色加分0-30分
报告的特色和难度系数，掌握程度予以评价

总分
报告得分+答辩/特色加分

ＰＩ控制器设计与讨论
1.引言
本文讨论的对象是智能交通/高速公路系统IVHS(习题7.3)，系统通过电子技术为驾驶者和控制系统提供实时路面情况，该系统还提供更方便服务，驾驶者可不用自己操纵驾驶，而系统自动控制车辆，且保持车相对速度，以实现交通的有序进行。
本文前半部分讨论了PI控制器参数的设计问题，根据性能指标推导控制器参数的约束条件，通过相关算法找到问题的解，后半部分讨论系统参数变化时对系统性能的影响，通过比较得出综合性能较优的控制器参数。
2.给定系统的控制器设计
2.1性能分析：
（1）阶跃响应零稳态误差。由原系统为零型系统，故必须通过控制器在原点加入至少一个开环极点，可选择的由PI控制器、PID控制器等，而单独使用比例控制器已达不到目标。

（2）要求系统对爬坡响应的跟踪误差小于25%，故Kv>4，从而有：
。
（3）阶跃响应的超调量小于5%，所以如果校正后系统近似为2阶系统，要求阻尼比>0.7。
（4） 调节时间Ts=4/<1.5s（2%），所以有>4/1.5=2.67。

就目前分析，单独使用比例控制器不能满足要求，而在考虑选择PID控制器之前我们选择相对简单的PI控制器：
。

图 1 原系统
做出根轨迹图并作出满足性能指标的区域：
rlocus(G);
hold on;
plot([-2.66 -2.66],[-20 20]);％指定性能指标在根轨迹图中所在区域
zeta=0.7;
plot([0 -20*zeta],[0 20*sqrt(1-zeta^2)],[0 -20*zeta],[0 -20*sqrt(1-zeta^2)]); ％指定性能指标在根轨迹图中所在区域。

图 2 原系统根轨迹图（黑线表示所要达到的指标边界）
系统波特图及其近似画法：
figure;
sys1=tf([1/16],[1]);%第一段直线近似
sys2=tf([1/8],[1 0]); %第二段直线近似
sys3=tf([1],[1 0 0]); %第三段直线近似
bodemag(sys1,{0.1,2})
hold on;
bodemag(sys2,{2,8})
bodemag(sys3,{8,100})
bode(num,den,{0.1,100})
grid;

图 3 原系统波特图
由波特图可看出原系统开环波特式增益小于1，低频区水平，转折频率为2rad/s和8rad/s，是一个稳定的二阶系统。
下图为其开环阶跃响应和斜坡响应：
subplot(2,1,1);
step(feedback(G,1))
subplot(2,1,2);
t=[0:0.01:5];
lsim(feedback(G,1),t,t)

图 4 原系统响应
阶跃响应稳态误差接近1，斜坡响应不能跟踪输入信号变化，可通过PI控制器增加开环系统类型以消除这两个缺陷。

图 5 校正后系统框图

2.2．参数KP、KI约束条件：
开环传函：

系统存在3个确定的极点在原点（两个自带的极点，一个PI控制器引入的极点），一个变化的零点，增益也是可变的，由性能指标要求即可确定主导极点所在区域。

（1）根轨迹渐近线与实轴的交点小于-2.66，即

得
（2）Kv>4，即

得出
（3）特征方程为
闭环传函

由劳斯判据，有

得

综上所述，对于Kp，Ki有4个约束条件（实为3个），做出取值区间如下图：
syms kp ki
f1='ki-64+0*kp';
f2='kp-ki/10+16';

f3='ki-14/3*kp';
figure
ezplot(f1,[0,160,-20,100]);hold on;
ezplot(f2,[0,160,-20,100]);hold on;
ezplot(f3,[0,160,-20,100]);hold on;
Ki和Kp取值区域：

图 6 Kp，Ki取值区域
图中箭头方向为解区域方向。
2.3考察边界值KI=64,KP=13.8
%通过linmod()函数将方框图转换成状态空间形式，再有状态空间得到传递函数
[a,b,c,d]=linmod('untitled');
sys=tf(ss(a,b,c,d))
得到传函：

13.8 s + 64
----------------------------
s^3 + 4 s^2 + 17.8 s + 64

图 7 系统各环节的比较和对性能的影响
由图7波特图知加入PI校正器增加了系统的类型。

图 8校正后根轨迹
由图8根轨迹看出，此时主导极点不能同时满足阻尼比和阻尼频率的条件。
阶跃和斜坡响应：

图 9 校正后阶跃响应（Ki=64,Kp=13.8）
由图7知超调量20.2%大于5%，稳定时间为2.68s也大于1.5s，故性能不满足要求。

图 10校正后斜坡响应（Ki=64,Kp=13.8）
由上图8知此时已达到斜坡响应的稳态误差要求要求
由代码得到性能指标数值PO、Ts和斜坡响应稳态误差
Kp=13.8;
Ki=64;
s=tf('s');
PI=Kp+Ki/s;
G_PI=series(G,PI);
G_PI_closed=feedback(G_PI,1);
figure
step(G_PI_closed);
[y_s,t_s]=step(G_PI_closed);
yss=dcgain(Gc);
po=(max(y_s)-yss)*100 %超调量PO
m=find(abs(y_s-yss)>0.02);
ts=t_s(length(m)) %调节时间

t_r=[0:0.001:5];
u=t_r;%斜坡信号
figure;
lsim(G_PI_closed,u,t_r);%斜坡响应
y_r=lsim(G_PI_closed,u,t_r);
err=u(length(u))-y_r(length(y_r))
由代码得
po =20.2016
ts =2.3831
err =0.2500
而由图得出的
po=20.2
Ts=2.68
err=0.25
误差分析：二者并不完全一致，其中Ts差别较大0.3/2=11%左右。注意到m数组如下，其中25——28缺省2个数，46——52缺省了6个数，所以将导致此方法计算的Ts有差异：
表格 1 m存储输出与稳态数值差大于3%的点在整个序列中对应的序号
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154

ts=t_s(length(m)+3+6)
>>ts=2.6852
此时就与图中值吻合了
对于一般情况，可用以下代码实现：
m=find(abs(y_s-1)>0.02);
dm=find(abs(y_s(1:(length(m)-1))-1)<0.02);%找到在稳态点之前所有与稳态值误差小于2%的离散点，即曲线穿越1的附近点
ts=t_s(length(m)+length(dm));%补足

2.4 校验系统性能
考察图6区域内Ki=100，Kp=40所对应PI控制器性能是否满足要求，下图为系统阶跃响应和斜坡响应：

图 11 KP=40，KI=100时系统响应
由图知超调量14%>5%，不满足性能要求，所以刚才界定的区域还不足以完全确定满足要求的解。下图得出此时闭环系统的零极点图：

图 12 KP=40，KI=100时，零极点图

由图12知系统主导极点在规定区域，但零点与极点的位置很接近，而零点有增大超调量的作用，下面为去掉该零点之后系统的响应：

图 13 系统去除零点后响应

由图13知此时超调量为0%<5%，满足要求，但此时稳定时间和斜坡响应稳态误差均不满足条件，由此可直观的看到零点对系统的影响。
下图为改变零点位置时阶跃响应，可反映零点位置对系统的影响（下一节通过穷举算法进一步讨论KPKI的取值范围对系统的影响）

图 14零点对性能的影响

2.5 基于穷举算法的控制器设计
下面通过穷举算法，找到一系列满足要求的Kp.Ki值，并将这些离散点与原来界定的解空间画在统一坐标系下：
ki=64;%穷举起点
kp=13.8;% 穷举起点
KPI=zeros(1000,2);%定义1000*2矩阵存储PI控制器参数
PO=zeros(1000,1);%存储PO
Ts=zeros(1000,1);%存储Ts
Err=zeros(1000,1);%存储斜坡响应稳态误差
i=1;%循环变量初值
G=zpk([],[-2,-8],[1]);%原系统
while ki<80 %定义一个边界
kp=3/14*ki; %边界
while kp<80
s=tf('s');
PI=kp+ki/s; %PI控制器
G_PI=series(G,PI);
G_PI_closed=feedback(G_PI,1);%闭环系统

[y_s,t_s]=step(G_PI_closed);%阶跃响应
po=(max(y_s)-1)*100;%PO
m=find(abs(y_s-1)>0.02);%Ts
dm=find(abs(y_s(1:(length(m)-1))-1)<0.02);
ts=t_s(length(m)+length(dm));
t_r=[0:0.001:5];
u=t_r;%
y_r=lsim(G_PI_closed,u,t_r);%
err=u(length(u))-y_r(length(y_r));%
if(po<5&&ts<1.5&&err<0.25)%
KPI(i,1)=ki;
KPI(i,2)=kp;
PO(i)=po;
Ts(i)=ts;
Err(i)=err;
i=i+1;
end;
kp=kp+2;
end;
ki=ki+2;
end;

ＫＰ，ＫＩ满足要求离散点取值（为方便显示，下表为KPKI取整后的值）
round(KPI);
表格 2KP，KI离散取值
KI
64
64
64
64
64
64
64
64
66
66
66
66
66
66
66
66
66
68
KP
34
36
38
40
42
44
46
48
34
36
38
40
42
44
46
48
50
35
KI
68
68
68
68
68
68
68
68
68
70
70
70
70
70
70
70
70
70
KP
37
39
41
43
45
47
49
51
53
37
39
41
43
45
47
49
51
53
KI
70
70
72
72
72
72
72
72
72
72
72
72
74
74
74
74
74
74
KP
55
57
41
43
45
47
49
51
53
55
57
59
46
48
50
52
54
56

由下述代码可做出解空间的分布图：
figure;
plot(KPI(:,1),KPI(:,2),’o’);
syms kp ki
f1='ki-64+0*kp';
f2='kp-ki/10+16'
f3='ki-14/3*kp';
hold on
ezplot(f1,[0,160,-20,100]);hold on;
ezplot(f2,[0,160,-20,100]);hold on;
ezplot(f3,[0,160,-20,100]);

图 15 KP，KI分布图
由图知满足约束条件并不足以规定KP,KI取值。
下图为其中几个PI控制器校正后所得的系统响应：
其中KPI_plot()为自定义函数，输入参数为PI控制器的参数，而做出校正后系统的阶跃响应和斜坡响应，参数1代表只做出阶跃响应曲线，若设置为2则表示只做出斜坡响应，详见附录所给代码：
figure
KPI_plot(KPI(1,:),1);%KPI_plot为自定义函数;
hold on;
KPI_plot(KPI(17,:),1);
hold on;
KPI_plot(KPI(38,:),1);
hold on;
legend('kp=33.7,ki=64','kp=50.1,ki=66','kp=57,ki=70')

图 16 性能校验
对这些点求平均值：
Kp=1000*mean(KPI(:,2))/54%矩阵大小为1000，而包含54个有意义的值，其余为0，下同
Ki=1000*mean(KPI(:,1))/54
得Kp = 45.6878,Ki 68.8889，取Kp=46，Ki=69。
首先得到校正后系统阶跃响应和斜坡响应曲线：

图 17 中心点性能校验

图 18中心点性能校验
性能指标符合要求，这一点不难想象（而在下一节我们将看到选择这一点作为控制器设计的参数所带来的好处，即PI控制器的鲁棒性好），做出零极点图：
pzmap(G_PI_closed);
hold on;
plot([-2.66 -2.66],[-20 20]);
zeta=0.7;
plot([0 -20*zeta],[0 20*sqrt(1-zeta^2)],[0 -20*zeta],[0 -20*sqrt(1-zeta^2)]);
hold off;

图 19系统零极点图
做出各环节波特图并计算相位裕度：
Kp=46;
Ki=69;
G=zpk([],[-2,-8],[1]);
s=tf('s');
PI=Kp+Ki/s;
G_PI=series(G,PI);
G_PI_closed=feedback(G_PI,1);
PI=Kp+Ki/s;
G_PI=series(G,PI);
figure;
margin(PI);
margin(G_PI);
margin(PI);
grid;
[p,m]=margin(G_PI)
得：p =In，fm =65.8188接近70°对于接近阻尼比0.7。

图 20 各环节波特图
由图知校正器增大了系统带宽，从而加快了响应速度，同时也增大了相位差。
2.6 性能指标的改进
如果系统对响应速度的要求高于超调量的要求，比如说调节时间要求小于1s（>4），超调量要求不大于10%（>0.6），仿上述设计流程，得出：

做出该区间：
syms kp ki
f1='ki-64+0*kp';
f2='kp-ki/10+16';
f3='ki-2*kp';
hold on
ezplot(f1,[0,160,-20,100]);hold on;
ezplot(f2,[0,160,-20,100]);hold on;
ezplot(f3,[0,160,-20,100]);hold on;

图 21
同样的方法可以得出一系列满足条件离散点，如下图

图 22
性能指标见下表：

表格 3
KI
64
64
64
69
69
69
74
74
74
79
79
KP
32
37
42
34.5
39.5
44.5
37
42
47
39.5
44.5
PO
4.492051
2.925229
2.383616
5.593106
4.236341
3.606199
6.853009
5.491437
5.07062
7.796231
6.694528
Ts
0.996713
0.726561
0.865214
0.970484
0.766934
0.942691
0.970484
0.781956
0.909712
0.93114
0.784087
Err
0.25
0.249985
0.249923
0.231884
0.231873
0.23183
0.216216
0.216208
0.216177
0.202532
0.202525
KI
79
84
84
84
84
84
89
89
89
94

KP
49.5
42
47
52
57
67
49.5
54.5
59.5
52

PO
6.202628
9.218092
7.833353
7.851369
8.035544
8.782046
8.934857
8.672828
9.27372
9.935864

Ts
0.950096
0.918025
0.773279
0.817075
0.924344
0.996271
0.762198
0.918862
0.915286
0.738163

Err
0.202503
0.190476
0.190471
0.190454
0.190415
0.190231
0.179771
0.179758
0.179728
0.17021

2．7系统参数变化时，系统性能的研究

图 23 参数ab可变系统
对于系统，其参数a,b可能发生改变，此时我们假定PI控制器参数是恒定的，给定Kp=33，Ki=69。其位于KI-KP坐标系中位置如下图，由图可知Kp=33，Ki=69时（Kp，Ki）（红点标记）与a=2,b=8时满足性能指标的解空间的相对位置，即处于解空间的临界区。

图 24 给定PI参数位置

图 25
校验此时系统性能，下图为阶跃响应和斜坡响应，由图可知阶跃响应超调量略大于5%，其余指标满足要求。

图 26 阶跃响应a=2 b=8

图 27 斜坡响应
类似于前面关于Kp、Ki取值范围的分析，分析KPKI不变时，为满足性能指标要求的a、b取值区间。
开环传函

系统存在1个确定的零点（-Ki/Kp）和在原点的极点（0），两个变化的极点（-a和-b），由性能指标要求可确定主导极点所在区域，且

得出
由Kv>4，得出abs a b
kp=33;
ki=69;
f1='a^2*b+a*b^2+46*(a+b)-69';
f2='a+b-5.8';
f3='a*b-17.25'
ezplot(f1,[-100,100,-100,100]);hold on;
ezplot(f2,[-100,100,-100,100]);hold on;
ezplot(f3,[-100,100,-100,100]);hold on;

图 29 约束条件界定区域

图 30 放大后
考虑到a,b之间是等价的，取区间的一半，如下图，箭头方向标注满足要求区域：

图 31 放大后
再次采用穷举法求a，b：
a=2;%初值
b=8;%初值
s=tf('s'); %
kp=33;
ki=69;
PI=kp+ki/s; %PI
AB=zeros(100,2);%存储a,b
PO=zeros(100,1);%
Ts=zeros(100,1);%
Err=zeros(100,1);%
i=1;
while a<4.4 %初值
b=8-a; %下界
while (b<11)&&(a*b-17.25<0) %边界
G=zpk([],[-a,-b],[1]);
G_PI=series(G,PI);
G_PI_closed=feedback(G_PI,1);%系统
[y_s,t_s]=step(G_PI_closed);%阶跃响应
po=(max(y_s)-1)*100;%³¬µ÷Á¿
m=find(abs(y_s-1)>0.02);
dm=find(abs(y_s(1:(length(m)-1))-1)<0.02);%
ts=t_s(length(m)+length(dm));%
t_r=[0:0.001:5];
u=t_r;%бÆÂÐźÅ
y_r=lsim(G_PI_closed,u,t_r);%
err=u(length(u))-y_r(length(y_r));%
if(po<5&&ts<1.5&&err<0.25)%
AB(i,1)=a;
AB(i,2)=b;
PO(i)=po;
Ts(i)=ts;
Err(i)=err;
i=i+1;
end;
b=b+0.1;
end;
a=a+0.1;
end;

由图显示了保持原性能指标的a,b可变化区域，为下图红色区域，对a,和b都采用0.1的步长穷举，只有212组满足条件，如下表。

图 32 a,b分布图
取各分段a,b取值的均值。
mean(AB(:,1))*100/21%数组大小100，有意义值21，其余为0
ans =
2.2762
mean(AB(:,2))*100/21
ans =
7.5286
得(a,b)=(2.2762,7.5286)，作出此时系统阶跃响应曲线：
G=zpk([],[-a,-b],[1]);
G_PI=series(G,PI);
G_PI_closed=feedback(G_PI,1);%系统闭环
step(G_PI_closed);
t_r=[0:0.001:5];
u=t_r;%斜坡信号
lsim(G_PI_closed,u,t_r);%斜坡响应

图 33阶跃响应

图 34 斜坡响应
下面考虑寻找一个使系统对于a,b参数变化更不敏感的PI控制器。
在上一节中（a=2,b=8保持不变时）通过穷举法找到了一系列满足性能指标的PI控制器参数，而且我们可以通过这些离散点确定连续解的取值区间（不慎严密），我们在最后对这些点求平均，得到了一个(Kp,Ki)组合，即Kp = 45.6878,Ki = 68.8889，近似为取为（46,69）。可以看到此点距区间边界近似为所有点中最远的，而这就意味着此时系统性能对Kp,Ki的变化不敏感，所以有理由认为此时闭环系统对原系统参数的变化不敏感。下面将分析，Kp=46,Ki=69时为满足性能要求，a,b的取值限制。
只需修改上述已有的穷举法程序中相关参数即可实现此Kp,Ki下的分析工作。进一步分析得， KpKi的变化只导致了（1）（2）（3）（4）的约束条件中个别参数，将程序中修改部分列出如下
figure;
syms a b
kp=46;
ki=69;
f1='a^2*b+a*b^2+46*(a+b)-69';
f2'a+b-8';
f3='a*b-17.25'
ezplot(f1,[-100,100,-100,100]);hold on;
ezplot(f2,[-100,100,-100,100]);hold on;
ezplot(f3,[-100,100,-100,100]);hold on;
此时，离散点的个数为1912个（求解时，根据a,b的对称性，求出一半，可关于第一象限角平分线镜像），下图显示其具体分布：

图 35 ab分布图

图 36 放大后
显然我们看到此时要求维持性能指标a,b的取值区间更大（离散点数191，而之前的PI控制器只有21）,所以说Kp=46,Ki=69时系统鲁棒性能要好与Kp=33,,Ki=69。实际上由图22知由于Kp=33,,Ki=69。
自此我们分析了控制器参数变化对系统性能的影响，以及原系统参数变化时校正后系统性能的影响，从而从两个方面考察系统，得出了一些结论。
４．结语
在大一时就知道Ｍａｔｌａｂ，是在数学建模中了解到的，当时觉得ｍａｔｌａｂ很神奇，似乎无所不能，那时候也没有真正接触ｍａｔｌａｂ，大二买了电脑后，开始使用ｍａｔｌａｂ，每天别人打游戏的时候我就把ｍａｔｌａｂ当做游戏，看ｄｅｍｏ，看函数，在平时做实验当中时不时地去用ｍａｔｌａｂ处理数据，以及作图，一直也没有系统地学习，只是经验的积累，大致上能熟练的进行操作。
大三选了这门课，学了一个学期学会了一些东西，但我发现并没有达到我预期的目标，我想这跟我不够踏实，而可以去追求“新”的东西，忽略了最基础也最重要的东西，这是很遗憾的。
关于这次课程设计，之前打算做模糊控制器，原因仅仅是：有意思，而且很新颖，而在这之前其实我是没有认真看过其他题目的。听了老师的建议之后回来仔细把题目看了一遍，然后从头到尾认真做了一遍才发现收获还是挺大的，而现在我想如果当初我真的选择模糊控制器作为课程设计对象，报告会是什么样子，在毫无理论基础的情况下，用ｍａｔｌａｂ提供的傻瓜式操作界面做出一个毫无实质内涵的东西是多么可怕，所以我觉得有必要反思一下一直以来的观念。对于一种编程语言一种理论（经典控制和现代控制），本来是没有必然的优劣之分的，而关键在于应用他的人，人的思想是有优劣的，对于学习或其他什么方面，我们最需要的其实不是一种工具，一种理论而是需要自己有自己的思想。
关于答辩的事儿，由于我买的是２２号下午两点半的票，所以我希望能尽早进行答辩，最好能在十一点之前,这是我的联系方式15972140096。
最后谢谢老师半年来的指导。

参考文献：
[1] Robert H.Bishop. 现代控制系统分析与设计. 北京：清华大学，2000
[2] MORRIES DRIEL. 线性控制系统工程. 北京：清华大学，2003
[3] Richard C.Dorf Robert H Shop. 现代控制系统. 北京：高等教育出版社，2001
[4] 胡寿松. 自动控制原理. 北京：科学出版社，2001
附录（程序清单）：
%原系统a,b(2,8)做出根轨迹图并作出满足性能指标的区域：
G=zpk([],[-2,-8],[1]);
rlocus(G);
hold on;
plot([-2.66 -2.66],[-20 20]);％指定性能指标在根轨迹图中所在区域
zeta=0.7;
plot([0 -20*zeta],[0 20*sqrt(1-zeta^2)],[0 -20*zeta],[0 -20*sqrt(1-zeta^2)]); ％指定性能指标在根轨迹图中所在区域。

%原系统波特图及其近似画法：
figure;
sys1=tf([1/16],[1]);%第一段直线近似
sys2=tf([1/8],[1 0]); %第二段直线近似
sys3=tf([1],[1 0 0]); %第三段直线近似
bodemag(sys1,{0.1,2})
hold on;
bodemag(sys2,{2,8})
bodemag(sys3,{8,100})
bode(num,den,{0.1,100})
grid;

%开环阶跃响应和斜坡响应：
subplot(2,1,1);
step(feedback(G,1))
subplot(2,1,2);
t=[0:0.01:5];
lsim(feedback(G,1),t,t)

%Kp，Ki有4个约束条件
syms kp ki
f1='ki-64+0*kp';
f2='kp-ki/10+16';
f3='ki-14/3*kp';
figure
ezplot(f1,[0,160,-20,100]);hold on;
ezplot(f2,[0,160,-20,100]);hold on;
ezplot(f3,[0160,-20,100]);hold on;

%通过linmod()函数将方框图转换成状态空间形式，再有状态空间得到传递函数
[a,b,c,d]=linmod('untitled');
sys=tf(ss(a,b,c,d))

%性能指标数值PO、Ts和斜坡响应稳态误差
Kp=13.8;
Ki=64;
s=tf('s');
PI=Kp+Ki/s;
G_PI=series(G,PI);
G_PI_closed=feedback(G_PI,1);
figure
step(G_PI_closed);
[y_s,t_s]=step(G_PI_closed);%УÕýºóϵͳ½×Ô¾ÏìÓ¦
yss=dcgain(G_PI_closed);
po=(max(y_s)-yss)*100 %³¬µ÷Á¿
m=find(abs(y_s-yss)>0.02);
ts=t_s(length(m)) %Îȶ¨Ê±¼ä
t_r=[0:0.001:5];
u=t_r;%бÆÂÐźÅ
figure;
lsim(G_PI_closed,u,t_r);%бÆÂÏìÓ¦
y_r=lsim(G_PI_closed,u,t_r);%бÆÂÏìÓ¦
err=u(length(u))-y_r(length(y_r))

%程序修改后计算Ts
m=find(abs(y_s-yss)>0.02);
dm=find(abs(y_s(1:(length(m)-1))-yss)<0.02);%找到在稳态点之前所有与稳态值误差小于2%的离散点，即曲线穿越1的附近点
ts=t_s(length(m)+length(dm));%补足

%通过穷举算法，找到一系列满足要求的Kp.Ki值
ki=64;%穷举起点
kp=13.8;% 穷举起点
KPI=zeros(1000,2);%定义1000*2矩阵存储PI控制器参数
PO=zeros(1000,1);%存储PO
Ts=zeros(1000,1);%存储Ts