基于 Matlab 的离散控制系统仿真

第一篇:基于 Matlab 的离散控制系统仿真

       2022 / 2022 学年第 1 学期

       计算机控制技术 实

       班 级 学 生 指 导 验 报 告

       学 号 1108030301 姓 名 蔡 梦 教 师 张 坤 鳌

       实验二 基于 Matlab 的离散控制系统仿真

       一、实验目的和要求:

       1、学习使用 Matlab 的命令对控制系统进行仿真的方法

       2、学习使用 Matlab 中的 Simulink 工具箱进行系统仿真的方法

       二、实验环境

       X86系列兼容型计算机,Matlab软件

       三、实验原理

       1、控制系统命令行仿真

       1)建立如图所示一阶系统控制模型并进行系统仿真:

       一阶系统闭环传递函数为G(S)=

       s1333s=s3,转换为离散系统脉冲传递函数并仿真。

       2)建立如图所示二阶系统控制模型并进行系统仿真:

       52s(s20.45)25251s(s20.45)=s220.45s52,二阶系统闭环传递函数为G(S)=转换为离散系统脉冲传递函数并仿真,改变参数,观察不同的系统的仿真结果。

       2、控制系统的 Simulink 仿真

       按图建立系统的 Simulink 模型,对不同的输入信号进行仿真,改变参数,观察不同的仿真结果。

       将上述系统离散化并仿真,观察仿真结果

       四、实验步骤

       1、根据实验原理对控制系统进行软件仿真

       2、观察记录输出的结果,与理论计算值相比较

       3、自行选择参数,练习仿真方法,观察不同的仿真结果

       5252s(s20.45)s(s20.45)525211s(s20.45)s(s20.45)进行软二阶系统闭环传递函数为G(S)=件仿真如下图:

       分别进行离散仿真:

       五、实验心得

       针对这次实验设计,我通过各种渠道,上课认真学习,请教老师、上网搜索,图书馆查阅,询问同学等学习到了很多知识,一步步了解最少拍控制系统设计,锻炼了自我学习能力。

       尽管学习上遇到了很多困难,结果也差强人意。但我们在不断处理困难的过程中磨练了处理事物的能力和耐心,也让同学间学会了互相学习,共享资源

第二篇:MATLAB与控制系统仿真实验报告

       《MATLAB与控制系统仿真》

       实验报告

       2022-2022学年 第 1 学期

       专业: 班级: 学号: 姓名:

       实验三 MATLAB图形系统一、实验目的:

       1.掌握绘制二维图形的常用函数。2.掌握绘制三维图形的常用函数。3.熟悉利用图形对象进行绘图操作的方法。4.掌握绘制图形的辅助操作。

       二、实验原理:

       1,二维数据曲线图

       (1)绘制单根二维曲线 plot(x,y);(2)绘制多根二维曲线 plot(x,y)当x是向量,y是有一维与x同维的矩阵时,则绘制多根不同颜色的曲线。当x,y是同维矩阵时,则以x,y对应列元素为横、纵坐标分别绘制曲线,曲线条数等于矩阵的列数。(3)含有多个输入参数的plot函数 plot(x1,y1,x2,y2,…,xn,yn)(4)具有两个纵坐标标度的图形 plotyy(x1,y1,x2,y2)2,图形标注与坐标控制 1)title(图形名称); 2)xlabel(x轴说明)3)ylabel(y轴说明)4)text(x,y图形说明)5)legend(图例1,图例2,…)

       6)axis([xmin xmax ymin ymax zmin zmax])3, 图形窗口的分割 subplot(m,n,p)4,三维曲线

       plot3(x1,y1,z1,选项1,x2,y2,选项2,…,xn,yn,zn,选项n)5,三维曲面

       mesh(x,y,z,c)与surf(x,y,z,c)。一般情况下,x,y,z是维数相同的矩阵。X,y是网格坐标矩阵,z是网格点上的高度矩阵,c用于指定在不同高度下的颜色范围。6,图像处理

       1)imread和imwrite函数 这两个函数分别用于将图象文件读入matlab工作空间,以及将图象数据和色图数据一起写入一定格式的图象文件。

       2)image和imagesc函数 这两个函数用于图象显示。为了保证图象的显示效果,一般还应使用colormap函数设置图象色图。

       三、实验仪器和设备:

       计算机一台(带有MATLAB6.5以上的软件环境)。

       四、预习要求:

       1.复习二维与三维图形的绘图函数。2.复习图形辅助操作。

       五、实验内容及步骤:

       1,设y[0.53sinx]cosx,在x=0~2π区间取101点,绘制函数曲线。21x

       2,已知y1=x2,y2=cos(2x),y3=y1*y2,完成下列操作:

       (1)在同一坐标系下用不同的颜色和线型绘制三条曲线;

       (2)分别用条形图、阶梯图、杆图和填充图绘制三条曲线。

       3,已知

       x,x02e y1In(x1x2),x02在-5<=x<=5区间绘制函数曲线。

       4,绘制函数的曲面图和等高线

       zcosxcosyex2y24

       其中x的21个值均匀分布在[-5,5]范围,y的31个值均匀分布在[0,10],要求使用subplot(2,1,1)和subplot(2,1,2)将产生的曲面图和登高图画在同一个窗口上。

       5.画出函数

       zx2y2sin(xy)的曲面及等高线图。

       x2y21绘制平面曲线,并分析参数a对其形状的影响。6.根据2a25a2

       四、心得体会:

       通过这次实验我能熟练掌握二维和三维图以及其他特殊图形的制作,弄清楚了基本的图形操作规则,大大加深了我对matlab的兴趣。

       实验二 MATLAB程序设计

       一、实验目的

       1.掌握利用if语句实现选择结构的方法。

       2.掌握利用switch语句实现多分支选择结构的方法。3.掌握利用for语句实现循环结构的方法。4.掌握利用while语句实现循环结构的方法。

       二、实验设备及条件

       计算机一台(带有MATLAB6.5以上的软件环境)。

       三、实验内容

       1.编写求解方程ax2bxc0的根的函数(这个方程不一定为一元二次方程,因a、b、c的不同取值而定),这里应根据a、b、c的不同取值分别处理,有输入参数提示,当a0,b0,c~0时应提示“为恒不等式!”。并输入几组典型值加以检验。

       clear,clc a=input('请输入一个数a=');b=input('请输入一个数b=');c=input('请输入一个数c=');m=b^2-4*a*c;if a==0

       if b==0

       '为恒不等式'

       end end

       m=b^2-4*a*c;if m>0

       x1=(-b sqrt(m))/(2*a)

       x2=(-b-sqrt(m))/(2*a)elseif m==0

       x=(-b)/(2*a)else

       '不存在正实根' end

       2.输入一个百分制成绩,要求输出成绩等级A 、A、B、C、D、E。其中100分为A ,90分~99分为A,80分~89分为B,70分~79分为C,60分~69分为D,60分以下为E。

       要求:(1)用switch语句实现。

       (2)输入百分制成绩后要判断该成绩的合理性,对不合理的成绩应输出出错信息。

       clear,clc for k=1:10

       a(k)={89 k};b(k)={79 k};

       c(k)={69 k};d(k)={59 k};end A=cell(3,6);A(1,:)={'a','b','c','d','e','f'};A(2,:)={85,76,95,100,40,65};for k=1:6

       switch A{2,k}

       case 100

       r='A ';

       case a

       r='A';

       case b

       r='B';

       case c

       r='C';

       case d

       r='D';

       otherwise

       r='E';

       end

       A(3,k)={r};end A A =

       'a'

       'b'

       'c'

       [85]

       [76]

       [95]

       'B'

       'C'

       'A'

       'd'

       'e'

       [100]

       [40]

       'A '

       'E'

       'f' [65] 'D'

       3.利用for循环语句编写计算n!的函数程序,取n分别为-89、0、3、5、10验证其正确性(输入n为负数时输出出错信息)。

       clear,clc n=input('请输入一个正数n=');if n<0

       '输入错误' elseif n==0

       'n!=0' elseif n==1

       'n!=1' else

       y=1;

       for i=1:1:n

       y=y*i;

       i=i 1;

       end

       y end 请输入一个正数n=-89

       ans =输入错误 请输入一个正数n=0

       ans =n!=0 请输入一个正数n=1

       ans =n!=1 请输入一个正数n=3

       y =6 请输入一个正数n=10

       y =3628800

       四、实验心得体会:

       通过本次实验课,我能熟练运用for循环语句,switch条件语句以及if条件语句的新用法,和在C中的区别。尽管如此,但是在实验中依然容易把for循环跟C语言中的for语句弄混,最后经过不懈努力下,终于弄明白了两者之间的差别,使我能更好的运用这些指令语句。

