理解数组运算与矩阵运算功能。
代码如下:
t=0:
0.1:
8;
y=1-2*exp
(1)-t.*sin(t);
figure('name','huizhishuxuemoxing');
plot(t,y);
title('y(t)=1-2e-tsin(t)');
xlabel('t');
ylabel('y');
运行结果:
4、理解函数文件与命令文件的区别,并自编编函数文件并调用。
自定义函数:
functions=sqrt(r)
s=pi*r^2+r^3
调用函数主代码:
r=2;
s=sqrt(r)
运行结果:
s=
1.4142
5、学会通过Help,熟悉MATLAB中为用户提供的功能各异的函数与命令。
实验二、数值积分算法练习与函数调用
一、实验目的
理解数值积分法,熟练掌握MATLAB的函数调用。
二、实验示例介绍:
1、用Euler法求初值问题的数值解:
设方程如下:
du/dt=u-2t/u(2-1)式
u(0)=1t=[0,1]
取步长h=0.1,名为FZSYZ1.M:
t0=0;
tf=1;
x0=1;
h1=0.1;
t=[t0:
h1:
tf];
n=length(t)
u=x0;
uu
(1)=u;
fori=2:
n
du=2*t(i-1)/u;
u=du*h1+u;
uu(i)=u;
end
figure
(1)
plot(t,uu)
运行结果:
2、在MATLAB中提供了现成的数值积分法的函数,如ode1,ode23,ode45求解微分方程组,下面介绍ode23它为二阶/三阶的RKF方法在MATLAB的ToolBox文件夹中MATLAB/funfun下的M文件,在此介绍其调用方法与应用例子如下:
[t,x]=ode23(‘系统函数名’,to,tf,x0,to1,trace)
其中,系统函数名为描述系统状态方程的M函数的名称,该函数名在调用时应该用引号括起来(文件字符串)
t0,tf为起始与终止时间
x0:
系统初始状态变量的值(列向量)
to1:
控制解的精度。
(缺省值为t01=10在ode23中)
trace:
输出形式控制变量,非零则程序运行的每步都显示出来。
Trace为缺省值——不显示中间结果。
t:
输出参数返回积分的时间离散值(列向量)
x:
输出参量,返回每个时间点值的解的列向量
注意:
系统函数的编写格式为固定的
Functionxdot=fun21(t,x)
Xdot=x-t^2
并以fun21.m存盘。
第二步:
编写如下程序并以fzsy22.m存盘。
t0=0;tf=3;t01=e-6;
x0=1;trace=1;
[t,x]=ode23(‘fun21’,t0,tf,x0,t01,trace)
Plot(t,x)
第三步:
在命令窗口下运行fzsy21即可求出x的解,并画出曲线。
3、实验具体内容、步骤、要求:
(1)运行交互式软件中函数调用,学习程序:
(2)试将(2-2)方程改为用Euler编程求解试比较用ode23求解结果。
先编写如下程序:
functiondx=fun22(t,x)
dx=x-t.^2;
代码如下:
t0=0;
tf=3;
options=odeset('RelTol',1e-6,'AbsTol',1e-6);
[T,X]=ode23('fun22',[t0,tf],1,options);
figure,plot(T,X);
title('ode23');
t0=0;
tf=3;
x0=1;
h1=0.1;
t=[t0:
h1:
tf];
n=length(t)
x=x0;
xx
(1)=x;
fori=2:
n
dx=x-t(i-1)^2;
x=dx*h1+x;
xx(i)=x;
end
figure
(2)
plot(t,xx)
运行结果:
(3)试将(2-1)方程改为用ode23算法调用函数求解,并试比较结果。
自定义函数fu233.m
functionudot=fun233(t,u)
udot=u-2*t/u;
end
代码如下:
t0=0;
tf=1;tol=1e-6;
u0=1;trace=1;
[t,u]=ode23('fun233',t0,tf,u0,tol,trace);
subplot(211)
title('ode23')
holdon
plot(t,u)
t0=0;
tf=1;
x0=1;
h1=0.1;
t=[t0:
h1:
tf];
n=length(t)
u=x0;
uu
(1)=u;
fori=2:
n
du=u-2*t(i-1)/u;
u=du*h1+u;
uu(i)=u;
end
subplot(212)
title('Euler')
holdon
plot(t,uu)
运行结果:
(4)利用ode23或ode45求解线性时不变系统微分方程
自定义函数fun34.m
functiondy=fun34(t,y)
A=[-0.5,1;-1,-0.5];
dy=A*y;
代码如下:
t0=0;tf=4;
y0=[01];
options=odeset('RelTol',1e-6,'AbsTol',[1e-61e-6]);
[T,Y1]=ode23('fun34',[t0,tf],y0,options);
figure('name','name');
plot(T,Y1(:
1),'g*',T,Y1(:
2),'m-');
title('ode23');
xlabel('t');
ylabel('y');
legend('y1','y2');
运行结果:
(5)求出G1(S)=2/(S^2+2S+1)与G2(S)=1/(2S^3+3S^2+1)的单位阶跃响应,并分别求出状态空间模型。
先作题一代码如下:
n1=[2];
d1=[121];
step(n1,d1);
[A1,B1,C1,D1]=tf2ss(n1,d1)
n2=1;d2=[2331];
figure
(2)
step(n2,d2);
[A2,B2,C2,D2]=tf2ss(n2,d2)
n=conv(200,[12]);
d=conv([11],[11042]);
[Z,P,K]=tf2zp(n,d)
figure(3)
step(n,d);
运行结果:
(6)设方程式为y=-40y,y(0)=2用各种数值积分方法与不同步长求方程式的数值解,并比较之。
代码如下:
figure('name','Eulerdy=-40y');
t0=0;
tf=1;
h=0.01;
t=t0:
h:
tf;
n=length(t);
y0=2;
y_Euler
(1)=y0;
fori=2:
n
Q(i-1)=-40*y_Euler(i-1)*h;
y_Euler(i)=y_Euler(i-1)+Q(i-1);
end
[T,Y]=ode23('fun36',[t0,tf],y0);
plot(T,Y,'m*');
holdon
plot(t,y_Euler);
title(['h=',num2str(h)]);
xlabel('t');
ylabel('y');
运行结果:
实验三、控制工具箱与SIMULINK软件应用
一、实验目的:
熟悉工具箱及其使用,进行系统仿真分析,通过仿真对系统进行校正校验。
二、实验预习要求:
必须先复习教材及上课介绍的有关控制工具箱命令与SIMULINK仿真工具的使用,并对实验题目作好准备。
三、学会调出、运行已由SIMULINK建立的仿真模型。
在打开MATLAB的窗口下,
(1)、输入SIMULINK或双击命令窗口的工具栏左起第二个按钮(NEWSIMULINKModel)将会打开,LibrarySIMULINK然后指向FILE菜单下拉菜单open调出fzsy31.mdl文件。
然后在fzy31.mdl文件的菜单观察并记录有关设置参数,然后指向Start下拉菜单单击一次则观察输出图形。
该仿真模型(SIMULINK系统结构图)如下:
(2)按照下图修改上述仿真模型,观察分析系统输出。
四、实验设计题目与要求
解:
对于前3问的实验程序代码为:
clc;clearall;
figure('name','根轨迹')
num1=[1];
den1=[100];
figure
(1),rlocus(num1,den1);
title('增加比例补偿器');
num2=[11];
den2=conv([100],[15]);
figure
(2),rlocus(num2,den2);
title('增加超前补偿器');
num3=[15];
den3=conv([100],[11]);
figure(3),rlocus(num3,den3);
title('增加滞后补偿器');
运行结果:
分析:
通过以上的图像可知,
(1)增加比例补偿器后,系统仍处在临界稳定状态,属于不稳定系统,增加Kc的作用仅是增加闭环纯虚数极点的大小,不能影响系统的稳定。
(2)增加超前补偿器后系统的根轨迹都在左半平面内,当Kc≠0时,系统处稳定,在一定的范围内增加Kc的作用是增强系统的稳定性。
(3)增加滞后补偿器后,系统的三条根轨迹有两条都在右半平面,系统处在不稳定,增加Kc就进一步增强系统的不稳定性。
对于第4问,仿真系统图如下:
模块I的仿真详细图
系统总图
分析,模块I、II、III调节器分别对应前三问中的三种情况,并且在每个模块中分别设置Kc=0.1、0.5、1,通过示波器观察Kc值不同时对同一个系统的影响。
仿真图如下:
由上图可知,随着Kc的增大系统的快速性逐渐的增大,但是系统的震荡性越大。
解:
实验才程序如下,图像如下:
clc;clearall;
z1=[];p1=[0-2];k1=4;
G=zpk(z1,p1,k1);
%PI补偿器K1(s)=(s+1)/s,K2(s)=(s+0.5)/s
K1=tf([11],[10]);
K2=tf([10.5],[10]);
A=K1*G;
sysA=feedback(A,-1)
sysB=feedback(G,K1)
C=K2*G;
sysC=feedback(G,-1)
subplot(231),step(sysA)%;title('系统A阶跃响应');
subplot(232),step(sysB);title('系统B阶跃响应');
subplot(233),step(sysC);title('系统C阶跃响应');
t=0:
0.001:
10;
u=t;%斜坡输入信号
subplot(234),lsim(sysA,u,t),title('系统A斜坡响应');
subplot(235),lsim(sysB,u,t);title('系统B斜坡响应');
subplot(236),lsim(sysC,u,t);title('系统C斜坡响应');
运行结果:
分析:
只有系统B是稳定的,只有针对系统B而言才有稳态误差;阶跃响应终值为0,稳态误差为1;斜坡响应终值约为1,稳态误差为无穷大。
取Kp=9,Ki=7,凑试Kd如下:
Kp=9;Ki=7;
Kd=[5:
13];
figure(3)
fori=1:
length(Kd)
G3=tf([Kd(i),Kp,Ki],[1,0]);
sysG3=G1*G3;
sysG33=feedback(sysG3,1);
step(sysG33);
axis([0,13,0.2,1.6])
holdon;
解:
实验程序如下:
取Kp=9,凑试Ki如下:
Kp=9;
Ki=[1:
2:
10];
figure
(2)
fori=1:
length(Ki)
G2=tf([Kp,Ki(i)],[1,0]);
sysG2=G2*G1;
sysG22=feedback(sysG2,1);
step(sysG22);
holdon;
首先凑试Kp方法如下:
num=[1];
den=conv([1,-1],[1,-1]);
G1=tf(num,den);
Kp=[1:
3:
15];
fori=1:
length(Kp)
G=Kp(i)*G1;
sysG=tf(G);
sysG11=feedback(sysG,1);
step(sysG11);
holdon
凑试Kp如上图凑试Ki如上图凑试Kd如上图
综合上面的PID参数分别凑试可得出当Kp=9,Ki=7,Kd=12~13之间,校正后的系统阶跃输出如下:
实验四、数字控制系统仿真与综合应用
解,将上面的代码输入后图像如下:
解:
(1)建立了仿真系统如下:
仿真系统的仿真如上
仿真示波器的输出如下:
(2)先设计连续系统PID调节器,然后将其离散化,建立simulink仿真系统,观察加入零阶保持器后系统的输出,根据指标再将微调一下离散PID参数,使其达到理想指标。
连续的PID凑试法程序如下:
Kd凑试如下:
Kp=0.4;Ki=0.001;
Kd=[0:
0.02:
0.1];
figure(3)
fori=1:
length(Kd)
G3=tf([Kd(i),Kp,Ki],[1,0]);
sysG3=G1*G3;
sysG33=feedback(sysG3,1);
step(sysG33);
holdon
axis([0,10,0.2,1.6])
end
Ki凑试如下:
p=0.4;
Ki=[0.001:
0.1:
1];
figure
(2)
fori=1:
length(Ki)
G2=tf([Kp,Ki(i)],[1,0]);
sysG2=G2*G1;
sysG22=feedback(sysG2,1);
step(sysG22);
holdon
axis([0,13,0,2])
end
Kp凑试如下:
clc;closeall;clearall;
num=[2000];
den=conv(conv([1,0],[1,10]),[1,20]);
G1=tf(num,den);
Kp=[0.1:
0.1:
1];
fori=1:
length(Kp)
G=Kp(i)*G1;
sysG=tf(G);
sysG11=feedback(sysG,1);
step(sysG11);
holdon
axis([0,13,0.2,1.6])
end
有Kp凑试图像可知,Kp从0.1~1逐渐变大过程中,系统反应速度逐渐增大,但考虑Kp越大导致的超调也会越大,所以取Kp=0.4,超调为20%;
由图像可知随着Ki的增大系统的超调增大速度增快,但系统震荡越大,取Ki=0.001,超调为小于20%;
由Kd凑试图像可知,当Kd=0.02时,系统效果较理想。
建立simulink仿真系统如下:
仿真系统总图
调节器subsystem详细图
综合simulink仿真中加入零阶保持器的仿真结果与该PID凑试m文件的仿真结果可知当存在D时系统经过长时间的有文波调节过程后文波越来越大变得不稳定,
所以只用PI调节,取Kp=0.3,Ki=0.001,Kd=0。
最终的系统仿真图如下:
解:
控制系统仿真图如下:
仿真输出图:
分析:
离散环节的采样周期T1根据以下因素确定:
(1)系统的频带宽度;
(2)实际采样开关硬件的性能;(3)实现数字控制器的计算程序的执行时间的长短。
对于种种原因采样周期T1比较大,但对系统中德连续部分按照采样周期选择仿真步长T2,将会出现较大的误差,因此T1>T2。