【自动控制原理仿真实验】 控制系统仿真实验（实验二）
控制系统仿真实验（实验二）
- 一、实验目的
- 二、实验简介
- 1.欧拉法
- 2.梯形法
- 3.四阶Runge-kutta法
- 4.离散相似法
- 保持器类型
- 三、实验过程
- 仿真时间及仿真步距的估计
- 1.整体离散法
- 2.分环节离散方法
- 3.数值积分法
- 四.仿真分析
- 相对作用强度
- 仿真步长的影响
一、实验目的
①构建与实验1相同的控制系统；应用欧拉法，梯形法，Runge-kutta法，分环节离散法等方法进行仿真。研究不同仿真方法的优缺点以及仿真步距对仿真精度的影响。
②（选做1）应用整体离散方法进行仿真，并加入对比研究。
③（选做2）对比讲过的仿真方法，自行设计一种仿真算法并和已有的仿真算法进行对比研究
二、实验简介
1.欧拉法
欧拉法原理:
逐次替代，最后求出所要求的解，并达到一定的精度。
将区间 [ a , b ] [a, b] [a,b] 分成 n n n 段, 那么方程在第 x i x_{i} xi 点有 y ′ ( x i ) = f ( X i , y ( x i ) ) y^{\prime}\left(x_{i}\right)=f\left(X_{i}, y\left(x_{i}\right)\right) y′(xi)=f(Xi,y(xi)),再用向前差商近似代替导数则为:
( y ( x i + 1 ) − y ( x i ) ) h = f ( x i , y ( x i ) ) \frac{\left(y\left(x_{i}+1\right)-y\left(x_{i}\right)\right)}{h}=f\left(x_{i},y\left(x_{i}\right)\right) \quad h(y(xi+1)−y(xi))=f(xi,y(xi))
在这里, h \mathrm{h} h 是步长, 即相邻两个结点间的距离。因此可以根据 x i \mathrm{xi} xi 点和 y i \mathrm{yi} yi
点的数值计算出 y i + 1 y_{i}+1 yi+1 来:
y i + 1 = y i + h × f ( x i , y i ) , i = 0 , 1 , 2 , n y_{i+1}=y_{i}+h \times f\left(x_{i}, y_{i}\right), \quad i=0,1,2, n yi+1=yi+h×f(xi,yi),i=0,1,2,n
这就是欧拉格式, 若初值 y i + 1 y_{i}+1 yi+1 是已知的, 则可依据上式逐步算出数
值解 y 1 , y 2 ⋯ ⋯ … y n y_{1}, y_{2} \cdots \cdots \ldots y n y1,y2⋯⋯…yn
X ( k + 1 ) = X ( k ) + Δ X ( k + 1 ) Δ X ( k + 1 ) = e h ( k ) ⋅ T = X ˙ ( k ) ⋅ T = [ A X ( k ) + B U ( k ) ] T \begin{aligned} &X(k+1)=X(k)+\Delta X(k+1)\\ &\Delta X(k+1)=e_{h}(k) \cdot T \\ &=\dot{X}(k) \cdot T=[A X(k)+B U(k)] T \\ \end{aligned} X(k+1)=X(k)+ΔX(k+1)ΔX(k+1)=eh(k)⋅T=X˙(k)⋅T=[AX(k)+BU(k)]T
特点：单步，显式，一阶求导精度，截断误差为二阶。
缺点：简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大
2.梯形法
梯形法原理: 将向前欧拉公式中的导数 f ( x i , y i ) 改为微元两端 f ( x i , y i ) 和 f ( x i + 1 , y i + 1 ) 的平均, 即梯形公式。 \begin{aligned} &\text { 梯形法原理: }\\ &\text { 将向前欧拉公式中的导数 } \mathrm{f}(\mathrm{xi}, \mathrm{yi}) \text { 改为微元两端 } \mathrm{f}(\mathrm{xi}, \mathrm{yi}) \text { 和 }\\ &\mathrm{f}(\mathrm{xi}+1, \mathrm{y} \mathrm{i}+1) \text { 的平均, 即梯形公式。 } \end{aligned} 梯形法原理: 将向前欧拉公式中的导数 f(xi,yi) 改为微元两端 f(xi,yi) 和 f(xi+1,yi+1) 的平均, 即梯形公式。
将被积函数近似为直线函数，被积的部分近似为梯形，要求得较准确的数值，可以将要求积的区间分成多个小区间
e h = 1 2 ⋅ { e ( k T ) + e [ ( k + 1 ) T ] } e ( k ) = X ˙ ( k ) = A X ( k ) + B U ( k ) \begin{aligned} &e_{h}=\frac{1}{2} \cdot\{e(k T)+e[(k+1) T]\} \\ &e(k)=\dot{X}(k)=A X(k)+B U(k) \end{aligned} eh=21⋅{e(kT)+e[(k+1)T]}e(k)=X˙(k)=AX(k)+BU(k)
估计K+1时刻的状态变量 X 0 ( k + 1 ) = X ( k ) + T e ( k ) e ( k + 1 ) = A [ X ( k ) + T e ( k ) ] + B U ( k + 1 ) X ( k + 1 ) = X ( k ) + T 2 [ e ( k ) + e ( k + 1 ) ] \begin{aligned} &\text { 估计K+1时刻的状态变量 }\\ &X_{0}(k+1)=X(k)+T e(k)\\ &e(k+1)=A[X(k)+T e(k)]+B U(k+1)\\ &X(k+1)=X(k)+\frac{T}{2}[e(k)+e(k+1)] \end{aligned} 估计K+1时刻的状态变量 X0(k+1)=X(k)+Te(k)e(k+1)=A[X(k)+Te(k)]+BU(k+1)X(k+1)=X(k)+2T[e(k)+e(k+1)]
3.四阶Runge-kutta法
在区间 [ x n , x n + 1 [\mathrm{xn}, \mathrm{xn}+1 [xn,xn+1 ] 内多取几个点, 将他们的斜率加权平均，作为导数
的近似。
令初值问题表述如下。
y ′ = f ( t , y ) , y ( t 0 ) = y 0 y^{\prime}=f(t, y), \quad y\left(t_{0}\right)=y_{0} y′=f(t,y),y(t0)=y0
则, 对于该问题的 RK4 由如下方程给出:
y n + 1 = y n + h 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) y_{n+1}=y_{n}+\frac{h}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right) yn+1=yn+6h(k1+2k2+2k3+k4)
其中
k 1 = f ( t n , y n ) k 2 = f ( t n + h 2 , y n + h 2 k 1 ) k 3 = f ( t n + h 2 , y n + h 2 k 2 ) k 4 = f ( t n + h , y n + h k 3 ) \begin{aligned} &k_{1}=f\left(t_{n}, y_{n}\right) \\ &k_{2}=f\left(t_{n}+\frac{h}{2}, y_{n}+\frac{h}{2} k_{1}\right) \\ &k_{3}=f\left(t_{n}+\frac{h}{2}, y_{n}+\frac{h}{2} k_{2}\right) \\ &k_{4}=f\left(t_{n}+h, y_{n}+h k_{3}\right) \end{aligned} k1=f(tn,yn)k2=f(tn+2h,yn+2hk1)k3=f(tn+2h,yn+2hk2)k4=f(tn+h,yn+hk3)
这样, 下一个值 ( y n + 1 ) \left(y_{n+1}\right) (yn+1) 由现在的值 ( y 2 ) \left(y_{2}\right) (y2) 加上时间间隔 ( h ) (h) (h) 和一个估算的
斜率的乘积决定。该斜率是以下斜率的加权平均:
k 1 k1 k1 是时间段开始时的斜率;
k2 是时间段中点的斜率，通过欧拉法采用斜率 k l \mathrm{kl} kl 来决定 y \mathrm{y} y
在点 t n + h / 2 \mathrm{tn}+\mathrm{h} / 2 tn+h/2 的值:
k 2 \mathrm{k} 2 k2 是时间段中点的斜率, 通过欧拉法采用斜率 k 1 \mathrm{k} 1 k1 来决定 y \mathrm{y} y
在点 t n + h / 2 \mathrm{tn}+\mathrm{h} / 2 tn+h/2 的值:
k3 也是中点的斜率, 但是这次采用斜率 k 2 \mathrm{k} 2 k2 决定 y \mathrm{y} y 值;
k 4 \mathrm{k} 4 k4 是时间段终点的斜率，其 y \mathrm{y} y 值用 k 3 \mathrm{k} 3 k3 决定。
当四个斜率取平均时，中点的斜率有更大的权值:
slope = k 1 + 2 k 2 + 2 k 3 + k 4 6 \text { slope }=\frac{k_{1}+2 k_{2}+2 k_{3}+k_{4}}{6} slope =6k1+2k2+2k3+k4
R K 4 \mathrm{RK} 4 RK4 法是四阶方法, 也就是说每步的误差是 h 5 h^{5} h5 阶, 而总积係误差为 h 4 h^{4} h4
阶。
e 1 = A X ( k ) + B U ( k ) e 2 = A X ( k + 1 2 ) + B U ( k + 1 2 ) e 3 = A X ( k + 1 2 ) + B U ( k + 1 2 ) x 0 ( k + 1 2 ) = x ( k ) + T 2 e 1 e 4 = A X ( k + 1 ) + B U ( k + 1 ) x 0 ( k + 1 2 ) = x ( k ) + 1 2 e 2 x 0 ( k + 1 ) = x ( k ) + T e 3 \begin{aligned} &e_{1}=A X(k)+B U(k) \\ &e_{2}=A X\left(k+\frac{1}{2}\right)+B U\left(k+\frac{1}{2}\right) \\ &e_{3}=A X\left(k+\frac{1}{2}\right)+B U\left(k+\frac{1}{2}\right) \quad x_{0}\left(k+\frac{1}{2}\right)=x(k)+\frac{T}{2} e_{1} \\ &e_{4}=A X(k+1)+B U(k+1) & \begin{array}{l} x_{0}\left(k+\frac{1}{2}\right)=x(k)+\frac{1}{2} e_{2} \\ x_{0}(k+1)=x(k)+T e_{3} \end{array} \\ &\qquad \begin{aligned} & \end{aligned} \end{aligned} e1=AX(k)+BU(k)e2=AX(k+21)+BU(k+21)e3=AX(k+21)+BU(k+21)x0(k+21)=x(k)+2Te1e4=AX(k+1)+BU(k+1)x0(k+21)=x(k)+21e2x0(k+1)=x(k)+Te3
4.离散相似法
离散相似法的原理
用周期为的虚拟采样开关将连续模型的输入、输出分别离散化，要求离散化后的输出 在采样时刻的值等同于原输出在同一时刻的值，以后的每一步计算均在这个模型基础上进行，而原来的连续模型不再参与计算，如图1所示。对比数值积分法虽然也进行了离散化处理，但在离散化过程中每一步都用到连续系统的模型(导函数 )，离散一步计算一步。
在系统输入输出处加上采样开关
开关的频度为仿真步距
为了重现信号应再加上保持器与补偿器
保持器类型
零阶保持器：
在信号传递过程中,把第nT时刻的采样信号值一直保持到第(n+1)T时刻的前一瞬时,把第(n+1)T时刻的采样值一直保持到(n+2)T时刻,依次类推,从而把一个脉冲序列变成一个连续的阶梯信号，因为在每一个采样区间内连续的阶梯信号的值均为常值,亦即其一阶导数为零,故称为零阶保持器
响应状态
u ~ ( k T + Δ t ) = a 0 + a 1 Δ t + a 2 Δ t 2 + … + a n Δ t n \tilde{u}(k T+\Delta t)=a_{0}+a_{1} \Delta t+a_{2} \Delta t^{2}+\ldots+a_{n} \Delta t^{n} u~(kT+Δt)=a0+a1Δt+a2Δt2+…+anΔtn
上式称为n阶外推公式，代表的是n阶保持器。
三、实验过程
仿真时间及仿真步距的估计
S T = ( 5 ∼ 20 ) n T S T=(5 \sim 20) n T ST=(5∼20)nT D T = n T 50 ∼ n T 10 D T=\frac{n T}{50} \sim \frac{n T}{10} DT=50nT∼10nT
ST——仿真时间;
DT------仿真步距;
n———被控对象传递函数的阶次;
T———被控对象传递函数的时间常数。
如果被控对象有若干个，ST中的nT则应以其中最大的为准,DT中的nT以最小则为准
1.整体离散法
我们把采样开关和零阶保持器加在系统入口处，得到的离散相似系统。
其输入矩阵，输入矩阵，输出矩阵。
使用MATLAB绘制其仿真曲线，得图三.
clear all;
DT=1;ST=80;LP=ST/DT;R=1;
p=0;x=[0;0]
%近似求取F(T)和Fm (T)
A=[-0.3 -0.2;1 0];B=[1;0];
n=6;m=2;
df1=DT*A;
df=eye(m);
fai=eye(m);
for i=1:ndf=df*df1/i;fai=fai+df;
end
dfm=eye(m);
faim=eye(m);
for i=1:n-1dfm=dfm*df1/(i+1);faim=faim+dfm;
end
faim=DT*faim*B;for i=1:LP%近似求取F(T)和Fm (T) 到 T的6次项x=fai*x+faim*R;y(i)=0.2*x(1)+0.2*x(2);t(i)=DT*i;
endplot(t,y,'k')
grid on
2.分环节离散方法
close all;
clear all;
clc;
K1 = 2;
T1 = 10;kp = 1;
ki = 1;
u = 0;
x = zeros(2,1);
dt = 1;
st = 50;
t = [];y = [];
obj1a = exp(-dt/T1);
obj1b = K1*(1 - obj1a);
for i = 0:dt:ste = 1 - x(2);x(1) = x(1) + ki*dt*e;u = kp*e + x(1);x(2) = obj1a * x(2) + obj1b * u;t = [t i];y = [y x(2)];
end
plot(t,y,'k-','Linewidth',1);
grid on;
3.数值积分法
我们分别用欧拉法、梯形法、4阶Runge-kutta用下列方法进行仿真，得到图五三条曲线。
clear all;
DT=1;ST=80;LP=ST/DT;
p0=1
p=0;p1=0;p2=0;p4=0;
x=[0 0]';x1=x;x2=x;x4=x;
A=[-0.3,-0.2;1,0]
B=[1;0];
%计算精确解差分方程系数F(T)和Fm (T)
n=6;m=2;
df1=DT*A;
df=eye(m);
fai=eye(m);
for i=1:ndf=df*df1/i;fai=fai+df;
end
dfm=eye(m);
faim=eye(m);
for i=1:n-1dfm=dfm*df1/(i+1);faim=faim+dfm;
end
faim=DT*faim*B;
for i=1:LP%计算精确解x=fai*x+faim*p0;y(i)=0.2*x(1)+0.2*x(2);t(i)=i*DT;p=p+abs(x(2))*DT;%计算欧拉公式x1=x1+DT*(A*x1+B*p0);y1(i)=0.2*x1(1)+0.2*x1(2);p1=p1+abs(x1(2))*DT;%计算梯形公式e21=A*x2+B*p0;x20=x2+DT*e21;e22=A*x20+B*p0;x2=x2+DT/2*(e21+e22);y2(i)=0.2*x2(1)+0.2*x2(2);p2=p2+abs(x2(2))*DT;%计算四阶龙格－库塔公式e41=A*x4+B*p0;x42=x4+DT/2*e41;e42=A*x42+B*p0;x43=x4+DT/2*e42;e43=A*x43+B*p0;x44=x4+DT*e43;e44=A*x44+B*p0;x4=x4+DT/6*(e41+2*e42+2*e43+e44);y4(i)=0.2*x4(1)+0.2*x4(2);p4=p4+abs(x4(2))*DT;end
ep1=abs(p-p1)/p*100 %计算欧拉公式的相对作用强度
ep2=abs(p-p2)/p*100 %计算梯形公式的相对作用强度
ep4=abs(p-p4)/p*100 %计算龙格－库塔公式的相对作用强度
%输出仿真结果
plot(t,y,'r',t,y1,'g',t,y2,'b',t,y4,'k')
legend('精确值','欧拉法','梯形法','四阶龙格－库塔法')
四.仿真分析
相对作用强度
由于每个离散-再现环节都对原来的信号进行了近似，所以只使用一次的离散-再现环节的整体离散系统是较为准确的。所以我们将在系统的入口处使用零阶保持器时得到的差分方程，取其系数到 T 的 6 次项作为系统的精确解。与其他仿真方法进行对比。
我们延用函数的相对作用强度概念来定量分析系统的仿真精度。
把系统在过渡过程时间段内，系统输出曲线下的绝对面积称为这个系统的作用强度。即，系统输出 y(t) 的作用强度为
式中： p——作用强度； ts ——考虑的时间段。
假设系统在精确解下的作用强度为 p0 ，数值解下的作用强度为 p ，则系统数值解的相对作用强度为
从仿真运行结果可知：
欧拉公式的相对作用强度： ep1 = 0.6708%
梯形公式的相对作用强度： ep2 = 0.0353%
四阶龙格－库塔公式的相对作用强度：ep3 = 3.9009e-04%
分环节相似离散的相对作用强度：ep4 =0.0793%
欧拉公式仿真结果较差、而梯形公式、四阶龙格－库塔公式以及分环节离散化都有很高的仿真精度，特别是4 阶龙格－库塔公式已经非常接近精确解。
仿真步长的影响
我们将步长分别扩大到2和缩小到0.5，计算此时的相对作用强度并画出对应的曲线。
从仿真运行结果可知：
欧拉公式的相对作用强度： ep1 = 664.7106%
梯形公式的相对作用强度： ep2 = 0.2815%
四阶龙格－库塔公式的相对作用强度：ep3 = 7.8156e-04%
分环节相似离散的相对作用强度：ep4 =0.8642%
从仿真运行结果可知：
欧拉公式的相对作用强度： ep1 = 0.3179%
梯形公式的相对作用强度： ep2 = 0.0042%
四阶龙格－库塔公式的相对作用强度：ep3 = 3.2835e-05%
分环节相似离散的相对作用强度：ep4 =0.0573%
【自动控制原理仿真实验】 控制系统仿真实验（实验二）相关推荐
- matlab仿真实训要求,南昌大学《MATLAB与控制系统仿真》实验报告
<南昌大学<MATLAB与控制系统仿真>实验报告>由会员分享,可在线阅读,更多相关<南昌大学<MATLAB与控制系统仿真>实验报告(36页珍藏版)>请在 ...
- 基于51单片机宠物自动投料喂食器控制系统仿真设计( proteus仿真+程序+讲解视频）
基于51单片机宠物自动投料喂食器控制系统仿真设计( proteus仿真+程序+讲解视频) 仿真图proteus 7.8及以上 程序编译器:keil 4/keil 5 编程语言:C语言 设计编号:S00 ...
- 控制系统仿真技术（二）-连续系统的数字仿真二
太原理工大学控制系统仿真技术实验报告 连续系统的数字仿真 1.分别利用欧拉法和预估-校正法求下图所示系统的阶跃响应,并对其结果进行比较. %欧拉法求阶跃响应 r=2;num0=8;den0=[1 3 ...
- 【自动控制原理仿真实验】 稳定性及稳态误差实验（实验三）
稳定性及稳态误差实验(实验三) 一.实验要求 二.实验简介 稳态误差的分类 1.原理性误差 2.实际性误差 三.实验过程 1.控制系统特征方程与稳定性的关系 2.分析给定输入的系统稳态误差影响 0型二 ...
- 神经网络控制simulink仿真,神经网络控制系统仿真
关于神经网络仿真的一些概念问题 40 1.常用的有sigmoid型函数.tansig函数.logsig函数等.采用不同函数,神经网络的运算效果不同.实际问题中,函数的选择是根据试验结果决定的,也就是试 ...
- lorenz系统simulink仿真_simulink控制系统仿真之控制系统的分析方法（2）（频域分析法）...
3.频域分析法1.频率响应曲线 频率特性曲线包括三种常用形式:Nyquist曲线(极坐标曲线),Bode图(对数坐标图)和Nichols图(对数幅相图). (1)Nyquist曲线 极坐标图,又称奈奎 ...
- matlab求临界稳定时的k,MATLAB自动控制原理仿真
一.某控制系统结构图如图所示, (1) 试用SIMULINK 建立系统仿真模型,且该系统中K=1保存路径为:E :\lsfz : (2) 利用所建立的SIMULINK 仿真模型求该系统闭环传递函数及开 ...
- 南昌大学matlab实验2,南昌大学MATLAB与控制系统仿真实验报告-资源下载人人文库网...
南昌大学<MATLAB与控制系统仿真>实验报告 实 验 报 告实验课程: MATLAB与控制系统仿真姓 名:学 号:专业班级: 2016年 6月目 录实验一 MATLAB的环境与基本运算( ...
- 传感器实验——控制电机
传感器实验--控制电机 电机实验 所选设备 12V直流电机 使用说明 小伙伴们,玩过4驱车吗?4驱车上动力是谁?没玩过也不要紧,电机听过没有~今天我们来电机. 直流电机(direct currentm ...
最新文章
- Invoking Page() in async task.
- 嘉兴新型智慧城市建设带来的三个问号
- ABAP 程序中退出操作(CHECK, EXIT, RETURN, LEAVE PROGRAM...)
- 为什么薄膜干涉的厚度要很小_薄膜干涉的薄膜为什么不能太厚？1.这里的厚是用什么来衡量的？2.假如一个厚为1mm的薄膜，是否薄？...
- PHP学级与年级的转换函数_PHP addslashes()和stripslashes()：字符串转义与还原
- 只读域控制器在Server Core中的部署
- 使用TargetSources
- [Unity3D]脚本中Start()和Awake()的区别
- mapinfo在线地图插件_官方插件“战争游戏”使用指南 公测同步开启
- Spartan-6 FPGA 如何使用ISE下载程序
- IO基础操作（文件）
- wav转mp3，wav怎么转换成mp3？
- python动态仪表图_matplotlib仪表动态更新
- SSH框架 Bean property * is not writable or has an invalid setter method错误分析与解决方法
- MySQL#在指定的列中添加数据
- 秒杀项目05-页面优化技术
- 教你如何在线播放FLV格式的文件
- 御龙在天怎么找回服务器,御龙在天手游人物找回 误删的角色如何找回
- 毕业季基于spring的基于安卓APP的基于ssm框架的基于微信小程序的管理系统设计与开发（开题+源码+讲解+论文）
- VS2019 C语言，在一个项目中添加多个包含main函数的源文件并分别调试运行