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

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%，故Kv4，从而有：。（3）阶跃响应的超调量小于5%，所以如果校正后系统近似为2阶系统，要求阻尼比0.7。（4） 调节时间Ts=4/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-zeta2),0 -20*zeta,0 -20*sqrt(1-zeta2); 指定性能指标在根轨迹图中所在区域。图 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）Kv4，即得出（3）特征方程为闭环传函由劳斯判据，有 得综上所述，对于Kp，Ki有4个约束条件（实为3个），做出取值区间如下图：syms kp kif1=ki-64+0*kp;f2=kp-ki/10+16;f3=ki-14/3*kp;figureezplot(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-s3 + 4 s2 + 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 %超调量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.2016ts =2.3831err =0.2500而由图得出的po=20.2 Ts=2.68 err=0.25误差分析：二者并不完全一致，其中Ts差别较大0.3/2=11%左右。注意到m数组如下，其中2528缺省2个数，4652缺省了6个数，所以将导致此方法计算的Ts有差异：表格 1 m存储输出与稳态数值差大于3%的点在整个序列中对应的序号1234567891011121314151617181920212223242528293031323334353637383940414243444546525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154ts=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%，不满足性能要求，所以刚才界定的区域还不足以完全确定满足要求的解。下图得出此时闭环系统的零极点图：图 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);%存储POTs=zeros(1000,1);%存储TsErr=zeros(1000,1);%存储斜坡响应稳态误差i=1;%循环变量初值G=zpk(,-2,-8,1);%原系统while ki80 %定义一个边界 kp=3/14*ki; %边界 while kp0.02);%Tsdm=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(po5&ts1.5&err4），超调量要求不大于10%（0.6），仿上述设计流程，得出：做出该区间：syms kp kif1=ki-64+0*kp;f2=kp-ki/10+16;f3=ki-2*kp;hold onezplot(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性能指标见下表：表格 3KI6464646969697474747979KP32374234.539.544.537424739.544.5PO4.4920512.9252292.3836165.5931064.2363413.6061996.8530095.4914375.070627.7962316.694528Ts0.9967130.7265610.8652140.9704840.7669340.9426910.9704840.7819560.9097120.931140.784087Err0.250.2499850.2499230.2318840.2318730.231830.2162160.2162080.2161770.2025320.202525KI79848484848489898994KP49.5424752576749.554.559.552PO6.2026289.2180927.8333537.8513698.0355448.7820468.9348578.6728289.273729.935864Ts0.9500960.9180250.7732790.8170750.9243440.9962710.7621980.9188620.9152860.738163Err0.2025030.1904760.1904710.1904540.1904150.1902310.1797710.1797580.1797280.1702127系统参数变化时，系统性能的研究图 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），由性能指标要求可确定主导极点所在区域，且得出由Kv4，得出abKi/4=17.25图 28 根轨迹特征方程：闭环传函：由劳斯判据得：综上所述，对于a,b有四个约束条件，做出a,b取值区间如下图：figure;syms a bkp=33;ki=69;f1=a2*b+a*b2+46*(a+b)-69;f2=a+b-5.8;f3=a*b-17.25ezplot(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; %PIAB=zeros(100,2);%存储a,bPO=zeros(100,1);%Ts=zeros(100,1);%Err=zeros(100,1);%i=1;while a4.4 %初值 b=8-a; %下界 while (b11)&(a*b-17.250.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(po5&ts1.5&err0.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)%程序修改后计算Tsm=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);%存储POTs=zeros(1000,1);%存储TsErr=zeros(1000,1);%存储斜坡响应稳态误差i=1;%循环变量初值G=zpk(,-2,-8,1);%原系统while ki80 %定义一个边界 kp=3/14*ki; %边界 while kp0.02);%Tsdm=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(po5&ts1.5&err0.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解空间的分布图：figure;plot(KPI(:,1),KPI(:,2),o);syms kp kif1=ki-64+0*kp;f2=kp-ki/10+16;f3=ki-14/3*kp;hold onezplot(f1,0,160,-20,100);hold on;ezplot(f2,0,160,-20,100);hold on;ezplot(f3,0,160,-20,100);%其中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)%对这些点求平均值：Kp=1000*mean(KPI(:,2)/54%矩阵大小为1000，而包含54个有意义的值，其余为0，下同Ki=1000*mean(KPI(:,1)/54%零极点图与界定区域pzmap(G_PI_closed);hold on;plot(-2.66 -2.66,-20 20);zeta=0.7;plot(0 -20*zeta,0 20*sqrt(1-zeta2),0 -20*zeta,0 -20*sqrt(1-zeta2);hold off;%做出各环节波特图并计算相位裕度：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)%对于a,b有四个约束条件，做出a,b取值区间figure;syms a bkp=33;