第三篇:控制系统的Matlab仿真与设计课后答案

       MATLAB课后习题答案 2.1 x=[15 22 33 94 85 77 60] x(6)x([1 3 5])x(4:end)x(find(x>70))2.3 A=zeros(2,5);

       A(:)=-4:5

       L=abs(A)>3 islogical(L)

       X=A(L)2.4 A=[4,15,-45,10,6;56,0,17,-45,0] find(A>=10&A<=20)2.5 p1=conv([1,0,2],conv([1,4],[1,1]));p2=[1 0 1 1];[q,r]=deconv(p1,p2);cq='商多项式为

       ';cr='余多项式为

       ';disp([cq,poly2str(q,'s')]),disp([cr,poly2str(r,'s')])2.6 A=[11 12 13;14 15 16;17 18 19];PA=poly(A)

       PPA=poly2str(PA,'s')3.1 n=(-10:10)';y=abs(n);plot(n,y,'r.','MarkerSize',20)axis equal grid on xlabel('n')3.2 x=0:pi/100:2*pi;y=2*exp(-0.5*x).*sin(2*pi*x);plot(x,y),grid on;3.3 t=0:pi/50:2*pi;x=8*cos(t);y=4*sqrt(2)*sin(t);z=-4*sqrt(2)*sin(t);plot3(x,y,z,'p');

       title('Line in 3-D Space');text(0,0,0,'origin');

       xlabel('X'),ylable('Y'),zlable('Z');grid;3.4

       theta=0:0.01:2*pi;

       rho=sin(2*theta).*cos(2*theta);polar(theta,rho,'k');3.5

       [x,y,z]=sphere(20);z1=z;

       z1(:,1:4)=NaN;c1=ones(size(z1));surf(3*x,3*y,3*z1,c1);hold on z2=z;

       c2=2*ones(size(z2));

       c2(:,1:4)=3*ones(size(c2(:,1:4)));surf(1.5*x,1.5*y,1.5*z2,c2);colormap([0,1,0;0.5,0,0;1,0,0]);grid on hold off 第四章

       function f=factor(n)if n<=1 f=1;else

       f=factor(n-1)*n;end

       function[s,p]=fcircle(r)s=pi*r*r;p=2*pi*r;

       function k=jcsum1(n)k=0;i=0;while i<=n k=k 2^i;i=i 1;end

       function k=jcsum(n)k=0;

       for i=0:n k=k 2^i;end

       4.1for m=100:999

       m1=fix(m/100);m2=rem(fix(m/10),10);m3=rem(m,10);

       if

       m==m1*m1*m1 m2*m2*m2 m3*m3*m3 disp(m)

       end end

       4.2[s,p]=fcircle(10)4.3y=0;n=100;for i=1:n

       y=y 1/i/i;end y

       4.4s=0;for i=1:5

       s=s factor(i);end s

       4.5sum=0;i=1;while sum<2000 sum=sum i;i=i 1;end;n=i-2

       4.6jcsum(63)jcsum1(63)4.1 for m=100:999

       m1=fix(m/100);

       m2=rem(fix(m/10),10);

       m3=rem(m,10);

       if m==m1*m1*m1 m2*m2*m2 m3*m3*m3

       disp(m)

       end end 4.3 y=0;n=100;

       for i=1:n

       y=y 1/i/i;end y 4.4 s=0;for i=1:5

       s=s factor(i);end s 4.5

       sum=0;i=1;

       while sum<2000

       sum=sum i;

       i=i 1;end;n=i-2 4.6

       i=0;k=0;while i<=63

       k=k 2^i;

       i=i 1;end k i

       i=0;k=0;for i=0:63

       k=k 2^i;end i k

       第五章

       function f=fxyz(u)x=u(1);y=u(2);z=u(3);f=x y.^2./x/4 z.^2./y 2./z;

       5.1A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];

       b=[13,-9,6,0]';x=Ab

       5.2[U,fmin]=fminsearch('fxyz',[0.5,0.5,0.5])

       5.3X=linspace(0,2*pi,50);Y=sin(X);

       P=polyfit(X,Y,3)AX=linspace(0,2*pi,50);Y=sin(X);Y1=polyval(P,X)

       plot(X,Y,':O',X,Y1,'-*')

       5.4x=0:2.5:10;h=[0:30:60]';

       T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];xi=[0:0.5:10];hi=[0:10:60]';

       temps=interp2(x,h,T,xi,hi,'cubic');

       mesh(xi,hi,temps);

       5.1 A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';x=Ab 5.3 X=linspace(0,2*pi,50);Y=sin(X);P=polyfit(X,Y,3)AX=linspace(0,2*pi,50);Y=sin(X);Y1=polyval(P,X)plot(X,Y,':O',X,Y1,'-*')6.1syms x

       y=finverse(1/tan(x))6.2syms x y

       f=1/(1 x^2);g=sin(y);fg=compose(f,g)6.3syms x

       g=(exp(x) x*sin(x))^(1/2);dg=diff(g)

       6.4F=int(int('x*exp(-x*y)','x'),'y')

       6.5syms x

       F=ztrans(x*exp(-x*10))6.6a=[0 1;-2-3];syms s

       inv(s*eye(2)-a);

       6.7 f=solve('a*x^2 b*x c')6.8 f=solve('x y z=1','x-y z=2','2*x-y-z=1')

       6.9y=dsolve('D2y 2*Dy 2*y=0','y(0)=1','Dy(0)=0')ezplot(y),grid on

       6.10a=maple('simplify(sin(x)^2 cos(x)^2);')

       6.11f=maple('laplace(exp(-3*t)*sin(t),t,s);')6.12 syms t x

       F=sin(x*t 2*t);L=laplace(F)第七章

       function

       [sys,x0,str,ts]=ww(t,x,u,flag)%¶¨ÒåÁ¬ÐøϵͳµÄSº¯Êý A=[0,1;-0.4,-0.2];B=[0;0.2];C=[1,0];D=0;

       switch flag, case 0,[sys,x0,str,ts]=mdlInitializeSizes(A,B,C,D);case 1,sys=mdlDerivatives(t,x,u,A,B,C,D);case 2,sys=mdlUpdate(t,x,u);case 3,sys=mdlOutputs(t,x,u,A,B,C,D);case 4,sys=mdlGetTimeOfNextVarHit(t,x,u);case 9,sys=mdlTerminate(t,x,u);otherwise

       error(['Unhandled flag = ',num2str(flag)]);end

       %=============================== function

       [sys,x0,str,ts]=mdlInitializeSizes(A,B,C,D)sizes = simsizes;sizes.NumContStates = 2;sizes.NumDiscStates = 0;sizes.NumOutputs = 1;sizes.NumInputs = 1;sizes.DirFeedthrough = 1;sizes.NumSampleTimes = 1;sys = simsizes(sizes);x0 = [0;0];str = [];ts = [0 0];

       %=============================== function

       sys=mdlDerivatives(t,x,u,A,B,C,D)

       sys = A*x B*u;

       %===============================

       function sys=mdlUpdate(t,x,u)sys = [];

       %=============================== function

       sys=mdlOutputs(t,x,u,A,B,C,D)sys = C*x D*u;

       %=============================== function

       sys=mdlGetTimeOfNextVarHit(t,x,u)

       sampleTime = 1;sys = t sampleTime;

       %===============================

       function sys=mdlTerminate(t,x,u)sys = [];

       7.1

       7.2

       7.3

       7.4

       7.5

       7.6

       7.7

       第八章

       8.1num=[5];den=[1,2,2];sys=tf(num,den)8.1.2s = tf('s');

       H = [5/(s^2 2*s 2)];H.inputdelay =2

       8.1.3h=tf([0.5,0],[1,-0.5,0.5],0.1)

       8.2num=2*[1,0.5];den=[1,0.2,1.01];

       sys=tf(num,den)

       [z,p,k]=tf2zp(num,den);zpk(z,p,k)

       [A,B,C,D]=tf2ss(num,den);ss(A,B,C,D)

       8.3 num=[1,5];den=[1,6,5,1];ts=0.1;

       sysc=tf(num,den);sysd=c2d(sysc,ts,'tustin')8.4.0

       8.4.1 %¶Ôϵͳ·½¿òͼÿ¸ö»·½Ú½øÐбàºÅ,ÓÐ8¸öͨµÀ,ÁÐдÿ¸öͨµÀ´«µÝº¯Êý r1=1;r2=2;c1=3;c2=4;G1=r1;G2=tf(1,[c1,0]);

       G3=1;%ÊÇ·ÖÀëµãºÍ»ãºÏµãµÄÁ¬Ïß,²»Äܺϲ¢,´«º¯Îª1 G4=-1;G5=1/r2;

       G6=tf(1,[c2,0]);G7=-1;G8=-1;%½¨Á¢ÎÞÁ¬½ÓµÄ״̬¿Õ¼äÄ£ÐÍ G=append(G1,G2,G3,G4,G5,G6,G7,G8)

       %д³öϵͳµÄÁ¬½Ó¾ØÕó

       Q=[1 4 0 %ͨµÀ1µÄÊäÈëÊÇͨµÀ4 2 1 7 %ͨµÀ2µÄÊäÈëÊÇͨµÀ1,7 3 2 0

       2 0 5 3 8 6 5 0 7 5 0

       6 0];%¸ººÅÔÚ´«º¯ÖÐÌåÏÖ %ÁгöϵͳµÄ×ܵÄÊäÈëºÍÊä³ö¶ËµÄ±àºÅ

       inputs=1;outputs=6;

       %Éú³É×éºÏºóϵͳµÄ״̬¿Õ¼äÄ£ÐÍ sys=connect(G,Q,inputs,outputs)

       8.4.2r1=1;r2=2;c1=3;c2=4;[A,B,C,D]=linmod('x84');[num,den]=ss2tf(A,B,C,D);sys=tf(num,den)

       8.5A=[1,1,0;0,1,0;0,0,2];B=[0,0;1,0;0,-2];n=size(A)

       Tc=ctrb(A,B);if n==rank(Tc)

       disp('ϵͳÍêÈ«ÄÜ¿Ø');else

       disp('ϵͳ²»ÍêÈ«ÄÜ¿Ø');end 第九章

       function [rtab,info]=routh(den)info=[];

       vec1=den(1:2:length(den));

       nrT=length(vec1);

       vec2=den(2:2:length(den)-1);rtab=[vec1;vec2,zeros(1,nrT-length(vec2))];for k=1:length(den)-2, alpha(k)=vec1(1)/vec2(1);

       for i=1:length(vec2),a3(i)=rtab(k,i 1)-alpha(k)*rtab(k 1,i 1);

       end

       if sum(abs(a3))==0 a3=polyder(vec2);

       info=[info,'All elements in row ',...int2str(k 2)' are zeros;'];

       elseif abs(a3(1))

       info=[info,'Replaced first element;'];

       end

       rtab=[rtab;a3, zeros(1,nrT-length(a3))];vec1=vec2;vec2=a3;end

       9.1num=[2,5,1];den=[1,2,3];bode(num,den);grid on;figure;

       nyquist(num,den);

       9.2num=5*[1,5,6];den=[1,6,10,8];step(num,den);grid on;figure;

       impulse(num,den);grid on;9.3kosi=0.7;wn=6;

       num=wn^2;den=[1,2*kosi*wn,wn^2];step(num,den);grid on;figure;

       impulse(num,den);grid on;9.4den=[1,2,8,12,20,16,16];[rtab,info]=routh(den)a=rtab(:,1)

       if all(a>0)

       disp('ϵͳÊÇÎȶ¨µÄ');

       else

       disp('ϵͳÊDz»Îȶ¨µÄ');

       end

       9.5num=7*[1,5];den=conv([1,0,0],conv([1,10],[1,1]));

       [gm,pm,wg,wc]=margin(num,den)

       9.1 >> sys=tf([2,5,1],[1,2,3])Transfer function: 2 s^2 5 s 1---------------s^2 2 s 3

       >> rlocus(sys)>> nyquist(sys)>> bode(sys)9.2

       >> G=tf(conv([5],[1,5,6]),[1,6,10,8]);>> step(G)>> impulse(G)

       sys=tf([5,25,30],[1,6,10,8]);>> step(sys)>> impulse(sys)9.4>>

       GH=tf(conv([7],[1,5]),conv([1,0,0],conv([1,10],[1,1])));

       >> [Gm,Pm,Wcg,Wcp]=margin(GH)Gm =

       0 Pm =

       -47.2870 Wcg =

       0 Wcp = 1.4354 >>

       GH=tf(conv([7],[1,5]),conv([1,0,0],conv([1,10],[1,1])));

       >> [Gm,Pm,Wcg,Wcp]=margin(GH)GH_close=feedback(GH,1)step(GH_close),grid on Gm =

       0 Pm =

       -47.2870 Wcg =

       0 Wcp =

       1.4354 第十章

       function s=bpts2s(bp,ts,delta)kosi=sqrt(1-1./(1 ((1./pi).*log(1./bp)).^2));

       wn=-log(delta.*sqrt(1-kosi.^2))/(kosi.*ts);

       s=-kosi.*wn j.*wn.*sqrt(1-kosi.^2);

       function

       [ngc,dgc]=fa_lead(ng0,dg0,Pm,wc,w)

       ngv=polyval(ng0,j*wc);dgv=polyval(dg0,j*wc);g=ngv/dgv;

       thetag=angle(g);mg=abs(g);thetar=Pm*pi/180;

       tz=(1 mg*cos(thetar-thetag))/(-wc*mg*sin(thetar-thetag));tp=(cos(thetar-thetag) mg)/(wc*sin(thetar-thetag));ngc=[tz,1];dgc=[tp,1];

       function

       [ngc,dgc]=fg_lag_pm(ng0,dg0,w,Pm)

       [mu,pu]=bode(ng0,dg0,w);wgc=spline(pu,w,Pm 5-180);%²åÖµÇóÈ¡Âú×ãÏà½ÇÔ£¶ÈµÄ½ÇƵÂÊ×÷ΪÆÚÍûµÄ¼ôÇÐƵÂÊ

       ngv=polyval(ng0,j*wgc);dgv=polyval(dg0,j*wgc);g=ngv/dgv;

       alph=abs(1/g);T=10/alph*wgc, ngc=[alph*T,1];dgc=[T,1];

       function

       [ngc,dgc]=fg_lag_wc(ng0,dg0,w,wc)

       ngv=polyval(ng0,j*wc);dgv=polyval(dg0,j*wc);g=ngv/dgv;

       alph=abs(1/g);T=10/(alph*wc);ngc=[alph*T,1];dgc=[T,1];

       function

       [ngc,dgc]=fg_lead_pd(ng0,dg0,wc)ngv=polyval(ng0,j*wc);dgv=polyval(dg0,j*wc);g=ngv/dgv;mg0=abs(g);

       t=sqrt(((1/mg0)^2-1)/(wc^2));%·ùÖµÏà¼ÓΪÁã

       ngc=[t,1];dgc=[1];%Gc(s)=1 Ts

       function

       [ngc,dgc]=fg_lead_pm(ng0,dg0,Pm,w)

       [mu,pu]=bode(ng0,dg0,w);%¼ÆËãÔ-ϵͳµÄ¶ÔÊýƵÂÊÏìÓ¦Êý¾Ý

       [gm,pm,wcg,wcp]=margin(mu,pu,w);%ÇóÈ¡Ô-ϵͳµÄÏà½ÇÔ£¶ÈºÍ¼ôÇÐƵÂÊ

       alf=ceil(Pm-pm 5);%¼ÆËã¿ØÖÆÆ÷ÌṩµÄ×î´ó³¬Ç°½Ç¶È£¬

       phi=(alf)*pi/180;%½«×î´ó³¬Ç°½Çת»»Îª»¡¶Èµ¥Î» a=(1 sin(phi))/(1-sin(phi));%¼ÆËãaÖµ

       dbmu=20*log10(mu);%ϵͳµÄ¶ÔÊý·ùÖµ

       mm=-10*log10(a);%wm´¦µÄ¿ØÖÆÆ÷¶ÔÊý·ùÖµ wgc=spline(dbmu,w,mm);%²îÖµÇóÈ¡wm£¬ÈÏΪwm£½wc T=1/(wgc*sqrt(a));%¼ÆËãT

       ngc=[a*T,1];dgc=[T,1];

       function

       [ngc,dgc]=fg_lead_pm_wc(ng0,dg0,Pm,wc,w)

       [mu,pu]=bode(ng0,dg0,w);ngv=polyval(ng0,j*wc);dgv=polyval(dg0,j*wc);

       g=ngv/dgv;%ÇóÔ-ϵͳÔÚÆÚÍûµÄ¼ôÇÐƵÂÊ´¦µÄƵÂÊÏìÓ¦Êý¾ÝG0(jwc)

       theta=180*angle(g)/pi;%Ô-ϵͳÔÚÆÚÍûµÄ¼ôÇÐƵÂÊ´¦µÄÏà½ÇÔ£¶È£¬µ¥Î»Îª¶È

       alf=ceil(Pm-(theta 180) 5);%

       ×î´ó³¬Ç°½Ç phi=(alf)*pi/180;

       a=(1 sin(phi))/(1-sin(phi));dbmu=20*log10(mu);mm=-10*log10(a);wgc=spline(dbmu,w,mm);T=1/(wgc*sqrt(a));

       KK=128;s1=-2 i*2*sqrt(3);a=2 ng0=[10];dg0=conv([1,0],conv([1,2],[1,8]));g0=tf(ng0,dg0);

       [ngc,dgc,k]=rg_lag(ng0,dg0,KK,s1,a);

       gc=tf(ngc,dgc)function s=kw2s(kosi,wn)s=-kosi.*wn j*wn.*sqrt(1-kosi.^2);

       10.1ng0=[1];dg0=10000*[1 0-1.1772];

       g0=tf(ng0,dg0);%Âú×㿪»·ÔöÒæµÄΪУÕýϵͳµÄ´«µÝº¯Êý

       s=kw2s(0.7,0.5)%ÆÚÍûµÄ±Õ»·Ö÷µ¼¼«µã

       ngc=rg_lead(ng0,dg0,s);gc=tf(ngc,1)g0c=tf(g0*gc);rlocus(g0,g0c);

       b1=feedback(g0,1);%δУÕýϵͳµÄ±Õ»·´«µÝº¯Êý

       b2=feedback(g0c,1);%УÕýºóϵͳµÄ±Õ»·´«µÝº¯Êý

       figure,step(b1,'r--',b2,'b');grid on %»æÖÆУÕýÇ°ºóϵͳµÄµ¥Î»½×Ô¾

       KK=20;s1=-2 i*sqrt(6);a=1 ng0=[10];dg0=conv([1,0],[1,4]);g0=tf(ng0,dg0);

       [ngc,dgc,k]=rg_lag(ng0,dg0,KK,s1,a);gc=tf(ngc,dgc)g0c=tf(KK*g0*gc);

       b1=feedback(k*g0,1);b2=feedback(g0c,1);step(b1,'r--',b2,'b');grid on

       g0c=tf(KK*g0*gc);rlocus(g0,g0c);

       b1=feedback(k*g0,1);

       b2=feedback(g0c,1);figure,step(b1,'r--',b2,'b');grid on

       ng0=[1];dg0=conv([1,0,0],[1,5]);g0=tf(ng0,dg0);w=logspace(-3,3);KK=1;Pm=50;

       [ngc,dgc]=lead4(ng0,dg0,KK,Pm,w);

       gc=tf(ngc,dgc);g0c=tf(KK*g0*gc);bode(KK*g0,w);hold

       on,bode(g0c,w);grid on,hold off [gm,pm,wcg,wcp]=margin(g0c)Kg=20*log10(gm)g1=feedback(g0c,1);bode(g1),grid on, [mag,phase,w]=bode(g1);a=find(mag<=0.707*mag(1));wb=w(a(1))max(mag)

       b=find(mag==max(mag))wr=w(b)

       KK=40;Pm=50;ng0= KK *[1];

       dg0=conv([1,0],conv([1,1],[1,4]));

       g0=tf(ng0,dg0);

       w=logspace(-2,4);

       [ngc,dgc]=fg_lead_pm(ng0,dg0,Pm, w)gc=tf(ngc,dgc),g0c=tf(g0*gc);

       b1=feedback(g0,1);b2=feedback(g0c,1);

       step(b1,'r--', b2,'b');grid on figure, bode(g0,'r--',g0c,'b',w), grid on,[gm,pm,wcg,wcp]=margin(g0c), Km=20*log10(gm)

       KK=200;bp=0.3;ts=0.7;delta=0.05;

       ng0=[1];dg0=conv([1,0],conv([0.1,1],conv([0.02 1],conv([0.01,1],[0.005 1]))));g0=tf(ng0,dg0);

       w=logspace(-4,3);t=[0:0.1:3];[mag,phase]=bode(KK*g0,w);[gm0,pm0,wg0,wc0]=margin(mag,phase,w),gm0=20*log10(gm0)%gm0 =-15.6769

       %2¡¢È·¶¨ÆÚÍûµÄ¿ª»·´«µÝº¯Êý mr=0.6 2.5*bp;

       wc=ceil((2 1.5*(mr-1) 2.5*(mr-1)^2)*pi/ts), h=(mr 1)/(mr-1)w1=2*wc/(h 1), w2=h*w1 w1=wc/10;w2=25;ng1=[1/w1,1];dg1=conv([1/w2,1],conv([1,0],[1,0]));

       g1=tf(ng1,dg1);

       g=polyval(ng1,j*wc)/polyval(dg1,j*wc);K=abs(1/g);%¼ôÇÐƵÂÊ´¦·ùֵΪ1£¬ÇóKÖµ g1=tf(K*g1)

       %3¡¢È·¶¨·´À¡»·½Ú´«µÝº¯Êý h=tf(dg1,ng1);Kh=1/K;h=tf(Kh*h)%ÆÚÍûƵÂÊÌØÐԵĵ¹ÌØÐÔ

       %4¡¢ÑéËãÐÔÄÜÖ¸±ê

       g2=feedback(KK*g0,h);%УÕýºó£¬ÏµÍ³µÄ¿ª»·´«µÝº¯Êý b1=feedback(KK*g0,1);b2=feedback(g2,1);

       bode(KK*g0,'r--',g2,'b',h,'g',w);grid on

       figure,step(b1, 'r--',b2, 'b',t);grid on,[pos,tr,ts,tp]=stepchar(b2,delta)

       function

       [ngc,dgc]=lag2(ng0,dg0,w,KK,Pm)[mu,pu]=bode(KK*ng0,dg0,w);wgc=spline(pu,w,Pm 5-180), ngv=polyval(KK*ng0,j*wgc);dgv=polyval(dg0,j*wgc);g=ngv/dgv;

       alph=abs(1/g), T=10/alph*wgc, ngc=[alph*T,1];dgc=[T,1];

       function

       [ngc,dgc]=lead4(ng0,dg0,KK,Pm,w)[mu,pu]=bode(KK*ng0,dg0,w);[gm,pm,wcg,wcp]=margin(mu,pu,w);alf=ceil(Pm-pm 5);phi=(alf)*pi/180;

       a=(1 sin(phi))/(1-sin(phi)), dbmu=20*log10(mu);mm=-10*log10(a);

       wgc=spline(dbmu,w,mm), T=1/(wgc*sqrt(a)),ngc=[a*T,1];dgc=[T,1];

       function

       [ngc,dgc]=ra_lead(ng0,dg0,s1)ngv=polyval(ng0,s1);dgv=polyval(dg0,s1);g=ngv/dgv;thetag=angle(g);mg=abs(g);thetas=angle(s1);ms=abs(s1);

       tz=(sin(thetas)-mg*sin(thetag-thetas))/(mg*ms*sin(thetag));tp=-(mg*sin(thetas) sin(thetag thetas))/(ms*sin(thetag));ngc=[tz,1];dgc=[tp,1];

       function

       [ngc,dgc,k]=rg_lag(ng0,dg0,KK,s1,a)

       ngv=polyval(ng0,s1);dgv=polyval(dg0,s1);g=dgv/ngv;

       k=abs(g);%ÆÚÍûÖ÷µ¼¼«µã´¦µÄ¸ù¹ì¼£ÔöÒæ beta=k/KK;

       [kosi1,wn1]=s2kw(s1);

       zc=-wn1*sin(a*pi/180)/sin(pi-atan(sqrt(1-kosi1^2)/kosi1)-(a*pi/180));%ÀûÓÃÕýÏÒ¶¨Àí pc=beta*zc;

       ngc=beta*[1,-zc];dgc=[1,-pc];

       function

       varargout=rg_lead(ng0,dg0,s1)if nargout==1

       ngv=polyval(ng0,s1);dgv=polyval(dg0,s1);g=ngv/dgv;

       thetal=pi-angle(g);

       zc=real(s1)-imag(s1)/tan(thetal);

       t=-1/zc;

       varargout{1}=[t,1];elseif nargout==2 ngv=polyval(ng0,s1);dgv=polyval(dg0,s1);

       g=ngv/dgv;theta=angle(g);phi=angle(s1);

       if theta>0 phi_c=pi-theta;

       end

       if theta<0;phi_c=-theta

       end

       theta_z=(phi phi_c)/2;theta_p=(phi-phi_c)/2;

       z_c=real(s1)-imag(s1)/tan(theta_z);

       p_c=real(s1)-imag(s1)/tan(theta_p);

       nk=[1-z_c];varargout{2}=[1-p_c];kc=abs(p_c/z_c);

       if theta<0 kc=-kc

       end

       varargout{1}=kc*nk;else

       error('Êä³ö±äÁ¿ÊýÄ¿²»ÕýÈ·£¡');end

       function [bp,ts]=s2bpts(s,delta)[kosi,wn]=s2kw(s);

       bp=exp(-kosi.*pi./sqrt(1-kosi.^2));

       ts=-1./(kosi.*wn)*log(delta.*sqrt(1-kosi.^2));

       function [kosi,wn]=s2kw(s)kosi=1./sqrt(1 (imag(s)/real(s)).^2);

       wn=-real(s)./kosi;

       %Èç¹ûwnΪ¸ºÖµ£¬ÔòwnÈ¡Õý£¬²¢ÇÒkosiÈ¡·´ iwn=(wn<0);wn(iwn)=-wn(iwn);kosi(iwn)=-kosi(iwn);

       function

       [pos,tr,ts,tp]=stepchar(g0,delta)

       [y,t]=step(g0);[mp,ind]=max(y);dimt=length(t);yss=y(dimt);

       pos=100*(mp-yss)/yss;tp=t(ind);for i=1:dimt

       if y(i)>=1 tr=t(i);

       break;

       end end;

       for i=1:length(y)

       if

       y(i)<=(1-delta)*yss|y(i)>=(1 delta)*yss ts=t(i);

       end end

       第十一章

       11.1a=[0 1 0;0 0 1;-1-5-6];b=[0 0 1]';

       p=[-2 4j;-2-4j;-10];K=acker(a,b,p)eig(a-b*K)

       11.2a=[0 1 0;0 0 1;-6-11-6];b=[1,0,0]';

       p=[-2 2*sqrt(3)*j;-2-2*sqrt(3)*j;-10];

       K=acker(a,b,p)eig(a-b*K)

       11.6A=[-1 0 0;0-2-3;0 0-3];B=[1 0;2 3;-3-3];C=[1 0 0;1 1 1 ];

       [G,K,L]=decoupling(A,B,C)

       11.8A=[0 20.6;1 0];b=[0 1]';c=[0 1];d=0;

       G=ss(A,b,c,d);

       Q=diag([1,0,0,0,0]);R=1;

       p=[-1.8 2.4j;-1.8-2.4j];[k,P]=lqr(A,b,Q,R);l=(acker(A',c',p))' Gc=-reg(G,k,l);zpk(Gc), eig(Gc.a), t=0:0.05:2;

       G_1=feedback(G*Gc,1);a1=eig(G_1.a), y_1=step(G_1,t);

       第十二章

       function

       [t,xx]=diffstate(G,H,x0,u0,N,T)xk=x0;u=u0;t=0 for k=1:N xk=G*xk H*u;x(:,k)=xk;t=[t,k*T];end;xx=[x0,x];

       12.1

       function sys=M601(t,x)u=1;

       sys=[x(2);x(3);-800*x(1)-80*x(2)-24*x(3) u];

       function

       [t,y]=ode4(A,B,C,D,x0,h,r,v,t0,tf)

       Ab=A-B*v*C;B=B;C=C;x=x0';y=0;t=t0;

       N=round((tf-t0)/h);for i=1:N k1=Ab*x B*r;

       k2=Ab*(x h*k1/2) B*r;k3=Ab*(x h*k2/2) B*r;k4=Ab*(x h*k3) B*r;x=x h*(k1 2*k2 2*k3 k4)/6;y=[y,C*x];t=[t,t(i) h];end

       12.1

       tspan=[0,10];x0=[0,0,0]';

       [t,y]=ode45('M601',tspan,x0);y1=800*y(:,1);plot(t,y1);

       12.2 num=10;den=conv([1,0],conv([1,2],[1,3]));

       [A,B,C,D]=tf2ss(num,den);x0=[0,0,0];v=1;t0=0;tf=10;h=0.01;r=1;

       [t,y]=ode4(A,B,C,D,x0,h,r,v,t0,tf);

       plot(t,y),grid

       12.3 12.4 g=[-2.8-1.4 0 0;1.4 0 0 0;-1.8-0.3-1.4-0.6;0 0 0.6 0];h=[1 0 1 0]';c=[0 0 0 1];d=0;

       x0=[0 0 0 0]';u=1;N=30;T=0.1;

       [t,xx]=diffstate(g,h,x0,u,N,T);plot(t,xx);y=c*xx;figure stairs(t,y)grid on

       12.6 第十四章

       14.1

       clear all;load optcar.mat;

       t=signals(1,:);p=signals(2,:);v=signals(3,:);a=signals(4,:);theta=signals(5,:);

       subplot(4,1,1);plot(t,p);grid on;ylabel('λÖÃ(m)');subplot(4,1,2);plot(t,v);grid on;ylabel('ËÙ¶È(m/s)');subplot(4,1,3);plot(t,a);grid on;ylabel('¼ÓËÙ¶È(m/s2)');subplot(4,1,4);plot(t,theta);grid on;ylabel('½Ç¶È(¶È)');

       14.1

       clear all load

       car.mat %½«µ¼Èëµ½car.matÖеķÂÕæʵÑéÊý¾Ý¶Á³ö t=signals(1,:);x=signals(2,:);theta=signals(3,:);x1=signals(4,:);theta1=signals(5,:);

       plot(t,x,t,x1);ylabel('С³µÎ»ÖÃ(m)'),grid on;%

       »æÖÆ¿ØÖÆÁ¦×÷ÓÃϽüËÆÄ£Ðͺ;«È·Ä£ÐÍxµÄµ¥Î»½×Ô¾ÏìÓ¦ÇúÏß figure % »æÖÆ¿ØÖÆÁ¦×÷ÓÃϽüËÆÄ£Ðͺ;«È·Ä£

       ÐÍthetaµÄµ¥Î»½×Ô¾ÏìÓ¦ÇúÏß plot(t,theta,t,theta1);ylabel('°Ú½ÇÖµ(rad)'),grid on;

