计算机辅助设计实验报告-2

PAGEPAGE22控制系统计算机辅助设计综合实验指导实验名称：连续与离散系统校正实验，系统可控性与可观性实验，系统的simulink仿真实验陈茜编实验人：王奔学号：200830810215班级：08电气工程及其自动化(2)班信息工程系实验任务书1.有一个单位负反馈控制系统，如果控制对象的传递函数为设计要求：①相角裕度≥45°；②当系统的输入信号是单位斜坡信号时，稳态误差ess≤0.04。③要求绘制出校正后系统和未校正系统的Bode图及其闭环系统的单位阶跃响应曲线，并进行对比。2.有一个单位负反馈控制系统，如果控制对象的传递函数为：试设计一个串联滞后校正装置。设计要求：①相角裕度≥45°；②当系统的输入信号是单位斜坡信号时，稳态误差ess≤0.04。③要求绘制出校正后系统和未校正系统的Bode图及其闭环系统的单位阶跃响应曲线，并进行对比。3.有一个单位负反馈控制系统，如果控制对象的传递函数为试设计一个串联超前滞后校正装置，设计要求：①相角裕度≥45°；②当系统的输入信号是单位斜坡信号时，稳态误差ess≤0.04。③要求绘制出校正后系统和未校正系统的Bode图及其闭环系统的单位阶跃响应曲线，并进行对比。4.系统结构图如图所示，其中，采样周期Ts=0.01s，被控对象，为零阶保持器。用W变换法设计一超前校正装置D(z)，使系统相位裕度γ≥50°，校验设计后系统的性能指标。5.系统结构图如图所示，其中，采样周期Ts=0.01s，被控对象，为零阶保持器。用对数频率法设计D(z)，使系统开环增益k≥30(1/s)，截止频率ωc≥15(1/s),相位裕度γ≥50°1使

，求出未校正系统的开环系统的开环传递函数

,的传递函数模型参数。

2作未校正系统的Bode图，求出截止频率

和相位裕度

，并与系统性能指标要求比较

3用串联超前校正，选择Wc=20(1/s)，确定D（s）的参数。

4求出校正后系统的相位裕度γ，并检验是否满足要求，若不满足要求，则返回到3

5用双线性变换算法求出D（z）.

