仿真技术实验程序及思考题解答完整版

实验一 连续系统的数字仿真

一、实验目的

1． 熟悉Matlab 中m 文件的编写；

2． 掌握龙格－库塔法的基本原理。

二、实验设备

计算机、MATLAB 软件

三、实验内容

假设单变量系统如图所示。试根据四阶龙格－库塔法，求系统输出y 的动态响应。

1．首先把原系统转化为状态空间表达式：?????=+=?CX

y bu AX X ，根据四阶龙格－库塔公式，

可得到： ??

???=++++=+++1143211)22(6k k k k CX y K K K K h X X （1） 其中： ?????????+++=+++=+++=+=)

()()2()2()2()2()(3423121h t bu hK X A K h t bu K h X A K h t bu K h X A K t bu AX K k k k k k k k k （2） 根据（1）、（2）式编写仿真程序。

2．在Simulink 环境下重新对上述系统进行仿真，并和1中结果进行比较。

四、实验结果及分析

要求给出系统输出响应曲线，并分析计算步长对龙格－库塔法的影响。

计算步长对龙格－库塔法的影响：单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加，不但引起计算量的增大，而且可能导致舍入误差严重积累，因此同积分的数值计算一样，微分方程的解法也有选择步长的问题。

源程序：

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;

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

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

K4=A*(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

plot(t,y)

Tf=5 h=0.02

五、思考题

1．试说明四阶龙格－库塔法与计算步长关系，它与欧拉法有何区别。

计算步长对龙格－库塔法的影响：单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步数就增加，不但引起计算量的增大，而且可能导致舍入误差严重积累，因此同积分的数值计算一样，微分方程的解法也有选择步长的问题。

区别：四阶龙格库塔法与欧拉法都是基于在初值附近展开成泰勒级数的原理，所不同的是取泰勒级数的项数，欧拉公式仅取到h项，四阶龙格库塔法取到h4项。

实验二面向结构图的仿真

一、实验目的

1．掌握连接矩阵及系统状态方程的确定方法；

2．掌握面向结构图的仿真方法。

二、实验设备

计算机、MATLAB软件

三、实验内容

仿真基本步骤：

1．给定输入信号，确定典型环节及环节参数；

2．确定连接矩阵；

3．输入仿真时间和计算步长；

4．求H,H-1和Q阵，确定A、B阵；

5．根据龙格－库塔法求状态方程的解；

6．根据1～5编写仿真程序。

四、实验结果及分析

源程序：

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+h*K1/2)+B*r;

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

K4=A*(x+h*K3)+B*r;

x=x+h*(K1+2*K2+2*K3+K4)/6;

y=[y;Wc*x];t=[t;t(i)+h];

end

plot(t,y) 仿真时间Tf=10 计算步长h=0.05

图一 仿真曲线

五、思考题

1． 典型环节的确定必须满足什么条件？

答： G(S)=bs

a ds ++=c U(S)X(S)式中u 为典型环节的输入，x 为典型环节的输出。 1.为了保证H 的逆1-H 存在，应严格按照i

b ≠0的原则确定每个典型环节，既避免以纯比例、

纯微分环节作为典型环节。

2.在输入向量不全为阶跃函数的情况下，只要在确定典型环节时，注意使含有微分项系数即

d≠0）的环节不直接与参考输入连接。

（

i

实验三连续系统的快速仿真

一、实验目的

1．熟悉增广矩阵的构建方法；

2．掌握连续系统的快速仿真基本原理。

二、实验设备

计算机、MATLAB软件

三、实验内容

假设某系统结构图如下，要求采用快速仿真方法求系统输出响应。Array

仿真基本步骤：

1．给定输入信号，确定典型环节及环节参数；

2．确定连接矩阵；

3．输入仿真时间和计算步长；

4．求H,H-1和Q阵，确定A、B阵；

5．构建增广矩阵；

6．采用增广矩阵法求齐次方程的解和系统输出响应。

根据1～6编写仿真程序。

四、实验结果和分析

源程序：

r=10;

[num,den]=series(1,[1,0],2,[1,2]);

[num,den]=series(num,den,10,[1,10]);

[num,den]=cloop(num,den);

[A1,B1,C1,D1]=tf2ss(num,den);

A=[A1,B1;0,0,0,0];

C=[C1,0];

Tf=10;

h=0.05;

k1=eye(size(A));

k2=A*h;

k3=k2*k2/2;

k4=k3*k2/3;

k=k1+k2+k3+k4;

x=[zeros(size(A1,1),1);10];

y=0;

t=0;

for i=1:Tf/h

x=k*x;

y=[y,C*x];

t=[t,i*h];

end

plot(t,y)

五、思考题

1．增广矩阵法和四阶龙格－库塔法有何区别？

答： ...t A !

31t A 21At I e 3322A ++++=！t 如果取t A e 的泰勒级数前五项，则增广矩阵的计算与四阶龙格库塔法相同。

四阶龙格库塔法具有很高的精度，但运行速度很慢，而增广矩阵法既满足了精度又满足了速度。

实验四 离散相似法仿真

一、实验目的

1． 掌握离散化的基本原理；

2． 掌握非线性系统的离散化仿真方法。

二、实验设备

计算机、MATLAB 软件

三、实验内容

已知非线性系统结构图如下，求系统输出响应。

1． 给定参考输入，输入仿真时间，采样周期（T=0.1s ）；

2． 将被控对象离散化；

3． 给定参数初始值；

4． 计算误差和非线性环节输出；

5．计算系统输出；

6．参数值更新；

7．若仿真时间到，则结束，否则转1。

编写仿真程序，求解系统输出响应。

四、实验结果及分析

源程序：

%backlash.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;

主程序：

r=10;

[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,'zoh');

x=zeros(size(G,1),1);

y=0;

t=0;

s=1;

e1=0;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];

