CAD仿真实验程序及波形图

自动化101班 周强 33号 实验一 程序： r=5;

numo=[1];deno=[1 4 8 5]; numh=[1];denh=1;

[num,den]=feedback(numo,deno,numh,denh); [A,b,C,d]=tf2ss(num,den);

Tf=input('仿真实验Tf= ');h=input('计算步长h= '); x=[zeros(length(A),1)];y=0;t=0; for i=1:Tf/h; K1=A*x+b*r;

K2=A*(x+K1*h/2)+b*r; K3=A*(x+K2*h/2)+b*r; K4=A*(x+K3*h)+b*r;

x=x+(K1+2*K2+2*K3+K4)*h/6; y=[y;C*x];t=[t;t(i)+h]; end

plot(t,y);

取仿真时间Tf= 5 计算步长h= 0.02 波形图：

0.90.80.70.60.50.40.30.20.1000.511.522.533.544.55系统方框图：

波形图：

实验二 程序： r=10;

P=[0 1 1 0;2 1 2 0;10 1 10 0]; W=[0 0 -1;1 0 0;0 1 0]; W0=[1;0;0];Wc=[0 0 1];

Tf=input('仿真实验 Tf = ');h=input('计算步长 h = ');

A1=diag(P(:,1)); B1=diag(P(:,2)); C1=diag(P(:,3)); D1=diag(P(:,4)); H=B1-D1*W;Q=C1*W-A1;

A=inv(H)*Q;B=inv(H)*C1*W0;

x=[zeros(length(A),1)];y=[zeros(length(Wc(:,1)),1)]; t=0;

for i=1:Tf/h K1=A*x+B*r;

K2=A*(x+K1*h/2)+B*r; K3=A*(x+K2*h/2)+B*r; K4=A*(x+K3*h)+B*r;

x=x+(K1+2*K2+2*K3+K4)*h/6; y=[y;Wc*x];t=[t;t(i)+h]; end

plot(t,y)

取仿真时间 Tf = 10 计算步长 h = 0.01 波形图：

121086420012345678910实验三 程序：

[num1,den1]=series(1,[1,0],2,[1,2]); [num2,den2]=series(num1,den1,10,[1,10]); [num,den]=cloop(num2,den2); [A,B,C,D]=tf2ss(num,den);

Tf=input('仿真实验 Tf = ');h=input('计算步长 h = ');

A1=[A,B;zeros(1,size(A,2)),zeros(1,size(B,2))]; C1=[C,0]; y=0;

x=[zeros(size(A,1),1);10]; t=0;

for i=1:tf/h

K1=eye(size(A1)); K2=A1*h; K3=K2*K2/2; K4=K3*K2/3; K=K1+K2+K3+K4; x=K*x;y=[y,C1*x]; t=[t,i*h]; end

plot(t,y)

取仿真时间 Tf = 10 计算步长 h = 0.01 波形图：

121086420012345678910实验四 子程序： ocklash.m

function[x,u1]=backlash(u1,u,x1,s) if(u>u1)

if((u-s)>=x1) x=u-s; else x=x1; end

else if(u

if((u+s)<=x1)x=u+s; else x=x1; end

else x=x1; end end

u1=u; end 主程序：

[num,den]=series([1,0.5],[1,0.1],[1],[1,0]); [num,den]=series(num,den,2,[1,2]); [num,den]=series(num,den,10,[1,10]); [A,B,C,D]=tf2ss(num,den); tf=10; T=0.025;

[G,H,C,D]=c2dm(A,B,C,D,T,'z0h'); x=zeros(size(G,1),1);

r=10;y=0;t=0;e1=0;s=1;u1=0; for i=1:tf/T e=r-y(end);

u=backlash(e1,e,u1,s); x=G*x+H*u; y=[y,C*x+D*u]; t=[t,i*T]; end

plot(t,y)

波形图：

181614121086420024681012实验五 程序： r=1; T=0.02; Tf=20; kp=1; ki=0; kd=0;

[A,B,C,D]=tf2ss(1,[1 1 0]);

[G,H,C,D]=c2dm(A,B,C,D,T,'zoh'); u1=0; e2=0; e1=0;

x=zeros(2,1); y=0; t=0;

for i=1:Tf/T e=r-y(end);

u2=u1+kp*(e-e1)+ki*e+kd*(e-2e1+e2); x=G*x+H*u2; y=[y,C*x+D*u2]; t=[t,i*T];

e2=e1;e1=e;u1=u2; end

plot(t,y); 波形图：

1.41.210.80.60.40.2002468101214161820

实验六 程序：

num0=40;den0=conv([1,0],[1,2]);

[Gm1,Pm1,Wcg1,Wcp1]=margin(num0,den0); r=50;r0=Pm1;

w=logspace(-1,3);

[mag1,phasel]=bode(num0,den0,w); for epsilon=5:15

phic=(r-r0+epsilon)*pi/180;

alpha=(1+sin(phic))/(1-sin(phic)); [i1,ii]=min(abs(mag1-1/sqrt(alpha))); wc=w(ii);

T=1/(wc*sqrt(alpha));

numc=[alpha*T,1];denc=[T,1];

[num,den]=series(num0,den0,numc,denc); [Gm,Pm,Wcg,Wcp]=margin(num,den); if(Pm>=r); break; end end

printsys(numc,denc) printsys(num,den)

[mag2,phase2]=bode(numc,denc,w); [mag,phase]=bode(num,den,w);

subplot(2,1,1);semilogx(w,20*log10(mag),w,20*log10(mag1),'--',w,20*log10(mag2),'-.'); grid;ylabel('幅值（dB）

');title('--Go,-.Gc,-GoGc'); subplot(2,1,2);

semilogx(w,phase,w,phasel,'--',w,phase2,'-.',w,(w-180-w),':');

grid;ylabel('相位（dB）');xlabel('频率（rad/sec）')

title(['校正后：幅值裕量

=',num2str(20*log10(Gm)),'dB， ','相位裕量

=',num2str(Pm),'度']); disp(['校正前：幅值裕量

=',num2str(20*log10(Gm1)),'dB,','相位裕量=',num2str(Pm1),'度']); disp(['校正后：幅值裕量

=',num2str(20*log10(Gm)),'dB,','相位裕量=',num2str(Pm),'度']); 得出结果： num/den =

0.22541 s + 1 -------------- 0.053537 s + 1

num/den =

9.0165 s + 40

------------------------------- 0.053537 s^3 + 1.1071 s^2 + 2 s

校正前：幅值裕量=InfdB,相位裕量=17.9642度 校正后：幅值裕量=InfdB,相位裕量=50.7196度 波形图：

--Go,-.Gc,-GoGc50幅值（dB）0-50-100-110100100101102103校正后：幅值裕量=InfdB，相位裕量=50.7196度相位（度）0-100-200-11010010频率（rad/sec）1102103