6检查系统的性能6.已知系统状态空间方程：，。对系统进行可控性|、可观性分析以及极点配置控制器的设计。实验要求：判别系统的可控性，求解系统的变换矩阵。如果系统完全可控，导出系统的第一可控标准型。判断系统的可观测性，求解系统的变换矩阵如果系统完全可观测，导出系统的第一可观测标准型。确定状态反馈矩阵K，使系统的闭环极点配置在位置上，并绘制出配置后系统单位阶跃时间响应曲线。7.建立如图所示双闭环计算机控制直流调速系统的simulink模型。仿真观测转速和电流跟踪特性等性能指标，调整数字PI调节器参数，以满足系统的性能指标要求。8.有限拍随动系统的设计：1）某一离散控制系统的被控队形的传递函数为：，保持器为零阶保持器，采样周期试设计单位速度输入时的有限拍控制器。2）某一离散控制系统的被控队形的传递函数为：，保持器为零阶保持器，采样周期试设计单位速度输入时的有限拍无纹波控制器。9．用函数调用的方法，求出，要求写出脚本式M文件和函数式M文件10.找出101——n之间的满足（1）它是完全平方数；（2）在满足完全平方数的前提下，有两位数字相同。统计满足条件的数的个数。11通过键盘输入一个整数，将该数变成按相反顺序的数，例如12345变成54321答案:1.主程序:源程序：num=8000;den=conv([1,0],conv([1,4],[1,80]));G=tf(num,den);%未校正系统的开环传递函数[Gm,Pm,Wcg,Wcp]=margin(G);%1)未校正系统的频域响应参数，计算需要相角欲度Pm(γ)w=0.1:0.1:10000;%确定频率的取值范围和频率采样的间隔值[mag,phase]=bode(G,w);magdb=20*log10(mag);%对数幅频、相频特性向量，计算需要magdb%(对数幅频响应值)phim1=45;deta=11;%2）(4.1.9)式设置系统设计参数phim1()；deta()phim=phim1-Pm+deta;%(4.1.9)式相位超前角φbita=(1-sin(phim*pi/180))/(1+sin(phim*pi/180));%求出β值n=find(magdb+10*log10(1/bita)<=0.0001);%找出满足magdb+10*log10(1/bita)<=0.0001%式的magdb向量所有的下标值wc=n(1);%通常magdb(1)+10*log10(1/bita)>0.0001;取的第1项为(比实际频率值大%10倍)，这是因为w=0.1:0.1:10000，而下标向量的值=1:1:100000w1=(wc/10)*sqrt(bita);w2=(wc/10)/sqrt(bita);%3)(4.1.11)和(4.1.12)式numc=[1/w1,1];denc=[1/w2,1];%(4.1.13)式，取K=1Gc=tf(numc,denc);%校正器的传递函数GmdB=20*log10(Gm);%4)下面程序校验系统校正后的系统指标GcG=Gc*G;[Gmc,Pmc,wcgc,wcpc]=margin(GcG);%GcG是校正后系统的开环传递函数GmcdB=20*log10(Gmc);disp('未校正系统的开环传递函数和频域响应参数：h，γ，wc')G,[GmdB,Pm,Wcp],disp('校正装置传递函数和校正后系统开环传递函数')Gc,GcG,disp('校正后系统的频域响应参数：h，γ，wc')[GmcdB,Pmc,wcpc],disp('校正装置的参数T和β值：T，β')T=1/w1;[T,bita],bode(G,GcG);figure(2);margin(GcG)[G,P,Wcgg,Wcpp]=margin(GcG);Psystemnew=feedback(GcG,1);figure(3);step(systemnew)运行结果：未校正系统的开环传递函数和频域响应参数：h，γ，wcTransferfunction:8000s^3+84s^2+320sans=10.526815.85789.5715校正装置传递函数和校正后系统开环传递函数Transferfunction:0.1458s+10.03602s+1Transferfunction:1166s+80000.03602s^4+4.025s^3+95.53s^2+320s校正后系统的频域响应参数：h，γ，wcans=16.046243.520813.7993校正装置的参数T和β值：T，βans=0.14580.2470截图：2.主程序:num=8000;den=conv([1,0],conv([1,4],[1,80]));G=tf(num,den);%未校正系统的传递函数gamma_cas=45;delta=10;%设计要求的相角裕度，(4.1.21)式中的设置为gamma_l=gamma_cas+delta;%根据(4.1.21)式=+△w=0.01:0.01:1000;%设置计算频率点[mag,phase]=bode(G,w);n=find(180+phase-(gamma_l)<=0.1);%求出对应的频率，wgamma_l=n(1)/100;%n(1)的值与校正后实际频率相差100倍[mag,phase]=bode(G,wgamma_l);%求出在频率点的相位值rr=-20*log10(mag);beta=10^(rr/20);%根据(4.1.22)式：w2=wgamma_l/10;w1=beta*w2;%根据第④步，求出和numc=[1/w2,1];denc=[1/w1,1];%根据(4.1.23)式，求出系统校正网络的传递函数Gc=tf(numc,denc)GcG=Gc*G%校正后系统的开环传递函数sysnew=feedback(GcG,1);figure(1);step(sysnew);figure(2);bode(G,GcG);figure(3);margin(GcG);beta[G,P,Wcgg,Wcpp]=margin(GcG);P运行结果：Transferfunction:3.344s+122.38s+1Transferfunction:2.676e004s+800022.38s^4+1881s^3+7246s^2+320sbeta=0.1494截图：3.主程序:num=1600;den=conv([1,0],([1,4]));G=tf(num,den);%未校正系统的传递函数[h,gamma,wg,wc]=margin(G);h=20*log10(h);%未校正系统的频域响应参数w=0.001:0.001:100;%设置计算频率点[mag,phase]=bode(G,w);%未校正系统的幅频与相频特性值向量disp('未校正系统的参数：h,wc,γ');[h,wc,gamma],gamma1=45;delta=6;%设计要求的相角裕度，(4.1.32)式中的设置为phim=gamma1-gamma+delta;%(4.1.32)式alpha=(1+sin(phim*pi/180))/(1-sin(phim*pi/180));%(4.1.33)式magdb=20*log10(mag);%用分贝值表示的未校正系统的幅值裕度n=find(magdb+10*log10(alpha)<=0.0001);%找出满足magdb+10*log10(alpha)<=0.0001式%的magdb向量所有的下标值wc=n(1);wcc=wc/1000;%的值与wc相差1000倍w3=wcc/sqrt(alpha);w4=sqrt(alpha)*wcc;%(4.1.34)式numc1=[1/w3,1];denc1=[1/w4,1];%(4.1.35)式Gc1=tf(numc1,denc1);%超前校正部分的传递函数w1=wcc/10;w2=w1/alpha;%取=；=numc2=[1/w1,1];denc2=[1/w2,1];%(4.1.36)式Gc2=tf(numc2,denc2);%滞后校正部分的传递函数Gc12=Gc1*Gc2;%串联超前-滞后校正网络的传递函数GcG=Gc12*G;%校正后系统的传递函数[Gmc,Pmc,wcgc,wcpc]=margin(GcG);GmcdB=20*log10(Gmc);disp('超前校正部分的传递函数'),Gc1,disp('滞后校正部分的传递函数'),Gc2,disp('串联超前-滞后校正网络的传递函数'),Gc12,disp('校正后系统的开环传递函数'),GcG,disp('校正后系统的性能参数：h,wc,γ及α值'),[GmcdB,wcpc,Pmc,alpha],bode(G,GcG),figure(2),margin(GcG),beta;GcGc=feedback(GcG,1);step(numc,denc);运行结果：未校正系统的参数：h,wc,γans=Inf39.89915.7249超前校正部分的传递函数Transferfunction:0.03902s+10.006604s+1滞后校正部分的传递函数Transferfunction:0.1605s+10.9484s+1串联超前-滞后校正网络的传递函数Transferfunction:0.006263s^2+0.1995s+10.006263s^2+0.955s+1校正后系统的开环传递函数Transferfunction:10.02s^2+319.3s+16000.006263s^4+0.98s^3+4.82s^2+4s校正后系统的性能参数：h,wc,γ及α值ans=Inf18.472025.67435.9083截图：4.主程序:margin(sysw);[Gm,Pm,Wcg,Wcp]=margin(sysw);w=0.1:0.1:10000;%确定频率的取值范围和频率采样的间隔值[mag,phase]=bode(sysw,w);magdb=20*log10(mag);%对数幅频、相频特性向量，计算需要magdb%(对数幅频响应值)phim1=45;deta=8;%2）(4.1.9)式设置系统设计参数phim1()；deta()phim=phim1-Pm+deta;%(4.1.9)式相位超前角φbita=(1-sin(phim*pi/180))/(1+sin(phim*pi/180));%求出β值n=find(magdb+10*log10(1/bita)<=0.0001);%找出满足magdb+10*log10(1/bita)<=0.0001%式的magdb向量所有的下标值wc=n(1);%通常magdb(1)+10*log10(1/bita)>0.0001;取的第1项为(比实际频率值大%10倍)，这是因为w=0.1:0.1:10000，而下标向量的值=1:1:100000w1=(wc/10)*sqrt(bita);w2=(wc/10)/sqrt(bita);%3)(4.1.11)和(4.1.12)式numc=[1/w1,1];denc=[1/w2,1];%(4.1.13)式，取K=1disp('控制器W传递函数');Gc=tf(numc,denc);%校正器的传递函数disp('校正后系统的开环W传递函数');GcG=Gc*sysw%校正后系统的开环传递函数margin(GcG)sysc=c2d(Gc,0.01,'tustin')syszk=sysc*dsysg;margin(syszk)运行结果：Transferfunction:4.983e-005z+4.967e-005z^2-1.99z+0.99Samplingtime:0.01控制器W传递函数校正后系统的开环W传递函数Transferfunction:-4.744e-008s^3-0.005683s^2+1.134s+11.084s^3+2.084s^2+s+1.232e-012Transferfunction:1.05z-1.041z-0.9908Samplingtime:0.01截图：5.主程序:s=0.01;[numh,denh]=pade(Ts/2,1);sysh=tf(numh,denh);numo=1;deno1=[0.210];sysgo1=tf(numo,deno1);sysgho=sysh*sysgo1;[kp,kv,ka,jienci]=wucha(sysgho);k=30/kv;sysgho=k*sysgho;margin(sysgho);[Gm,Pm,Wcg,Wcp]=margin(sysgho),h=20*log10(Gm);w=0.1:0.1:10000;%确定频率的取值范围和频率采样的间隔值[mag,phase]=bode(sysgho,w);magdb=20*log10(mag);%对数幅频、相频特性向量，计算需要magdb%(对数幅频响应值)phim1=50;deta=10;%2）(4.1.9)式设置系统设计参数phim1()；deta()phim=phim1-Pm+deta;%(4.1.9)式相位超前角φbita=(1-sin(phim*pi/180))/(1+sin(phim*pi/180));%求出β值n=find(magdb+10*log10(1/bita)<=0.0001);%找出满足magdb+10*log10(1/bita)<=0.0001%式的magdb向量所有的下标值wc=n(1);%通常magdb(1)+10*log10(1/bita)>0.0001;取的第1项为(比实际频率值大%10倍)，这是因为w=0.1:0.1:10000，而下标向量的值=1:1:100000w1=(wc/10)*sqrt(bita);w2=(wc/10)/sqrt(bita);%3)(4.1.11)和(4.1.12)式numc=[1/w1,1];denc=[1/w2,1];%(4.1.13)式，取K=1Gc=tf(numc,denc);disp('系统校正后的开环传递函数');GcG=Gc*sysghomargin(GcG);Gcd=c2d(Gc,0.01,'tustin')误差子程序:function[Kp,Kv,Ka,jienci]=wucha(sys)%sys：TF模型[num,den]=tfdata(sys,'v');%求出模型的分子和分母的多项表达式[Kp,jienci]=jixian(num,den,0);%Kp：位置误差系数；jienci：系统类型[Kv,jienci]=jixian(num,den,1);%Kv：速度误差系数[Ka,jienci]=jixian(num,den,2);%Ka：加速度误差系数disp('系统误差系数')Kp,Kv,Ka,disp('系统类型是'),jienci极限子程序:function[lim,jienci]=jixian(knum,kden,nn)%求Kp、Kv、Ka，n=0，单位位置输入，n=1：单位速度输入，nn=2：单位加速度输入，lim：Kp、Kv、Ka，jienn：系统类型jien=length(kden);jienn=0;%jien：系统分母多项式的长度fori=1:jienifkden(jien-i+1)<=0.00001%首先从多项式常数项开始判断，若为0，则系统类型值jienn=1+jienn;%jienn+1jienn；然后依s升幂次序判断各项的系数elsebreak%若系数不为0，则退出循环，得到系统类型值endendjiemn=length(knum);knumk=knum(jiemn);jienci=jienn;%knumk：分子的常数项ifnn==0&jienn==0%单位位置输入，0型系统lim=knumk/kden(jien);%kden(jien)：分母的常数项elseifnn==0lim='inf';%单位位置输入，1、2型系统endendifnn==1&jienn==0%单位速度输入，0型系统lim=0;elseifnn==1&jienn==1%单位速度输入，1型系统lim=knumk/kden(jien-1);elseifnn==1%单位速度输入，2型以上系统lim='inf';endendendifnn==2&jienn<2%单位加速度输入，0，1型系统lim=0;lim=0;elseifnn==2&jienn==2%单位加速度输入，2型系统lim=knumk/kden(jien-2);elseifnn==2%单位加速度输入，3型系统以上lim='inf';endendend运行结果：系统误差系数Kp=infKv=1Ka=0系统类型是jienci=1Gm=6.7079Pm=19.6906Wcg=31.5241Wcp=11.7478系统校正后的开环传递函数Transferfunction:-3.661s^2+1434s+120000.005232s^4+2.319s^3+91.46s^2+400sTransferfunction:4.076z-3.755z-0.6791截图：6.主程序:(1)A=[100;02-2;-1-10];B=[2;0;1];C=[120];%系统状态方程模型n=length(A);%求系统阶次nQ=ctrb(A,B);%求解系统可变换矩阵m=rank(Q);%求系统可控性矩阵的秩mifm==n%if-else-end程序判断系统是否完全可控，满足m=n系统状态完全可控Ac1=inv(Q)*A*Q;%Bc1=inv(Q)*B;Cc1=C*Q;%；disp('Systemiscontrollable.');disp('SystemFirstControllableCanonnicalFormis:');Ac1,Bc1,Cc1disp('TheTransformationMartrixis:');Qelse%m<n系统状态不完全可控disp('systemstateVariablecannotbetotallycontrolled');disp('TherankofSystemControllableMartixis:');m%可控的状态变量数end(2)A=[100;02-2;-1-10];B=[2;0;1];C=[120];%系统状态方程模型n=length(A);%求系统的秩nQ=obsv(A,C);%求解系统状态可观性矩阵,（6.2.3）式m=rank(Q);%求系统可观性矩阵的秩mifm==n%若m＝n则系统完全可观P=inv(Q);Ao1=inv(P)*A*P;Bo1=inv(P)*B;Co1=C*P;%（6.2.5）式disp('SystemisObservable.');%系统是可观测的disp('SystemFirstObservableCanonnicalFormis:');Ao1,Bo1,Co1disp('TheTransformationMartrixis:');Pelsedisp('systemstateVariablecannotbetotallyObserved');disp('TherankofSystemObservableMartixis:');mend(3)A=[100;02-2;-1-10];B=[2;0;1];C=[120];%系统状态方程模型D=[0];disp('原系统的极点为');p=eig(A)'P=[-2;-3.4+sqrt(-5.29);-3.4-sqrt(-5.29)];%P为所期望的极点K=place(A,B,P)%求状态反馈增益；或用acker(A,B,P)命令disp('配置后系统的极点为');p=eig(A-B*K)'disp('极点配置后的闭环系统为')sysnew=ss(A-B*K,B,C,D)%（6.3.1）式step(sysnew/dcgain(sysnew))%绘制时间响应曲线运行结果：(1)Systemiscontrollable.SystemFirstControllableCanonnicalFormis:Ac1=00-2100013Bc1=100Cc1=2-22TheTransformationMartrixis:Q=2220-201-20(2)SystemisObservable.SystemFirstObservableCanonnicalFormis:Ao1=0.00001.00000-0.0000-0.00001.0000-2.0000-0.00003.0000Bo1=2.0000-2.00002.0000Co1=1.000000TheTransformationMartrixis:P=-2.0000-2.00001.00001.50001.0000-0.50001.00000.2500-0.2500(3)原系统的极点为p=-0.73212.73211.0000K=114.6250299.1500-217.4500配置后系统的极点为p=-3.4000-2.3000i-3.4000+2.3000i-2.0000极点配置后的闭环系统为a=x1x2x3x1-228.3-598.3434.9x202-2x3-115.6-300.2217.5b=u1x12x20x31c=x1x2x3y1120d=u1y10Continuous-timemodel.截图：7.截图：8.主程序:（1）symssTKztkGs=(1/s)*K/s/(s+3);ft=ilaplace(Gs);ftt=subs(ft,t,k*T);GoGpZ=(1-z^-1)*ztrans(ftt);GoGpZ=simplify(GoGpZ);Ez=(1-z^-1)^2;Dcz=(1-Ez)/(Ez*GoGpZ);Dcz=simplify(Dcz);Dcz=subs(Dcz,T,0.5);Dcz=subs(Dcz,K,5);[dnum,dden]=numden(Dcz);dnum=sym2poly(dnum);dden=sym2poly(dden);disp('单位速度输入时的有限拍控制器Dc(z)=');[zpk]=tf2zp(dnum,dden);Dcz=zpk(z,p,k)disp('系统脉冲传递函数D(z)=C(z)/R(z)');Dz=1-Ez;