u1=u;

e1=e;

end;

plot(t,y)

五、思考题

1．试比较零阶保持器与双线性变换法区别。

答：零阶保持器与双线性变换都能够实现将连续系统的离散化，只是所使用的方法不同，零阶保持器是在系统中加入了s e Ts --1实现离散化，而双线性变换是直接将1

12+-=z z h s 代入传递函数实现离散化。

实验五 采样控制系统的数字仿真

一、实验目的

3． 掌握采样控制系统基本原理；

4． 掌握常用的离散化方法；

5． 掌握采样控制系统的仿真方法。

二、实验设备

计算机、MATLAB 软件

三、实验内容

已知采样系统结构如下图所示，求系统的输出响应。

将系统中连续部分采用零阶保持器离散化，系统采用离散化方法进行仿真（同步采样），其具体步骤如下：

1、 给定参考输入，输入仿真时间，采样周期（T=0.1s ），PID 控制器参数；

2、将被控对象按零阶保持器方法离散化；

3、给定参数初始值；

4、计算误差和控制器输出；

5、计算系统输出；

6、参数值更新

7、若仿真时间到，则结束，否则转1。

编写仿真程序，并调整控制器参数，使得系统输出响应较为合理。

四、实验结果及分析

源程序：

r=1;

Tf=20;

Tm=1;

T=0.01;

kp=0.4;

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

五、思考题

2． 同步采样与异步采样在程序实现上有何区别？

答：同步采样：T=h ，异步采样T=N*h

同步采样时，在程序中只要一个for 循环即可实现，而异步采样需要两个for 循环实现。

实验六 串联超前校正

一、实验目的

1． 掌握超前校正装置的基本特性；

2． 熟悉串联超前校正方法。

二、实验设备

计算机、MATLAB 软件

三、实验内容 设有一单位负反馈系统，其开环传递函数为：)

