系统建模与仿真作业

作业一：统计试验法解决报童问题

1.问题提出

报童问题是决策方面的一个简单例子。对一个报童来说，他每天由报社买来报纸，然后到街上去卖，当然他希望获得最大的利润。如果把报童问题看作一个系统，那么报纸、顾客、报童和利润就成了该系统的重要组成部分。问题在于，根据市场需求，报童寻求一个什么样的定货法则或策略才能使他所获得的利润最大。因为定货多了，如果市场需求小，卖不了将导致利润下降，甚至亏本；定货少了，如果市场需求量大，又失去了赚钱的机会。这样一来，定货法则或定货策略就成了报童能否赚取最大利润的关键。当然，最好是每天需要多少份就定多少份，但市场需求量是一个随机变量，这就需要每天用统计方法作出决策。

报童以每份0.5元的价格买进报纸，以0.8元的价格出售。根据长期统计，报童根据以前的卖报记录知道，每天的需求量有几种可能及出现的概率。对现有数据分析，得出报童每天最佳买进报纸量，使报童的平均总收入最大。

2.模型建立与求解

假定，B n—当天买进或订购的报纸数量

S n—当天社会需要报纸的数量，显然，它是一个随机变量

D n—当天卖掉的报纸数量

再假定，报童每天买进和卖出每份报纸的价格分别用P B和P S表示，且P S>P B，即卖出价大于买入价，则第n 天的利润为

P n=S n P S-B n P B

报童决定：当天订购报纸的数量等于前一天的市场需求量，即

B n=D n-1

而当天卖掉的报纸的数量n S 则由以下两个条件来决定

D n≤ B n时，S n=D n

D n> B n时，S n=B n

即是说，如果当天订购的报纸的数量大于需求量，当天卖掉的报纸的数量只能等于需求量；如果当天买的报纸的数量小于需求量，即当天卖掉的报纸的数量

等于当天订购量和当天需求量二者中的小者。

现在的问题是，前一天的需求量如何决定。报童决定利用过去一年的统计数字确定D n-1。报童根据以前的卖报记录知道，每天的需求量有以下几种可能：40份，41份，42份，43份，44份，45份，46份；并且统计出了相对频数，得到如右表的一组数据。

其需求量的平均值D n=∑ i P i=43. 10

最后，报童做了一个轮盘，并将其分成了七份，每份的大小分别等于每个需求量对应的频数，即需求量分别为40，41，…,46份报纸时,其轮盘上对应的面积分别为0.05,0.1,…,0.1类推。

这样，报童每天去订货之前转一次轮盘，指针所指的数量就作为前一天的需求量D n-1。假定，第一天转了一次，D n-1=45，即D0=45。作为第二天买报的依据；第二天又转了一次，D1 =44表示当天需求量，说明这一天订购44份，由于B n>D n，所以只卖掉44份，即S n=D n。……第六天又转了一次，D5 =41表示当天的需求量，说明这一天订购40份，由于B n

3.实验结果

报童问题的10天的仿真结果

表采用不同的决策准则时报童所获得的利润

附录程序

#include

#include

#define N 17 //天数#define Ps 0.8 //卖出价#define Pb 0.5 //进货价

int random()

{

int Dn;

float temp;

temp=rand()*100/(float)RAND_M AX;

if(temp<5)

Dn=40;

else if(temp<15)

Dn=41;

else if(temp<35)

Dn=42;

else if(temp<65)

Dn=43;

else if(temp<80)

Dn=44;

else if(temp<90)

Dn=45;

else

Dn=46;

return Dn;

}

void main()

