华中科技大学电气学院matlab选修课大作业-PI控制器

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

器

2008级

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

大作业

姓名赖智鹏

学号u200811806

所在院系电气与电子工程学院

班级电气0809

日期1月16日

同组人员

作业评分评阅人

设计报告评分表

ＰＩ控制器设计与讨论

1.引言

本文讨论的对象是智能交通/高速公路系统IVHS(习题7.3)，系统通过电子技术为驾驶者和控制系统提供实时路面情况，该系统还提供更方便服务，驾驶者可不用自己操纵驾驶，而系统自动控制车辆，且保持车相对速度，以实现交通的有序进行。

本文前半部分讨论了PI控制器参数的设计问题，根据性能指标推导控制器参数的约束条件，通过相关算法找到问题的解，后半部分讨论系统参数变化时对系统性能的影响，通过比较得出综合性能较优的控制器参数。

2.给定系统的控制器设计

2.1性能分析：

（1）阶跃响应零稳态误差。由原系统为零型系统，故必须通过控制器在原点加入至少一个开环极点，可选择的由PI控制器、PID控制器等，而单独使用比例控制器已达不到目标。（2）要求系统对爬坡响应的跟踪误差小于25%，故Kv>4，从而有：

001

4

16ssKvsGcGsGc->->==>∑∑。（3）阶跃响应的超调量小于5%，所以如果校正后系统近似为2阶系统，要求阻尼比ξ>0.7。

（4）调节时间Ts=4/nξω4/1.5=2.67。

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

i

PK

GcKS=+。

图1原系统

做出根轨迹图并作出满足性能指标的区域：

rlocus(G);

holdon;

plot([-2.66-2.66],[-2020]);％指定性能指标在根轨迹图中所在区域zeta=0.7;

plot([0-20*zeta],[020*sqrt(1-zeta^2)],[0-20*zeta],[0-20*sqrt(1-zeta^2)]);％指定性能指标在根轨迹图中所在区域。

-20-15

-10

-5

510

15

20

RootLocus

RealAxis

ImaginaryAxis

图2原系统根轨迹图（黑线表示所要达到的指标边界）

系统波特图及其近似画法：figure;

sys1=tf([1/16],[1]);%第一段直线近似sys2=tf([1/8],[10]);%第二段直线近似sys3=tf([1],[100]);%第三段直线近似bodemag(sys1,{0.1,2})holdon;

bodemag(sys2,{2,8})bodemag(sys3,{8,100})bode(num,den,{0.1,100})grid;

-100-80

-60

-40

-20

Magnitude(dB)10

10

10

10

-180

-135-90-450Phase(deg)

BodeDiagram

Frequency(rad/sec)

图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)

0.511.52

2.53

3.5

4

4.55

2

4

6

System:untitled1Time(sec):3.79Amplitude:0.183

Input:In(1)

Time(sec):3.79Amplitude:3.79

Input:In(1)Time(sec):4.6Amplitude:4.6

System:untitled1Time(sec):4.6Amplitude:0.2280

0.02

0.04

0.06

StepResponse

Time(sec)

Amplitude

LinearSimulationResults

Time(sec)

Amplitude

图4原系统响应

阶跃响应稳态误差接近1，斜坡响应不能跟踪输入信号变化，可通过PI控制器增加开环系统

类型以消除这两个缺陷。

图5校正后系统框图

2.2．参数KP、KI约束条件：

开环传函：

(2)(8)

KpsKi

GHsss?+=

++

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

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

280(/)5/22.662

pzKpKiKpKimnα+--=

==-+4，即

1

416

s

KvsGcGKi->=

=

>∑得出Ki64>

（3）特征方程为32

10(16)

0ssKpsKi++++=闭环传函

3210(16)KpsKi

ssKpsKi?+++++

由劳斯判据，有

3210

116010

0160100010

0sKpsKi

KpKi

ssKi

++-

得

160100KpKi+->0Ki>

综上所述，对于Kp，Ki有4个约束条件（实为3个），做出取值区间如下图：

symskpki

f1='ki-64+0*kp';f2='kp-ki/10+16';f3='ki-14/3*kp';figure

ezplot(f1,[0,160,-20,100]);holdon;ezplot(f2,[0,160,-20,100]);holdon;ezplot(f3,[0,160,-20,100]);holdon;Ki和Kp取值区域：

-200

20

40

60

80

100

ki

kp

ki-14/3kp=0

图6Kp，Ki取值区域

图中箭头方向为解区域方向。

2.3考察边界值KI=64,KP=1

3.8

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

13.8s+64

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

BodeDiagramGm=Inf,Pm=Inf

Frequency(rad/sec)

-150-100-50050

100

Magnitude(dB)

10

10

10

10

10

Phase(deg)

图7系统各环节的比较和对性能的影响

由图7波特图知加入PI校正器增加了系统的类型。

-25-20

-15-10-5

0510

152025RootLocus

RealAxis

ImaginaryAxis

图8校正后根轨迹

由图8根轨迹看出，此时主导极点不能同时满足阻尼比和阻尼频率的条件。

阶跃和斜坡响应：

Time(sec)

Amplitude

0.2

0.4

0.6

0.8

1

1.2

1.4

图9校正后阶跃响应（Ki=64,Kp=13.8）

