控制系统仿真实验（实验二）
一、实验目的
①构建与实验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:n
df=df*df1/i;
fai=fai+df;
end
dfm=eye(m);
faim=eye(m);
for i=1:n-1
dfm=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;
end
plot(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:st
e = 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:n
df=df*df1/i;
fai=fai+df;
end
dfm=eye(m);
faim=eye(m);
for i=1:n-1
dfm=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%
4098
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