第四篇:控制系统的MATLAB仿真与设计课后答案

       第二章

       1>>x=[15 22 33 94 85 77 60] >>x(6)>>x([1 3 5])>>x(4:end)>>x(find(x>70))2>>T=[1-2 3-4 2-3];>>n=length(T);>>TT=T';>>for k=n-1:-1:0 >>B(:,n-k)=TT.^k;>>end >>B >>test=vander(T)3>>A=zeros(2,5);>>A(:)=-4:5 >>L=abs(A)>3 >>islogical(L)>>X=A(L)4>>A=[4,15,-45,10,6;56,0,17,-45,0] >>find(A>=10&A<=20)5>>p1=conv([1,0,2],conv([1,4],[1,1]));>>p2=[1 0 1 1];>>[q,r]=deconv(p1,p2);>>cq='商多项式为 ';cr='余多项式为 ';>>disp([cq,poly2str(q,'s')]),disp([cr,poly2str(r,'s')])6>>A=[11 12 13;14 15 16;17 18 19];>>PA=poly(A)>>PPA=poly2str(PA,'s')第三章

       1>>n=(-10:10)';>>y=abs(n);>>plot(n,y,'r.','MarkerSize',20)>>axis equal >>grid on

       >>xlabel('n')2>>x=0:pi/100:2*pi;>>y=2*exp(-0.5*x).*sin(2*pi*x);>>plot(x,y),grid on;3>>t=0:pi/50:2*pi;>>x=8*cos(t);>>y=4*sqrt(2)*sin(t);

       >>z=-4*sqrt(2)*sin(t);>>plot3(x,y,z,'p');>>title('Line in 3-D Space');>>text(0,0,0,'origin');>>xlabel('X'),ylable('Y'),zlable('Z');grid;4>>theta=0:0.01:2*pi;>>rho=sin(2*theta).*cos(2*theta);>>polar(theta,rho,'k');5>>[x,y,z]=sphere(20);>>z1=z;>>z1(:,1:4)=NaN;>>c1=ones(size(z1));>>surf(3*x,3*y,3*z1,c1);>>hold on >>z2=z;>>c2=2*ones(size(z2));>>c2(:,1:4)=3*ones(size(c2(:,1:4)));>>surf(1.5*x,1.5*y,1.5*z2,c2);>>colormap([0,1,0;0.5,0,0;1,0,0]);>>grid on >>hold off

       第四章

       1>>for m=100:999 m1=fix(m/100);m2=rem(fix(m/10),10);m3=rem(m,10);if m==m1*m1*m1 m2*m2*m2 m3*m3*m3 disp(m)end end M文件:

       function[s,p]=fcircle(r)s=pi*r*r;p=2*pi*r;主程序:

       [s,p]=fcircle(10)3>>y=0;n=100;for i=1:n y=y 1/i/i;end >>y M文件:

       function f=factor(n)if n<=1 f=1;else

       f=factor(n-1)*n;end

       主程序: >>s=0;for i=1:5 s=s factor(i);end >>s 5>>sum=0;i=1;while sum<2000 sum=sum i;i=i 1;end;>>n=i-2 6 for循环M文件: function k=jcsum(n)k=0;for i=0:n k=k 2^i;end

       主程序: >>jcsum(63)

       While循环M文件: function k=jcsum1(n)k=0;i=0;while i<=n k=k 2^i;i=i 1;end

       主程序:

       >>jcsum1(63)第五章

       1>>A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];>>b=[13,-9,6,0]';>>x=Ab M文件:

       function f=fxyz(u)x=u(1);y=u(2);z=u(3);f=x y.^2./x/4 z.^2./y 2./z;主程序:

       [U,fmin]=fminsearch('fxyz',[0.5,0.5,0.5])3>>X=linspace(0,2*pi,50);>>Y=sin(X);>>P=polyfit(X,Y,3)>>AX=linspace(0,2*pi,50);>>Y=sin(X);>>Y1=polyval(P,X)>>plot(X,Y,':O',X,Y1,'-*')4>>x=0:2.5:10;>>h=[0:30:60]';>>T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];>>xi=[0:0.5:10];>>hi=[0:10:60]';>>temps=interp2(x,h,T,xi,hi,'cubic');>>mesh(xi,hi,temps);第六章

       1>>syms x

       >>y=finverse(1/tan(x))2>>syms x y

       >>f=1/(1 x^2);g=sin(y);>>fg=compose(f,g)3>>syms x

       >>g=(exp(x) x*sin(x))^(1/2);>>dg=diff(g)4>>F=int(int('x*exp(-x*y)','x'),'y')5>>syms x

       >>F=ztrans(x*exp(-x*10))6>>a=[0 1;-2-3];>>syms s

       >>inv(s*eye(2)-a);7>>f=solve('a*x^2 b*x c')8>>f=solve('x y z=1','x-y z=2','2*x-y-z=1')9>>y=dsolve('D2y 2*Dy 2*y=0','y(0)=1','Dy(0)=0')>>ezplot(y),grid on

       10>>a=maple('simplify(sin(x)^2 cos(x)^2);')11>>f=maple('laplace(exp(-3*t)*sin(t),t,s);')

       12>>syms t x

       >>F=sin(x*t 2*t);>>L=laplace(F)第七章

       第八章

       1-1>>h=tf([5,0],[1,2,2])1-2>>s = tf('s');>>H = [5/(s^2 2*s 2)];>>H.inputdelay =2 1-3>>h=tf([0.5,0],[1,-0.5,0.5],0.1)2>>num=2*[1,0.5];den=[1,0.2,1.01];>>sys=tf(num,den)>>[z,p,k]=tf2zp(num,den);>>zpk(z,p,k)>>[A,B,C,D]=tf2ss(num,den);>>ss(A,B,C,D)3 >>num=[1,5];den=[1,6,5,1];ts=0.1;>>sysc=tf(num,den);>>sysd=c2d(sysc,ts,'tustin')

       >>r1=1;r2=2;c1=3;c2=4;>>[A,B,C,D]=linmod('x84');>>[num,den]=ss2tf(A,B,C,D);>>sys=tf(num,den)5>>A=[1,1,0;0,1,0;0,0,2];B=[0,0;1,0;0,-2];>>n=size(A)>>Tc=ctrb(A,B);if n==rank(Tc)disp('系统完全能控');else

       disp('系统不完全能控');end

       第九章

       1>>num=[2,5,1];den=[1,2,3];>>bode(num,den);grid on;>>figure;>>nyquist(num,den);2>>num=5*[1,5,6];den=[1,6,10,8];>>step(num,den);grid on;>>figure;>>impulse(num,den);grid on;3>>kosi=0.7;wn=6;>>num=wn^2;den=[1,2*kosi*wn,wn^2];>>step(num,den);grid on;>>figure;>>impulse(num,den);grid on;4 M文件:

       function [rtab,info]=routh(den)info=[];vec1=den(1:2:length(den));nrT=length(vec1);vec2=den(2:2:length(den)-1);rtab=[vec1;vec2, zeros(1,nrT-length(vec2))];for k=1:length(den)-2, alpha(k)=vec1(1)/vec2(1);for i=1:length(vec2), a3(i)=rtab(k,i 1)-alpha(k)*rtab(k 1,i 1);

       end

       if sum(abs(a3))==0 a3=polyder(vec2);info=[info,'All elements in row ',...int2str(k 2)' are zeros;'];elseif abs(a3(1))

       rtab=[rtab;a3, zeros(1,nrT-length(a3))];vec1=vec2;vec2=a3;end

       主程序:

       >>den=[1,2,8,12,20,16,16];>>[rtab,info]=routh(den)>>a=rtab(:,1)if all(a>0)disp('系统是稳定的');else

       disp('系统是不稳定的');end

       5>>num=7*[1,5];den=conv([1,0,0],conv([1,10],[1,1]));>>[gm,pm,wg,wc]=margin(num,den)第十章 M文件:

       function varargout=rg_lead(ng0,dg0,s1)if nargout==1 ngv=polyval(ng0,s1);dgv=polyval(dg0,s1);g=ngv/dgv;thetal=pi-angle(g);zc=real(s1)-imag(s1)/tan(thetal);t=-1/zc;varargout{1}=[t,1];elseif nargout==2 ngv=polyval(ng0,s1);dgv=polyval(dg0,s1);g=ngv/dgv;theta=angle(g);phi=angle(s1);if theta>0 phi_c=pi-theta;end

       if theta<0;phi_c=-theta end

       theta_z=(phi phi_c)/2;theta_p=(phi-phi_c)/2;z_c=real(s1)-imag(s1)/tan(theta_z);p_c=real(s1)-imag(s1)/tan(theta_p);nk=[1-z_c];varargout{2}=[1-p_c];kc=abs(p_c/z_c);if theta<0 kc=-kc end

       varargout{1}=kc*nk;else

       error('输出变量数目不正确!');end

       主程序:

       >> ng0=[1];dg0=10000*[1 0-1.1772];>>g0=tf(ng0,dg0);%满足开环增益的为校正系统的传递函数 >>s=kw2s(0.7,0.5)%期望的闭环主导极点

       >>ngc=rg_lead(ng0,dg0,s);>>gc=tf(ngc,1)>>g0c=tf(g0*gc);>>rlocus(g0,g0c);>>b1=feedback(g0,1);%未校正系统的闭环传递函数 >>b2=feedback(g0c,1);%校正后系统的闭环传递函数 >>figure,step(b1,'r--',b2,'b');grid on %绘 2 M文件:

       function [ngc,dgc,k]=rg_lag(ng0,dg0,KK,s1,a)ngv=polyval(ng0,s1);dgv=polyval(dg0,s1);g=dgv/ngv;k=abs(g);%期望主导极点处的根轨迹增益 beta=k/KK;[kosi1,wn1]=s2kw(s1);zc=-wn1*sin(a*pi/180)/sin(pi-atan(sqrt(1-kosi1^2)/kosi1)-(a*pi/180));%利用正弦定理 pc=beta*zc;ngc=beta*[1,-zc];dgc=[1,-pc];主程序:

       >>KK=20;s1=-2 i*sqrt(6);a=1 >>ng0=[10];dg0=conv([1,0],[1,4]);>>g0=tf(ng0,dg0);>>[ngc,dgc,k]=rg_lag(ng0,dg0,KK,s1,a);>>gc=tf(ngc,dgc)>>g0c=tf(KK*g0*gc);>>b1=feedback(k*g0,1);>>b2=feedback(g0c,1);>>step(b1,'r--',b2,'b');grid on M文件:

       function [ngc,dgc,k]=rg_lag(ng0,dg0,KK,s1,a)ngv=polyval(ng0,s1);dgv=polyval(dg0,s1);g=dgv/ngv;k=abs(g);%期望主导极点处的根轨迹增益 beta=k/KK;[kosi1,wn1]=s2kw(s1);zc=-wn1*sin(a*pi/180)/sin(pi-atan(sqrt(1-kosi1^2)/kosi1)-(a*pi/180));%利用正弦定理 pc=beta*zc;ngc=beta*[1,-zc];dgc=[1,-pc];主程序:

       >>KK=128;s1=-2 i*2*sqrt(3);a=2 >>ng0=[10];dg0=conv([1,0],conv([1,2],[1,8]));>>g0=tf(ng0,dg0);>>[ngc,dgc,k]=rg_lag(ng0,dg0,KK,s1,a);>>gc=tf(ngc,dgc)>>g0c=tf(KK*g0*gc);>>rlocus(g0,g0c);>>b1=feedback(k*g0,1);>>b2=feedback(g0c,1);>>figure,step(b1,'r--',b2,'b');grid on 4 M文件:

       function [ngc,dgc]=lead4(ng0,dg0,KK,Pm,w)[mu,pu]=bode(KK*ng0,dg0,w);[gm,pm,wcg,wcp]=margin(mu,pu,w);alf=ceil(Pm-pm 5);phi=(alf)*pi/180;a=(1 sin(phi))/(1-sin(phi)), dbmu=20*log10(mu);mm=-10*log10(a);wgc=spline(dbmu,w,mm), T=1/(wgc*sqrt(a)), ngc=[a*T,1];dgc=[T,1];主程序:

       >>ng0=[1];dg0=conv([1,0,0],[1,5]);>>g0=tf(ng0,dg0);>>w=logspace(-3,3);>>KK=1;Pm=50;>>[ngc,dgc]=lead4(ng0,dg0,KK,Pm,w);>>gc=tf(ngc,dgc);g0c=tf(KK*g0*gc);>>bode(KK*g0,w);hold on,bode(g0c,w);grid on,hold off >>[gm,pm,wcg,wcp]=margin(g0c)>>Kg=20*log10(gm)>>g1=feedback(g0c,1);>>bode(g1),grid on, >>[mag,phase,w]=bode(g1);>>a=find(mag<=0.707*mag(1));>>wb=w(a(1))

       >>max(mag)>>b=find(mag==max(mag))>>wr=w(b)5 M文件:

       function [ngc,dgc]=fg_lead_pm(ng0,dg0,Pm,w)[mu,pu]=bode(ng0,dg0,w);%计算原系统的对数频率响应数据 [gm,pm,wcg,wcp]=margin(mu,pu,w);%求取原系统的相角裕度和剪切频率

       alf=ceil(Pm-pm 5);%计算控制器提供的最大超前角度,phi=(alf)*pi/180;%将最大超前角转换为弧度单位 a=(1 sin(phi))/(1-sin(phi));%计算a值

       dbmu=20*log10(mu);%系统的对数幅值

       mm=-10*log10(a);%wm处的控制器对数幅值

       wgc=spline(dbmu,w,mm);%差值求取wm,认为wm=wc T=1/(wgc*sqrt(a));%计算T ngc=[a*T,1];dgc=[T,1];主程序:

       >>KK=40;Pm=50;>>ng0= KK *[1];dg0=conv([1,0],conv([1,1],[1,4]));>>g0=tf(ng0,dg0);>>w=logspace(-2,4);>>[ngc,dgc]=fg_lead_pm(ng0,dg0,Pm,w)>>gc=tf(ngc,dgc),g0c=tf(g0*gc);>>b1=feedback(g0,1);b2=feedback(g0c,1);>>step(b1,'r--', b2,'b');grid on

       >>figure, bode(g0,'r--',g0c,'b',w), grid on, >>[gm,pm,wcg,wcp]=margin(g0c), Km=20*log10(gm)

第五篇:Matlab 控制系统 传递函数模型

       MATLAB及控制系统

       仿真实验

       班

       级:

       智能0702 姓

       名:

       刘保卫

       学

       号:

       06074053(18)实验四 控制系统数学模型转换及MATLAB实现

       一、实验目的

       熟悉MATLAB 的实验环境。

       掌握MATLAB 建立系统数学模型的方法。

       二、实验内容

       (注:实验报告只提交第2 题)

       1、复习并验证相关示例。(1)系统数学模型的建立 包括多项式模型(Transfer Function,TF),零极点增益模型(Zero-Pole,ZP),状态空间模型

       (State-space,SS);(2)模型间的相互转换 系统多项式模型到零极点模型(tf2zp),零极点增益模型到多项式模型(zp2tf),状态空间模

       型与多项式模型和零极点模型之间的转换(tf2ss,ss2tf,zp2ss…);(3)模型的连接

       模型串联(series),模型并联(parallel),反馈连接(feedback)

       2、用MATLAB 做如下练习。(1)用2 种方法建立系统 程序如下:

       %建立系统的多项式模型(传递函数)%方法一,直接写表达式 s=tf('s')Gs1=(s 2)/(s^2 5*s 10)%方法二,由分子分母构造 num=[1 2];den=[1 5 10];Gs2=tf(num,den)figure pzmap(Gs1)figure pzmap(Gs1)grid on

       运行结果:

       易知两种方法结果一样 的多项式模型。

       Transfer function: s Transfer function: s 2--------------s^2 5 s 10

       Transfer function: s 2--------------s^2 5 s 10(2)用2 种方法建立系统程序如下: %方法一 s=tf('s')Gs1=10*(s 1)/((s 1)*(s 5)*(s 10))% zpk模型 ZPK=zpk(Gs1)

       %方法二 % tf模型

       num=[10 10];den=conv([1 1],conv([1 5],[1 10]));Gs2=tf(num,den)% zpk模型 ZPK=zpk(Gs2)figure pzmap(Gs1)figure pzmap(Gs1)grid on

       运行结果:

       易知两种方法结果一样

       的零极点模型和多项式模型。

       Transfer function: s

       Transfer function: 10 s 10------------------------s^3 16 s^2 65 s 50 Zero/pole/gain: 10(s 1)------------------(s 10)(s 5)(s 1)

       Transfer function: 10 s 10------------------------s^3 16 s^2 65 s 50

       Zero/pole/gain: 10(s 1)------------------(s 10)(s 5)(s 1)

       (3)如图,已知G(s)和H(s)两方框对应的微分方程是:

       且初始条件为零。试求传递函数C(s)/R(s)及E(s)/R(s)。

       程序如下:

       %求微分方程的传递函数C(s)/R(s),E(s)/R(s)%求Gs=Cs/Rs n1=[20];d1=[6 10];Gs=tf(n1,d1)

       %求Hs=Bs/Cs n2=[10];d2=[20 5];Hs=tf(n2,d2)

       % C(s)/R(s)sys=feedback(Gs,Hs)

       % E(s)/R(s)=(Es/Cs)*(Cs/Rs)ER=sys/Gs

       运行结果:

       Transfer function:

       20--------

       % Gs=Cs/Rs 6 s 10

       Transfer function:

       10--------

       % Hs=Bs/Cs 20 s 5

       Transfer function:

       400 s 100---------------------

       s^2 230 s 250

       Transfer function: 2400 s^2 4600 s 1000------------------------

       2400 s^2 4600 s 5000

       % C(s)/R(s)% E(s)/R(s)=(Es/Cs)*(Cs/Rs)