2()(0+=s s k s G ，要求系统的稳态速度误差系数k v =20(1/s)，相位裕量r>50，幅值裕量k g >10db ，试确定串联校正装置。

设计基本步骤：

1． 根据性能指标确定开环增益k ；

2． 利用确定的开环增益k 画出未校正系统的Bode 图，并求其相位裕量r0和幅值裕量

k g ；

3． 确定所需增加的超前相位角ε+-=Φ0r r c ，ε一般取5～15；

4． 令超前校正装置最大超前角c m Φ=Φ，则c

c Φ+Φ-=sin 1sin 1α；

5． 将校正装置的最大超前相位角处的频率m w 作为校正后系统的剪切频率c w ，则：

α1

|)(|0=c jw G ；

6． 根据c m w w =，可求得：αc w T 1

=；

7． 画出校正后系统Bode 图，检验指标是否满足，如不满足，增大ε重新设计。

四、实验结果及分析

源程序见课本P220

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*lo g10(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('相位（度）');xlabel('频率（radc ）')

title(['校正后：幅值裕量=',num2str(20*log10(Gm)),'dB ，','相位裕量=',num2str(Pm),'度']);

disp(['校正前：幅值裕量=',num2str(20*log10(Gm1)),'dB,','相位裕量

=',num2str(Pm1),'度']);

disp(['校正后：幅值裕量=',num2str(20*log10(Gm)),'dB,','相位裕量

=',num2str(Pm),'度']);

五、思考题

1．简述串联超前校正的基本思想。

答：串联校正装置的主要作用是通过其相位超前效应来改变频率响应曲线的形状，产生足够大的相位超前脚，以补偿原来系统中元件造成的过大相位滞后。因此校正时应使校正装置的最大超前相位出现在校正后系统的开环剪切频率。

实验七 控制系统的计算机辅助分析

一、实验目的

1．掌握控制系统的稳定性分析；

2．掌握控制系统的时域分析；

3．学会利用MATLAB 软件求解控制系统性能指标。

二、实验设备

计算机、MATLAB 软件

三、实验内容

1．已知闭环系统的开环传递函数为

1

22532423)(2345234+++++++++=s s s s s s s s s s G 判断系统的稳定性，并给出不稳定极点。

num=[3 2 1 4 2];den=[3 5 1 2 2 1];

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

ii=find(real(p)>0);n1=length(ii);

if (n1>0)

disp('The Unstable Poles are:');

disp(p(ii));

else disp('System is Stable');

end

The Unstable Poles are:

0.4103 + 0.6801i

0.4103 - 0.6801i

2．对于典型二阶系统

2222)(n

n n s s s G ωξωω++= 试绘制出自然振荡频率ωn =6，阻尼比ξ＝0.2时系统的单位阶跃响应，并求出稳态误差、上升时间，最大超调量及调节时间。

wn=6;zeta=0.2;

num=wn^2;den=[1,2*zeta*wn,wn^2];

t=0:0.01:10;

[y,x,t]=step(num,den,t);

plot(t,y)

%求稳态误差

num0=wn^2;den0=[1,2*zeta*wn,0];

sys=tf(num0,den0);

sys1=1+sys;

sys2=tf(sys1.den,sys1.num);

num1=[1,0];den1=[1];sys3=tf(num1,den1);

num2=1;den2=[1 0];r=tf(num2,den2);

dcg=dcgain(sys2*sys3*r);

disp(['稳态误差 dcg=' num2str(dcg)]);

%求超调量

[M,lab]=max(y);

M=((M-1)/1)*100;

disp(['最大超调量 M=' num2str(M) '%'])

%求上升时间

t1=t(lab);

[Val,lab1]=min(abs(y(1:lab)-1));

Tr=t(lab1);

disp(['上升时间 Tr=' num2str(Tr)]);

%求调节时间

s=1;

L=length(t);

while s

lab2=find(abs(y(L:end)-1)>0.05);

if length(lab2)>0;

s=0;

else L=L-1;

end

end

L=L+1;

Ts=t(L);

disp(['调节时间 Ts=' num2str(Ts)]);

四、实验结果及分析

实验八积分分离PID控制算法及仿真

一、实验目的

1．掌握积分分离PID控制算法原理；

2．掌握积分分离PID控制器设计。

二、实验设备

计算机、MATLAB软件

三、实验内容

已知某被控对象传递函数为1

60)(80+=-s e s G s

，采样时间为20s ，延迟时间为4个采样时间，试设计积分分离PID 控制器，并求出系统输出响应。

积分分离算法：

∑=--+

+=k j d i p T k e k e k T j e k k e k k u 0))1()(()()()(β

式中，β为积分项开关系数，

???>≤=εεβ)(0)(1

k e k e

设计基本步骤：

1、给定参考输入，输入仿真时间，采样周期（T=0.1s ），积分分离PID 控制器参数；

2、被控对象离散化：y(k)=-den(2)y(k-1)+num(2)u(k-5)；

3、给定输入信号r(k)=40；

4、设定阈值ε>0，计算当前误差e(k)，若ε>)(k e 时，采用PD 控制，避免超调过大，同时使系统响应较快；当ε≤)(k e 时，采用PID 控制，保证精度。

5、计算系统响应y(k)。

四、实验结果及分析

%Integration Separation PID Controller

clear all ;

close all ;

ts=20;

%Delay plant

sys=tf([1],[60,1],'inputdelay',80);%构造exp(80)，inputdelay 不能少 dsys=c2d(sys,ts,'zoh');

[num,den]=tfdata(dsys,'v');

u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;

y_1=0;y_2=0;y_3=0;

error_1=0;

ei=0;

for k=1:1:200

time(k)=k*ts;

rin(k)=40;

kp=0.80;

ki=0.005;

kd=3.0;

%Delay plant

yout(k)=-den(2)*y_1+num(2)*u_5;

%I separation

error(k)=rin(k)-yout(k);

M=2;

if M==1 %Using integration separation

if abs(error(k))>=30&abs(error(k))<=40

beta=0.3;

elseif abs(error(k))>=20&abs(error(k))<=30

beta=0.6;

elseif abs(error(k))>=10&abs(error(k))<=20

beta=0.9;

else

beta=1.0;

end

elseif M==2

beta=1.0; %Not using integration separation

end

kp=0.80;

ki=0.005;

kd=3.0;

ei=ei+error(k)*ts; %积分作用

u(k)=kp*error(k)+kd*(error(k)-error_1)/ts+beta*ki*ei;

if u(k)>=110 % 限制PID输出值

u(k)=110;

end

if u(k)<=-110

u(k)=-110;

end

u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=u(k);

y_3=y_2;y_2=y_1;y_1=yout(k);

error_2=error_1;

error_1=error(k);

end

figure(1);

plot(time,rin,'b',time,yout,'r');

xlabel('time(s)');ylabel('rin,yout');

figure(2);

plot(time,u,'r');

xlabel('time(s)');ylabel('u');

五、思考题

1、试比较PID与积分分离PID控制算法。