由图7知超调量20.2%大于5%，稳定时间为2.68s也大于1.5s，故性能不满足要求。

0.511.522.53

3.54

4.55Time(sec)

Amplitude

图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%超调量POm=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有差异：

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)5%，不满足性能要求，所以刚才界定的区域还不足以完全确定满足要求的解。下图得出此时闭环系统的零极点图：

Pole-ZeroMap

RealAxis

ImaginaryAxis

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

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

Time(sec)

Amplitude

LinearSimulationResults

Time(sec)

Amplitude

图13系统去除零点后响应

由图13知此时超调量为0%4），超调

量要求不大于10%（ξ>0.6），仿上述设计流程，得出：

Ki/Kp2

160100KpKi+->

做出该区间：symskpki

f1='ki-64+0*kp';f2='kp-ki/10+16';f3='ki-2*kp';holdon

ezplot(f1,[0,160,-20,100]);holdon;ezplot(f2,[0,160,-20,100]);holdon;ezplot(f3,[0,160,-20,100]);holdon;

ki

kp

ki-2kp=0

图21

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

20

40

60

80

100120

140

160

180

200

020406080100120

140160180200ki

kpfigure

ki-2kpfigure=0

图22

性能指标见下表：

r

K

I

79848484848489898994

K

P

49.5424752576749.554.559.552

P

O

6.202

628

9.218

092

7.833

353

7.851

369

8.035

544

8.782

046

8.934

857

8.672

828

9.273

72

9.935

864

T

s

0.950

096

0.918

025

0.773

279

0.817

075

0.924

344

0.996

271

0.762

198

0.918

862

0.915

286

0.738

163

E

r

r

0.202

503

0.190

476

0.190

471

0.190

454

0.190

415

0.190

231

0.179

771

0.179

758

0.179

728

0.170

21

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

图23参数ab可变系统

对于系统

1

()()

G

sasb

=

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

-20

20

40

60

80

100

ki

k

p

ki-14/3kp=0

图24给定PI参数位置

图25

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

StepResponse

Time(sec)

Amplitude

0.2

0.4

0.6

0.8

1

1.2

1.4

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

LinearSimulationResults

Time(sec)

Amplitude

0.511.522.53

3.54

4.55

图27斜坡响应

类似于前面关于Kp、Ki取值范围的分析，分析KPKI不变时，为满足性能指标要求的a、b取值区间。

开环传函

,33,69

()()KpsKi

GHKpKissasb?+===++

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

0(/)0.24()/22.662pzabKpKiabmnα+--=

==-+

由Kv>4，得出

ab

0ab+>

综上所述，对于a,b有四个约束条件，做出a,b取值区间如下图：figure;symsabkp=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]);holdon;ezplot(f2,[-100,100,-100,100]);holdon;ezplot(f3,[-100,100,-100,100]);holdon;

-100

-80

-60

-40

-20

020

40

60

80

100

a

b

ab-17.25=0

图29约束条件界定区域

a

b

ab-17.25=0

图30放大后

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

a

b

ab-17.25=0

图31放大后

再次采用穷举法求a，b：a=2;%初值b=8;%初值s=tf('s');%kp=33;ki=69;

PI=kp+ki/s;%PI

AB=zeros(100,2);%存储a,bPO=zeros(100,1);%Ts=zeros(100,1);%Err=zeros(100,1);%i=1;

whilea0.02);

dm=find(abs(y_s(1:(length(m)-1))-1)0.02);

ts=t_s(length(m))%?��?

t_r=[0:0.001:5];

u=t_r;%D±??D?o?

figure;

lsim(G_PI_closed,u,t_r);%D±???ìó|

y_r=lsim(G_PI_closed,u,t_r);%D±???ìó|

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);%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

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解空间的分布图：

figure;

plot(KPI(:,1),KPI(:,2),'o');

symskpki

f1='ki-64+0*kp';

f2='kp-ki/10+16';

f3='ki-14/3*kp';

holdon

ezplot(f1,[0,160,-20,100]);holdon;

ezplot(f2,[0,160,-20,100]);holdon;

ezplot(f3,[0,160,-20,100]);

%其中KPI_plot()为自定义函数，输入参数为PI控制器的参数，而做出校正后系统的阶跃响应和斜坡响应，参数1代表只做出阶跃响应曲线，若设置为2则表示只做出斜坡响应，详见附录所给代码：

figure;

KPI_plot(KPI(1,:),1);%KPI_plot为自定义函数;

holdon;

KPI_plot(KPI(17,:),1);

holdon;

KPI_plot(KPI(38,:),1);

holdon;

legend('kp=33.7,ki=64','kp=50.1,ki=66','kp=57,ki=70')

%对这些点求平均值：

Kp=1000*mean(KPI(:,2))/54%矩阵大小为1000，而包含54个有意义的值，其余为0，下同Ki=1000*mean(KPI(:,1))/54

%零极点图与界定区域

pzmap(G_PI_closed);

holdon;

plot([-2.66-2.66],[-2020]);

zeta=0.7;

plot([0-20*zeta],[020*sqrt(1-zeta^2)],[0-20*zeta],[0-20*sqrt(1-zeta^2)]);

holdoff;

%做出各环节波特图并计算相位裕度：

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;