{

int D[N+1],B[N+1],S[N+1];

int i;

double P[N+1];

double pn=0;

int j;

int a=0;

printf("D B S P

pn\n");

for(i=0; i<=N; i++)

{

D[i]=random();

}

for(i=0; i<=N; i++)

{

D[i]=random();

/***1*****B n=D n-1************ ***/

// if(i!=N) B[i+1]=D[i];

/***2***B n=

（D n-1+D n-2)/2*********/

/* if(i==0)B[i+1]=D[i];

if(i!=N)

B[i+2]=(D[i]+D[i+1])/2;*/

/***3****** B n=

43**************/

// B[i]=43;

/***4****B n=(D0+D1+…+D n-1)]/n ***/

for(j=0,a=0;j

{

a=a+D[j];

}

if(i!=0) B[i]=a/i;

///////////////////////////////////////////////// if(D[i]

S[i]=D[i];

else

S[i]=B[i];

if(i!=0)

{

P[i]=S[i]*Ps-B[i]*Pb;

pn+=P[i];

printf("%d %d %d %.1f %. 1f\n", D[i],B[i],S[i],P[i],pn);

}

}

作业二：产生独立分布随机数和相关分布随机数

1．（0-1）之间的均匀分布：

概率密度函数： ??

?≤≤=其他01

01)(x x P ； 产生思想：采用乘同余法产生； 具体实现方法：x n+1=ux n (mod n)；

参数：u 一般取23

a ±3或5，a 为正整数；x 0为初始值一般取2b+1，

b 取正整数；

m 一般取计算机的字长，其是控制所产生随机数的精度（即：小数点后的位数）；

程序实现中取u=11，m=2^32，x0的取值是3，计算10000点随机数。 程序如下：

c=3; u =11; x(1)=3; M=2^32 for i=2：1：10000;

x(i)=mod(u *x(i-1)+c,M); end ; x=x./M; e=mean(x); d=var(x); hist(x,20);

title('（0-1）均匀分布直方图');

text(0,-60,strcat('均值是',num2str(e)));

text(0.5,-60,strcat('方差是',num2str(d)));

运行结果如下：

运行结果为：均值为0.49937，方差为0.083626表示所产生数据的方差，而（0-1）之间的均匀分布的随机数的数学期望为0.5，与上面所求出的0. 49937很接近，方差0. 083626近似与0，于是这种产生方法已经符合要求。

2．高斯分布：

高斯分布的概率密度函数如下： 2

2

2)

(21

)(σ

σπu x e x p --=

其产生方法是在均匀分布随机数的基础上通过函数变换法来产生。产生步骤是①产生均匀分布的随机数。②产生服从标准正态分布的随机数。③由标准正态分布产生一般正态分布。

)

2cos(ln 2i i i v u z π-=

1）参数：μ=0，σ=0 ，计算10000点随机数。 程序如下：

x=rand(1,10000); y=rand(1,10000); i=1：1：10000;

z(i)=sqrt(-2*log(x(i))).*cos(2.*pi.*y(i)); w(i)=sqrt(-2*log(x(i))).*sin(2.*pi.*y(i));

e=mean(z)

d=var(z)

n=-5：0.05：5;

subplot(1,2,1);

hist(z,n);

title('正态分布直方图-1');

text(-6,-20,strcat('均值是',num2str(e)));

text(2,-20,strcat('方差是',num2str(d)));

subplot(1,2,2);

hist(w,n);

e=mean(w)

d=var(w)

title('正态分布直方图-2');

text(-6,-20,strcat('均值是',num2str(e)));

text(2,-20,strcat('方差是',num2str(d)));

运行结果如下：

运行结果为：（1）均值为-0.0069，方差为 1.0014，标准高斯分布数学期望为0，与上面所求出的-0.0069很接近；方差1，与1.0014近似，于是这种产生方法已经符合要求。

2）参数：μ=2，σ=2，计算10000点随机数。

程序如下：

一般高斯分布

a=2;b=2;

x=rand(1,10000);

y=rand(1,10000);

i=1：1：10000;

z(i)=sqrt(-2*log(x(i))).*cos(2.*pi.*y(i)); v=b*u+a;

e=mean(v)

d=var(v)

n=-10：0.2：10;

hist(v,n);

title('一般正态分布直方图');

text(-6,-20,strcat('均值是',num2str(e))); text(2,-20,strcat('方差是',num2str(d)));

运行结果如下：

运行结果为：在横坐标为2附近得到峰值；均值为-2.0202，方差为3.7686，一般高斯分布数学期望为2，与上面所求出的-2.0202很接近；方差4，与3.7686接近，于是这种产生方法已经符合要求。

3．瑞利分布随机数的产生

概率密度函数：

??

?

??<≥=-000)(2

2

22x x e

x x P x σσ；

产生思想：利用直接抽样法产生； 具体实现方法：

a ．先调用产生（0-1）之间的均匀分布的函数（y==rand(n,1)）；

b ．然后作2ln z y =-；

c ．令z y σ=，于是向量y 就是要产生的瑞利分布的随机数； 计算10000点随机数。 程序如下：

x=rand(10000,1); for i=1：10000 if (x(i)==0) x(i)=0.0001; end end

u=(-2)*log(x); y=1*sqrt(u); e=mean(y); d=var(y); hist(y,20);

title('瑞利分布直方图');

text(0,-120,strcat('均值是',num2str(e)));

text(3,-120,strcat('方差是',num2str(d))); 运行结果如下：

均值和方差为：1.2448和0.42292，

而瑞利分布的数学期望计算式为：

σ

其中σ=1，代入计算得：1.253，与上面所求出的随机数的平均值1. 2448相当

接近，瑞利分布方差的计算公式为：2

24σ

π-当1σ=时代入计算得0.42920与

0. 42292相当接近。

4．指数分布随机数的产生

概率密度函数：

?????<≥=-0

01)(x x e

x P x

λλ；

产生思想：利用直接抽样法； 具体实现方法：

先调用产生（0-1）之间的均匀分布的函数（x=rand(n,1)）；

然后作x

y ln 1

λ

-

=（λ为参数）于是向量y 就是所要产生的指数分布的随机

向量；

参数估计：λ=1，计算10000点随机数。

程序如下：

n=10000;

m=1;

x=rand(n,1);

for i=1;n

if (x(i)==0)

x(i)=0.0001;

end

end

u=log(x);

y=-(1/m)*u;

e=mean(y);

d=var(y);

hist(y,50);

title('指数分布直方图');

text(0,-170,strcat('均值是',num2str(e))); text(3,-170,strcat('方差是',num2str(d))); 运行结果如下：

均值和方差为：0.99154和0.97058，而指数分布的数学期望计算式为：1/m ＝1，与上面所求出的随机数的平均值0. 99154相当接近，指数分布方差的计算公式为：(1/m)的平方，当m=1时代入计算得1与0. 97058相当符合。

作业三：雷达系统仿真

雷达系统仿真流程为：

雷达目标回波：

这里的雷达回波由三部分组成即：目标回波、地物杂波、系统噪声。其中目标回波这里采用线性调频信号模拟，地物杂波采用相关对数-正态分布的随机数，系统噪声采用高斯白噪声。回波信号由目标回波（动目标），地物杂波（静目标），及系统高斯白噪声组成。

线性调频信号：x=rect(t/mk)exp(j πkt 2) (k=B/mk)

目标回波：y ＝rect(t/mk)*exp(j*2*pi*((f1+k*t/2).*t+fd*i*T)) 地物杂波（静目标）：y ＝rect(t/mk)*exp(j*2*pi*((f1+k*t/2).*t)) 系统噪声（高斯白噪声）： z ＝0.2*randn(1,N)。

1．雷达的回波模型：

回波

加噪声的目标回波的图形：

2．放大器：

其中高频放大器的增益为k ，带通滤波器的设计采用的是FIR 滤波器，其中通频带的带宽取为信号的带宽B 的5/3倍,过渡带的宽度为B ，采用汉明窗进行加权，数字滤波器的通频带和阻带的值利用T w Ω=（w 为数字角频率，Ω为模拟角频率，T 为采样周期）计算。

高放图：

3．混频器：

X （t

])(2cos[21

])(2cos[212121t f f t f f -++ππ，这里取带通滤波器的通带的中心频率为2

1f f -于是通过带通滤波器以后的输出为：])(2cos[21

21t f f -π，于是目标回波

经高频放大器以后通过混频器后的时域图形频域以及滤波器的频率响应分别如下图所示：

4．中频放大器的设计：

中频放大器的设计与高频放大器的设计除滤波器的通带范围不同以外其余的均相同，模型框图与高频放大器完全相同，回波经混频以后通过中频放大器后的时域图形频域以及滤波器的频率响应分别如下图所示：

5．相位检波器：

低通滤波器的设计与前面的滤波器设计相类似，输出的两路信号就是I、Q 两路信号，通过检波后信号的时域图形、包络以及频率域的图形如下：

6．A/D变换器

采样频率为5MHz，x0=(Vmax/2a)*int{xi*2a / Vmax }，其中a为AD位数

7．脉冲压缩

脉冲压缩就是把输入的I、Q两路信号中I路信号作为实部Q路信号作为虚部和起来对其压缩处理。脉冲压缩后的图形：

8．对消器（MTI）：

MTI采用一次对消：y(n)=x(n)-x(n-1)；n=1,2,3…N

对消后的图形如下：

9．取模：

()

2

/12

2Q

I x x R +=，取模后的图形如下：

10．积累（视频积累，非相参）：单极点：()()()()11-+-=n ky n x k n y ；双极点：

()()()()21121---+-=n y k n y k n x n y ，积累后的图形如下：

11．恒虚警（CFAR）：

其程序如下：

function [t,s,g,f0,fs,f1]=huibo%产生目标回波信号x,系统噪声y，地物杂波z以及回波p

f0=3*10^7; %发射信号频率

w=0; %发射信号初始相位

c=3*10^8; %光速

l=c/f0; %雷达信号波长(载波波长)

R=40000; %目标范围

Vd=200; %雷达与目标之间的径向速度

fd=2*Vd/l; %多普勒频率

Tr=600/f0; %脉冲重复周期

N=10; %雷达脉冲串长度

f1=f0/5; %调频带宽是发射信号频率的1/5

k=10*f1/Tr;

fs=3*f0; %仿真采样频率

Ts=1/fs;

Btar=4*pi*R/l;

M=floor(Tr*fs); %一个脉冲重复周期内的采样点数M=600

mt1=floor(2*Tr*fs/5);%mt1=720

mt2=floor(3*Tr*fs/10);%mt2=540

mt3=floor(7*Tr*fs/10);%mt3=1260

mt4=floor(3*Tr*fs/5);%mt4=1080

mt5=mt1-mt2;%mt5=180

Vgain=6;

s=zeros(M,N); %回波幅度起伏

for m=1：M

for n=1：N

v(m,n)=(u(mt1-m)-u(mt2-m))*Vgain*cos(2*pi*f0*(m-mt2)*Ts+2*pi*k*Ts*(m-mt2)/2*(m-mt2)*Ts+2*pi*fd*n*Tr);

g(m,n)=(u(mt3-m)-u(mt4-m))*Vgain*cos(2*pi*(f0*(m-mt4)*Ts+k*Ts*(m-mt4) /2*(m-mt4)*Ts));

end

end

for i=1：M/10

t(i)=Vgain*cos(2*pi*(k*Ts*(M/10-i)/2*(M/10-i)*Ts));