下载文档原格式
1名词解释：（1）系统：按照某种规律组合起来，互相作用、互相依存的所有实体的集合或总和（2）连续系统：系统状态量随时间连续变化，可以通过微分方程或者偏微分方程来描述。
（3）离散事件系统：系统状态是在离散的随机时点上发生变化，且状态在一段时间内保持不变（4）系统仿真过程：建立模型并通过模型在计算机上的运行对模型进行检验、修正和分析的过程2、什么是系统建模与仿真技术？系统建模与仿真技术是以相似原理、模型理论、系统技术、信息技术以及建模与仿真应用领域的有关专业技术为基础，以计算机系统、与应用相关的设备及仿真器为工具，利用模型参与已有或设想的系统进行研究、分析、设计、加工、生产、试验、运行、评估、维护和报废（全生命周期）活动的一门多学科的综合技术。
3、画图说明计算机仿真的三要素及三个基本活动。
系统仿真有三个基本的活动，即系统建模、仿真建模和仿真实验，联系这三个活动的是仿真三要素：系统、模型、计算机（包括硬件和软件）。
它们关系如图所示。
4、什么是数学模型的有效性？解释复制有效、预测有效和结构有效的含义。
数学模型所产生的行为数据与实际过程系统数据源的相似程度称为模型的有效性。
通常数学模型的有效性按复制有效、预测有效和结构有效分为三级，后面的相似程度高于前面的相似程度（1）若数学模型产生的数据与过程系统数据源相匹配，称为复制有效。
（2）在过程系统数据源取得之前，可以得到数学模型产生的数据与过程系统数据源的匹配情况，称为预测有效。
（3）数学模型不仅具有预测有效特性，而且可以反映出产生这些行为数据的内在原因，称为结构有效。
5、动态数学模型求解的实时性要求是什么？常用哪些方法提高实时性？动态数学模型运行特点是按选定的积分时间步长，每跨进一个步长，需将全部数据模型求解一遍，一直运行到收到停止命令。
经验证明：积分步长选1s 可以达到实时要求。
提高模型实时性常用的方法有：（1）通过预先试算找出规律，尽量避开非线性代数方程组的迭代计算；（2）使用回归或者辨识的方法获取简化降阶模型；（3）使用欧拉法求解高阶微分方程；（4）偏微分方程简化为常微分方程；（5）采用稳态加动态补偿方法获取动态响应。
实验1 Witness仿真软件认识一、实验目的熟悉Witness 的启动；熟悉Witness2006用户界面；熟悉Witness 建模元素；熟悉Witness 建模与仿真过程。
二、实验内容1、运行witness软件，了解软件界面及组成；2、以一个简单流水线实例进行操作。
小部件（widget）要经过称重、冲洗、加工和检测等操作。
执行完每一步操作后小部件通过充当运输工具和缓存器的传送带（conveyer）传送至下一个操作单元。
小部件在经过最后一道工序“检测”以后，脱离本模型系统。
三、实验步骤仿真实例操作:模型元素说明：widget 为加工的小部件名称；weigh、wash、produce、inspect 为四种加工机器，每种机器只有一台；C1、C2、C3 为三条输送链；ship 是系统提供的特殊区域，表示本仿真系统之外的某个地方；操作步骤：1：将所需元素布置在界面：2：更改各元素名称：如;3:编辑各个元素的输入输出规则：4: 运行一周（5 天*8 小时*60 分钟=2400 分钟），得到统计结果。
5:仿真结果及分析：Widget：各机器工作状态统计表：分析：第一台机器效率最高位100%，第二台机器效率次之为79%，第三台和第四台机器效率低下，且空闲时间较多，可考虑加快传送带C2、C3的传送速度以及提高第二台机器的工作效率，以此来提高第三台和第四台机器的工作效率。
6：实验小结：通过本次实验，我对Witness的操作界面及基本操作有了一个初步的掌握，同学会了对于一个简单的流水线生产线进行建模仿真，总体而言，实验非常成功。
实验2 单品种流水线生产计划设计一、实验目的1.理解系统元素route的用法。
2.了解优化器optimization的用法。
3.了解单品种流水线生产计划的设计。
4.找出高生产效率、低临时库存的方案。
二、实验内容某一个车间有5台不同机器，加工一种产品。
该种产品都要求完成7道工序，而每道工序必须在指定的机器上按照事先规定好的工艺顺序进行。
实验二 动态系统的Simulink 仿真一、实验目的：1、掌握Simulink 使用的基本方法；2、熟悉连续系统仿真设计的基本方法；二、实验内容：1、编写M 脚本文件编写一个M 脚本文件，绘制函数⎪⎩⎪⎨⎧>+-≤<≤=3,630,0,sin )(x x x x x x x y在区间[-5，5]中的图形。
x=-5:0.1:5; % 设定系统输入范围与仿真步长leng=length(x); % 计算系统输入序列长度for i=1:leng % 计算系统输出序列if x(i)<=0 % 逻辑判断y(i)=sin(x(i));else if (x(i)>0&&x(i)<=3)y(i)=x(i);elsey(i)=-x(i)+6;endendendplot(x,y);grid;2、编写和调用M 函数编写一个M 函数，表示出如下函数关系t=0:0.1:3;leng=length(t);for i=1:lengif t(i)<=1;y(i)=t(i).^2;elsey(i)=t(i).^(1/2);endendplot(t,y);grid;⎪⎩⎪⎨⎧>∈=1,]1 ,0[,212t u t uy并用M脚本文件调用该函数，绘制其在[0，3]区间内的图像。
3.一个生长在罐中的细菌的简单模型。
要求给各模块和信号线改名称、改颜色或增加阴影。
假定细菌的出生率和当前细菌的总数成正比，死亡率和当前的总数的平方成正比。
若以x 代表当前细菌的总数，则细菌的出生率可表示为：_birth=ratebx细菌的死亡率可表示为：2death=rate_px细菌总数的总变化率可表示为出生率与死亡率之差。
因此系统可表示为如下的微分方程形式：2px=x-bx假定h;/5==，当前细菌的总数为1000，建立其simulink模型，.0phb/05并绘制细菌总数变化图。
4.根据种群增长曲线的数学方程进行simulink仿真，并正确设置参数，绘制出种群增长的“J”型曲线和“S”型曲线。
系统建模与仿真实验报告院系：管理科学与工程学院专业：质量与可靠性工程班级：1005104学号：100510432姓名：谢纪伟实验目录一.问题描述.二.系统数据.三. 建立过程的简单流程图.四.模型实体设计.五. 建立模型.六.运行模型.七.实验改进.八.结果分析.实验报告一.问题描述.电路板生产商要引入一个新产品，需要适当扩大现有生产线的产能，因此对现有生产线进行研究，经提前分析，发现生产过程存在瓶颈，现在对此生产线进行建模，并通过用extendsim建立的模型所得到的数据对现有生产线进行分析，并通过分析得到解决问题的办法。
二.系统数据.1.根据确定的时间表，5种型号电路板按照固定批量送入生产线中，时间表每隔120min重复一次，如下表所示：电路板种类在...min进入批量电路板种类在...min进入批量1 0 20 5 80 252 20 30 1 120 203 40 25 2 140 304 60 30 ………………进料时间表2.第一步操作是通过一台清洁工作站，每一个电路板需要至少36s，至多54s 的时间，一般情况需要48s。
3.清洁后的电路板装入自动插件机中，这台机器最多能同是处理6个电路板，每个板耗时5min。
4.当完成大部分标准插件的工作，电路板被置于一个10m的传送带上，通过波峰焊接机。
传送带上能放下30个电路板，每分钟移动1米。
5.此外，有三个工作站，用来插件机无法完成的非标准元件。
这个操作的耗时量根据板的种类而不同，如下表：电路板种类处理时间（min）电路板种类处理时间（min）1 2.5 4 3.02 2.0 5 2.03 2.5非标准元件的处理时间6.最后一步是高温加速老化试验，在这个过程中，电路板被组合成24个一组，放入烤箱中，循环通电20min。
三.建立过程的简单流程图电路板清洁自动插件波峰焊非标准插件非标准插件非标准插件高温老化离开四.模型实体设计.模拟电路板到达模拟缓冲器模拟插件机模拟convey item模拟非标准插件机三个物体汇合在一个通道将24个电路板组成一个批量对成批的电路板进行高温老化将成批的电路板还原成单独的电路板将加工后的电路板输出五.建立模型.1.定义全局单位时间.搭建模型从选择合适的全局时间单位开始。
系统建模与仿真报告姓名：葛海军学号：0411420841系统建模与仿真作业一． 产生十种随机分布的数：1．（0-1）之间的均匀分布：概率密度函数：⎩⎨⎧≤≤=其他0101)(x x P ； 产生思想：采用乘同余法产生；具体实现方法：n n ux x =+1 （mod m ）；参数：取正整数，为初始值一般取为正整数；，或一般取b b x a a u 1253203+±；m 一般取计算机的字长，其是控制所产生随机数的精度（即：小数点后的位数）； 程序（具体程序见附录）实现中取u=11，m=100000，0x 的取值是随机赋的；参数估计：在matlab 命令窗口键入y=junyun(10240)；就可以产生10240个随机数保存在向量y 中，然后再键入zhifangtu （y ，100）（调用直方图来对其进行检验），运行结果如下：然后在计算这10240个数的均值和方差在命令窗口键入z=canshu （y ），运行结果为：z=[0.50038 0.083263]其中0.50038表示所产生的数据的均值，0.083263表示所产生数据的方差，而（0-1）之间的均匀分布的随机数的数学期望为0.5，与上面所求出的0.50038很接近，方差0.083263近似与0，于是这种产生方法已经符合要求。
2．瑞利分布随机数的产生概率密度函数：⎪⎩⎪⎨⎧<≥=-000)(2222x x e x x P xσσ； 产生思想：利用直接抽样法产生；具体实现方法：a ．先调用产生（0-1）之间的均匀分布的函数（y=junyun(n)）产生一组（0-1）之间均匀分布的随机数保存在向量x 里；b ．然后作2ln z y =-；c ．另z y σ=，于是向量y 就是要产生的瑞利分布的随机数；参数估计：在matlab 命令窗口键入y=ruili(1，10240)；就可以产生10240个随机数保存在向量y 中，然后再键入zhifangtu （y ，100）（调用直方图来对其进行检验），运行结果如下：然后在计算这10240个数的均值和方差在命令窗口键入z=canshu （y ），运行结果为：z=[1.255 0.43138]其中1.255表示所产生的数据的均值，0.43138表示所产生数据的方差，而瑞利分布的数学期望计算式为：1σσ=，代入计算得： 1.253，与上面所求出的随机数的平均值 1.2555相当接近，瑞利分布方差的计算公式为：224σπ-当1σ=时代入计算得0.42920与0.43138相当接近，于是这种产生方法已经符合要求。
无穷大功率电源供电系统仿真假设无穷大功率电源供电系统，在0.02s时刻变压器低压母线发生三相短路故障，仿真其短路电流周期分量幅值和冲击电流的大小。
线路参数L=50km，x1=0.4Ω/km，r1=0.17Ω/km；变压器Sn=20MV·A，短路电压Us%=10.5，短路损耗ΔPs=135kw，空载损耗ΔP0=22kw，空载电流I0%=0.8，变比kT=110/11，高低压绕组均为Y行联接；并设供电点电压为110KV。
其对应的Simulink仿真模型如图1-1所示。
图1-1 无穷大功率电源供电系统的Simulink仿真图表1-1 图1-1仿真电路中各模块名称及提取路径模块名提取路径无穷大功率电源Three-Phase Source SimPowerSystems/Eletrical Sources三相并联RLC负荷模块5MW SimPowerSystems/Elements串联RLC支路Three-phaseParallelRLCBranch SimPowerSystems/Elements三相故障模块Three-phase-Fault SimPowerSystems/Elements三相电压电流测量模块V-I-M SimPowerSystems/Measurements示波器模块Scope Simulink/Sinks电力系统图形用户界面Poweigui SimPowerSystems双绕组变压器模块Three-PhaseTransformer SimPowerSystems/Elements图1-2 电源模块的参数设置变压器T 采用“Three-PhaseTransformer （Two Windings ）”模型。
根据给定的数据，计算折算到110kv 侧的参数如下：变压器的电阻为2233221351101010 4.0820000s N T N PU R S ∆⨯=⨯=⨯Ω=Ω 变压器的电抗为22332%10.5110101063.5310010020000s N T N U U X S ⨯=⨯=⨯Ω=Ω⨯ 变压器的漏感：63.53/(2)0.2022 3.1450T T L X f H H π===⨯⨯变压器的励磁电阻为2233301101010 5.51022N m U R P =⨯=⨯Ω=⨯Ω∆ 变压器的励磁电抗为22330100100110101075625%0.820000N m N U X I S ⨯=⨯=⨯Ω=Ω⨯ 变压器的励磁电感为75625/(2)240.82 3.1450m m L X f H H π===⨯⨯变压器模块中的参数采用有名值则设置如图1-3所示图1-3采用有名值时变压器模块的参数设置如果要采用标幺值，则在Similink 的三相变压器模型中，一次、二次绕组漏感和电阻的标幺值以额定功率和一次、二次侧各自的额定线电压为基准值，励磁电阻和励磁电感以额定功率和一次额定线电压为基准值。
《系统建模与仿真》课程作业要求鉴于本重修课程上课对象主要为工业工程毕业班的同学，很多同学在企业进行实习，上课学生太少，因此后续课堂授课不再进行，而改由同学们自学。
自学参考资料主要有：（1）课程讲义（需要到曹旭东处去拿底稿进行复印）及其模型（在上下载）；（2）WITNESS软件操作指南，下载链接：/s/1gd43COV。
课程作业相关说明：（1）课程最终考核依据：为提交报告（和模型）+前期出席情况。
（2）提交截止时间：2014-5-23（3）提交方式：纸质打印报告一份（交由曹旭东同一收齐后，在截止时间之前联系并交给我），报告（和模型，若有）电子档发送至jiannywang@ （4）报告要求：文档结构和排版参考毕业论文的格式（5）课程作业题型(一)（A类课题）根据课程讲义独立章节的模型进行流程扩展或条件变更，进行模型的设计、仿真结果的统计分析和改善设计，并最终形成仿真设计分析报告。
(二)（B类课题）根据实习企业的生产、物流或供应链系统、也可以是服务系统，选择系统中相对独立的部分（子系统）进行系统仿真方案设计，需要包括：（i）系统现状分析，主要是描述清楚系统运作过程相关的内容，可能包括系统设施布局、组成要素（机器设备、仓库、车辆、人员、订单等）、组成要素的运作流程和数据（加工时间、运输时间、订单量等）、存在的问题（例如：设备利用率不高、系统库存过高、订单满足率过低等），以及其他相关内容。
（ii）仿真目标和可行方案，主要说明该系统进行仿真可能达成的目标，以及为了达成该目标可能采取的方案（这些方案的效果需要通过仿真方能够进行评价，在报告中可以通过其他方法进行这些方案效果进行简要分析）。
(三)（C类课题）系统建模与仿真研究报告，主要是进行系统建模与仿真研究文献综述，或系统建模与仿真的实际应用介绍等。
注：（1）三类题型优劣排序为A—〉B—〉C（2）报告成文不少于8页；（3）报告格式排版需要工整清楚，图表清晰；（4）报告需要封面页：写上课程名称、课题名称、姓名、班级、学号、日期等。
安徽工业大学
《系统建模与仿真实验》
结课作业
作业题目：多产品生产系统建模与仿真
学院：
专业：
班级：
学号：
姓名：
2016年11月1日
目录
1.问题背景 (3)
1.1背景 (3)
1.2相关数据 (3)
1.3要求 (3)
2.建模仿真步骤 (4)
2.1模型建立 (4)
2.2参数设置 (4)
2.3模型运行 (6)
2.4系统问题 (7)
3.改善方案及效果 (7)
1.问题背景
1.1背景：一个工厂有5个不同的车间（普通车间，钻床车间，铣床车间，磨床车间，检测车间），加工3种类型产品。
每种产品都在5个不同的车间完成5道工序。
1.2相关数据：
1.3要求：
（1）建立上述多产品生产系统的仿真模型；
（2）查看系统运行情况报告；
（3）如存在瓶颈，请给出改善方案。
2.建模仿真步骤2.1模型建立
2.2参数设置
（1）发生器设置
运行结果
（2）处理器设置
2.3模型运行
2.4系统问题
由上图可知在暂存区1和暂存区5都有产品堆积现象。
堆积的下一工序就是瓶颈工序。
3.改善方案及效果
对生产产品工序分析，应运并行工程原理进行相似工艺合并，可以适当加大产品周转占用量，同时根据自身生产能力，制定相应的批量。
分别取0.3，0.5，0.8，2时，系统的bode图绘制：num=[1];den=[1 8 27 38 26 0];rlocus(num,den[r,k]=rlocfind(num,denGridxlabel('Real Axis',ylabel('Imaginary Axis'title('Root Locus'选定图中根轨迹与虚轴的交点，单击鼠标左键得：selected_point =0.0021 + 0.9627ik =28.7425r =-2.8199 + 2.1667i-2.8199 - 2.1667i-2.3313-0.0145 + 0.9873i-0.0145 - 0.9873inum=[1 12]；den=[1 23 242 1220 1000];rlocus(num,den[k,r]=rlocfind(num,dengridxlabel('Real Axis',ylabel('Imaginary Axis'title('Root Locus' 选定图中根轨迹与虚轴的交点，单击鼠标左键得：selected_point =0.0059 + 9.8758ik =1.0652e+003r=-11.4165 + 2.9641i-11.4165 - 2.9641i-0.0835 + 9.9528i-0.0835 - 9.9528i3.（3）num=[0.05 1]；den=[0.0008568 0.01914 0.1714 1 0];rlocus(num,den[k,r]=rlocfind(num,denGridxlabel('Real Axis',ylabel('Imaginary Axis'title('Root Locus'选定图中根轨迹与虚轴的交点，单击鼠标左键得：selected_point =0.0237 + 8.3230ik =7.6385r =-0.0916 + 8.4713i-0.0916 - 8.4713i-11.0779 + 1.2238i-11.0779 - 1.2238i4.2先令G(s=1/s,则可得其单位阶跃响应波形图为然后逐步添加如下：第一步、添加共轭极点-1+j1和-1-j1得到G(s=1/[s(s2+2s+2]，运行可得其单位阶跃响应波形为第二步、添加共轭极点-3+j2和-3-j2得到G(s=1/[s(s2+2s+2( s2+6s+13]，运行后可得其单位阶跃响应波形为（2）通过添加零、极点凑系统：先令G(s=1/(s+1,则可得其单位阶跃响应波形然后逐步添加如下：第一步、添加共轭极点-6+j8和-6-j8得到G(s=1/[(s+1(s2+12s+100]，运行后可得其单位阶跃响应波形为第二步、添加极点-10得到G(s=1/[(s+1(s2+12s+100(s+10]，运行后可得其单位阶跃响应波形为第三步、添加零点-12得到G(s=(s+12/[(s+1(s2+12s+100(s+10], 运行后可得其单位阶跃响应波形为（3）通过添加零、极点凑系统：先令G(s=1/s,则可得其单位阶跃响应波形图为然后逐步添加如下：第一步、添加极点-1/0.0714得到G(s=1/[s(0.0714s+1], 运行后可得其单位阶跃响应波形为第二步、添加一对共轭极点，即分子添加项（0.012s2+0.1s+1）后可得到G(s=1/[s(0.0714s+1( 0.012s2+0.1s+1]运行后可得其单位阶跃响应波形为第三步、添加极点-20得到G(s=1/[s(0.0714s+1( 0.012s2+0.1s+1(0.05s+1]，运行后可得其单位阶跃响应波形为4.11 分别取，0.3，0.5，0.8，2时，系统的bode图绘制：2.（1）Bode图的绘制：num=[0 0 0 0 10 ]；den=[5 24 -5 0 0]；bode(num,den gridM a g n i t u d e (d B )10-210-110101102P h a s e (d e g )Bode DiagramFrequency (rad/sec)②Nyquist 图的绘制： num=[ 0 0 0 0 10]; den=[ 5 24 -5 0 0]; [z,p,k]=tf2zp(num,den; p p = 0 0 -5.0000 0.2000③Nichols图的绘制num=[ 0 0 0 0 10];den=[ 5 24 -5 0 0]；[mag,phase]=nichols(num,den; plot(phase,20*log10(mag ngrid④Step曲线的绘制num=[ 0 0 0 0 10];den=[ 5 24 -5 0 0];step(num,dengridnum=[0 0 0 0 8 8 ];1 bode曲线的绘制：num=[0 0 0 0 8 8 ];den=[1 21 100 150 0 0 ]; bode(num,dengrid②Nyquist曲线的绘制：num=[ 0 0 0 0 8 8];den=[ 1 21 100 150 0 0]; [z,p,k]=tf2zp(num,den; nyquist(num,denp =-15.0000-3.0000 + 1.0000i-3.0000 - 1.0000i③Nichols图的绘制：num=[ 0 0 0 0 8 8];den=[ 1 21 100 150 0 0]; [mag,phase]=nichols(num,den; plot(phase,20*log10(mag ngrid④Step曲线的绘制num=[ 0 0 0 0 8 8];den=[ 1 21 100 150 0 0]; step(num,dengrid①bode曲线的绘制：num=[ 0 0 0 1.333 4];den=[0.0001 0.008 0.17 1 0]; bode(num,dengrid②Nyquist曲线的绘制：num=[ 0 0 0 1.333 4];den=[0.0001 0.008 0.17 1 0]; [z,p,k]=tf2zp(num,den; p nyquist(num,denp =-50.0000-20.0000-10.0000③Nichols图的绘制：num=[ 0 0 0 1.333 4 ];den=[0.0001 0.008 0.17 1 0]; [mag,phase]=nichols(num,den; plot(phase,20*log10(mag; ngrid④Step曲线的绘制num=[ 0 0 0 1.333 4];den=[0.0001 0.008 0.17 1 0];step(num,dengridnum=[ 0 0 1 1];den=[0.1 1 0 0];[gm,pm,wcg,wcp]=margin(num,den; gm,pm,wcg,wcpgm =Infpm =44.4594wcg =Infwcp =1.26475.3num0=10;>> den0=conv([1,0],conv([1,1],[1,2];>> w=logspace(-1,1.2;>> [gm1,pm1,wcg1,wcp1]=margin(num0,den0; Warning: The closed-loop system is unstable. > In lti.margin at 89In margin at 92>> [mag1,phase1]=bode(num0,den0,w;>> [gm1,pm1,wcg1,wcp1]ans =0.6000 -12.9919 1.4142 1.8020>> margin(num0,den0>> grid;>> wc=1.41;>> beit=10;>> T2=10/wc;>> lw=20*log10(w/1.41-4.44;>> [il,ii]=min(abs(lw+20; w1=w(ii;>> numc1=[1/w1,1];denc1=[1/ (beit*w1,1];>> numc2=[ T2,1];denc2=[ beit*T2,1];>> [numc,denc]=series(numc1,denc1,numc2,denc2;>> [num,den]=series(num0,den0,numc,denc;printsys(numc,dencnum/den =31.0168 s^2 + 11.4656 s + 1---------------------------31.0168 s^2 + 71.3593 s + 1>> disp('校正之后的系统开环传递函数为:';printsys(num,den校正之后的系统开环传递函数为:num/den =310.1682 s^2 + 114.6557 s + 10--------------------------------------------------------------31.0168 s^5 + 164.4098 s^4 + 277.1116 s^3 + 145.7186 s^2 + 2 s>> [mag2,phase2]=bode(numc,denc,w;>> [mag,phase]=bode(num,den,w;>> [gm,pm,wcg,wcp]=margin(num,den;>> subplot(2,1,1;semilogx(w,20*log10(mag,w,20*log10(mag1,'--',w,20*log10(mag2,'-.'; >> grid; ylabel('幅值(db'; title('--Go,-Gc,GoGc';>> subplot(2,1,2; semilogx(w,phase,w,phase1,'--',w,phase2,'-',w,(w-180-w,':';>> grid; ylabel('相位(0'; xlabel('频率(rad/sec';>> title(['校正后:幅值裕量=',num2str(20*log10(gm,'db','相位裕量=',num2str(pm,'0'];。
作业一：统计试验法解决报童问题 1.问题提出 报童问题是决策方面的一个简单例子。对一个报童来说，他每天由报社买来报纸，然后到街上去卖，当然他希望获得最大的利润。如果把报童问题看作一个系统，那么报纸、顾客、报童和利润就成了该系统的重要组成部分。问题在于，根据市场需求，报童寻求一个什么样的定货法则或策略才能使他所获得的利润最大。因为定货多了，如果市场需求小，卖不了将导致利润下降，甚至亏本；定货少了，如果市场需求量大，又失去了赚钱的机会。这样一来，定货法则或定货策略就成了报童能否赚取最大利润的关键。当然，最好是每天需要多少份就定多少份，但市场需求量是一个随机变量，这就需要每天用统计方法作出决策。 报童以每份0.5元的价格买进报纸，以0.8元的价格出售。根据长期统计，报童根据以前的卖报记录知道，每天的需求量有几种可能及出现的概率。对现有数据分析，得出报童每天最佳买进报纸量，使报童的平均总收入最大。 2.模型建立与求解 假定， Bn— 当天买进或订购的报纸数量
Sn— 当天社会需要报纸的数量，显然，它是一个随机变量
Dn— 当天卖掉的报纸数量
再假定，报童每天买进和卖出每份报纸的价格分别用 PB和PS表示，且 PS>PB，即卖出价大于买入价，则第 n 天的利润为 Pn=SnPS-BnPB 报童决定：当天订购报纸的数量等于前一天的市场需求量，即 Bn=Dn-1
而当天卖掉的报纸的数量n S 则由以下两个条件来决定
Dn ≤ Bn时，Sn=Dn Dn> Bn时，Sn=Bn 即是说，如果当天订购的报纸的数量大于需求量，当天卖掉的报纸的数量只能等于需求量；如果当天买的报纸的数量小于需求量，即当天卖掉的报纸的数量等于当天订购量和当天需求量二者中的小者。 现在的问题是，前一天的需求量如何决定。报童决定利用过去一年的统计数字确定Dn-1。报童根据以前的卖报记录知道，每天的需求量有以下几种可能：40份，41份，42份，43份，44份， 45份，46份；并且统计出了相对频数，得到如右表的一组数据。 需求量Dn 40 41 42 43 44 45 46 相对频数 Pn 0.05 0.1 0.2 0.3 0.15 0.1 0.1 其需求量的平均值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份，由于Bn3.实验结果 报童问题的10天的仿真结果
Sn Pn ∑Pn 44 12.7 12.7 43 12.4 25.1 44 12.1 37.2 40 11.0 48.2 40 12.0 60.2 41 12.3 72.5 41 11.3 83.8 41 12.3 96.1 42 12.6 108.7 43 11.4 120.1 表采用不同的决策准则时报童所获得的利润
N0 17天的利润 100天的利润 500天的利润 报童使用的决策规则
1 208.3元 1215.3元 6122.1元 Bn=Dn-1 2 207.0元 1223.8元 6159.8元 Bn=（Dn-1+Dn-2)/2
3 209.7元 1242.8元 6234.8元 Bn= 43
4 209.8元 1241.1元 6235.4元 B1=D0 B2=(D0+D1)/2 Bn=(D0+D1+…+Dn-1)]/n
说明： Bn-每天买入的报纸数量 Dn-每天社会的需求量 PB-每份报纸的买入价 Ps-每份报纸的卖出价 过去累计的需求量及相应的概率： 40 41 42 43 44 45 46 0.05 0.10 0.20 0.30 0.15 0.10 0.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_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; 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*****Bn=Dn-1***************/ // if(i!=N) B[i+1]=D[i]; /***2***Bn=（Dn-1+Dn-2)/2*********/ /* if(i==0)B[i+1]=D[i]; if(i!=N) B[i+2]=(D[i]+D[i+1])/2;*/ /***3****** Bn= 43**************/ // B[i]=43; /***4****Bn=(D0+D1+…+Dn-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）之间的均匀分布：
概率密度函数： 其他0101)(x
xP
；
产生思想：采用乘同余法产生； 具体实现方法：xn+1=uxn(mod n)； 参数：u一般取23a±3或5，a为正整数；x0为初始值一般取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．高斯分布： 高斯分布的概率密度函数如下： 22
2)(21)(uxexp
其产生方法是在均匀分布随机数的基础上通过函数变换法来产生。产生步骤是①产生均匀分布的随机数。②产生服从标准正态分布的随机数。③由标准正态分布产生一般正态分布。 )2cos(ln2iiivuz
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近似，于是这种产生方法已经符合要求。
