系统建模与仿真作业

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

2、份0.5元的价格买进报纸，以0.8元的价格出售。根据长期统计，报童根据以前的卖报记录知道，每天的需求量有几种可能及出现的概率。对现有数据分析，得出报童每天最佳买进报纸量，使报童的平均总收入最大。2.模型建立与求解假定， Bn 当天买进或订购的报纸数量Sn 当天社会需要报纸的数量，显然，它是一个随机变量 Dn 当天卖掉的报纸数量 再假定，报童每天买进和卖出每份报纸的价格分别用 PB 和PS表示，且 PS >PB ，即卖出价大于买入价，则第 n 天的利润为Pn=SnPS-BnPB报童决定：当天订购报纸的数量等于前一天的市场需求量，即 Bn=Dn-1而当天卖掉的报纸的数量 n S 则由以下两个

3、条件来决定Dn Bn 时， Sn=DnDn > Bn 时， Sn=Bn即是说，如果当天订购的报纸的数量大于需求量，当天卖掉的报纸的数量只能等于需求量；如果当天买的报纸的数量小于需求量，即当天卖掉的报纸的数量等于当天订购量和当天需求量二者中的小者。 现在的问题是，前一天的需求量如何决定。报童决定利用过去一年的统计数字确定Dn-1。报童根据以前的卖报记录知道，每天的需求量有以下几种可能：40份，41份，42份，43份，44份， 45份，46份；并且统计出了相对频数，得到如右表的一组数据。需求量Dn40414243444546 相对频数 Pn 0.050.10.20.30.150.10.1其需

4、求量的平均值Dn= i Pi=43. 10最后，报童做了一个轮盘，并将其分成了七份，每份的大小分别等于每个需求量对应的频数，即需求量分别为40，41，,46份报纸时,其轮盘上对应的面积分别为0.05,0.1,0.1类推。这样，报童每天去订货之前转一次轮盘，指针所指的数量就作为前一天的需求量Dn-1。假定，第一天转了一次，Dn-1=45，即D0=45。作为第二天买报的依据；第二天又转了一次，D1 =44 表示当天需求量，说明这一天订购44份，由于Bn>Dn，所以只卖掉44份，即Sn=Dn 。第六天又转了一次， D5 =41表示当天的需求量，说明这一天订购 40份，由于Bn<Dn所以当

5、天只卖掉40份，即 Sn=Bn ，依此类推。 3.实验结果报童问题的10天的仿真结果nDnBnDnBnDn>BnSnPnPn 04514445Y4412.712.724344Y 4312.425.134243Y4412.137.244042Y4011.048.254140Y4012.060.264341Y4112.372.574143Y4111.383.884241Y4112.396.194642Y4212.6108.7104346Y4311.4120.1表 采用不同的决策准则时报童所获得的利润N0 17天的利润100天的利润500天的利润报童使用的决策规则1208.3元1215.3元6

6、122.1元Bn=Dn-12207.0元1223.8元6159.8元Bn=（Dn-1+Dn-2)/23209.7元1242.8元6234.8元Bn= 43 4209.8元1241.1元6235.4元B1=D0B2=(D0+D1)/2Bn=(D0+D1+Dn-1)/n说明：Bn-每天买入的报纸数量Dn-每天社会的需求量PB-每份报纸的买入价Ps-每份报纸的卖出价过去累计的需求量及相应的概率：40 41 42 43 44 45 460.05 0.10 0.20 0.30 0.15 0.10 0.10附录程序#include<stdio.h>#include<stdlib.h>

7、;#define N 17 /天数#define Ps 0.8 /卖出价#define Pb 0.5 /进货价int random()int Dn;float temp;temp=rand()*100/(float)RAND_MAX;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;elseDn=46;return Dn;void main()intDN+1,

8、BN+1,SN+1;inti;doublePN+1;doublepn=0; int j;inta=0;printf("D B S P pnn");for(i=0; i<=N; i+)Di=random();for(i=0; i<=N; i+)Di=random();/*1*Bn=Dn-1*/if(i!=N) Bi+1=Di;/*2*Bn=（Dn-1+Dn-2)/2*/*if(i=0)Bi+1=Di;if(i!=N)Bi+2=(Di+Di+1)/2;*/*3* Bn= 43*/Bi=43;/*4*Bn=(D0+D1+Dn-1)/n*/for(j=0,a=0;j&l

9、t;i;j+)a=a+Dj;if(i!=0)Bi=a/i;/if(Di<Bi)Si=Di;elseSi=Bi;if(i!=0)Pi=Si*Ps-Bi*Pb;pn+=Pi;printf("%d %d %d %.1f %.1fn", Di,Bi,Si,Pi,pn);作业二：产生独立分布随机数和相关分布随机数1（0-1）之间的均匀分布：概率密度函数：；产生思想：采用乘同余法产生；具体实现方法：xn+1=uxn(mod n)； 参数：u一般取23 a±3或5，a为正整数；x0为初始值一般取2b+1，b取正整数；m一般取计算机的字长，其是控制所产生随机数的精度（即：小

10、数点后的位数）；程序实现中取u=11，m=232，x0的取值是3，计算10000点随机数。程序如下：c=3; u=11; x(1)=3; M=232for 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，方

11、差为0.083626表示所产生数据的方差，而（0-1）之间的均匀分布的随机数的数学期望为0.5，与上面所求出的0. 49937很接近，方差0. 083626近似与0，于是这种产生方法已经符合要求。2高斯分布：高斯分布的概率密度函数如下：其产生方法是在均匀分布随机数的基础上通过函数变换法来产生。产生步骤是产生均匀分布的随机数。产生服从标准正态分布的随机数。由标准正态分布产生一般正态分布。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.

12、*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

13、('均值是',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

14、=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瑞利分布随机数的产生概率密度函数：；产生思想：利用直接抽样法产生；具体实现方法

15、：a先调用产生（0-1）之间的均匀分布的函数（y=rand(n,1)）；b然后作；c令，于是向量就是要产生的瑞利分布的随机数；计算10000点随机数。程序如下：x=rand(10000,1);for i=1：10000 if (x(i)=0) x(i)=0.0001; endend 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('方差是'

16、,num2str(d);运行结果如下：均值和方差为：1.2448和0.42292，而瑞利分布的数学期望计算式为： 其中=1，代入计算得：1.253，与上面所求出的随机数的平均值1. 2448相当接近，瑞利分布方差的计算公式为：当时代入计算得0.42920与0. 42292相当接近。4指数分布随机数的产生概率密度函数：；产生思想：利用直接抽样法；具体实现方法：先调用产生（0-1）之间的均匀分布的函数（x=rand(n,1)）；然后作（为参数）于是向量就是所要产生的指数分布的随机向量；参数估计：=1，计算10000点随机数。程序如下：n=10000;m=1; x=rand(n,1);for i=1

17、;n if (x(i)=0) x(i)=0.0001; endend 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/m1，与上面所求出的随机数的平均值0. 99154相当接近，指数分布方差的计算公式为：(1/m)的平方，当m=

18、1时代入计算得1与0. 97058相当符合。作业三：雷达系统仿真雷达系统仿真流程为：目标回波地物杂波系统噪声天线高 放混 频中 放相干检波A/DA/D脉冲压缩缓存MTI取模积累CFAR雷达目标回波：这里的雷达回波由三部分组成即：目标回波、地物杂波、系统噪声。其中目标回波这里采用线性调频信号模拟，地物杂波采用相关对数-正态分布的随机数，系统噪声采用高斯白噪声。回波信号由目标回波（动目标），地物杂波（静目标），及系统高斯白噪声组成。线性调频信号：x=rect(t/mk)exp(jkt²) (k=B/mk)目标回波：yrect(t/mk)*exp(j*2*pi*(f1+k*t/2).*t+

19、fd*i*T)地物杂波（静目标）：yrect(t/mk)*exp(j*2*pi*(f1+k*t/2).*t)系统噪声（高斯白噪声）： z0.2*randn(1,N)。1雷达的回波模型：目标回波 十杂波回波 回波噪声 加噪声的目标回波的图形： 2放大器：其中高频放大器的增益为k，带通滤波器的设计采用的是FIR滤波器，其中通频带的带宽取为信号的带宽B的5/3倍,过渡带的宽度为B，采用汉明窗进行加权，数字滤波器的通频带和阻带的值利用（为数字角频率，为模拟角频率，为采样周期）计算。高放图：3混频器： X X（t）带通滤波器 参考源其中，参考源为：，两路信号相乘以后的输出为：，这里取带通滤波器的通带的中

20、心频率为于是通过带通滤波器以后的输出为：，于是目标回波经高频放大器以后通过混频器后的时域图形频域以及滤波器的频率响应分别如下图所示：4中频放大器的设计：中频放大器的设计与高频放大器的设计除滤波器的通带范围不同以外其余的均相同，模型框图与高频放大器完全相同，回波经混频以后通过中频放大器后的时域图形频域以及滤波器的频率响应分别如下图所示：5相位检波器：低通滤波器的设计与前面的滤波器设计相类似，输出的两路信号就是I、Q两路信号，通过检波后信号的时域图形、包络以及频率域的图形如下：6A/D变换器 采样频率为5MHz，x0=(Vmax/2ª)*intxi*2ª / Vmax ，其中a

21、为AD位数7脉冲压缩 脉冲压缩就是把输入的I、Q两路信号中I路信号作为实部Q路信号作为虚部和起来对其压缩处理。脉冲压缩后的图形：8对消器（MTI）：MTI采用一次对消：y(n)=x(n)-x(n-1)； n=1,2,3N对消后的图形如下：9取模：，取模后的图形如下：10积累（视频积累，非相参）：单极点：；双极点：，积累后的图形如下：11恒虚警（CFAR）：其程序如下：function t,s,g,f0,fs,f1=huibo%产生目标回波信号x,系统噪声y，地物杂波z以及回波pf0=3*107; %发射信号频率w=0; %发射信号初始相位c=3*108; %光速l=c/f0; %雷达信号波长(

22、载波波长)R=40000; %目标范围Vd=200; %雷达与目标之间的径向速度fd=2*Vd/l; %多普勒频率Tr=600/f0; %脉冲重复周期N=10; %雷达脉冲串长度f1=f0/5; %调频带宽是发射信号频率的1/5k=10*f1/Tr;fs=3*f0; %仿真采样频率Ts=1/fs;Btar=4*pi*R/l; M=floor(Tr*fs); %一个脉冲重复周期内的采样点数M=600mt1=floor(2*Tr*fs/5);%mt1=720mt2=floor(3*Tr*fs/10);%mt2=540mt3=floor(7*Tr*fs/10);%mt3=1260mt4=floor(

23、3*Tr*fs/5);%mt4=1080mt5=mt1-mt2;%mt5=180Vgain=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); endendfor i=1：M/10 t(i)=V

24、gain*cos(2*pi*(k*Ts*(M/10-i)/2*(M/10-i)*Ts);endfor n=1：N y(1：M,n)=gaussian(M)' %系统噪声服从高斯分布endfor n=1：N z(1：M,n)=swerling2pu(M)' %地物杂波服从swerling2型分布ends=v+g+y+z; m,n=size(s);x=zeros(1,m*n);q=zeros(1,m*n);x=s(：)'q=g(：)'p=x;i=0：length(p)-1;subplot(2,1,1);plot(i,x),title('目标回波信号'

25、); %目标回波信号xsubplot(2,1,2); l=canshu(x);b=x-l(1);plotpu(b),title('目标回波频谱'); function y=gaofang(s,f0,fs,f1)%高频放大器 M,N=size(s);Vgain=3; %高放增益w1=(f0-4*f1/3)*2*pi/fs; %w1=26pi/45w3=(f0-f1/3)*2*pi/fs; %w3=29pi/45w2=(f0+7*f1/3)*2*pi/fs; %w2=37pi/45w4=(f0+4*f1/3)*2*pi/fs; %w4=34pi/45th=min(w3-w1),(w

26、2-w4);%w3-w1=pi/15,w2-w4=pi/15M1=ceil(6.6*pi/th)+1; %M1=99w5=(w1+w3)/2;w6=(w4+w2)/2;%w5=11pi/18,w6=71pi/90h=wide(w6,M1)-wide(w5,M1);w=(hamming(M1)'h=h.*w;L=length(h);s=Vgain*s;z=zeros(M+L-1,N);y=zeros(M,N);K=ceil(L/2);for i=1：N for j=1：10 z(1+(j-1)*M/10：j*M/10+L-1,i)=conv(h,s(1+(j-1)*M/10：j*M/10

27、,i)')' y(1+(j-1)*M/10：j*M/10,i)=z(1+(j-1)*M/10+K：j*M/10+K,i); endendw=y(：)'subplot(3,1,1);i=0：length(w)-1;plot(i,w),title('高放');subplot(3,1,2);plotpu(w),title('频谱');subplot(3,1,3);plotpu(h),title('幅频特性');function y=hunpin(s,f0,fs,f1)%进行混频，输出为中频信号 M,N=size(s);i=0：M

28、-1;z1=cos(2*pi*2*f0*i/3/fs);y1=zeros(M,N);for i=1：N y1(1：M,i)=(z1.*s(1：M,i)')'endw1=(f0-5*f1)*2*pi/(3*fs);%w1=pi/9w3=(f0-2*f1)*2*pi/(3*fs);%w3=1.6pi/9w2=(f0+8*f1)*2*pi/(3*fs);%w2=0.4piw4=(f0+5*f1)*2*pi/(3*fs);%w4=pi/3th=min(w3-w1),(w2-w4); %w3-w1=0.2pi/3,w2-w4=0.2pi/3M1=ceil(6.6*pi/th)+1; %M

29、1=99w5=(w1+w3)/2;w6=(w4+w2)/2;%w5=1.3pi/9,w6=1.1pi/3h=wide(w6,M1)-wide(w5,M1);wth=(hamming(M1)'h=h.*wth;L=length(h);z=zeros(M+L-1,N);y=zeros(M,N);K=ceil(L/2);for i=1：N for j=1：10 z(1+(j-1)*M/10：j*M/10+L-1,i)=conv(h,s(1+(j-1)*M/10：j*M/10,i)')' y(1+(j-1)*M/10：j*M/10,i)=z(1+(j-1)*M/10+K：j*M

30、/10+K,i); endendw=y(：)'v=y1(：)'subplot(4,1,1);i=0：length(w)-1;plot(i,w),title('混频');subplot(4,1,2);plotpu(v),title('频谱');subplot(4,1,3);plotpu(w);subplot(4,1,4);plotpu(h),title('幅频特性'); function y=zhongfang(s,f0,fs,f1)%中频放大器 M,N=size(s);Vgain=5; %中放增益w1=(f0-5*f1)*2*pi

31、/(3*fs);%w1=pi/9w3=(f0-2*f1)*2*pi/(3*fs);%w3=1.6pi/9w2=(f0+8*f1)*2*pi/(3*fs);%w2=0.4piw4=(f0+5*f1)*2*pi/(3*fs);%w4=pi/3th=min(w3-w1),(w2-w4); %w3-w1=0.2pi/3,w2-w4=0.2pi/3M1=ceil(6.6*pi/th)+1; %M1=99w5=(w1+w3)/2;w6=(w4+w2)/2;%w5=1.3pi/9,w6=1.1pi/3h=wide(w6,M1)-wide(w5,M1);w=(hamming(M1)'h=h.*w;L=

32、length(h);s=Vgain*s;z=zeros(M+L-1,N);y=zeros(M,N);K=ceil(L/2);for i=1：N for j=1：10 z(1+(j-1)*M/10：j*M/10+L-1,i)=conv(h,s(1+(j-1)*M/10：j*M/10,i)')' y(1+(j-1)*M/10：j*M/10,i)=z(1+(j-1)*M/10+K：j*M/10+K,i); endendw=y(：)'subplot(3,1,1);i=0：length(w)-1;plot(i,w),title('中放');subplot(3,1,

33、2);plotpu(w),title('频谱');subplot(3,1,3);plotpu(h),title('幅频特性');function I,Q=xiangganjianbo(s,fs,f0,f1)%相位相干检波 M,N=size(s);i=0：M-1;z1=cos(2*pi*f0*i/(3*fs);z2=sin(2*pi*f0*i/(3*fs);y1=zeros(M,N);y2=zeros(M,N);for i=1：N y1(1：M,i)=(z1.*s(1：M,i)')' y2(1：M,i)=(z2.*s(1：M,i)')

34、9;endwt1=2*pi*(4*f1)/(3*fs);%wt1=4pi/45wt2=2*pi*(7*f1)/(3*fs);%wt2=7pi/45th=wt2-wt1;M1=ceil(6.6*pi/th)+1;%M1=70wt3=(wt1+wt2)/2;%wt3=11pi/90hd=wide(wt3,M1);w=(hamming(M1)'h=hd.*w;L=length(h);z1=zeros(M+L-1,N);z2=zeros(M+L-1,N);I=zeros(M,N);Q=zeros(M,N);K=ceil(L/2);for i=1：N for j=1：10 z1(1+(j-1)*

35、M/10：j*M/10+L-1,i)=conv(h,y1(1+(j-1)*M/10：j*M/10,i)')' z2(1+(j-1)*M/10：j*M/10+L-1,i)=conv(h,y2(1+(j-1)*M/10：j*M/10,i)')' I(1+(j-1)*M/10：j*M/10,i)=z1(1+(j-1)*M/10+K：j*M/10+K,i); Q(1+(j-1)*M/10：j*M/10,i)=z2(1+(j-1)*M/10+K：j*M/10+K,i); endendw1=I(：)'w2=Q(：)'subplot(4,1,1);j=0：le

36、ngth(w1)-1;plot(j,w1),title('I路信号');subplot(4,1,2);plot(j,w2),title('Q路信号');subplot(4,1,3);w=w1+i*w2;plotpu(w),title('频谱');subplot(4,1,4);plotpu(h),title('幅频特性');function x,y=AD(z1,z2)%模数转换m1=min(min(z1);m2=min(min(z2);Vmax=6;N=12;x=Vmax/(2N)*floor(z1-m1)*2N/Vmax);y=V

37、max/(2N)*floor(z2-m2)*2N/Vmax);M,N=size(z1);for i=1：N x(1：M,i)=Vmax/(2N)*floor(z1(1：M,i)-m1)*2N/Vmax); y(1：M,i)=Vmax/(2N)*floor(z2(1：M,i)-m2)*2N/Vmax);endw1=x(：)'w2=y(：)'subplot(2,1,1);plot(w1);subplot(2,1,2);plot(w2);function I1,Q1=maichongyasuo(I,Q,h)%对正交两路信号进行脉冲压缩 M,N=size(I);K=M/10;f0=3*

38、107;fs=3*f0; Tr=600/f0;mt2=floor(3*Tr*fs/10);mt4=floor(3*Tr*fs/5);i1=I(mt2+1：mt2+K,1：N);q1=Q(mt2+1：mt2+K,1：N);for i=1：N S1(1：K,i)=fft(i1(1：K,i); S2(1：K,i)=fft(q1(1：K,i); w_ham=(hamming(K); h1=h.*w_ham' H1=fft(h); s(1：K,i)=(S1(1：K,i)+j*S2(1：K,i).*H1' S3(1：K,i)=ifft(s(1：K,i);endi2=I(mt4+1：mt4+

39、K,1：N);q2=Q(mt4+1：mt4+K,1：N);for i=1：N P1(1：K,i)=fft(i2(1：K,i); P2(1：K,i)=fft(q2(1：K,i); w_ham=(hamming(K); h1=h.*w_ham' H1=fft(h); p(1：K,i)=(P1(1：K,i)+j*P2(1：K,i).*H1' P3(1：K,i)=ifft(p(1：K,i);endI1=I(1：mt2,1：N)',real(S3)',I(mt2+K+1：mt4,1：N)',real(P3)',I(mt4+K+1：M,1：N)'&#

40、39;Q1=I(1：mt2,1：N)',imag(S3)',I(mt2+K+1：mt4,1：N)',imag(P3)',I(mt4+K+1：M,1：N)''w1=I1(：)'w2=Q1(：)'w3=sqrt(w1.2+w2.2);i=0：N*M-1;subplot(3,1,1);plot(i,w1),title('I路信号');subplot(3,1,2);plot(i,w2),title('Q路信号');subplot(3,1,3);plot(i,w3),title('两路合成一路');function s1,s2=MTI(x1,x2)%做动目标检测 M,N=size(x1);s1=zeros(M,N);for i=1：M for j=1：N if(i=1) s1(i,j)=x1(i,j); else s1(i,j)=x1(i,j)-x1(i-1