1.蒙特卡洛仿真方法的统计特性
假设一个 m m m个系统输出数据 { y i } i = 1 m \{y_i\}_{i=1}^m {yi}i=1m, N N N次循环得到 N N N组数据 { { y i } i = 1 m } j = 1 N \{\{y_i\}_{i=1}^m\}_{j=1}^N {{yi}i=1m}j=1N。那么实际上会有以下两组统计特性指标:
1.2数值统计特性:
(1)数学期望: E ( y ) ≈ 1 N ∑ k = 1 N y k E(y) \approx \frac{1}{N}\sum_{k=1}^Ny_k E(y)≈N1∑k=1Nyk
(2)方差: σ y 2 ≈ 1 N ∑ K = 1 N ( y k − E y ) 2 \sigma_y^2 \approx \frac{1}{N}\sum_{K=1}^N{(y_k-Ey)}^2 σy2≈N1∑K=1N(yk−Ey)2
(3)概率密度函数: p ( y ) ≈ N y N Δ y p(y)\approx\frac{N_y}{N\Delta y} p(y)≈NΔyNy,其中若将区间 ( a , b ) (a,b) (a,b)划分为 L L L段,则 Δ y = b − a L \Delta y=\frac{b-a}{L} Δy=Lb−a,其中 N y N_y Ny是数值满足不等式:
( y − 0.5 Δ y ) < y < ( y + 0.5 Δ y ) (y-0.5\Delta y)<y<(y+0.5\Delta y) (y−0.5Δy)<y<(y+0.5Δy)的个数,用曲线拟合可以得到概率密度函数 p ( y ) p(y) p(y)
1.3 时间统计特性计算(Y(t)是平稳随机过程):
(1)时间均值 E { y K ( t ) } = l i m T → ∞ 1 T ∫ 0 T y K ( t ) d t E\{y_K(t)\}=lim_{T\rightarrow \infty}\frac{1}{T}\int_0^Ty_K(t)dt E{yK(t)}=limT→∞T1∫0TyK(t)dt
(2)时间均方差 σ y K 2 = l i m T → ∞ 1 T ∫ 0 T { y K ( t ) − E { y K ( t ) } } 2 d t \sigma_{y_K}^2=lim_{T\rightarrow \infty}\frac{1}{T}\int_0^T\{y_K(t)-E\{y_K(t)\}\}^2dt σyK2=limT→∞T1∫0T{yK(t)−E{yK(t)}}2dt
(3)自相关函数 R y K ( t 1 , t 2 ) = l i m T → ∞ 1 T ∫ 0 T y K ( t ) y K ( t + τ ) d t R_{y_K}(t_1,t_2)=lim_{T\rightarrow \infty}\frac{1}{T}\int_0^Ty_K(t)y_K(t+\tau)dt RyK(t1,t2)=limT→∞T1∫0TyK(t)yK(t+τ)dt
2.仿真实例
以导弹末端制导系统仿真为例,说明该仿真方法的应用
2.1仿真原理及算法
1.系统组成
一个导弹末端制导系统由控制系统,导引头,目标运动,导弹运动和相对运动等若干环节组成,这些环节构成了末端制导系统(如下图所示)。符号说明:
R
:
R:
R:视线距,
λ
:
\lambda:
λ:视线角,
a
c
a_c
ac:控制指令,
a
L
a_L
aL:导弹控制系统对制导指令的响应(实际视线法向加速度)。采用比例制导法,那么
a
c
=
K
v
c
λ
˙
a_c=Kv_c\dot{\lambda}
ac=Kvcλ˙,其中
K
K
K为比例导引系数,
v
c
v_c
vc为导弹实际接近目标的速度,
λ
˙
\dot{\lambda}
λ˙为目标的实际角速度。
2.工作原理
由视线角速度 λ ˙ \dot\lambda λ˙进过制导计算机形成制导指令,导弹控制系统执行制导指令,改变导弹的姿态运动并且给出导弹视线法向加速度,影响导弹同目标运动的相对运动,使得 R ( t F ) R(t_F) R(tF)(即在一段飞行时间内导弹的视线距尽可能小)。
仿真目的在于通过仿真活得脱靶量 y ( t F ) y(t_F) y(tF)在不同的 m m m个飞行时间内由于随机目标机动运动造成的均方根 σ y ( 1 ) , σ y ( 2 ) , σ y ( 3 ) . . . , σ y ( m ) \sigma_{y(1)},\sigma_{y(2)},\sigma_{y(3)}...,\sigma_{y(m)} σy(1),σy(2),σy(3)...,σy(m),以及相关的数值和时间统计特性。
3.制导系统的线性化模型
仿真过程中:
a
T
a_T
aT是目标的运动加速度是随机量,其代表的是一类“电报机动”的随机干扰模型,即目标要么正的最大加速度机动,要么负的最大加速度机动,机动开始时间和正负方向都是随机的。当运动规律确定的时候,开始时间随机的随机扰动可以描述为:
x
(
t
)
=
h
(
t
−
T
)
p
T
(
t
)
=
1
t
F
(
0
≤
t
≤
t
F
)
,
0
(
e
l
s
e
)
x(t)=h(t-T)\\ p_T(t)=\frac{1}{t_F}(0\leq t\leq t_F),0(else)
x(t)=h(t−T)pT(t)=tF1(0≤t≤tF),0(else)
在本例中,目标做随机机动,其视线法向加速度
a
T
a_T
aT的幅值为29.4m/s
,加速度的正负符号和开始时间
T
T
T为在区间
(
0
,
t
F
)
(0,t_F)
(0,tF)内的等概率的随机量。其描述如下:
a
T
(
t
)
=
C
o
e
f
∗
a
T
∗
s
i
g
n
(
t
−
T
)
C
o
e
f
=
1
(
R
a
n
d
o
m
>
0.5
)
,
−
1
(
r
a
n
d
o
m
<
0.5
)
,
p
T
(
t
)
=
1
t
F
(
0
≤
t
≤
t
F
)
,
0
(
e
l
s
e
)
a_T(t)=C_{oef}*a_T*sign(t-T)\\ C_{oef}=1(Random>0.5),-1(random<0.5),p_T(t)=\frac{1}{t_F}(0\leq t\leq t_F),0(else)
aT(t)=Coef∗aT∗sign(t−T)Coef=1(Random>0.5),−1(random<0.5),pT(t)=tF1(0≤t≤tF),0(else)
其中
R
(
0
)
R(0)
R(0)为视线距初值,
t
F
=
R
(
0
)
v
c
t_F=\frac{R(0)}{v_c}
tF=vcR(0),
y
(
t
F
)
y(t_F)
y(tF)定义为脱靶量,
t
g
o
=
t
F
−
t
t_{go}=t_F-t
tgo=tF−t为剩余飞行时间,
v
c
t
g
o
v_ct_{go}
vctgo为剩余飞行距离,视线角为
λ
=
y
v
c
t
g
o
\lambda=\frac{y}{v_ct_{go}}
λ=vctgoy。一些其他的参数标注为v_c=1219m/s
,K=3
,系统时间常数T=1s
,导弹速度v_m=914.4m/s
。
设
X
=
(
y
,
y
˙
,
a
c
,
a
L
)
T
X=(y,\dot{y},a_c,a_L)^T
X=(y,y˙,ac,aL)T系统的状态方程可以描述为:
X
˙
=
f
(
t
,
X
,
a
T
,
t
F
)
\dot{X}=f(t,X,a_T,t_F)
X˙=f(t,X,aT,tF)
其中:
X
˙
=
(
X
(
2
)
,
a
T
−
X
(
4
)
,
K
(
X
(
2
)
t
F
−
t
+
X
(
1
)
(
t
F
−
t
)
2
)
,
X
(
3
)
−
X
(
4
)
T
)
T
\dot X = (X(2),a_T-X(4),K(\frac{X(2)}{t_F-t}+\frac{X(1)}{(t_F-t)^2}),\frac{X(3)-X(4)}{T})^T
X˙=(X(2),aT−X(4),K(tF−tX(2)+(tF−t)2X(1)),TX(3)−X(4))T。
其流程图可以表述如下:
2.2 仿真代码及结果
clc,clear;close all;
v_c = 1219;K = 3;T = 1;
v_m = 914.4;a_T = 29.4;
coffient.K = K;
coffient.a_T = a_T;
t_f = 10;h = 0.1;
%重要参数区间
E = [];D = [];time_point = [];
for time = h:h:t_f
T = time*rand;%产生常数T
x0 = [0 0 0 0]';
%下面进行参数结构体赋值
coffient.t_F = time;
coffient.T = T;
[tout,yout] = Adams_4v(@(t,x)missle_model(t,x,coffient),[0 time],x0,(time - 0)/50);
%下面计算具体参数
time_point = [time_point;time];
E = [E;mean(yout(:,1))];
D = [D;sqrt(mean((yout(:,1) - mean(yout(:,1))).^2))];
end
%下面进行绘图操作
subplot(121)
hold on;box on;grid on;
plot(time_point,E,'bo');
xlabel('t_F/s');ylabel('E(y_P{t_F})/m');
title('数学期望仿真结果');
subplot(122)
hold on;box on;grid on;
plot(time_point,D,'r-+');
xlabel('t_F/s');ylabel('\sigma(y_P{t_F})/m');
title('方差仿真结果');
%导弹末端制导模型
function f = missle_model(t,x,coffient)
K = coffient.K;
t_F = coffient.t_F;
aT = coffient.a_T;
T = coffient.T;
f = zeros(4,1);
%下面给出随机a_T
if rand >= 0.5
C_oef = 1;
else
C_oef = -1;
end
if t>T
flag = 1;
else
flag = 0;
end
a_T = C_oef*aT*flag;
%下面是函数导数
f(1) = x(2);
f(2) = a_T - x(4);
f(3) = K*(x(2)/(t_F - t) + x(1)/(t_F - t)^2);
f(4) = (x(3) - x(4))/T;
end
%下面是4阶亚当姆斯预估-校正法
function [t_,y] = Adams_4v(f,tspan,x0,h)
t_min = tspan(1);t_max = tspan(2);
t_f1 = [t_min,t_min+h,t_min+2*h,t_min+3*h];
t_f2 = [t_min+h,t_min+2*h,t_min+3*h,t_min+4*h];
fn1 = zeros(1,4);fn2 = zeros(1,4);
y = [];t_ = [];x = x0;t = t_min;
%下面是用RK-4初始化提供初值
while(t<t_f1(end))
y = [y;x'];
t_ = [t_;t];
K1 = h*f(t,x);
K2 = h*f(t + h/2,x + 1/2*K1);
K3 = h*f(t + h/2,x + 1/2*K2);
K4 = h*f(t + h,x + K3);
x = x + (K1 + 2*K2 + 2*K3 + K4)/6;
t = t + h;
end
for t = (t_min+3*h):h:t_max
y = [y;x'];t_ = [t_;t];
x_ = x + h/24*(55*f(t,x)-59*f(t-h,x)+37*f(t-2*h,x)-9*f(t-3*h,x));
f_ = f(t+h,x_);
x = x + h/24*(9*f_ + 19*f(t,x) - 5*f(t-h,x) + f(t-2*h,x));
end
end
结果: